Annotation of sys/arch/m68k/fpsp/ssin.sa, Revision 1.1
1.1 ! nbrk 1: * $OpenBSD: ssin.sa,v 1.3 2003/11/07 10:36:10 miod Exp $
! 2: * $NetBSD: ssin.sa,v 1.3 1994/10/26 07:50:01 cgd Exp $
! 3:
! 4: * MOTOROLA MICROPROCESSOR & MEMORY TECHNOLOGY GROUP
! 5: * M68000 Hi-Performance Microprocessor Division
! 6: * M68040 Software Package
! 7: *
! 8: * M68040 Software Package Copyright (c) 1993, 1994 Motorola Inc.
! 9: * All rights reserved.
! 10: *
! 11: * THE SOFTWARE is provided on an "AS IS" basis and without warranty.
! 12: * To the maximum extent permitted by applicable law,
! 13: * MOTOROLA DISCLAIMS ALL WARRANTIES WHETHER EXPRESS OR IMPLIED,
! 14: * INCLUDING IMPLIED WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A
! 15: * PARTICULAR PURPOSE and any warranty against infringement with
! 16: * regard to the SOFTWARE (INCLUDING ANY MODIFIED VERSIONS THEREOF)
! 17: * and any accompanying written materials.
! 18: *
! 19: * To the maximum extent permitted by applicable law,
! 20: * IN NO EVENT SHALL MOTOROLA BE LIABLE FOR ANY DAMAGES WHATSOEVER
! 21: * (INCLUDING WITHOUT LIMITATION, DAMAGES FOR LOSS OF BUSINESS
! 22: * PROFITS, BUSINESS INTERRUPTION, LOSS OF BUSINESS INFORMATION, OR
! 23: * OTHER PECUNIARY LOSS) ARISING OF THE USE OR INABILITY TO USE THE
! 24: * SOFTWARE. Motorola assumes no responsibility for the maintenance
! 25: * and support of the SOFTWARE.
! 26: *
! 27: * You are hereby granted a copyright license to use, modify, and
! 28: * distribute the SOFTWARE so long as this entire notice is retained
! 29: * without alteration in any modified and/or redistributed versions,
! 30: * and that such modified versions are clearly identified as such.
! 31: * No licenses are granted by implication, estoppel or otherwise
! 32: * under any patents or trademarks of Motorola, Inc.
! 33:
! 34: *
! 35: * ssin.sa 3.3 7/29/91
! 36: *
! 37: * The entry point sSIN computes the sine of an input argument
! 38: * sCOS computes the cosine, and sSINCOS computes both. The
! 39: * corresponding entry points with a "d" computes the same
! 40: * corresponding function values for denormalized inputs.
! 41: *
! 42: * Input: Double-extended number X in location pointed to
! 43: * by address register a0.
! 44: *
! 45: * Output: The funtion value sin(X) or cos(X) returned in Fp0 if SIN or
! 46: * COS is requested. Otherwise, for SINCOS, sin(X) is returned
! 47: * in Fp0, and cos(X) is returned in Fp1.
! 48: *
! 49: * Modifies: Fp0 for SIN or COS; both Fp0 and Fp1 for SINCOS.
! 50: *
! 51: * Accuracy and Monotonicity: The returned result is within 1 ulp in
! 52: * 64 significant bit, i.e. within 0.5001 ulp to 53 bits if the
! 53: * result is subsequently rounded to double precision. The
! 54: * result is provably monotonic in double precision.
! 55: *
! 56: * Speed: The programs sSIN and sCOS take approximately 150 cycles for
! 57: * input argument X such that |X| < 15Pi, which is the usual
! 58: * situation. The speed for sSINCOS is approximately 190 cycles.
! 59: *
! 60: * Algorithm:
! 61: *
! 62: * SIN and COS:
! 63: * 1. If SIN is invoked, set AdjN := 0; otherwise, set AdjN := 1.
! 64: *
! 65: * 2. If |X| >= 15Pi or |X| < 2**(-40), go to 7.
! 66: *
! 67: * 3. Decompose X as X = N(Pi/2) + r where |r| <= Pi/4. Let
! 68: * k = N mod 4, so in particular, k = 0,1,2,or 3. Overwirte
! 69: * k by k := k + AdjN.
! 70: *
! 71: * 4. If k is even, go to 6.
! 72: *
! 73: * 5. (k is odd) Set j := (k-1)/2, sgn := (-1)**j. Return sgn*cos(r)
! 74: * where cos(r) is approximated by an even polynomial in r,
! 75: * 1 + r*r*(B1+s*(B2+ ... + s*B8)), s = r*r.
! 76: * Exit.
! 77: *
! 78: * 6. (k is even) Set j := k/2, sgn := (-1)**j. Return sgn*sin(r)
! 79: * where sin(r) is approximated by an odd polynomial in r
! 80: * r + r*s*(A1+s*(A2+ ... + s*A7)), s = r*r.
! 81: * Exit.
! 82: *
! 83: * 7. If |X| > 1, go to 9.
! 84: *
! 85: * 8. (|X|<2**(-40)) If SIN is invoked, return X; otherwise return 1.
! 86: *
! 87: * 9. Overwrite X by X := X rem 2Pi. Now that |X| <= Pi, go back to 3.
! 88: *
! 89: * SINCOS:
! 90: * 1. If |X| >= 15Pi or |X| < 2**(-40), go to 6.
! 91: *
! 92: * 2. Decompose X as X = N(Pi/2) + r where |r| <= Pi/4. Let
! 93: * k = N mod 4, so in particular, k = 0,1,2,or 3.
! 94: *
! 95: * 3. If k is even, go to 5.
! 96: *
! 97: * 4. (k is odd) Set j1 := (k-1)/2, j2 := j1 (EOR) (k mod 2), i.e.
! 98: * j1 exclusive or with the l.s.b. of k.
! 99: * sgn1 := (-1)**j1, sgn2 := (-1)**j2.
! 100: * SIN(X) = sgn1 * cos(r) and COS(X) = sgn2*sin(r) where
! 101: * sin(r) and cos(r) are computed as odd and even polynomials
! 102: * in r, respectively. Exit
! 103: *
! 104: * 5. (k is even) Set j1 := k/2, sgn1 := (-1)**j1.
! 105: * SIN(X) = sgn1 * sin(r) and COS(X) = sgn1*cos(r) where
! 106: * sin(r) and cos(r) are computed as odd and even polynomials
! 107: * in r, respectively. Exit
! 108: *
! 109: * 6. If |X| > 1, go to 8.
! 110: *
! 111: * 7. (|X|<2**(-40)) SIN(X) = X and COS(X) = 1. Exit.
! 112: *
! 113: * 8. Overwrite X by X := X rem 2Pi. Now that |X| <= Pi, go back to 2.
! 114: *
! 115:
! 116: SSIN IDNT 2,1 Motorola 040 Floating Point Software Package
! 117:
! 118: section 8
! 119:
! 120: include fpsp.h
! 121:
! 122: BOUNDS1 DC.L $3FD78000,$4004BC7E
! 123: TWOBYPI DC.L $3FE45F30,$6DC9C883
! 124:
! 125: SINA7 DC.L $BD6AAA77,$CCC994F5
! 126: SINA6 DC.L $3DE61209,$7AAE8DA1
! 127:
! 128: SINA5 DC.L $BE5AE645,$2A118AE4
! 129: SINA4 DC.L $3EC71DE3,$A5341531
! 130:
! 131: SINA3 DC.L $BF2A01A0,$1A018B59,$00000000,$00000000
! 132:
! 133: SINA2 DC.L $3FF80000,$88888888,$888859AF,$00000000
! 134:
! 135: SINA1 DC.L $BFFC0000,$AAAAAAAA,$AAAAAA99,$00000000
! 136:
! 137: COSB8 DC.L $3D2AC4D0,$D6011EE3
! 138: COSB7 DC.L $BDA9396F,$9F45AC19
! 139:
! 140: COSB6 DC.L $3E21EED9,$0612C972
! 141: COSB5 DC.L $BE927E4F,$B79D9FCF
! 142:
! 143: COSB4 DC.L $3EFA01A0,$1A01D423,$00000000,$00000000
! 144:
! 145: COSB3 DC.L $BFF50000,$B60B60B6,$0B61D438,$00000000
! 146:
! 147: COSB2 DC.L $3FFA0000,$AAAAAAAA,$AAAAAB5E
! 148: COSB1 DC.L $BF000000
! 149:
! 150: INVTWOPI DC.L $3FFC0000,$A2F9836E,$4E44152A
! 151:
! 152: TWOPI1 DC.L $40010000,$C90FDAA2,$00000000,$00000000
! 153: TWOPI2 DC.L $3FDF0000,$85A308D4,$00000000,$00000000
! 154:
! 155: xref PITBL
! 156:
! 157: INARG equ FP_SCR4
! 158:
! 159: X equ FP_SCR5
! 160: XDCARE equ X+2
! 161: XFRAC equ X+4
! 162:
! 163: RPRIME equ FP_SCR1
! 164: SPRIME equ FP_SCR2
! 165:
! 166: POSNEG1 equ L_SCR1
! 167: TWOTO63 equ L_SCR1
! 168:
! 169: ENDFLAG equ L_SCR2
! 170: N equ L_SCR2
! 171:
! 172: ADJN equ L_SCR3
! 173:
! 174: xref t_frcinx
! 175: xref t_extdnrm
! 176: xref sto_cos
! 177:
! 178: xdef ssind
! 179: ssind:
! 180: *--SIN(X) = X FOR DENORMALIZED X
! 181: bra t_extdnrm
! 182:
! 183: xdef scosd
! 184: scosd:
! 185: *--COS(X) = 1 FOR DENORMALIZED X
! 186:
! 187: FMOVE.S #:3F800000,FP0
! 188: *
! 189: * 9D25B Fix: Sometimes the previous fmove.s sets fpsr bits
! 190: *
! 191: fmove.l #0,fpsr
! 192: *
! 193: bra t_frcinx
! 194:
! 195: xdef ssin
! 196: ssin:
! 197: *--SET ADJN TO 0
! 198: CLR.L ADJN(a6)
! 199: BRA.B SINBGN
! 200:
! 201: xdef scos
! 202: scos:
! 203: *--SET ADJN TO 1
! 204: MOVE.L #1,ADJN(a6)
! 205:
! 206: SINBGN:
! 207: *--SAVE FPCR, FP1. CHECK IF |X| IS TOO SMALL OR LARGE
! 208:
! 209: FMOVE.X (a0),FP0 ...LOAD INPUT
! 210:
! 211: MOVE.L (A0),D0
! 212: MOVE.W 4(A0),D0
! 213: FMOVE.X FP0,X(a6)
! 214: ANDI.L #$7FFFFFFF,D0 ...COMPACTIFY X
! 215:
! 216: CMPI.L #$3FD78000,D0 ...|X| >= 2**(-40)?
! 217: BGE.B SOK1
! 218: BRA.W SINSM
! 219:
! 220: SOK1:
! 221: CMPI.L #$4004BC7E,D0 ...|X| < 15 PI?
! 222: BLT.B SINMAIN
! 223: BRA.W REDUCEX
! 224:
! 225: SINMAIN:
! 226: *--THIS IS THE USUAL CASE, |X| <= 15 PI.
! 227: *--THE ARGUMENT REDUCTION IS DONE BY TABLE LOOK UP.
! 228: FMOVE.X FP0,FP1
! 229: FMUL.D TWOBYPI,FP1 ...X*2/PI
! 230:
! 231: *--HIDE THE NEXT THREE INSTRUCTIONS
! 232: LEA PITBL+$200,A1 ...TABLE OF N*PI/2, N = -32,...,32
! 233:
! 234:
! 235: *--FP1 IS NOW READY
! 236: FMOVE.L FP1,N(a6) ...CONVERT TO INTEGER
! 237:
! 238: MOVE.L N(a6),D0
! 239: ASL.L #4,D0
! 240: ADDA.L D0,A1 ...A1 IS THE ADDRESS OF N*PIBY2
! 241: * ...WHICH IS IN TWO PIECES Y1 & Y2
! 242:
! 243: FSUB.X (A1)+,FP0 ...X-Y1
! 244: *--HIDE THE NEXT ONE
! 245: FSUB.S (A1),FP0 ...FP0 IS R = (X-Y1)-Y2
! 246:
! 247: SINCONT:
! 248: *--continuation from REDUCEX
! 249:
! 250: *--GET N+ADJN AND SEE IF SIN(R) OR COS(R) IS NEEDED
! 251: MOVE.L N(a6),D0
! 252: ADD.L ADJN(a6),D0 ...SEE IF D0 IS ODD OR EVEN
! 253: ROR.L #1,D0 ...D0 WAS ODD IFF D0 IS NEGATIVE
! 254: TST.L D0
! 255: BLT.W COSPOLY
! 256:
! 257: SINPOLY:
! 258: *--LET J BE THE LEAST SIG. BIT OF D0, LET SGN := (-1)**J.
! 259: *--THEN WE RETURN SGN*SIN(R). SGN*SIN(R) IS COMPUTED BY
! 260: *--R' + R'*S*(A1 + S(A2 + S(A3 + S(A4 + ... + SA7)))), WHERE
! 261: *--R' = SGN*R, S=R*R. THIS CAN BE REWRITTEN AS
! 262: *--R' + R'*S*( [A1+T(A3+T(A5+TA7))] + [S(A2+T(A4+TA6))])
! 263: *--WHERE T=S*S.
! 264: *--NOTE THAT A3 THROUGH A7 ARE STORED IN DOUBLE PRECISION
! 265: *--WHILE A1 AND A2 ARE IN DOUBLE-EXTENDED FORMAT.
! 266: FMOVE.X FP0,X(a6) ...X IS R
! 267: FMUL.X FP0,FP0 ...FP0 IS S
! 268: *---HIDE THE NEXT TWO WHILE WAITING FOR FP0
! 269: FMOVE.D SINA7,FP3
! 270: FMOVE.D SINA6,FP2
! 271: *--FP0 IS NOW READY
! 272: FMOVE.X FP0,FP1
! 273: FMUL.X FP1,FP1 ...FP1 IS T
! 274: *--HIDE THE NEXT TWO WHILE WAITING FOR FP1
! 275:
! 276: ROR.L #1,D0
! 277: ANDI.L #$80000000,D0
! 278: * ...LEAST SIG. BIT OF D0 IN SIGN POSITION
! 279: EOR.L D0,X(a6) ...X IS NOW R'= SGN*R
! 280:
! 281: FMUL.X FP1,FP3 ...TA7
! 282: FMUL.X FP1,FP2 ...TA6
! 283:
! 284: FADD.D SINA5,FP3 ...A5+TA7
! 285: FADD.D SINA4,FP2 ...A4+TA6
! 286:
! 287: FMUL.X FP1,FP3 ...T(A5+TA7)
! 288: FMUL.X FP1,FP2 ...T(A4+TA6)
! 289:
! 290: FADD.D SINA3,FP3 ...A3+T(A5+TA7)
! 291: FADD.X SINA2,FP2 ...A2+T(A4+TA6)
! 292:
! 293: FMUL.X FP3,FP1 ...T(A3+T(A5+TA7))
! 294:
! 295: FMUL.X FP0,FP2 ...S(A2+T(A4+TA6))
! 296: FADD.X SINA1,FP1 ...A1+T(A3+T(A5+TA7))
! 297: FMUL.X X(a6),FP0 ...R'*S
! 298:
! 299: FADD.X FP2,FP1 ...[A1+T(A3+T(A5+TA7))]+[S(A2+T(A4+TA6))]
! 300: *--FP3 RELEASED, RESTORE NOW AND TAKE SOME ADVANTAGE OF HIDING
! 301: *--FP2 RELEASED, RESTORE NOW AND TAKE FULL ADVANTAGE OF HIDING
! 302:
! 303:
! 304: FMUL.X FP1,FP0 ...SIN(R')-R'
! 305: *--FP1 RELEASED.
! 306:
! 307: FMOVE.L d1,FPCR ;restore users exceptions
! 308: FADD.X X(a6),FP0 ;last inst - possible exception set
! 309: bra t_frcinx
! 310:
! 311:
! 312: COSPOLY:
! 313: *--LET J BE THE LEAST SIG. BIT OF D0, LET SGN := (-1)**J.
! 314: *--THEN WE RETURN SGN*COS(R). SGN*COS(R) IS COMPUTED BY
! 315: *--SGN + S'*(B1 + S(B2 + S(B3 + S(B4 + ... + SB8)))), WHERE
! 316: *--S=R*R AND S'=SGN*S. THIS CAN BE REWRITTEN AS
! 317: *--SGN + S'*([B1+T(B3+T(B5+TB7))] + [S(B2+T(B4+T(B6+TB8)))])
! 318: *--WHERE T=S*S.
! 319: *--NOTE THAT B4 THROUGH B8 ARE STORED IN DOUBLE PRECISION
! 320: *--WHILE B2 AND B3 ARE IN DOUBLE-EXTENDED FORMAT, B1 IS -1/2
! 321: *--AND IS THEREFORE STORED AS SINGLE PRECISION.
! 322:
! 323: FMUL.X FP0,FP0 ...FP0 IS S
! 324: *---HIDE THE NEXT TWO WHILE WAITING FOR FP0
! 325: FMOVE.D COSB8,FP2
! 326: FMOVE.D COSB7,FP3
! 327: *--FP0 IS NOW READY
! 328: FMOVE.X FP0,FP1
! 329: FMUL.X FP1,FP1 ...FP1 IS T
! 330: *--HIDE THE NEXT TWO WHILE WAITING FOR FP1
! 331: FMOVE.X FP0,X(a6) ...X IS S
! 332: ROR.L #1,D0
! 333: ANDI.L #$80000000,D0
! 334: * ...LEAST SIG. BIT OF D0 IN SIGN POSITION
! 335:
! 336: FMUL.X FP1,FP2 ...TB8
! 337: *--HIDE THE NEXT TWO WHILE WAITING FOR THE XU
! 338: EOR.L D0,X(a6) ...X IS NOW S'= SGN*S
! 339: ANDI.L #$80000000,D0
! 340:
! 341: FMUL.X FP1,FP3 ...TB7
! 342: *--HIDE THE NEXT TWO WHILE WAITING FOR THE XU
! 343: ORI.L #$3F800000,D0 ...D0 IS SGN IN SINGLE
! 344: MOVE.L D0,POSNEG1(a6)
! 345:
! 346: FADD.D COSB6,FP2 ...B6+TB8
! 347: FADD.D COSB5,FP3 ...B5+TB7
! 348:
! 349: FMUL.X FP1,FP2 ...T(B6+TB8)
! 350: FMUL.X FP1,FP3 ...T(B5+TB7)
! 351:
! 352: FADD.D COSB4,FP2 ...B4+T(B6+TB8)
! 353: FADD.X COSB3,FP3 ...B3+T(B5+TB7)
! 354:
! 355: FMUL.X FP1,FP2 ...T(B4+T(B6+TB8))
! 356: FMUL.X FP3,FP1 ...T(B3+T(B5+TB7))
! 357:
! 358: FADD.X COSB2,FP2 ...B2+T(B4+T(B6+TB8))
! 359: FADD.S COSB1,FP1 ...B1+T(B3+T(B5+TB7))
! 360:
! 361: FMUL.X FP2,FP0 ...S(B2+T(B4+T(B6+TB8)))
! 362: *--FP3 RELEASED, RESTORE NOW AND TAKE SOME ADVANTAGE OF HIDING
! 363: *--FP2 RELEASED.
! 364:
! 365:
! 366: FADD.X FP1,FP0
! 367: *--FP1 RELEASED
! 368:
! 369: FMUL.X X(a6),FP0
! 370:
! 371: FMOVE.L d1,FPCR ;restore users exceptions
! 372: FADD.S POSNEG1(a6),FP0 ;last inst - possible exception set
! 373: bra t_frcinx
! 374:
! 375:
! 376: SINBORS:
! 377: *--IF |X| > 15PI, WE USE THE GENERAL ARGUMENT REDUCTION.
! 378: *--IF |X| < 2**(-40), RETURN X OR 1.
! 379: CMPI.L #$3FFF8000,D0
! 380: BGT.B REDUCEX
! 381:
! 382:
! 383: SINSM:
! 384: MOVE.L ADJN(a6),D0
! 385: TST.L D0
! 386: BGT.B COSTINY
! 387:
! 388: SINTINY:
! 389: CLR.W XDCARE(a6) ...JUST IN CASE
! 390: FMOVE.L d1,FPCR ;restore users exceptions
! 391: FMOVE.X X(a6),FP0 ;last inst - possible exception set
! 392: bra t_frcinx
! 393:
! 394:
! 395: COSTINY:
! 396: FMOVE.S #:3F800000,FP0
! 397:
! 398: FMOVE.L d1,FPCR ;restore users exceptions
! 399: FSUB.S #:00800000,FP0 ;last inst - possible exception set
! 400: bra t_frcinx
! 401:
! 402:
! 403: REDUCEX:
! 404: *--WHEN REDUCEX IS USED, THE CODE WILL INEVITABLY BE SLOW.
! 405: *--THIS REDUCTION METHOD, HOWEVER, IS MUCH FASTER THAN USING
! 406: *--THE REMAINDER INSTRUCTION WHICH IS NOW IN SOFTWARE.
! 407:
! 408: FMOVEM.X FP2-FP5,-(A7) ...save FP2 through FP5
! 409: MOVE.L D2,-(A7)
! 410: FMOVE.S #:00000000,FP1
! 411: *--If compact form of abs(arg) in d0=$7ffeffff, argument is so large that
! 412: *--there is a danger of unwanted overflow in first LOOP iteration. In this
! 413: *--case, reduce argument by one remainder step to make subsequent reduction
! 414: *--safe.
! 415: cmpi.l #$7ffeffff,d0 ;is argument dangerously large?
! 416: bne.b LOOP
! 417: move.l #$7ffe0000,FP_SCR2(a6) ;yes
! 418: * ;create 2**16383*PI/2
! 419: move.l #$c90fdaa2,FP_SCR2+4(a6)
! 420: clr.l FP_SCR2+8(a6)
! 421: ftst.x fp0 ;test sign of argument
! 422: move.l #$7fdc0000,FP_SCR3(a6) ;create low half of 2**16383*
! 423: * ;PI/2 at FP_SCR3
! 424: move.l #$85a308d3,FP_SCR3+4(a6)
! 425: clr.l FP_SCR3+8(a6)
! 426: fblt.w red_neg
! 427: or.w #$8000,FP_SCR2(a6) ;positive arg
! 428: or.w #$8000,FP_SCR3(a6)
! 429: red_neg:
! 430: fadd.x FP_SCR2(a6),fp0 ;high part of reduction is exact
! 431: fmove.x fp0,fp1 ;save high result in fp1
! 432: fadd.x FP_SCR3(a6),fp0 ;low part of reduction
! 433: fsub.x fp0,fp1 ;determine low component of result
! 434: fadd.x FP_SCR3(a6),fp1 ;fp0/fp1 are reduced argument.
! 435:
! 436: *--ON ENTRY, FP0 IS X, ON RETURN, FP0 IS X REM PI/2, |X| <= PI/4.
! 437: *--integer quotient will be stored in N
! 438: *--Intermeditate remainder is 66-bit long; (R,r) in (FP0,FP1)
! 439:
! 440: LOOP:
! 441: FMOVE.X FP0,INARG(a6) ...+-2**K * F, 1 <= F < 2
! 442: MOVE.W INARG(a6),D0
! 443: MOVE.L D0,A1 ...save a copy of D0
! 444: ANDI.L #$00007FFF,D0
! 445: SUBI.L #$00003FFF,D0 ...D0 IS K
! 446: CMPI.L #28,D0
! 447: BLE.B LASTLOOP
! 448: CONTLOOP:
! 449: SUBI.L #27,D0 ...D0 IS L := K-27
! 450: CLR.L ENDFLAG(a6)
! 451: BRA.B WORK
! 452: LASTLOOP:
! 453: CLR.L D0 ...D0 IS L := 0
! 454: MOVE.L #1,ENDFLAG(a6)
! 455:
! 456: WORK:
! 457: *--FIND THE REMAINDER OF (R,r) W.R.T. 2**L * (PI/2). L IS SO CHOSEN
! 458: *--THAT INT( X * (2/PI) / 2**(L) ) < 2**29.
! 459:
! 460: *--CREATE 2**(-L) * (2/PI), SIGN(INARG)*2**(63),
! 461: *--2**L * (PIby2_1), 2**L * (PIby2_2)
! 462:
! 463: MOVE.L #$00003FFE,D2 ...BIASED EXPO OF 2/PI
! 464: SUB.L D0,D2 ...BIASED EXPO OF 2**(-L)*(2/PI)
! 465:
! 466: MOVE.L #$A2F9836E,FP_SCR1+4(a6)
! 467: MOVE.L #$4E44152A,FP_SCR1+8(a6)
! 468: MOVE.W D2,FP_SCR1(a6) ...FP_SCR1 is 2**(-L)*(2/PI)
! 469:
! 470: FMOVE.X FP0,FP2
! 471: FMUL.X FP_SCR1(a6),FP2
! 472: *--WE MUST NOW FIND INT(FP2). SINCE WE NEED THIS VALUE IN
! 473: *--FLOATING POINT FORMAT, THE TWO FMOVE'S FMOVE.L FP <--> N
! 474: *--WILL BE TOO INEFFICIENT. THE WAY AROUND IT IS THAT
! 475: *--(SIGN(INARG)*2**63 + FP2) - SIGN(INARG)*2**63 WILL GIVE
! 476: *--US THE DESIRED VALUE IN FLOATING POINT.
! 477:
! 478: *--HIDE SIX CYCLES OF INSTRUCTION
! 479: MOVE.L A1,D2
! 480: SWAP D2
! 481: ANDI.L #$80000000,D2
! 482: ORI.L #$5F000000,D2 ...D2 IS SIGN(INARG)*2**63 IN SGL
! 483: MOVE.L D2,TWOTO63(a6)
! 484:
! 485: MOVE.L D0,D2
! 486: ADDI.L #$00003FFF,D2 ...BIASED EXPO OF 2**L * (PI/2)
! 487:
! 488: *--FP2 IS READY
! 489: FADD.S TWOTO63(a6),FP2 ...THE FRACTIONAL PART OF FP1 IS ROUNDED
! 490:
! 491: *--HIDE 4 CYCLES OF INSTRUCTION; creating 2**(L)*Piby2_1 and 2**(L)*Piby2_2
! 492: MOVE.W D2,FP_SCR2(a6)
! 493: CLR.W FP_SCR2+2(a6)
! 494: MOVE.L #$C90FDAA2,FP_SCR2+4(a6)
! 495: CLR.L FP_SCR2+8(a6) ...FP_SCR2 is 2**(L) * Piby2_1
! 496:
! 497: *--FP2 IS READY
! 498: FSUB.S TWOTO63(a6),FP2 ...FP2 is N
! 499:
! 500: ADDI.L #$00003FDD,D0
! 501: MOVE.W D0,FP_SCR3(a6)
! 502: CLR.W FP_SCR3+2(a6)
! 503: MOVE.L #$85A308D3,FP_SCR3+4(a6)
! 504: CLR.L FP_SCR3+8(a6) ...FP_SCR3 is 2**(L) * Piby2_2
! 505:
! 506: MOVE.L ENDFLAG(a6),D0
! 507:
! 508: *--We are now ready to perform (R+r) - N*P1 - N*P2, P1 = 2**(L) * Piby2_1 and
! 509: *--P2 = 2**(L) * Piby2_2
! 510: FMOVE.X FP2,FP4
! 511: FMul.X FP_SCR2(a6),FP4 ...W = N*P1
! 512: FMove.X FP2,FP5
! 513: FMul.X FP_SCR3(a6),FP5 ...w = N*P2
! 514: FMove.X FP4,FP3
! 515: *--we want P+p = W+w but |p| <= half ulp of P
! 516: *--Then, we need to compute A := R-P and a := r-p
! 517: FAdd.X FP5,FP3 ...FP3 is P
! 518: FSub.X FP3,FP4 ...W-P
! 519:
! 520: FSub.X FP3,FP0 ...FP0 is A := R - P
! 521: FAdd.X FP5,FP4 ...FP4 is p = (W-P)+w
! 522:
! 523: FMove.X FP0,FP3 ...FP3 A
! 524: FSub.X FP4,FP1 ...FP1 is a := r - p
! 525:
! 526: *--Now we need to normalize (A,a) to "new (R,r)" where R+r = A+a but
! 527: *--|r| <= half ulp of R.
! 528: FAdd.X FP1,FP0 ...FP0 is R := A+a
! 529: *--No need to calculate r if this is the last loop
! 530: TST.L D0
! 531: BGT.W RESTORE
! 532:
! 533: *--Need to calculate r
! 534: FSub.X FP0,FP3 ...A-R
! 535: FAdd.X FP3,FP1 ...FP1 is r := (A-R)+a
! 536: BRA.W LOOP
! 537:
! 538: RESTORE:
! 539: FMOVE.L FP2,N(a6)
! 540: MOVE.L (A7)+,D2
! 541: FMOVEM.X (A7)+,FP2-FP5
! 542:
! 543:
! 544: MOVE.L ADJN(a6),D0
! 545: CMPI.L #4,D0
! 546:
! 547: BLT.W SINCONT
! 548: BRA.B SCCONT
! 549:
! 550: xdef ssincosd
! 551: ssincosd:
! 552: *--SIN AND COS OF X FOR DENORMALIZED X
! 553:
! 554: FMOVE.S #:3F800000,FP1
! 555: bsr sto_cos ;store cosine result
! 556: bra t_extdnrm
! 557:
! 558: xdef ssincos
! 559: ssincos:
! 560: *--SET ADJN TO 4
! 561: MOVE.L #4,ADJN(a6)
! 562:
! 563: FMOVE.X (a0),FP0 ...LOAD INPUT
! 564:
! 565: MOVE.L (A0),D0
! 566: MOVE.W 4(A0),D0
! 567: FMOVE.X FP0,X(a6)
! 568: ANDI.L #$7FFFFFFF,D0 ...COMPACTIFY X
! 569:
! 570: CMPI.L #$3FD78000,D0 ...|X| >= 2**(-40)?
! 571: BGE.B SCOK1
! 572: BRA.W SCSM
! 573:
! 574: SCOK1:
! 575: CMPI.L #$4004BC7E,D0 ...|X| < 15 PI?
! 576: BLT.B SCMAIN
! 577: BRA.W REDUCEX
! 578:
! 579:
! 580: SCMAIN:
! 581: *--THIS IS THE USUAL CASE, |X| <= 15 PI.
! 582: *--THE ARGUMENT REDUCTION IS DONE BY TABLE LOOK UP.
! 583: FMOVE.X FP0,FP1
! 584: FMUL.D TWOBYPI,FP1 ...X*2/PI
! 585:
! 586: *--HIDE THE NEXT THREE INSTRUCTIONS
! 587: LEA PITBL+$200,A1 ...TABLE OF N*PI/2, N = -32,...,32
! 588:
! 589:
! 590: *--FP1 IS NOW READY
! 591: FMOVE.L FP1,N(a6) ...CONVERT TO INTEGER
! 592:
! 593: MOVE.L N(a6),D0
! 594: ASL.L #4,D0
! 595: ADDA.L D0,A1 ...ADDRESS OF N*PIBY2, IN Y1, Y2
! 596:
! 597: FSUB.X (A1)+,FP0 ...X-Y1
! 598: FSUB.S (A1),FP0 ...FP0 IS R = (X-Y1)-Y2
! 599:
! 600: SCCONT:
! 601: *--continuation point from REDUCEX
! 602:
! 603: *--HIDE THE NEXT TWO
! 604: MOVE.L N(a6),D0
! 605: ROR.L #1,D0
! 606:
! 607: TST.L D0 ...D0 < 0 IFF N IS ODD
! 608: BGE.W NEVEN
! 609:
! 610: NODD:
! 611: *--REGISTERS SAVED SO FAR: D0, A0, FP2.
! 612:
! 613: FMOVE.X FP0,RPRIME(a6)
! 614: FMUL.X FP0,FP0 ...FP0 IS S = R*R
! 615: FMOVE.D SINA7,FP1 ...A7
! 616: FMOVE.D COSB8,FP2 ...B8
! 617: FMUL.X FP0,FP1 ...SA7
! 618: MOVE.L d2,-(A7)
! 619: MOVE.L D0,d2
! 620: FMUL.X FP0,FP2 ...SB8
! 621: ROR.L #1,d2
! 622: ANDI.L #$80000000,d2
! 623:
! 624: FADD.D SINA6,FP1 ...A6+SA7
! 625: EOR.L D0,d2
! 626: ANDI.L #$80000000,d2
! 627: FADD.D COSB7,FP2 ...B7+SB8
! 628:
! 629: FMUL.X FP0,FP1 ...S(A6+SA7)
! 630: EOR.L d2,RPRIME(a6)
! 631: MOVE.L (A7)+,d2
! 632: FMUL.X FP0,FP2 ...S(B7+SB8)
! 633: ROR.L #1,D0
! 634: ANDI.L #$80000000,D0
! 635:
! 636: FADD.D SINA5,FP1 ...A5+S(A6+SA7)
! 637: MOVE.L #$3F800000,POSNEG1(a6)
! 638: EOR.L D0,POSNEG1(a6)
! 639: FADD.D COSB6,FP2 ...B6+S(B7+SB8)
! 640:
! 641: FMUL.X FP0,FP1 ...S(A5+S(A6+SA7))
! 642: FMUL.X FP0,FP2 ...S(B6+S(B7+SB8))
! 643: FMOVE.X FP0,SPRIME(a6)
! 644:
! 645: FADD.D SINA4,FP1 ...A4+S(A5+S(A6+SA7))
! 646: EOR.L D0,SPRIME(a6)
! 647: FADD.D COSB5,FP2 ...B5+S(B6+S(B7+SB8))
! 648:
! 649: FMUL.X FP0,FP1 ...S(A4+...)
! 650: FMUL.X FP0,FP2 ...S(B5+...)
! 651:
! 652: FADD.D SINA3,FP1 ...A3+S(A4+...)
! 653: FADD.D COSB4,FP2 ...B4+S(B5+...)
! 654:
! 655: FMUL.X FP0,FP1 ...S(A3+...)
! 656: FMUL.X FP0,FP2 ...S(B4+...)
! 657:
! 658: FADD.X SINA2,FP1 ...A2+S(A3+...)
! 659: FADD.X COSB3,FP2 ...B3+S(B4+...)
! 660:
! 661: FMUL.X FP0,FP1 ...S(A2+...)
! 662: FMUL.X FP0,FP2 ...S(B3+...)
! 663:
! 664: FADD.X SINA1,FP1 ...A1+S(A2+...)
! 665: FADD.X COSB2,FP2 ...B2+S(B3+...)
! 666:
! 667: FMUL.X FP0,FP1 ...S(A1+...)
! 668: FMUL.X FP2,FP0 ...S(B2+...)
! 669:
! 670:
! 671:
! 672: FMUL.X RPRIME(a6),FP1 ...R'S(A1+...)
! 673: FADD.S COSB1,FP0 ...B1+S(B2...)
! 674: FMUL.X SPRIME(a6),FP0 ...S'(B1+S(B2+...))
! 675:
! 676: move.l d1,-(sp) ;restore users mode & precision
! 677: andi.l #$ff,d1 ;mask off all exceptions
! 678: fmove.l d1,FPCR
! 679: FADD.X RPRIME(a6),FP1 ...COS(X)
! 680: bsr sto_cos ;store cosine result
! 681: FMOVE.L (sp)+,FPCR ;restore users exceptions
! 682: FADD.S POSNEG1(a6),FP0 ...SIN(X)
! 683:
! 684: bra t_frcinx
! 685:
! 686:
! 687: NEVEN:
! 688: *--REGISTERS SAVED SO FAR: FP2.
! 689:
! 690: FMOVE.X FP0,RPRIME(a6)
! 691: FMUL.X FP0,FP0 ...FP0 IS S = R*R
! 692: FMOVE.D COSB8,FP1 ...B8
! 693: FMOVE.D SINA7,FP2 ...A7
! 694: FMUL.X FP0,FP1 ...SB8
! 695: FMOVE.X FP0,SPRIME(a6)
! 696: FMUL.X FP0,FP2 ...SA7
! 697: ROR.L #1,D0
! 698: ANDI.L #$80000000,D0
! 699: FADD.D COSB7,FP1 ...B7+SB8
! 700: FADD.D SINA6,FP2 ...A6+SA7
! 701: EOR.L D0,RPRIME(a6)
! 702: EOR.L D0,SPRIME(a6)
! 703: FMUL.X FP0,FP1 ...S(B7+SB8)
! 704: ORI.L #$3F800000,D0
! 705: MOVE.L D0,POSNEG1(a6)
! 706: FMUL.X FP0,FP2 ...S(A6+SA7)
! 707:
! 708: FADD.D COSB6,FP1 ...B6+S(B7+SB8)
! 709: FADD.D SINA5,FP2 ...A5+S(A6+SA7)
! 710:
! 711: FMUL.X FP0,FP1 ...S(B6+S(B7+SB8))
! 712: FMUL.X FP0,FP2 ...S(A5+S(A6+SA7))
! 713:
! 714: FADD.D COSB5,FP1 ...B5+S(B6+S(B7+SB8))
! 715: FADD.D SINA4,FP2 ...A4+S(A5+S(A6+SA7))
! 716:
! 717: FMUL.X FP0,FP1 ...S(B5+...)
! 718: FMUL.X FP0,FP2 ...S(A4+...)
! 719:
! 720: FADD.D COSB4,FP1 ...B4+S(B5+...)
! 721: FADD.D SINA3,FP2 ...A3+S(A4+...)
! 722:
! 723: FMUL.X FP0,FP1 ...S(B4+...)
! 724: FMUL.X FP0,FP2 ...S(A3+...)
! 725:
! 726: FADD.X COSB3,FP1 ...B3+S(B4+...)
! 727: FADD.X SINA2,FP2 ...A2+S(A3+...)
! 728:
! 729: FMUL.X FP0,FP1 ...S(B3+...)
! 730: FMUL.X FP0,FP2 ...S(A2+...)
! 731:
! 732: FADD.X COSB2,FP1 ...B2+S(B3+...)
! 733: FADD.X SINA1,FP2 ...A1+S(A2+...)
! 734:
! 735: FMUL.X FP0,FP1 ...S(B2+...)
! 736: fmul.x fp2,fp0 ...s(a1+...)
! 737:
! 738:
! 739:
! 740: FADD.S COSB1,FP1 ...B1+S(B2...)
! 741: FMUL.X RPRIME(a6),FP0 ...R'S(A1+...)
! 742: FMUL.X SPRIME(a6),FP1 ...S'(B1+S(B2+...))
! 743:
! 744: move.l d1,-(sp) ;save users mode & precision
! 745: andi.l #$ff,d1 ;mask off all exceptions
! 746: fmove.l d1,FPCR
! 747: FADD.S POSNEG1(a6),FP1 ...COS(X)
! 748: bsr sto_cos ;store cosine result
! 749: FMOVE.L (sp)+,FPCR ;restore users exceptions
! 750: FADD.X RPRIME(a6),FP0 ...SIN(X)
! 751:
! 752: bra t_frcinx
! 753:
! 754: SCBORS:
! 755: CMPI.L #$3FFF8000,D0
! 756: BGT.W REDUCEX
! 757:
! 758:
! 759: SCSM:
! 760: CLR.W XDCARE(a6)
! 761: FMOVE.S #:3F800000,FP1
! 762:
! 763: move.l d1,-(sp) ;save users mode & precision
! 764: andi.l #$ff,d1 ;mask off all exceptions
! 765: fmove.l d1,FPCR
! 766: FSUB.S #:00800000,FP1
! 767: bsr sto_cos ;store cosine result
! 768: FMOVE.L (sp)+,FPCR ;restore users exceptions
! 769: FMOVE.X X(a6),FP0
! 770: bra t_frcinx
! 771:
! 772: end
CVSweb
![[BACK]](/icons/back.gif){kind=icon}
Return to ssin.sa CVS log ![[TXT]](/icons/text.gif){kind=icon}
![[DIR]](/icons/dir.gif){kind=icon}
Up to [local] / sys / arch / m68k / fpsp