arch/m68k/fpsp040/stan.S

   1 |
   2 |       stan.sa 3.3 7/29/91
   3 |
   4 |       The entry point stan computes the tangent of
   5 |       an input argument;
   6 |       stand does the same except for denormalized input.
   7 |
   8 |       Input: Double-extended number X in location pointed to
   9 |               by address register a0.
  10 |
  11 |       Output: The value tan(X) returned in floating-point register Fp0.
  12 |
  13 |       Accuracy and Monotonicity: The returned result is within 3 ulp in
  14 |               64 significant bit, i.e. within 0.5001 ulp to 53 bits if the
  15 |               result is subsequently rounded to double precision. The
  16 |               result is provably monotonic in double precision.
  17 |
  18 |       Speed: The program sTAN takes approximately 170 cycles for
  19 |               input argument X such that |X| < 15Pi, which is the usual
  20 |               situation.
  21 |
  22 |       Algorithm:
  23 |
  24 |       1. If |X| >= 15Pi or |X| < 2**(-40), go to 6.
  25 |
  26 |       2. Decompose X as X = N(Pi/2) + r where |r| <= Pi/4. Let
  27 |               k = N mod 2, so in particular, k = 0 or 1.
  28 |
  29 |       3. If k is odd, go to 5.
  30 |
  31 |       4. (k is even) Tan(X) = tan(r) and tan(r) is approximated by a
  32 |               rational function U/V where
  33 |               U = r + r*s*(P1 + s*(P2 + s*P3)), and
  34 |               V = 1 + s*(Q1 + s*(Q2 + s*(Q3 + s*Q4))),  s = r*r.
  35 |               Exit.
  36 |
  37 |       4. (k is odd) Tan(X) = -cot(r). Since tan(r) is approximated by a
  38 |               rational function U/V where
  39 |               U = r + r*s*(P1 + s*(P2 + s*P3)), and
  40 |               V = 1 + s*(Q1 + s*(Q2 + s*(Q3 + s*Q4))), s = r*r,
  41 |               -Cot(r) = -V/U. Exit.
  42 |
  43 |       6. If |X| > 1, go to 8.
  44 |
  45 |       7. (|X|<2**(-40)) Tan(X) = X. Exit.
  46 |
  47 |       8. Overwrite X by X := X rem 2Pi. Now that |X| <= Pi, go back to 2.
  48 |
  49
  50 |               Copyright (C) Motorola, Inc. 1990
  51 |                       All Rights Reserved
  52 |
  53 |       THIS IS UNPUBLISHED PROPRIETARY SOURCE CODE OF MOTOROLA
  54 |       The copyright notice above does not evidence any
  55 |       actual or intended publication of such source code.
  56
  57 |STAN   idnt    2,1 | Motorola 040 Floating Point Software Package
  58
  59         |section        8
  60
  61 #include "fpsp.h"
  62
  63 BOUNDS1:        .long 0x3FD78000,0x4004BC7E
  64 TWOBYPI:        .long 0x3FE45F30,0x6DC9C883
  65
  66 TANQ4:  .long 0x3EA0B759,0xF50F8688
  67 TANP3:  .long 0xBEF2BAA5,0xA8924F04
  68
  69 TANQ3:  .long 0xBF346F59,0xB39BA65F,0x00000000,0x00000000
  70
  71 TANP2:  .long 0x3FF60000,0xE073D3FC,0x199C4A00,0x00000000
  72
  73 TANQ2:  .long 0x3FF90000,0xD23CD684,0x15D95FA1,0x00000000
  74
  75 TANP1:  .long 0xBFFC0000,0x8895A6C5,0xFB423BCA,0x00000000
  76
  77 TANQ1:  .long 0xBFFD0000,0xEEF57E0D,0xA84BC8CE,0x00000000
  78
  79 INVTWOPI: .long 0x3FFC0000,0xA2F9836E,0x4E44152A,0x00000000
  80
  81 TWOPI1: .long 0x40010000,0xC90FDAA2,0x00000000,0x00000000
  82 TWOPI2: .long 0x3FDF0000,0x85A308D4,0x00000000,0x00000000
  83
  84 |--N*PI/2, -32 <= N <= 32, IN A LEADING TERM IN EXT. AND TRAILING
  85 |--TERM IN SGL. NOTE THAT PI IS 64-BIT LONG, THUS N*PI/2 IS AT
  86 |--MOST 69 BITS LONG.
  87         .global PITBL
  88 PITBL:
  89   .long  0xC0040000,0xC90FDAA2,0x2168C235,0x21800000
  90   .long  0xC0040000,0xC2C75BCD,0x105D7C23,0xA0D00000
  91   .long  0xC0040000,0xBC7EDCF7,0xFF523611,0xA1E80000
  92   .long  0xC0040000,0xB6365E22,0xEE46F000,0x21480000
  93   .long  0xC0040000,0xAFEDDF4D,0xDD3BA9EE,0xA1200000
  94   .long  0xC0040000,0xA9A56078,0xCC3063DD,0x21FC0000
  95   .long  0xC0040000,0xA35CE1A3,0xBB251DCB,0x21100000
  96   .long  0xC0040000,0x9D1462CE,0xAA19D7B9,0xA1580000
  97   .long  0xC0040000,0x96CBE3F9,0x990E91A8,0x21E00000
  98   .long  0xC0040000,0x90836524,0x88034B96,0x20B00000
  99   .long  0xC0040000,0x8A3AE64F,0x76F80584,0xA1880000
 100   .long  0xC0040000,0x83F2677A,0x65ECBF73,0x21C40000
 101   .long  0xC0030000,0xFB53D14A,0xA9C2F2C2,0x20000000
 102   .long  0xC0030000,0xEEC2D3A0,0x87AC669F,0x21380000
 103   .long  0xC0030000,0xE231D5F6,0x6595DA7B,0xA1300000
 104   .long  0xC0030000,0xD5A0D84C,0x437F4E58,0x9FC00000
 105   .long  0xC0030000,0xC90FDAA2,0x2168C235,0x21000000
 106   .long  0xC0030000,0xBC7EDCF7,0xFF523611,0xA1680000
 107   .long  0xC0030000,0xAFEDDF4D,0xDD3BA9EE,0xA0A00000
 108   .long  0xC0030000,0xA35CE1A3,0xBB251DCB,0x20900000
 109   .long  0xC0030000,0x96CBE3F9,0x990E91A8,0x21600000
 110   .long  0xC0030000,0x8A3AE64F,0x76F80584,0xA1080000
 111   .long  0xC0020000,0xFB53D14A,0xA9C2F2C2,0x1F800000
 112   .long  0xC0020000,0xE231D5F6,0x6595DA7B,0xA0B00000
 113   .long  0xC0020000,0xC90FDAA2,0x2168C235,0x20800000
 114   .long  0xC0020000,0xAFEDDF4D,0xDD3BA9EE,0xA0200000
 115   .long  0xC0020000,0x96CBE3F9,0x990E91A8,0x20E00000
 116   .long  0xC0010000,0xFB53D14A,0xA9C2F2C2,0x1F000000
 117   .long  0xC0010000,0xC90FDAA2,0x2168C235,0x20000000
 118   .long  0xC0010000,0x96CBE3F9,0x990E91A8,0x20600000
 119   .long  0xC0000000,0xC90FDAA2,0x2168C235,0x1F800000
 120   .long  0xBFFF0000,0xC90FDAA2,0x2168C235,0x1F000000
 121   .long  0x00000000,0x00000000,0x00000000,0x00000000
 122   .long  0x3FFF0000,0xC90FDAA2,0x2168C235,0x9F000000
 123   .long  0x40000000,0xC90FDAA2,0x2168C235,0x9F800000
 124   .long  0x40010000,0x96CBE3F9,0x990E91A8,0xA0600000
 125   .long  0x40010000,0xC90FDAA2,0x2168C235,0xA0000000
 126   .long  0x40010000,0xFB53D14A,0xA9C2F2C2,0x9F000000
 127   .long  0x40020000,0x96CBE3F9,0x990E91A8,0xA0E00000
 128   .long  0x40020000,0xAFEDDF4D,0xDD3BA9EE,0x20200000
 129   .long  0x40020000,0xC90FDAA2,0x2168C235,0xA0800000
 130   .long  0x40020000,0xE231D5F6,0x6595DA7B,0x20B00000
 131   .long  0x40020000,0xFB53D14A,0xA9C2F2C2,0x9F800000
 132   .long  0x40030000,0x8A3AE64F,0x76F80584,0x21080000
 133   .long  0x40030000,0x96CBE3F9,0x990E91A8,0xA1600000
 134   .long  0x40030000,0xA35CE1A3,0xBB251DCB,0xA0900000
 135   .long  0x40030000,0xAFEDDF4D,0xDD3BA9EE,0x20A00000
 136   .long  0x40030000,0xBC7EDCF7,0xFF523611,0x21680000
 137   .long  0x40030000,0xC90FDAA2,0x2168C235,0xA1000000
 138   .long  0x40030000,0xD5A0D84C,0x437F4E58,0x1FC00000
 139   .long  0x40030000,0xE231D5F6,0x6595DA7B,0x21300000
 140   .long  0x40030000,0xEEC2D3A0,0x87AC669F,0xA1380000
 141   .long  0x40030000,0xFB53D14A,0xA9C2F2C2,0xA0000000
 142   .long  0x40040000,0x83F2677A,0x65ECBF73,0xA1C40000
 143   .long  0x40040000,0x8A3AE64F,0x76F80584,0x21880000
 144   .long  0x40040000,0x90836524,0x88034B96,0xA0B00000
 145   .long  0x40040000,0x96CBE3F9,0x990E91A8,0xA1E00000
 146   .long  0x40040000,0x9D1462CE,0xAA19D7B9,0x21580000
 147   .long  0x40040000,0xA35CE1A3,0xBB251DCB,0xA1100000
 148   .long  0x40040000,0xA9A56078,0xCC3063DD,0xA1FC0000
 149   .long  0x40040000,0xAFEDDF4D,0xDD3BA9EE,0x21200000
 150   .long  0x40040000,0xB6365E22,0xEE46F000,0xA1480000
 151   .long  0x40040000,0xBC7EDCF7,0xFF523611,0x21E80000
 152   .long  0x40040000,0xC2C75BCD,0x105D7C23,0x20D00000
 153   .long  0x40040000,0xC90FDAA2,0x2168C235,0xA1800000
 154
 155         .set    INARG,FP_SCR4
 156
 157         .set    TWOTO63,L_SCR1
 158         .set    ENDFLAG,L_SCR2
 159         .set    N,L_SCR3
 160
 161         | xref  t_frcinx
 162         |xref   t_extdnrm
 163
 164         .global stand
 165 stand:
 166 |--TAN(X) = X FOR DENORMALIZED X
 167
 168         bra             t_extdnrm
 169
 170         .global stan
 171 stan:
 172         fmovex          (%a0),%fp0      | ...LOAD INPUT
 173
 174         movel           (%a0),%d0
 175         movew           4(%a0),%d0
 176         andil           #0x7FFFFFFF,%d0
 177
 178         cmpil           #0x3FD78000,%d0         | ...|X| >= 2**(-40)?
 179         bges            TANOK1
 180         bra             TANSM
 181 TANOK1:
 182         cmpil           #0x4004BC7E,%d0         | ...|X| < 15 PI?
 183         blts            TANMAIN
 184         bra             REDUCEX
 185
 186
 187 TANMAIN:
 188 |--THIS IS THE USUAL CASE, |X| <= 15 PI.
 189 |--THE ARGUMENT REDUCTION IS DONE BY TABLE LOOK UP.
 190         fmovex          %fp0,%fp1
 191         fmuld           TWOBYPI,%fp1    | ...X*2/PI
 192
 193 |--HIDE THE NEXT TWO INSTRUCTIONS
 194         leal            PITBL+0x200,%a1 | ...TABLE OF N*PI/2, N = -32,...,32
 195
 196 |--FP1 IS NOW READY
 197         fmovel          %fp1,%d0                | ...CONVERT TO INTEGER
 198
 199         asll            #4,%d0
 200         addal           %d0,%a1         | ...ADDRESS N*PIBY2 IN Y1, Y2
 201
 202         fsubx           (%a1)+,%fp0     | ...X-Y1
 203 |--HIDE THE NEXT ONE
 204
 205         fsubs           (%a1),%fp0      | ...FP0 IS R = (X-Y1)-Y2
 206
 207         rorl            #5,%d0
 208         andil           #0x80000000,%d0 | ...D0 WAS ODD IFF D0 < 0
 209
 210 TANCONT:
 211
 212         cmpil           #0,%d0
 213         blt             NODD
 214
 215         fmovex          %fp0,%fp1
 216         fmulx           %fp1,%fp1               | ...S = R*R
 217
 218         fmoved          TANQ4,%fp3
 219         fmoved          TANP3,%fp2
 220
 221         fmulx           %fp1,%fp3               | ...SQ4
 222         fmulx           %fp1,%fp2               | ...SP3
 223
 224         faddd           TANQ3,%fp3      | ...Q3+SQ4
 225         faddx           TANP2,%fp2      | ...P2+SP3
 226
 227         fmulx           %fp1,%fp3               | ...S(Q3+SQ4)
 228         fmulx           %fp1,%fp2               | ...S(P2+SP3)
 229
 230         faddx           TANQ2,%fp3      | ...Q2+S(Q3+SQ4)
 231         faddx           TANP1,%fp2      | ...P1+S(P2+SP3)
 232
 233         fmulx           %fp1,%fp3               | ...S(Q2+S(Q3+SQ4))
 234         fmulx           %fp1,%fp2               | ...S(P1+S(P2+SP3))
 235
 236         faddx           TANQ1,%fp3      | ...Q1+S(Q2+S(Q3+SQ4))
 237         fmulx           %fp0,%fp2               | ...RS(P1+S(P2+SP3))
 238
 239         fmulx           %fp3,%fp1               | ...S(Q1+S(Q2+S(Q3+SQ4)))
 240
 241
 242         faddx           %fp2,%fp0               | ...R+RS(P1+S(P2+SP3))
 243
 244
 245         fadds           #0x3F800000,%fp1        | ...1+S(Q1+...)
 246
 247         fmovel          %d1,%fpcr               |restore users exceptions
 248         fdivx           %fp1,%fp0               |last inst - possible exception set
 249
 250         bra             t_frcinx
 251
 252 NODD:
 253         fmovex          %fp0,%fp1
 254         fmulx           %fp0,%fp0               | ...S = R*R
 255
 256         fmoved          TANQ4,%fp3
 257         fmoved          TANP3,%fp2
 258
 259         fmulx           %fp0,%fp3               | ...SQ4
 260         fmulx           %fp0,%fp2               | ...SP3
 261
 262         faddd           TANQ3,%fp3      | ...Q3+SQ4
 263         faddx           TANP2,%fp2      | ...P2+SP3
 264
 265         fmulx           %fp0,%fp3               | ...S(Q3+SQ4)
 266         fmulx           %fp0,%fp2               | ...S(P2+SP3)
 267
 268         faddx           TANQ2,%fp3      | ...Q2+S(Q3+SQ4)
 269         faddx           TANP1,%fp2      | ...P1+S(P2+SP3)
 270
 271         fmulx           %fp0,%fp3               | ...S(Q2+S(Q3+SQ4))
 272         fmulx           %fp0,%fp2               | ...S(P1+S(P2+SP3))
 273
 274         faddx           TANQ1,%fp3      | ...Q1+S(Q2+S(Q3+SQ4))
 275         fmulx           %fp1,%fp2               | ...RS(P1+S(P2+SP3))
 276
 277         fmulx           %fp3,%fp0               | ...S(Q1+S(Q2+S(Q3+SQ4)))
 278
 279
 280         faddx           %fp2,%fp1               | ...R+RS(P1+S(P2+SP3))
 281         fadds           #0x3F800000,%fp0        | ...1+S(Q1+...)
 282
 283
 284         fmovex          %fp1,-(%sp)
 285         eoril           #0x80000000,(%sp)
 286
 287         fmovel          %d1,%fpcr               |restore users exceptions
 288         fdivx           (%sp)+,%fp0     |last inst - possible exception set
 289
 290         bra             t_frcinx
 291
 292 TANBORS:
 293 |--IF |X| > 15PI, WE USE THE GENERAL ARGUMENT REDUCTION.
 294 |--IF |X| < 2**(-40), RETURN X OR 1.
 295         cmpil           #0x3FFF8000,%d0
 296         bgts            REDUCEX
 297
 298 TANSM:
 299
 300         fmovex          %fp0,-(%sp)
 301         fmovel          %d1,%fpcr                |restore users exceptions
 302         fmovex          (%sp)+,%fp0     |last inst - possible exception set
 303
 304         bra             t_frcinx
 305
 306
 307 REDUCEX:
 308 |--WHEN REDUCEX IS USED, THE CODE WILL INEVITABLY BE SLOW.
 309 |--THIS REDUCTION METHOD, HOWEVER, IS MUCH FASTER THAN USING
 310 |--THE REMAINDER INSTRUCTION WHICH IS NOW IN SOFTWARE.
 311
 312         fmovemx %fp2-%fp5,-(%a7)        | ...save FP2 through FP5
 313         movel           %d2,-(%a7)
 314         fmoves         #0x00000000,%fp1
 315
 316 |--If compact form of abs(arg) in d0=$7ffeffff, argument is so large that
 317 |--there is a danger of unwanted overflow in first LOOP iteration.  In this
 318 |--case, reduce argument by one remainder step to make subsequent reduction
 319 |--safe.
 320         cmpil   #0x7ffeffff,%d0         |is argument dangerously large?
 321         bnes    LOOP
 322         movel   #0x7ffe0000,FP_SCR2(%a6)        |yes
 323 |                                       ;create 2**16383*PI/2
 324         movel   #0xc90fdaa2,FP_SCR2+4(%a6)
 325         clrl    FP_SCR2+8(%a6)
 326         ftstx   %fp0                    |test sign of argument
 327         movel   #0x7fdc0000,FP_SCR3(%a6)        |create low half of 2**16383*
 328 |                                       ;PI/2 at FP_SCR3
 329         movel   #0x85a308d3,FP_SCR3+4(%a6)
 330         clrl   FP_SCR3+8(%a6)
 331         fblt    red_neg
 332         orw     #0x8000,FP_SCR2(%a6)    |positive arg
 333         orw     #0x8000,FP_SCR3(%a6)
 334 red_neg:
 335         faddx  FP_SCR2(%a6),%fp0                |high part of reduction is exact
 336         fmovex  %fp0,%fp1               |save high result in fp1
 337         faddx  FP_SCR3(%a6),%fp0                |low part of reduction
 338         fsubx  %fp0,%fp1                        |determine low component of result
 339         faddx  FP_SCR3(%a6),%fp1                |fp0/fp1 are reduced argument.
 340
 341 |--ON ENTRY, FP0 IS X, ON RETURN, FP0 IS X REM PI/2, |X| <= PI/4.
 342 |--integer quotient will be stored in N
 343 |--Intermediate remainder is 66-bit long; (R,r) in (FP0,FP1)
 344
 345 LOOP:
 346         fmovex          %fp0,INARG(%a6) | ...+-2**K * F, 1 <= F < 2
 347         movew           INARG(%a6),%d0
 348         movel          %d0,%a1          | ...save a copy of D0
 349         andil           #0x00007FFF,%d0
 350         subil           #0x00003FFF,%d0 | ...D0 IS K
 351         cmpil           #28,%d0
 352         bles            LASTLOOP
 353 CONTLOOP:
 354         subil           #27,%d0  | ...D0 IS L := K-27
 355         movel           #0,ENDFLAG(%a6)
 356         bras            WORK
 357 LASTLOOP:
 358         clrl            %d0             | ...D0 IS L := 0
 359         movel           #1,ENDFLAG(%a6)
 360
 361 WORK:
 362 |--FIND THE REMAINDER OF (R,r) W.R.T.   2**L * (PI/2). L IS SO CHOSEN
 363 |--THAT INT( X * (2/PI) / 2**(L) ) < 2**29.
 364
 365 |--CREATE 2**(-L) * (2/PI), SIGN(INARG)*2**(63),
 366 |--2**L * (PIby2_1), 2**L * (PIby2_2)
 367
 368         movel           #0x00003FFE,%d2 | ...BIASED EXPO OF 2/PI
 369         subl            %d0,%d2         | ...BIASED EXPO OF 2**(-L)*(2/PI)
 370
 371         movel           #0xA2F9836E,FP_SCR1+4(%a6)
 372         movel           #0x4E44152A,FP_SCR1+8(%a6)
 373         movew           %d2,FP_SCR1(%a6)        | ...FP_SCR1 is 2**(-L)*(2/PI)
 374
 375         fmovex          %fp0,%fp2
 376         fmulx           FP_SCR1(%a6),%fp2
 377 |--WE MUST NOW FIND INT(FP2). SINCE WE NEED THIS VALUE IN
 378 |--FLOATING POINT FORMAT, THE TWO FMOVE'S       FMOVE.L FP <--> N
 379 |--WILL BE TOO INEFFICIENT. THE WAY AROUND IT IS THAT
 380 |--(SIGN(INARG)*2**63   +       FP2) - SIGN(INARG)*2**63 WILL GIVE
 381 |--US THE DESIRED VALUE IN FLOATING POINT.
 382
 383 |--HIDE SIX CYCLES OF INSTRUCTION
 384         movel           %a1,%d2
 385         swap            %d2
 386         andil           #0x80000000,%d2
 387         oril            #0x5F000000,%d2 | ...D2 IS SIGN(INARG)*2**63 IN SGL
 388         movel           %d2,TWOTO63(%a6)
 389
 390         movel           %d0,%d2
 391         addil           #0x00003FFF,%d2 | ...BIASED EXPO OF 2**L * (PI/2)
 392
 393 |--FP2 IS READY
 394         fadds           TWOTO63(%a6),%fp2       | ...THE FRACTIONAL PART OF FP1 IS ROUNDED
 395
 396 |--HIDE 4 CYCLES OF INSTRUCTION; creating 2**(L)*Piby2_1  and  2**(L)*Piby2_2
 397         movew           %d2,FP_SCR2(%a6)
 398         clrw           FP_SCR2+2(%a6)
 399         movel           #0xC90FDAA2,FP_SCR2+4(%a6)
 400         clrl            FP_SCR2+8(%a6)          | ...FP_SCR2 is  2**(L) * Piby2_1
 401
 402 |--FP2 IS READY
 403         fsubs           TWOTO63(%a6),%fp2               | ...FP2 is N
 404
 405         addil           #0x00003FDD,%d0
 406         movew           %d0,FP_SCR3(%a6)
 407         clrw           FP_SCR3+2(%a6)
 408         movel           #0x85A308D3,FP_SCR3+4(%a6)
 409         clrl            FP_SCR3+8(%a6)          | ...FP_SCR3 is 2**(L) * Piby2_2
 410
 411         movel           ENDFLAG(%a6),%d0
 412
 413 |--We are now ready to perform (R+r) - N*P1 - N*P2, P1 = 2**(L) * Piby2_1 and
 414 |--P2 = 2**(L) * Piby2_2
 415         fmovex          %fp2,%fp4
 416         fmulx           FP_SCR2(%a6),%fp4               | ...W = N*P1
 417         fmovex          %fp2,%fp5
 418         fmulx           FP_SCR3(%a6),%fp5               | ...w = N*P2
 419         fmovex          %fp4,%fp3
 420 |--we want P+p = W+w  but  |p| <= half ulp of P
 421 |--Then, we need to compute  A := R-P   and  a := r-p
 422         faddx           %fp5,%fp3                       | ...FP3 is P
 423         fsubx           %fp3,%fp4                       | ...W-P
 424
 425         fsubx           %fp3,%fp0                       | ...FP0 is A := R - P
 426         faddx           %fp5,%fp4                       | ...FP4 is p = (W-P)+w
 427
 428         fmovex          %fp0,%fp3                       | ...FP3 A
 429         fsubx           %fp4,%fp1                       | ...FP1 is a := r - p
 430
 431 |--Now we need to normalize (A,a) to  "new (R,r)" where R+r = A+a but
 432 |--|r| <= half ulp of R.
 433         faddx           %fp1,%fp0                       | ...FP0 is R := A+a
 434 |--No need to calculate r if this is the last loop
 435         cmpil           #0,%d0
 436         bgt             RESTORE
 437
 438 |--Need to calculate r
 439         fsubx           %fp0,%fp3                       | ...A-R
 440         faddx           %fp3,%fp1                       | ...FP1 is r := (A-R)+a
 441         bra             LOOP
 442
 443 RESTORE:
 444         fmovel          %fp2,N(%a6)
 445         movel           (%a7)+,%d2
 446         fmovemx (%a7)+,%fp2-%fp5
 447
 448
 449         movel           N(%a6),%d0
 450         rorl            #1,%d0
 451
 452
 453         bra             TANCONT
 454
 455         |end