* * $VER: atan.s 33.1 (22.1.97) * * Calculates the arctangent of the source * * Version history: * * 33.1 22.1.97 (c) Motorola * * - snipped from M68060SP sources * machine 68040 fpu 1 XDEF _atan XDEF @atan ************************************************************************* * atan(): computes the arctangent of a normalized number * * * * INPUT *************************************************************** * * fp0 = extended precision input * * * * OUTPUT ************************************************************** * * fp0 = arctan(X) * * * * ACCURACY and MONOTONICITY ******************************************* * * The returned result is within 2 ulps in 64 significant bit, * * i.e. within 0.5001 ulp to 53 bits if the result is subsequently * * rounded to double precision. The result is provably monotonic * * in double precision. * * * * ALGORITHM *********************************************************** * * Step 1. If |X| >= 16 or |X| < 1/16, go to Step 5. * * * * Step 2. Let X = sgn * 2**k * 1.xxxxxxxx...x. * * Note that k = -4, -3,..., or 3. * * Define F = sgn * 2**k * 1.xxxx1, i.e. the first 5 * * significant bits of X with a bit-1 attached at the 6-th * * bit position. Define u to be u = (X-F) / (1 + X*F). * * * * Step 3. Approximate arctan(u) by a polynomial poly. * * * * Step 4. Return arctan(F) + poly, arctan(F) is fetched from a * * table of values calculated beforehand. Exit. * * * * Step 5. If |X| >= 16, go to Step 7. * * * * Step 6. Approximate arctan(X) by an odd polynomial in X. Exit. * * * * Step 7. Define X' = -1/X. Approximate arctan(X') by an odd * * polynomial in X'. * * Arctan(X) = sign(X)*Pi/2 + arctan(X'). Exit. * * * ************************************************************************* ATANA3 dc.l $BFF6687E,$314987D8 ATANA2 dc.l $4002AC69,$34A26DB3 ATANA1 dc.l $BFC2476F,$4E1DA28E ATANB6 dc.l $3FB34444,$7F876989 ATANB5 dc.l $BFB744EE,$7FAF45DB ATANB4 dc.l $3FBC71C6,$46940220 ATANB3 dc.l $BFC24924,$921872F9 ATANB2 dc.l $3FC99999,$99998FA9 ATANB1 dc.l $BFD55555,$55555555 ATANC5 dc.l $BFB70BF3,$98539E6A ATANC4 dc.l $3FBC7187,$962D1D7D ATANC3 dc.l $BFC24924,$827107B8 ATANC2 dc.l $3FC99999,$9996263E ATANC1 dc.l $BFD55555,$55555536 PPIBY2 dc.l $3FFF0000,$C90FDAA2,$2168C235,$00000000 NPIBY2 dc.l $BFFF0000,$C90FDAA2,$2168C235,$00000000 PTINY dc.l $00010000,$80000000,$00000000,$00000000 NTINY dc.l $80010000,$80000000,$00000000,$00000000 ATANTBL dc.l $3FFB0000,$83D152C5,$060B7A51,$00000000 dc.l $3FFB0000,$8BC85445,$65498B8B,$00000000 dc.l $3FFB0000,$93BE4060,$17626B0D,$00000000 dc.l $3FFB0000,$9BB3078D,$35AEC202,$00000000 dc.l $3FFB0000,$A3A69A52,$5DDCE7DE,$00000000 dc.l $3FFB0000,$AB98E943,$62765619,$00000000 dc.l $3FFB0000,$B389E502,$F9C59862,$00000000 dc.l $3FFB0000,$BB797E43,$6B09E6FB,$00000000 dc.l $3FFB0000,$C367A5C7,$39E5F446,$00000000 dc.l $3FFB0000,$CB544C61,$CFF7D5C6,$00000000 dc.l $3FFB0000,$D33F62F8,$2488533E,$00000000 dc.l $3FFB0000,$DB28DA81,$62404C77,$00000000 dc.l $3FFB0000,$E310A407,$8AD34F18,$00000000 dc.l $3FFB0000,$EAF6B0A8,$188EE1EB,$00000000 dc.l $3FFB0000,$F2DAF194,$9DBE79D5,$00000000 dc.l $3FFB0000,$FABD5813,$61D47E3E,$00000000 dc.l $3FFC0000,$8346AC21,$0959ECC4,$00000000 dc.l $3FFC0000,$8B232A08,$304282D8,$00000000 dc.l $3FFC0000,$92FB70B8,$D29AE2F9,$00000000 dc.l $3FFC0000,$9ACF476F,$5CCD1CB4,$00000000 dc.l $3FFC0000,$A29E7630,$4954F23F,$00000000 dc.l $3FFC0000,$AA68C5D0,$8AB85230,$00000000 dc.l $3FFC0000,$B22DFFFD,$9D539F83,$00000000 dc.l $3FFC0000,$B9EDEF45,$3E900EA5,$00000000 dc.l $3FFC0000,$C1A85F1C,$C75E3EA5,$00000000 dc.l $3FFC0000,$C95D1BE8,$28138DE6,$00000000 dc.l $3FFC0000,$D10BF300,$840D2DE4,$00000000 dc.l $3FFC0000,$D8B4B2BA,$6BC05E7A,$00000000 dc.l $3FFC0000,$E0572A6B,$B42335F6,$00000000 dc.l $3FFC0000,$E7F32A70,$EA9CAA8F,$00000000 dc.l $3FFC0000,$EF888432,$64ECEFAA,$00000000 dc.l $3FFC0000,$F7170A28,$ECC06666,$00000000 dc.l $3FFD0000,$812FD288,$332DAD32,$00000000 dc.l $3FFD0000,$88A8D1B1,$218E4D64,$00000000 dc.l $3FFD0000,$9012AB3F,$23E4AEE8,$00000000 dc.l $3FFD0000,$976CC3D4,$11E7F1B9,$00000000 dc.l $3FFD0000,$9EB68949,$3889A227,$00000000 dc.l $3FFD0000,$A5EF72C3,$4487361B,$00000000 dc.l $3FFD0000,$AD1700BA,$F07A7227,$00000000 dc.l $3FFD0000,$B42CBCFA,$FD37EFB7,$00000000 dc.l $3FFD0000,$BB303A94,$0BA80F89,$00000000 dc.l $3FFD0000,$C22115C6,$FCAEBBAF,$00000000 dc.l $3FFD0000,$C8FEF3E6,$86331221,$00000000 dc.l $3FFD0000,$CFC98330,$B4000C70,$00000000 dc.l $3FFD0000,$D6807AA1,$102C5BF9,$00000000 dc.l $3FFD0000,$DD2399BC,$31252AA3,$00000000 dc.l $3FFD0000,$E3B2A855,$6B8FC517,$00000000 dc.l $3FFD0000,$EA2D764F,$64315989,$00000000 dc.l $3FFD0000,$F3BF5BF8,$BAD1A21D,$00000000 dc.l $3FFE0000,$801CE39E,$0D205C9A,$00000000 dc.l $3FFE0000,$8630A2DA,$DA1ED066,$00000000 dc.l $3FFE0000,$8C1AD445,$F3E09B8C,$00000000 dc.l $3FFE0000,$91DB8F16,$64F350E2,$00000000 dc.l $3FFE0000,$97731420,$365E538C,$00000000 dc.l $3FFE0000,$9CE1C8E6,$A0B8CDBA,$00000000 dc.l $3FFE0000,$A22832DB,$CADAAE09,$00000000 dc.l $3FFE0000,$A746F2DD,$B7602294,$00000000 dc.l $3FFE0000,$AC3EC0FB,$997DD6A2,$00000000 dc.l $3FFE0000,$B110688A,$EBDC6F6A,$00000000 dc.l $3FFE0000,$B5BCC490,$59ECC4B0,$00000000 dc.l $3FFE0000,$BA44BC7D,$D470782F,$00000000 dc.l $3FFE0000,$BEA94144,$FD049AAC,$00000000 dc.l $3FFE0000,$C2EB4ABB,$661628B6,$00000000 dc.l $3FFE0000,$C70BD54C,$E602EE14,$00000000 dc.l $3FFE0000,$CD000549,$ADEC7159,$00000000 dc.l $3FFE0000,$D48457D2,$D8EA4EA3,$00000000 dc.l $3FFE0000,$DB948DA7,$12DECE3B,$00000000 dc.l $3FFE0000,$E23855F9,$69E8096A,$00000000 dc.l $3FFE0000,$E8771129,$C4353259,$00000000 dc.l $3FFE0000,$EE57C16E,$0D379C0D,$00000000 dc.l $3FFE0000,$F3E10211,$A87C3779,$00000000 dc.l $3FFE0000,$F919039D,$758B8D41,$00000000 dc.l $3FFE0000,$FE058B8F,$64935FB3,$00000000 dc.l $3FFF0000,$8155FB49,$7B685D04,$00000000 dc.l $3FFF0000,$83889E35,$49D108E1,$00000000 dc.l $3FFF0000,$859CFA76,$511D724B,$00000000 dc.l $3FFF0000,$87952ECF,$FF8131E7,$00000000 dc.l $3FFF0000,$89732FD1,$9557641B,$00000000 dc.l $3FFF0000,$8B38CAD1,$01932A35,$00000000 dc.l $3FFF0000,$8CE7A8D8,$301EE6B5,$00000000 dc.l $3FFF0000,$8F46A39E,$2EAE5281,$00000000 dc.l $3FFF0000,$922DA7D7,$91888487,$00000000 dc.l $3FFF0000,$94D19FCB,$DEDF5241,$00000000 dc.l $3FFF0000,$973AB944,$19D2A08B,$00000000 dc.l $3FFF0000,$996FF00E,$08E10B96,$00000000 dc.l $3FFF0000,$9B773F95,$12321DA7,$00000000 dc.l $3FFF0000,$9D55CC32,$0F935624,$00000000 dc.l $3FFF0000,$9F100575,$006CC571,$00000000 dc.l $3FFF0000,$A0A9C290,$D97CC06C,$00000000 dc.l $3FFF0000,$A22659EB,$EBC0630A,$00000000 dc.l $3FFF0000,$A388B4AF,$F6EF0EC9,$00000000 dc.l $3FFF0000,$A4D35F10,$61D292C4,$00000000 dc.l $3FFF0000,$A60895DC,$FBE3187E,$00000000 dc.l $3FFF0000,$A72A51DC,$7367BEAC,$00000000 dc.l $3FFF0000,$A83A5153,$0956168F,$00000000 dc.l $3FFF0000,$A93A2007,$7539546E,$00000000 dc.l $3FFF0000,$AA9E7245,$023B2605,$00000000 dc.l $3FFF0000,$AC4C84BA,$6FE4D58F,$00000000 dc.l $3FFF0000,$ADCE4A4A,$606B9712,$00000000 dc.l $3FFF0000,$AF2A2DCD,$8D263C9C,$00000000 dc.l $3FFF0000,$B0656F81,$F22265C7,$00000000 dc.l $3FFF0000,$B1846515,$0F71496A,$00000000 dc.l $3FFF0000,$B28AAA15,$6F9ADA35,$00000000 dc.l $3FFF0000,$B37B44FF,$3766B895,$00000000 dc.l $3FFF0000,$B458C3DC,$E9630433,$00000000 dc.l $3FFF0000,$B525529D,$562246BD,$00000000 dc.l $3FFF0000,$B5E2CCA9,$5F9D88CC,$00000000 dc.l $3FFF0000,$B692CADA,$7ACA1ADA,$00000000 dc.l $3FFF0000,$B736AEA7,$A6925838,$00000000 dc.l $3FFF0000,$B7CFAB28,$7E9F7B36,$00000000 dc.l $3FFF0000,$B85ECC66,$CB219835,$00000000 dc.l $3FFF0000,$B8E4FD5A,$20A593DA,$00000000 dc.l $3FFF0000,$B99F41F6,$4AFF9BB5,$00000000 dc.l $3FFF0000,$BA7F1E17,$842BBE7B,$00000000 dc.l $3FFF0000,$BB471285,$7637E17D,$00000000 dc.l $3FFF0000,$BBFABE8A,$4788DF6F,$00000000 dc.l $3FFF0000,$BC9D0FAD,$2B689D79,$00000000 dc.l $3FFF0000,$BD306A39,$471ECD86,$00000000 dc.l $3FFF0000,$BDB6C731,$856AF18A,$00000000 dc.l $3FFF0000,$BE31CAC5,$02E80D70,$00000000 dc.l $3FFF0000,$BEA2D55C,$E33194E2,$00000000 dc.l $3FFF0000,$BF0B10B7,$C03128F0,$00000000 dc.l $3FFF0000,$BF6B7A18,$DACB778D,$00000000 dc.l $3FFF0000,$BFC4EA46,$63FA18F6,$00000000 dc.l $3FFF0000,$C0181BDE,$8B89A454,$00000000 dc.l $3FFF0000,$C065B066,$CFBF6439,$00000000 dc.l $3FFF0000,$C0AE345F,$56340AE6,$00000000 dc.l $3FFF0000,$C0F22291,$9CB9E6A7,$00000000 X EQU -12 XDCARE EQU X+2 XFRAC EQU X+4 XFRACLO EQU X+8 ATANF EQU -24 ATANFHI EQU ATANF+4 ATANFLO EQU ATANF+8 TEMP_SIZE EQU 24 ;--ENTRY POINT FOR ATAN(X), HERE X IS FINITE, NON-ZERO, AND NOT NAN'S _atan fmove.d (4,sp),fp0 @atan link a0,#-TEMP_SIZE fmove.x fp0,(X,a0) move.l (X,a0),d1 move.l d1,d0 ; sign move.w (XFRAC,a0),d1 and.l #$7FFFFFFF,d1 cmp.l #$3FFB8000,d1 ; |X| >= 1/16? bge.b .ATANOK1 bra.w .ATANSM .ATANOK1 cmp.l #$4002FFFF,d1 ; |X| < 16 ? ble.b .ATANMAIN bra.w .ATANBIG ;--THE MOST LIKELY CASE, |X| IN [1/16, 16). WE USE TABLE TECHNIQUE ;--THE IDEA IS ATAN(X) = ATAN(F) + ATAN( [X-F] / [1+XF] ). ;--SO IF F IS CHOSEN TO BE CLOSE TO X AND ATAN(F) IS STORED IN ;--A TABLE, ALL WE NEED IS TO APPROXIMATE ATAN(U) WHERE ;--U = (X-F)/(1+XF) IS SMALL (REMEMBER F IS CLOSE TO X). IT IS ;--TRUE THAT A DIVIDE IS NOW NEEDED, BUT THE APPROXIMATION FOR ;--ATAN(U) IS A VERY SHORT POLYNOMIAL AND THE INDEXING TO ;--FETCH F AND SAVING OF REGISTERS CAN BE ALL HIDED UNDER THE ;--DIVIDE. IN THE END THIS METHOD IS MUCH FASTER THAN A TRADITIONAL ;--ONE. NOTE ALSO THAT THE TRADITIONAL SCHEME THAT APPROXIMATE ;--ATAN(X) DIRECTLY WILL NEED TO USE A RATIONAL APPROXIMATION ;--(DIVISION NEEDED) ANYWAY BECAUSE A POLYNOMIAL APPROXIMATION ;--WILL INVOLVE A VERY LONG POLYNOMIAL. ;--NOW WE SEE X AS +-2^K * 1.BBBBBBB....B <- 1. + 63 BITS ;--WE CHOSE F TO BE +-2^K * 1.BBBB1 ;--THAT IS IT MATCHES THE EXPONENT AND FIRST 5 BITS OF X, THE ;--SIXTH BITS IS SET TO BE 1. SINCE K = -4, -3, ..., 3, THERE ;--ARE ONLY 8 TIMES 16 = 2^7 = 128 |F|'S. SINCE ATAN(-|F|) IS ;-- -ATAN(|F|), WE NEED TO STORE ONLY ATAN(|F|). .ATANMAIN and.l #$F8000000,(XFRAC,a0) ; FIRST 5 BITS or.l #$04000000,(XFRAC,a0) ; SET 6-TH BIT TO 1 clr.l (XFRACLO,a0) ; LOCATION OF X IS NOW F fmove.x fp0,fp1 ; FP1 IS X fmul.x (X,a0),fp1 ; FP1 IS X*F, NOTE THAT X*F > 0 fsub.x (X,a0),fp0 ; FP0 IS X-F fadd.s #$3F800000,fp1 ; FP1 IS 1 + X*F fdiv.x fp1,fp0 ; FP0 IS U = (X-F)/(1+X*F) ;--WHILE THE DIVISION IS TAKING ITS TIME, WE FETCH ATAN(|F|) ;--CREATE ATAN(F) AND STORE IT IN ATANF, AND ;--SAVE REGISTERS FP2. move.l d2,-(sp) ; SAVE d2 TEMPORARILY move.l d1,d2 ; THE EXP AND 16 BITS OF X and.l #$00007800,d1 ; 4 VARYING BITS OF F'S FRACTION and.l #$7FFF0000,d2 ; EXPONENT OF F sub.l #$3FFB0000,d2 ; K+4 asr.l #1,d2 add.l d2,d1 ; THE 7 BITS IDENTIFYING F asr.l #7,d1 ; INDEX INTO TBL OF ATAN(|F|) lea (ATANTBL,pc),a1 add.l d1,a1 ; ADDRESS OF ATAN(|F|) move.l (a1)+,(ATANF,a0) move.l (a1)+,(ATANFHI,a0) move.l (a1)+,(ATANFLO,a0) ; ATANF IS NOW ATAN(|F|) move.l (X,a0),d1 ; LOAD SIGN AND EXPO. AGAIN and.l #$80000000,d1 ; SIGN(F) or.l d1,(ATANF,a0) ; ATANF IS NOW SIGN(F)*ATAN(|F|) move.l (sp)+,d2 ; RESTORE d2 ;--THAT'S ALL I HAVE TO DO FOR NOW, ;--BUT ALAS, THE DIVIDE IS STILL CRANKING! ;--U IN FP0, WE ARE NOW READY TO COMPUTE ATAN(U) AS ;--U + A1*U*V*(A2 + V*(A3 + V)), V = U*U ;--THE POLYNOMIAL MAY LOOK STRANGE, BUT IS NEVERTHELESS CORRECT. ;--THE NATURAL FORM IS U + U*V*(A1 + V*(A2 + V*A3)) ;--WHAT WE HAVE HERE IS MERELY A1 = A3, A2 = A1/A3, A3 = A2/A3. ;--THE REASON FOR THIS REARRANGEMENT IS TO MAKE THE INDEPENDENT ;--PARTS A1*U*V AND (A2 + ... STUFF) MORE LOAD-BALANCED fmove.x fp2,-(sp) ; save fp2 fmove.x fp0,fp1 fmul.x fp1,fp1 fmove.d (ATANA3,pc),fp2 fadd.x fp1,fp2 ; A3+V fmul.x fp1,fp2 ; V*(A3+V) fmul.x fp0,fp1 ; U*V fadd.d (ATANA2,pc),fp2 ; A2+V*(A3+V) fmul.d (ATANA1,pc),fp1 ; A1*U*V fmul.x fp2,fp1 ; A1*U*V*(A2+V*(A3+V)) fadd.x fp1,fp0 ; ATAN(U), FP1 RELEASED fmove.x (sp)+,fp2 ; restore fp2 fadd.x (ATANF,a0),fp0 ; ATAN(X) unlk a0 rts .ATANSM ;--|X| <= 1/16 ;--IF |X| < 2^(-40), RETURN X AS ANSWER. OTHERWISE, APPROXIMATE ;--ATAN(X) BY X + X*Y*(B1+Y*(B2+Y*(B3+Y*(B4+Y*(B5+Y*B6))))) ;--WHICH IS X + X*Y*( [B1+Z*(B3+Z*B5)] + [Y*(B2+Z*(B4+Z*B6)] ) ;--WHERE Y = X*X, AND Z = Y*Y. cmp.l #$3FD78000,d1 blt.w .ATANTINY ;--COMPUTE POLYNOMIAL fmovem.x fp2/fp3,-(sp) ; save fp2/fp3 fmul.x fp0,fp0 ; FPO IS Y = X*X fmove.x fp0,fp1 fmul.x fp1,fp1 ; FP1 IS Z = Y*Y fmove.d (ATANB6,pc),fp2 fmove.d (ATANB5,pc),fp3 fmul.x fp1,fp2 ; Z*B6 fmul.x fp1,fp3 ; Z*B5 fadd.d (ATANB4,pc),fp2 ; B4+Z*B6 fadd.d (ATANB3,pc),fp3 ; B3+Z*B5 fmul.x fp1,fp2 ; Z*(B4+Z*B6) fmul.x fp3,fp1 ; Z*(B3+Z*B5) fadd.d (ATANB2,pc),fp2 ; B2+Z*(B4+Z*B6) fadd.d (ATANB1,pc),fp1 ; B1+Z*(B3+Z*B5) fmul.x fp0,fp2 ; Y*(B2+Z*(B4+Z*B6)) fmul.x (X,a0),fp0 ; X*Y fadd.x fp2,fp1 ; [B1+Z*(B3+Z*B5)]+[Y*(B2+Z*(B4+Z*B6))] fmul.x fp1,fp0 ; X*Y*([B1+Z*(B3+Z*B5)]+[Y*(B2+Z*(B4+Z*B6))]) fmovem.x (sp)+,fp2/fp3 ; restore fp2/fp3 fadd.x (X,a0),fp0 .ATANTINY unlk a0 rts ;--|X| < 2^(-40), ATAN(X) = X .ATANBIG ;--IF |X| > 2^(100), RETURN SIGN(X)*(PI/2 - TINY). OTHERWISE, ;--RETURN SIGN(X)*PI/2 + ATAN(-1/X). cmp.l #$40638000,d1 bgt.w .ATANHUGE ;--APPROXIMATE ATAN(-1/X) BY ;--X'+X'*Y*(C1+Y*(C2+Y*(C3+Y*(C4+Y*C5)))), X' = -1/X, Y = X'*X' ;--THIS CAN BE RE-WRITTEN AS ;--X'+X'*Y*( [C1+Z*(C3+Z*C5)] + [Y*(C2+Z*C4)] ), Z = Y*Y. fmovem.x fp2/fp3,-(sp) ; save fp2/fp3 fmove.s #$BF800000,fp1 ; LOAD -1 fdiv.x fp0,fp1 ; FP1 IS -1/X ;--DIVIDE IS STILL CRANKING fmove.x fp1,fp0 ; FP0 IS X' fmul.x fp0,fp0 ; FP0 IS Y = X'*X' fmove.x fp1,(X,a0) ; X IS REALLY X' fmove.x fp0,fp1 fmul.x fp1,fp1 ; FP1 IS Z = Y*Y fmove.d (ATANC5,pc),fp3 fmove.d (ATANC4,pc),fp2 fmul.x fp1,fp3 ; Z*C5 fmul.x fp1,fp2 ; Z*B4 fadd.d (ATANC3,pc),fp3 ; C3+Z*C5 fadd.d (ATANC2,pc),fp2 ; C2+Z*C4 fmul.x fp3,fp1 ; Z*(C3+Z*C5), FP3 RELEASED fmul.x fp0,fp2 ; Y*(C2+Z*C4) fadd.d (ATANC1,pc),fp1 ; C1+Z*(C3+Z*C5) fmul.x (X,a0),fp0 ; X'*Y fadd.x fp2,fp1 ; [Y*(C2+Z*C4)]+[C1+Z*(C3+Z*C5)] fmul.x fp1,fp0 ; X'*Y*([B1+Z*(B3+Z*B5)] ; +[Y*(B2+Z*(B4+Z*B6))]) fadd.x (X,a0),fp0 fmovem.x (sp)+,fp2/fp3 ; restore fp2/fp3 tst.l d0 bpl.b .pos_big .neg_big fadd.x (NPIBY2,pc),fp0 unlk a0 rts .pos_big fadd.x (PPIBY2,pc),fp0 unlk a0 rts .ATANHUGE ;--RETURN SIGN(X)*(PIBY2 - TINY) = SIGN(X)*PIBY2 - SIGN(X)*TINY tst.l d0 bpl.b .pos_huge .neg_huge fmove.x (NPIBY2,pc),fp0 fadd.x (PTINY,pc),fp0 unlk a0 rts .pos_huge fmove.x (PPIBY2,pc),fp0 fadd.x (NTINY,pc),fp0 unlk a0 rts