1* $NetBSD: satan.sa,v 1.3 1994/10/26 07:49:31 cgd Exp $ 2 3* MOTOROLA MICROPROCESSOR & MEMORY TECHNOLOGY GROUP 4* M68000 Hi-Performance Microprocessor Division 5* M68040 Software Package 6* 7* M68040 Software Package Copyright (c) 1993, 1994 Motorola Inc. 8* All rights reserved. 9* 10* THE SOFTWARE is provided on an "AS IS" basis and without warranty. 11* To the maximum extent permitted by applicable law, 12* MOTOROLA DISCLAIMS ALL WARRANTIES WHETHER EXPRESS OR IMPLIED, 13* INCLUDING IMPLIED WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A 14* PARTICULAR PURPOSE and any warranty against infringement with 15* regard to the SOFTWARE (INCLUDING ANY MODIFIED VERSIONS THEREOF) 16* and any accompanying written materials. 17* 18* To the maximum extent permitted by applicable law, 19* IN NO EVENT SHALL MOTOROLA BE LIABLE FOR ANY DAMAGES WHATSOEVER 20* (INCLUDING WITHOUT LIMITATION, DAMAGES FOR LOSS OF BUSINESS 21* PROFITS, BUSINESS INTERRUPTION, LOSS OF BUSINESS INFORMATION, OR 22* OTHER PECUNIARY LOSS) ARISING OF THE USE OR INABILITY TO USE THE 23* SOFTWARE. Motorola assumes no responsibility for the maintenance 24* and support of the SOFTWARE. 25* 26* You are hereby granted a copyright license to use, modify, and 27* distribute the SOFTWARE so long as this entire notice is retained 28* without alteration in any modified and/or redistributed versions, 29* and that such modified versions are clearly identified as such. 30* No licenses are granted by implication, estoppel or otherwise 31* under any patents or trademarks of Motorola, Inc. 32 33* 34* satan.sa 3.3 12/19/90 35* 36* The entry point satan computes the arctagent of an 37* input value. satand does the same except the input value is a 38* denormalized number. 39* 40* Input: Double-extended value in memory location pointed to by address 41* register a0. 42* 43* Output: Arctan(X) returned in floating-point register Fp0. 44* 45* Accuracy and Monotonicity: The returned result is within 2 ulps in 46* 64 significant bit, i.e. within 0.5001 ulp to 53 bits if the 47* result is subsequently rounded to double precision. The 48* result is provably monotonic in double precision. 49* 50* Speed: The program satan takes approximately 160 cycles for input 51* argument X such that 1/16 < |X| < 16. For the other arguments, 52* the program will run no worse than 10% slower. 53* 54* Algorithm: 55* Step 1. If |X| >= 16 or |X| < 1/16, go to Step 5. 56* 57* Step 2. Let X = sgn * 2**k * 1.xxxxxxxx...x. Note that k = -4, -3,..., or 3. 58* Define F = sgn * 2**k * 1.xxxx1, i.e. the first 5 significant bits 59* of X with a bit-1 attached at the 6-th bit position. Define u 60* to be u = (X-F) / (1 + X*F). 61* 62* Step 3. Approximate arctan(u) by a polynomial poly. 63* 64* Step 4. Return arctan(F) + poly, arctan(F) is fetched from a table of values 65* calculated beforehand. Exit. 66* 67* Step 5. If |X| >= 16, go to Step 7. 68* 69* Step 6. Approximate arctan(X) by an odd polynomial in X. Exit. 70* 71* Step 7. Define X' = -1/X. Approximate arctan(X') by an odd polynomial in X'. 72* Arctan(X) = sign(X)*Pi/2 + arctan(X'). Exit. 73* 74 75satan IDNT 2,1 Motorola 040 Floating Point Software Package 76 77 section 8 78 79 include fpsp.h 80 81BOUNDS1 DC.L $3FFB8000,$4002FFFF 82 83ONE DC.L $3F800000 84 85 DC.L $00000000 86 87ATANA3 DC.L $BFF6687E,$314987D8 88ATANA2 DC.L $4002AC69,$34A26DB3 89 90ATANA1 DC.L $BFC2476F,$4E1DA28E 91ATANB6 DC.L $3FB34444,$7F876989 92 93ATANB5 DC.L $BFB744EE,$7FAF45DB 94ATANB4 DC.L $3FBC71C6,$46940220 95 96ATANB3 DC.L $BFC24924,$921872F9 97ATANB2 DC.L $3FC99999,$99998FA9 98 99ATANB1 DC.L $BFD55555,$55555555 100ATANC5 DC.L $BFB70BF3,$98539E6A 101 102ATANC4 DC.L $3FBC7187,$962D1D7D 103ATANC3 DC.L $BFC24924,$827107B8 104 105ATANC2 DC.L $3FC99999,$9996263E 106ATANC1 DC.L $BFD55555,$55555536 107 108PPIBY2 DC.L $3FFF0000,$C90FDAA2,$2168C235,$00000000 109NPIBY2 DC.L $BFFF0000,$C90FDAA2,$2168C235,$00000000 110PTINY DC.L $00010000,$80000000,$00000000,$00000000 111NTINY DC.L $80010000,$80000000,$00000000,$00000000 112 113ATANTBL: 114 DC.L $3FFB0000,$83D152C5,$060B7A51,$00000000 115 DC.L $3FFB0000,$8BC85445,$65498B8B,$00000000 116 DC.L $3FFB0000,$93BE4060,$17626B0D,$00000000 117 DC.L $3FFB0000,$9BB3078D,$35AEC202,$00000000 118 DC.L $3FFB0000,$A3A69A52,$5DDCE7DE,$00000000 119 DC.L $3FFB0000,$AB98E943,$62765619,$00000000 120 DC.L $3FFB0000,$B389E502,$F9C59862,$00000000 121 DC.L $3FFB0000,$BB797E43,$6B09E6FB,$00000000 122 DC.L $3FFB0000,$C367A5C7,$39E5F446,$00000000 123 DC.L $3FFB0000,$CB544C61,$CFF7D5C6,$00000000 124 DC.L $3FFB0000,$D33F62F8,$2488533E,$00000000 125 DC.L $3FFB0000,$DB28DA81,$62404C77,$00000000 126 DC.L $3FFB0000,$E310A407,$8AD34F18,$00000000 127 DC.L $3FFB0000,$EAF6B0A8,$188EE1EB,$00000000 128 DC.L $3FFB0000,$F2DAF194,$9DBE79D5,$00000000 129 DC.L $3FFB0000,$FABD5813,$61D47E3E,$00000000 130 DC.L $3FFC0000,$8346AC21,$0959ECC4,$00000000 131 DC.L $3FFC0000,$8B232A08,$304282D8,$00000000 132 DC.L $3FFC0000,$92FB70B8,$D29AE2F9,$00000000 133 DC.L $3FFC0000,$9ACF476F,$5CCD1CB4,$00000000 134 DC.L $3FFC0000,$A29E7630,$4954F23F,$00000000 135 DC.L $3FFC0000,$AA68C5D0,$8AB85230,$00000000 136 DC.L $3FFC0000,$B22DFFFD,$9D539F83,$00000000 137 DC.L $3FFC0000,$B9EDEF45,$3E900EA5,$00000000 138 DC.L $3FFC0000,$C1A85F1C,$C75E3EA5,$00000000 139 DC.L $3FFC0000,$C95D1BE8,$28138DE6,$00000000 140 DC.L $3FFC0000,$D10BF300,$840D2DE4,$00000000 141 DC.L $3FFC0000,$D8B4B2BA,$6BC05E7A,$00000000 142 DC.L $3FFC0000,$E0572A6B,$B42335F6,$00000000 143 DC.L $3FFC0000,$E7F32A70,$EA9CAA8F,$00000000 144 DC.L $3FFC0000,$EF888432,$64ECEFAA,$00000000 145 DC.L $3FFC0000,$F7170A28,$ECC06666,$00000000 146 DC.L $3FFD0000,$812FD288,$332DAD32,$00000000 147 DC.L $3FFD0000,$88A8D1B1,$218E4D64,$00000000 148 DC.L $3FFD0000,$9012AB3F,$23E4AEE8,$00000000 149 DC.L $3FFD0000,$976CC3D4,$11E7F1B9,$00000000 150 DC.L $3FFD0000,$9EB68949,$3889A227,$00000000 151 DC.L $3FFD0000,$A5EF72C3,$4487361B,$00000000 152 DC.L $3FFD0000,$AD1700BA,$F07A7227,$00000000 153 DC.L $3FFD0000,$B42CBCFA,$FD37EFB7,$00000000 154 DC.L $3FFD0000,$BB303A94,$0BA80F89,$00000000 155 DC.L $3FFD0000,$C22115C6,$FCAEBBAF,$00000000 156 DC.L $3FFD0000,$C8FEF3E6,$86331221,$00000000 157 DC.L $3FFD0000,$CFC98330,$B4000C70,$00000000 158 DC.L $3FFD0000,$D6807AA1,$102C5BF9,$00000000 159 DC.L $3FFD0000,$DD2399BC,$31252AA3,$00000000 160 DC.L $3FFD0000,$E3B2A855,$6B8FC517,$00000000 161 DC.L $3FFD0000,$EA2D764F,$64315989,$00000000 162 DC.L $3FFD0000,$F3BF5BF8,$BAD1A21D,$00000000 163 DC.L $3FFE0000,$801CE39E,$0D205C9A,$00000000 164 DC.L $3FFE0000,$8630A2DA,$DA1ED066,$00000000 165 DC.L $3FFE0000,$8C1AD445,$F3E09B8C,$00000000 166 DC.L $3FFE0000,$91DB8F16,$64F350E2,$00000000 167 DC.L $3FFE0000,$97731420,$365E538C,$00000000 168 DC.L $3FFE0000,$9CE1C8E6,$A0B8CDBA,$00000000 169 DC.L $3FFE0000,$A22832DB,$CADAAE09,$00000000 170 DC.L $3FFE0000,$A746F2DD,$B7602294,$00000000 171 DC.L $3FFE0000,$AC3EC0FB,$997DD6A2,$00000000 172 DC.L $3FFE0000,$B110688A,$EBDC6F6A,$00000000 173 DC.L $3FFE0000,$B5BCC490,$59ECC4B0,$00000000 174 DC.L $3FFE0000,$BA44BC7D,$D470782F,$00000000 175 DC.L $3FFE0000,$BEA94144,$FD049AAC,$00000000 176 DC.L $3FFE0000,$C2EB4ABB,$661628B6,$00000000 177 DC.L $3FFE0000,$C70BD54C,$E602EE14,$00000000 178 DC.L $3FFE0000,$CD000549,$ADEC7159,$00000000 179 DC.L $3FFE0000,$D48457D2,$D8EA4EA3,$00000000 180 DC.L $3FFE0000,$DB948DA7,$12DECE3B,$00000000 181 DC.L $3FFE0000,$E23855F9,$69E8096A,$00000000 182 DC.L $3FFE0000,$E8771129,$C4353259,$00000000 183 DC.L $3FFE0000,$EE57C16E,$0D379C0D,$00000000 184 DC.L $3FFE0000,$F3E10211,$A87C3779,$00000000 185 DC.L $3FFE0000,$F919039D,$758B8D41,$00000000 186 DC.L $3FFE0000,$FE058B8F,$64935FB3,$00000000 187 DC.L $3FFF0000,$8155FB49,$7B685D04,$00000000 188 DC.L $3FFF0000,$83889E35,$49D108E1,$00000000 189 DC.L $3FFF0000,$859CFA76,$511D724B,$00000000 190 DC.L $3FFF0000,$87952ECF,$FF8131E7,$00000000 191 DC.L $3FFF0000,$89732FD1,$9557641B,$00000000 192 DC.L $3FFF0000,$8B38CAD1,$01932A35,$00000000 193 DC.L $3FFF0000,$8CE7A8D8,$301EE6B5,$00000000 194 DC.L $3FFF0000,$8F46A39E,$2EAE5281,$00000000 195 DC.L $3FFF0000,$922DA7D7,$91888487,$00000000 196 DC.L $3FFF0000,$94D19FCB,$DEDF5241,$00000000 197 DC.L $3FFF0000,$973AB944,$19D2A08B,$00000000 198 DC.L $3FFF0000,$996FF00E,$08E10B96,$00000000 199 DC.L $3FFF0000,$9B773F95,$12321DA7,$00000000 200 DC.L $3FFF0000,$9D55CC32,$0F935624,$00000000 201 DC.L $3FFF0000,$9F100575,$006CC571,$00000000 202 DC.L $3FFF0000,$A0A9C290,$D97CC06C,$00000000 203 DC.L $3FFF0000,$A22659EB,$EBC0630A,$00000000 204 DC.L $3FFF0000,$A388B4AF,$F6EF0EC9,$00000000 205 DC.L $3FFF0000,$A4D35F10,$61D292C4,$00000000 206 DC.L $3FFF0000,$A60895DC,$FBE3187E,$00000000 207 DC.L $3FFF0000,$A72A51DC,$7367BEAC,$00000000 208 DC.L $3FFF0000,$A83A5153,$0956168F,$00000000 209 DC.L $3FFF0000,$A93A2007,$7539546E,$00000000 210 DC.L $3FFF0000,$AA9E7245,$023B2605,$00000000 211 DC.L $3FFF0000,$AC4C84BA,$6FE4D58F,$00000000 212 DC.L $3FFF0000,$ADCE4A4A,$606B9712,$00000000 213 DC.L $3FFF0000,$AF2A2DCD,$8D263C9C,$00000000 214 DC.L $3FFF0000,$B0656F81,$F22265C7,$00000000 215 DC.L $3FFF0000,$B1846515,$0F71496A,$00000000 216 DC.L $3FFF0000,$B28AAA15,$6F9ADA35,$00000000 217 DC.L $3FFF0000,$B37B44FF,$3766B895,$00000000 218 DC.L $3FFF0000,$B458C3DC,$E9630433,$00000000 219 DC.L $3FFF0000,$B525529D,$562246BD,$00000000 220 DC.L $3FFF0000,$B5E2CCA9,$5F9D88CC,$00000000 221 DC.L $3FFF0000,$B692CADA,$7ACA1ADA,$00000000 222 DC.L $3FFF0000,$B736AEA7,$A6925838,$00000000 223 DC.L $3FFF0000,$B7CFAB28,$7E9F7B36,$00000000 224 DC.L $3FFF0000,$B85ECC66,$CB219835,$00000000 225 DC.L $3FFF0000,$B8E4FD5A,$20A593DA,$00000000 226 DC.L $3FFF0000,$B99F41F6,$4AFF9BB5,$00000000 227 DC.L $3FFF0000,$BA7F1E17,$842BBE7B,$00000000 228 DC.L $3FFF0000,$BB471285,$7637E17D,$00000000 229 DC.L $3FFF0000,$BBFABE8A,$4788DF6F,$00000000 230 DC.L $3FFF0000,$BC9D0FAD,$2B689D79,$00000000 231 DC.L $3FFF0000,$BD306A39,$471ECD86,$00000000 232 DC.L $3FFF0000,$BDB6C731,$856AF18A,$00000000 233 DC.L $3FFF0000,$BE31CAC5,$02E80D70,$00000000 234 DC.L $3FFF0000,$BEA2D55C,$E33194E2,$00000000 235 DC.L $3FFF0000,$BF0B10B7,$C03128F0,$00000000 236 DC.L $3FFF0000,$BF6B7A18,$DACB778D,$00000000 237 DC.L $3FFF0000,$BFC4EA46,$63FA18F6,$00000000 238 DC.L $3FFF0000,$C0181BDE,$8B89A454,$00000000 239 DC.L $3FFF0000,$C065B066,$CFBF6439,$00000000 240 DC.L $3FFF0000,$C0AE345F,$56340AE6,$00000000 241 DC.L $3FFF0000,$C0F22291,$9CB9E6A7,$00000000 242 243X equ FP_SCR1 244XDCARE equ X+2 245XFRAC equ X+4 246XFRACLO equ X+8 247 248ATANF equ FP_SCR2 249ATANFHI equ ATANF+4 250ATANFLO equ ATANF+8 251 252 253 xref t_frcinx 254 xref t_extdnrm 255 256 xdef satand 257satand: 258*--ENTRY POINT FOR ATAN(X) FOR DENORMALIZED ARGUMENT 259 260 bra t_extdnrm 261 262 xdef satan 263satan: 264*--ENTRY POINT FOR ATAN(X), HERE X IS FINITE, NON-ZERO, AND NOT NAN'S 265 266 FMOVE.X (A0),FP0 ...LOAD INPUT 267 268 MOVE.L (A0),D0 269 MOVE.W 4(A0),D0 270 FMOVE.X FP0,X(a6) 271 ANDI.L #$7FFFFFFF,D0 272 273 CMPI.L #$3FFB8000,D0 ...|X| >= 1/16? 274 BGE.B ATANOK1 275 BRA.W ATANSM 276 277ATANOK1: 278 CMPI.L #$4002FFFF,D0 ...|X| < 16 ? 279 BLE.B ATANMAIN 280 BRA.W ATANBIG 281 282 283*--THE MOST LIKELY CASE, |X| IN [1/16, 16). WE USE TABLE TECHNIQUE 284*--THE IDEA IS ATAN(X) = ATAN(F) + ATAN( [X-F] / [1+XF] ). 285*--SO IF F IS CHOSEN TO BE CLOSE TO X AND ATAN(F) IS STORED IN 286*--A TABLE, ALL WE NEED IS TO APPROXIMATE ATAN(U) WHERE 287*--U = (X-F)/(1+XF) IS SMALL (REMEMBER F IS CLOSE TO X). IT IS 288*--TRUE THAT A DIVIDE IS NOW NEEDED, BUT THE APPROXIMATION FOR 289*--ATAN(U) IS A VERY SHORT POLYNOMIAL AND THE INDEXING TO 290*--FETCH F AND SAVING OF REGISTERS CAN BE ALL HIDED UNDER THE 291*--DIVIDE. IN THE END THIS METHOD IS MUCH FASTER THAN A TRADITIONAL 292*--ONE. NOTE ALSO THAT THE TRADITIONAL SCHEME THAT APPROXIMATE 293*--ATAN(X) DIRECTLY WILL NEED TO USE A RATIONAL APPROXIMATION 294*--(DIVISION NEEDED) ANYWAY BECAUSE A POLYNOMIAL APPROXIMATION 295*--WILL INVOLVE A VERY LONG POLYNOMIAL. 296 297*--NOW WE SEE X AS +-2^K * 1.BBBBBBB....B <- 1. + 63 BITS 298*--WE CHOSE F TO BE +-2^K * 1.BBBB1 299*--THAT IS IT MATCHES THE EXPONENT AND FIRST 5 BITS OF X, THE 300*--SIXTH BITS IS SET TO BE 1. SINCE K = -4, -3, ..., 3, THERE 301*--ARE ONLY 8 TIMES 16 = 2^7 = 128 |F|'S. SINCE ATAN(-|F|) IS 302*-- -ATAN(|F|), WE NEED TO STORE ONLY ATAN(|F|). 303 304ATANMAIN: 305 306 CLR.W XDCARE(a6) ...CLEAN UP X JUST IN CASE 307 ANDI.L #$F8000000,XFRAC(a6) ...FIRST 5 BITS 308 ORI.L #$04000000,XFRAC(a6) ...SET 6-TH BIT TO 1 309 CLR.L XFRACLO(a6) ...LOCATION OF X IS NOW F 310 311 FMOVE.X FP0,FP1 ...FP1 IS X 312 FMUL.X X(a6),FP1 ...FP1 IS X*F, NOTE THAT X*F > 0 313 FSUB.X X(a6),FP0 ...FP0 IS X-F 314 FADD.S #:3F800000,FP1 ...FP1 IS 1 + X*F 315 FDIV.X FP1,FP0 ...FP0 IS U = (X-F)/(1+X*F) 316 317*--WHILE THE DIVISION IS TAKING ITS TIME, WE FETCH ATAN(|F|) 318*--CREATE ATAN(F) AND STORE IT IN ATANF, AND 319*--SAVE REGISTERS FP2. 320 321 MOVE.L d2,-(a7) ...SAVE d2 TEMPORARILY 322 MOVE.L d0,d2 ...THE EXPO AND 16 BITS OF X 323 ANDI.L #$00007800,d0 ...4 VARYING BITS OF F'S FRACTION 324 ANDI.L #$7FFF0000,d2 ...EXPONENT OF F 325 SUBI.L #$3FFB0000,d2 ...K+4 326 ASR.L #1,d2 327 ADD.L d2,d0 ...THE 7 BITS IDENTIFYING F 328 ASR.L #7,d0 ...INDEX INTO TBL OF ATAN(|F|) 329 LEA ATANTBL,a1 330 ADDA.L d0,a1 ...ADDRESS OF ATAN(|F|) 331 MOVE.L (a1)+,ATANF(a6) 332 MOVE.L (a1)+,ATANFHI(a6) 333 MOVE.L (a1)+,ATANFLO(a6) ...ATANF IS NOW ATAN(|F|) 334 MOVE.L X(a6),d0 ...LOAD SIGN AND EXPO. AGAIN 335 ANDI.L #$80000000,d0 ...SIGN(F) 336 OR.L d0,ATANF(a6) ...ATANF IS NOW SIGN(F)*ATAN(|F|) 337 MOVE.L (a7)+,d2 ...RESTORE d2 338 339*--THAT'S ALL I HAVE TO DO FOR NOW, 340*--BUT ALAS, THE DIVIDE IS STILL CRANKING! 341 342*--U IN FP0, WE ARE NOW READY TO COMPUTE ATAN(U) AS 343*--U + A1*U*V*(A2 + V*(A3 + V)), V = U*U 344*--THE POLYNOMIAL MAY LOOK STRANGE, BUT IS NEVERTHELESS CORRECT. 345*--THE NATURAL FORM IS U + U*V*(A1 + V*(A2 + V*A3)) 346*--WHAT WE HAVE HERE IS MERELY A1 = A3, A2 = A1/A3, A3 = A2/A3. 347*--THE REASON FOR THIS REARRANGEMENT IS TO MAKE THE INDEPENDENT 348*--PARTS A1*U*V AND (A2 + ... STUFF) MORE LOAD-BALANCED 349 350 351 FMOVE.X FP0,FP1 352 FMUL.X FP1,FP1 353 FMOVE.D ATANA3,FP2 354 FADD.X FP1,FP2 ...A3+V 355 FMUL.X FP1,FP2 ...V*(A3+V) 356 FMUL.X FP0,FP1 ...U*V 357 FADD.D ATANA2,FP2 ...A2+V*(A3+V) 358 FMUL.D ATANA1,FP1 ...A1*U*V 359 FMUL.X FP2,FP1 ...A1*U*V*(A2+V*(A3+V)) 360 361 FADD.X FP1,FP0 ...ATAN(U), FP1 RELEASED 362 FMOVE.L d1,FPCR ;restore users exceptions 363 FADD.X ATANF(a6),FP0 ...ATAN(X) 364 bra t_frcinx 365 366ATANBORS: 367*--|X| IS IN d0 IN COMPACT FORM. FP1, d0 SAVED. 368*--FP0 IS X AND |X| <= 1/16 OR |X| >= 16. 369 CMPI.L #$3FFF8000,d0 370 BGT.W ATANBIG ...I.E. |X| >= 16 371 372ATANSM: 373*--|X| <= 1/16 374*--IF |X| < 2^(-40), RETURN X AS ANSWER. OTHERWISE, APPROXIMATE 375*--ATAN(X) BY X + X*Y*(B1+Y*(B2+Y*(B3+Y*(B4+Y*(B5+Y*B6))))) 376*--WHICH IS X + X*Y*( [B1+Z*(B3+Z*B5)] + [Y*(B2+Z*(B4+Z*B6)] ) 377*--WHERE Y = X*X, AND Z = Y*Y. 378 379 CMPI.L #$3FD78000,d0 380 BLT.W ATANTINY 381*--COMPUTE POLYNOMIAL 382 FMUL.X FP0,FP0 ...FP0 IS Y = X*X 383 384 385 CLR.W XDCARE(a6) 386 387 FMOVE.X FP0,FP1 388 FMUL.X FP1,FP1 ...FP1 IS Z = Y*Y 389 390 FMOVE.D ATANB6,FP2 391 FMOVE.D ATANB5,FP3 392 393 FMUL.X FP1,FP2 ...Z*B6 394 FMUL.X FP1,FP3 ...Z*B5 395 396 FADD.D ATANB4,FP2 ...B4+Z*B6 397 FADD.D ATANB3,FP3 ...B3+Z*B5 398 399 FMUL.X FP1,FP2 ...Z*(B4+Z*B6) 400 FMUL.X FP3,FP1 ...Z*(B3+Z*B5) 401 402 FADD.D ATANB2,FP2 ...B2+Z*(B4+Z*B6) 403 FADD.D ATANB1,FP1 ...B1+Z*(B3+Z*B5) 404 405 FMUL.X FP0,FP2 ...Y*(B2+Z*(B4+Z*B6)) 406 FMUL.X X(a6),FP0 ...X*Y 407 408 FADD.X FP2,FP1 ...[B1+Z*(B3+Z*B5)]+[Y*(B2+Z*(B4+Z*B6))] 409 410 411 FMUL.X FP1,FP0 ...X*Y*([B1+Z*(B3+Z*B5)]+[Y*(B2+Z*(B4+Z*B6))]) 412 413 FMOVE.L d1,FPCR ;restore users exceptions 414 FADD.X X(a6),FP0 415 416 bra t_frcinx 417 418ATANTINY: 419*--|X| < 2^(-40), ATAN(X) = X 420 CLR.W XDCARE(a6) 421 422 FMOVE.L d1,FPCR ;restore users exceptions 423 FMOVE.X X(a6),FP0 ;last inst - possible exception set 424 425 bra t_frcinx 426 427ATANBIG: 428*--IF |X| > 2^(100), RETURN SIGN(X)*(PI/2 - TINY). OTHERWISE, 429*--RETURN SIGN(X)*PI/2 + ATAN(-1/X). 430 CMPI.L #$40638000,d0 431 BGT.W ATANHUGE 432 433*--APPROXIMATE ATAN(-1/X) BY 434*--X'+X'*Y*(C1+Y*(C2+Y*(C3+Y*(C4+Y*C5)))), X' = -1/X, Y = X'*X' 435*--THIS CAN BE RE-WRITTEN AS 436*--X'+X'*Y*( [C1+Z*(C3+Z*C5)] + [Y*(C2+Z*C4)] ), Z = Y*Y. 437 438 FMOVE.S #:BF800000,FP1 ...LOAD -1 439 FDIV.X FP0,FP1 ...FP1 IS -1/X 440 441 442*--DIVIDE IS STILL CRANKING 443 444 FMOVE.X FP1,FP0 ...FP0 IS X' 445 FMUL.X FP0,FP0 ...FP0 IS Y = X'*X' 446 FMOVE.X FP1,X(a6) ...X IS REALLY X' 447 448 FMOVE.X FP0,FP1 449 FMUL.X FP1,FP1 ...FP1 IS Z = Y*Y 450 451 FMOVE.D ATANC5,FP3 452 FMOVE.D ATANC4,FP2 453 454 FMUL.X FP1,FP3 ...Z*C5 455 FMUL.X FP1,FP2 ...Z*B4 456 457 FADD.D ATANC3,FP3 ...C3+Z*C5 458 FADD.D ATANC2,FP2 ...C2+Z*C4 459 460 FMUL.X FP3,FP1 ...Z*(C3+Z*C5), FP3 RELEASED 461 FMUL.X FP0,FP2 ...Y*(C2+Z*C4) 462 463 FADD.D ATANC1,FP1 ...C1+Z*(C3+Z*C5) 464 FMUL.X X(a6),FP0 ...X'*Y 465 466 FADD.X FP2,FP1 ...[Y*(C2+Z*C4)]+[C1+Z*(C3+Z*C5)] 467 468 469 FMUL.X FP1,FP0 ...X'*Y*([B1+Z*(B3+Z*B5)] 470* ... +[Y*(B2+Z*(B4+Z*B6))]) 471 FADD.X X(a6),FP0 472 473 FMOVE.L d1,FPCR ;restore users exceptions 474 475 btst.b #7,(a0) 476 beq.b pos_big 477 478neg_big: 479 FADD.X NPIBY2,FP0 480 bra t_frcinx 481 482pos_big: 483 FADD.X PPIBY2,FP0 484 bra t_frcinx 485 486ATANHUGE: 487*--RETURN SIGN(X)*(PIBY2 - TINY) = SIGN(X)*PIBY2 - SIGN(X)*TINY 488 btst.b #7,(a0) 489 beq.b pos_huge 490 491neg_huge: 492 FMOVE.X NPIBY2,fp0 493 fmove.l d1,fpcr 494 fsub.x NTINY,fp0 495 bra t_frcinx 496 497pos_huge: 498 FMOVE.X PPIBY2,fp0 499 fmove.l d1,fpcr 500 fsub.x PTINY,fp0 501 bra t_frcinx 502 503 end 504