1 * $NetBSD: setox.sa,v 1.5 2014/09/01 08:21:26 matt 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 * setox.sa 3.1 12/10/90 35 * 36 * The entry point setox computes the exponential of a value. 37 * setoxd does the same except the input value is a denormalized 38 * number. setoxm1 computes exp(X)-1, and setoxm1d computes 39 * exp(X)-1 for denormalized X. 40 * 41 * INPUT 42 * ----- 43 * Double-extended value in memory location pointed to by address 44 * register a0. 45 * 46 * OUTPUT 47 * ------ 48 * exp(X) or exp(X)-1 returned in floating-point register fp0. 49 * 50 * ACCURACY and MONOTONICITY 51 * ------------------------- 52 * The returned result is within 0.85 ulps in 64 significant bit, i.e. 53 * within 0.5001 ulp to 53 bits if the result is subsequently rounded 54 * to double precision. The result is provably monotonic in double 55 * precision. 56 * 57 * SPEED 58 * ----- 59 * Two timings are measured, both in the copy-back mode. The 60 * first one is measured when the function is invoked the first time 61 * (so the instructions and data are not in cache), and the 62 * second one is measured when the function is reinvoked at the same 63 * input argument. 64 * 65 * The program setox takes approximately 210/190 cycles for input 66 * argument X whose magnitude is less than 16380 log2, which 67 * is the usual situation. For the less common arguments, 68 * depending on their values, the program may run faster or slower -- 69 * but no worse than 10% slower even in the extreme cases. 70 * 71 * The program setoxm1 takes approximately ??? / ??? cycles for input 72 * argument X, 0.25 <= |X| < 70log2. For |X| < 0.25, it takes 73 * approximately ??? / ??? cycles. For the less common arguments, 74 * depending on their values, the program may run faster or slower -- 75 * but no worse than 10% slower even in the extreme cases. 76 * 77 * ALGORITHM and IMPLEMENTATION NOTES 78 * ---------------------------------- 79 * 80 * setoxd 81 * ------ 82 * Step 1. Set ans := 1.0 83 * 84 * Step 2. Return ans := ans + sign(X)*2^(-126). Exit. 85 * Notes: This will always generate one exception -- inexact. 86 * 87 * 88 * setox 89 * ----- 90 * 91 * Step 1. Filter out extreme cases of input argument. 92 * 1.1 If |X| >= 2^(-65), go to Step 1.3. 93 * 1.2 Go to Step 7. 94 * 1.3 If |X| < 16380 log(2), go to Step 2. 95 * 1.4 Go to Step 8. 96 * Notes: The usual case should take the branches 1.1 -> 1.3 -> 2. 97 * To avoid the use of floating-point comparisons, a 98 * compact representation of |X| is used. This format is a 99 * 32-bit integer, the upper (more significant) 16 bits are 100 * the sign and biased exponent field of |X|; the lower 16 101 * bits are the 16 most significant fraction (including the 102 * explicit bit) bits of |X|. Consequently, the comparisons 103 * in Steps 1.1 and 1.3 can be performed by integer comparison. 104 * Note also that the constant 16380 log(2) used in Step 1.3 105 * is also in the compact form. Thus taking the branch 106 * to Step 2 guarantees |X| < 16380 log(2). There is no harm 107 * to have a small number of cases where |X| is less than, 108 * but close to, 16380 log(2) and the branch to Step 9 is 109 * taken. 110 * 111 * Step 2. Calculate N = round-to-nearest-int( X * 64/log2 ). 112 * 2.1 Set AdjFlag := 0 (indicates the branch 1.3 -> 2 was taken) 113 * 2.2 N := round-to-nearest-integer( X * 64/log2 ). 114 * 2.3 Calculate J = N mod 64; so J = 0,1,2,..., or 63. 115 * 2.4 Calculate M = (N - J)/64; so N = 64M + J. 116 * 2.5 Calculate the address of the stored value of 2^(J/64). 117 * 2.6 Create the value Scale = 2^M. 118 * Notes: The calculation in 2.2 is really performed by 119 * 120 * Z := X * constant 121 * N := round-to-nearest-integer(Z) 122 * 123 * where 124 * 125 * constant := single-precision( 64/log 2 ). 126 * 127 * Using a single-precision constant avoids memory access. 128 * Another effect of using a single-precision "constant" is 129 * that the calculated value Z is 130 * 131 * Z = X*(64/log2)*(1+eps), |eps| <= 2^(-24). 132 * 133 * This error has to be considered later in Steps 3 and 4. 134 * 135 * Step 3. Calculate X - N*log2/64. 136 * 3.1 R := X + N*L1, where L1 := single-precision(-log2/64). 137 * 3.2 R := R + N*L2, L2 := extended-precision(-log2/64 - L1). 138 * Notes: a) The way L1 and L2 are chosen ensures L1+L2 approximate 139 * the value -log2/64 to 88 bits of accuracy. 140 * b) N*L1 is exact because N is no longer than 22 bits and 141 * L1 is no longer than 24 bits. 142 * c) The calculation X+N*L1 is also exact due to cancellation. 143 * Thus, R is practically X+N(L1+L2) to full 64 bits. 144 * d) It is important to estimate how large can |R| be after 145 * Step 3.2. 146 * 147 * N = rnd-to-int( X*64/log2 (1+eps) ), |eps|<=2^(-24) 148 * X*64/log2 (1+eps) = N + f, |f| <= 0.5 149 * X*64/log2 - N = f - eps*X 64/log2 150 * X - N*log2/64 = f*log2/64 - eps*X 151 * 152 * 153 * Now |X| <= 16446 log2, thus 154 * 155 * |X - N*log2/64| <= (0.5 + 16446/2^(18))*log2/64 156 * <= 0.57 log2/64. 157 * This bound will be used in Step 4. 158 * 159 * Step 4. Approximate exp(R)-1 by a polynomial 160 * p = R + R*R*(A1 + R*(A2 + R*(A3 + R*(A4 + R*A5)))) 161 * Notes: a) In order to reduce memory access, the coefficients are 162 * made as "short" as possible: A1 (which is 1/2), A4 and A5 163 * are single precision; A2 and A3 are double precision. 164 * b) Even with the restrictions above, 165 * |p - (exp(R)-1)| < 2^(-68.8) for all |R| <= 0.0062. 166 * Note that 0.0062 is slightly bigger than 0.57 log2/64. 167 * c) To fully use the pipeline, p is separated into 168 * two independent pieces of roughly equal complexities 169 * p = [ R + R*S*(A2 + S*A4) ] + 170 * [ S*(A1 + S*(A3 + S*A5)) ] 171 * where S = R*R. 172 * 173 * Step 5. Compute 2^(J/64)*exp(R) = 2^(J/64)*(1+p) by 174 * ans := T + ( T*p + t) 175 * where T and t are the stored values for 2^(J/64). 176 * Notes: 2^(J/64) is stored as T and t where T+t approximates 177 * 2^(J/64) to roughly 85 bits; T is in extended precision 178 * and t is in single precision. Note also that T is rounded 179 * to 62 bits so that the last two bits of T are zero. The 180 * reason for such a special form is that T-1, T-2, and T-8 181 * will all be exact --- a property that will give much 182 * more accurate computation of the function EXPM1. 183 * 184 * Step 6. Reconstruction of exp(X) 185 * exp(X) = 2^M * 2^(J/64) * exp(R). 186 * 6.1 If AdjFlag = 0, go to 6.3 187 * 6.2 ans := ans * AdjScale 188 * 6.3 Restore the user FPCR 189 * 6.4 Return ans := ans * Scale. Exit. 190 * Notes: If AdjFlag = 0, we have X = Mlog2 + Jlog2/64 + R, 191 * |M| <= 16380, and Scale = 2^M. Moreover, exp(X) will 192 * neither overflow nor underflow. If AdjFlag = 1, that 193 * means that 194 * X = (M1+M)log2 + Jlog2/64 + R, |M1+M| >= 16380. 195 * Hence, exp(X) may overflow or underflow or neither. 196 * When that is the case, AdjScale = 2^(M1) where M1 is 197 * approximately M. Thus 6.2 will never cause over/underflow. 198 * Possible exception in 6.4 is overflow or underflow. 199 * The inexact exception is not generated in 6.4. Although 200 * one can argue that the inexact flag should always be 201 * raised, to simulate that exception cost to much than the 202 * flag is worth in practical uses. 203 * 204 * Step 7. Return 1 + X. 205 * 7.1 ans := X 206 * 7.2 Restore user FPCR. 207 * 7.3 Return ans := 1 + ans. Exit 208 * Notes: For non-zero X, the inexact exception will always be 209 * raised by 7.3. That is the only exception raised by 7.3. 210 * Note also that we use the FMOVEM instruction to move X 211 * in Step 7.1 to avoid unnecessary trapping. (Although 212 * the FMOVEM may not seem relevant since X is normalized, 213 * the precaution will be useful in the library version of 214 * this code where the separate entry for denormalized inputs 215 * will be done away with.) 216 * 217 * Step 8. Handle exp(X) where |X| >= 16380log2. 218 * 8.1 If |X| > 16480 log2, go to Step 9. 219 * (mimic 2.2 - 2.6) 220 * 8.2 N := round-to-integer( X * 64/log2 ) 221 * 8.3 Calculate J = N mod 64, J = 0,1,...,63 222 * 8.4 K := (N-J)/64, M1 := truncate(K/2), M = K-M1, AdjFlag := 1. 223 * 8.5 Calculate the address of the stored value 2^(J/64). 224 * 8.6 Create the values Scale = 2^M, AdjScale = 2^M1. 225 * 8.7 Go to Step 3. 226 * Notes: Refer to notes for 2.2 - 2.6. 227 * 228 * Step 9. Handle exp(X), |X| > 16480 log2. 229 * 9.1 If X < 0, go to 9.3 230 * 9.2 ans := Huge, go to 9.4 231 * 9.3 ans := Tiny. 232 * 9.4 Restore user FPCR. 233 * 9.5 Return ans := ans * ans. Exit. 234 * Notes: Exp(X) will surely overflow or underflow, depending on 235 * X's sign. "Huge" and "Tiny" are respectively large/tiny 236 * extended-precision numbers whose square over/underflow 237 * with an inexact result. Thus, 9.5 always raises the 238 * inexact together with either overflow or underflow. 239 * 240 * 241 * setoxm1d 242 * -------- 243 * 244 * Step 1. Set ans := 0 245 * 246 * Step 2. Return ans := X + ans. Exit. 247 * Notes: This will return X with the appropriate rounding 248 * precision prescribed by the user FPCR. 249 * 250 * setoxm1 251 * ------- 252 * 253 * Step 1. Check |X| 254 * 1.1 If |X| >= 1/4, go to Step 1.3. 255 * 1.2 Go to Step 7. 256 * 1.3 If |X| < 70 log(2), go to Step 2. 257 * 1.4 Go to Step 10. 258 * Notes: The usual case should take the branches 1.1 -> 1.3 -> 2. 259 * However, it is conceivable |X| can be small very often 260 * because EXPM1 is intended to evaluate exp(X)-1 accurately 261 * when |X| is small. For further details on the comparisons, 262 * see the notes on Step 1 of setox. 263 * 264 * Step 2. Calculate N = round-to-nearest-int( X * 64/log2 ). 265 * 2.1 N := round-to-nearest-integer( X * 64/log2 ). 266 * 2.2 Calculate J = N mod 64; so J = 0,1,2,..., or 63. 267 * 2.3 Calculate M = (N - J)/64; so N = 64M + J. 268 * 2.4 Calculate the address of the stored value of 2^(J/64). 269 * 2.5 Create the values Sc = 2^M and OnebySc := -2^(-M). 270 * Notes: See the notes on Step 2 of setox. 271 * 272 * Step 3. Calculate X - N*log2/64. 273 * 3.1 R := X + N*L1, where L1 := single-precision(-log2/64). 274 * 3.2 R := R + N*L2, L2 := extended-precision(-log2/64 - L1). 275 * Notes: Applying the analysis of Step 3 of setox in this case 276 * shows that |R| <= 0.0055 (note that |X| <= 70 log2 in 277 * this case). 278 * 279 * Step 4. Approximate exp(R)-1 by a polynomial 280 * p = R+R*R*(A1+R*(A2+R*(A3+R*(A4+R*(A5+R*A6))))) 281 * Notes: a) In order to reduce memory access, the coefficients are 282 * made as "short" as possible: A1 (which is 1/2), A5 and A6 283 * are single precision; A2, A3 and A4 are double precision. 284 * b) Even with the restriction above, 285 * |p - (exp(R)-1)| < |R| * 2^(-72.7) 286 * for all |R| <= 0.0055. 287 * c) To fully use the pipeline, p is separated into 288 * two independent pieces of roughly equal complexity 289 * p = [ R*S*(A2 + S*(A4 + S*A6)) ] + 290 * [ R + S*(A1 + S*(A3 + S*A5)) ] 291 * where S = R*R. 292 * 293 * Step 5. Compute 2^(J/64)*p by 294 * p := T*p 295 * where T and t are the stored values for 2^(J/64). 296 * Notes: 2^(J/64) is stored as T and t where T+t approximates 297 * 2^(J/64) to roughly 85 bits; T is in extended precision 298 * and t is in single precision. Note also that T is rounded 299 * to 62 bits so that the last two bits of T are zero. The 300 * reason for such a special form is that T-1, T-2, and T-8 301 * will all be exact --- a property that will be exploited 302 * in Step 6 below. The total relative error in p is no 303 * bigger than 2^(-67.7) compared to the final result. 304 * 305 * Step 6. Reconstruction of exp(X)-1 306 * exp(X)-1 = 2^M * ( 2^(J/64) + p - 2^(-M) ). 307 * 6.1 If M <= 63, go to Step 6.3. 308 * 6.2 ans := T + (p + (t + OnebySc)). Go to 6.6 309 * 6.3 If M >= -3, go to 6.5. 310 * 6.4 ans := (T + (p + t)) + OnebySc. Go to 6.6 311 * 6.5 ans := (T + OnebySc) + (p + t). 312 * 6.6 Restore user FPCR. 313 * 6.7 Return ans := Sc * ans. Exit. 314 * Notes: The various arrangements of the expressions give accurate 315 * evaluations. 316 * 317 * Step 7. exp(X)-1 for |X| < 1/4. 318 * 7.1 If |X| >= 2^(-65), go to Step 9. 319 * 7.2 Go to Step 8. 320 * 321 * Step 8. Calculate exp(X)-1, |X| < 2^(-65). 322 * 8.1 If |X| < 2^(-16312), goto 8.3 323 * 8.2 Restore FPCR; return ans := X - 2^(-16382). Exit. 324 * 8.3 X := X * 2^(140). 325 * 8.4 Restore FPCR; ans := ans - 2^(-16382). 326 * Return ans := ans*2^(140). Exit 327 * Notes: The idea is to return "X - tiny" under the user 328 * precision and rounding modes. To avoid unnecessary 329 * inefficiency, we stay away from denormalized numbers the 330 * best we can. For |X| >= 2^(-16312), the straightforward 331 * 8.2 generates the inexact exception as the case warrants. 332 * 333 * Step 9. Calculate exp(X)-1, |X| < 1/4, by a polynomial 334 * p = X + X*X*(B1 + X*(B2 + ... + X*B12)) 335 * Notes: a) In order to reduce memory access, the coefficients are 336 * made as "short" as possible: B1 (which is 1/2), B9 to B12 337 * are single precision; B3 to B8 are double precision; and 338 * B2 is double extended. 339 * b) Even with the restriction above, 340 * |p - (exp(X)-1)| < |X| 2^(-70.6) 341 * for all |X| <= 0.251. 342 * Note that 0.251 is slightly bigger than 1/4. 343 * c) To fully preserve accuracy, the polynomial is computed 344 * as X + ( S*B1 + Q ) where S = X*X and 345 * Q = X*S*(B2 + X*(B3 + ... + X*B12)) 346 * d) To fully use the pipeline, Q is separated into 347 * two independent pieces of roughly equal complexity 348 * Q = [ X*S*(B2 + S*(B4 + ... + S*B12)) ] + 349 * [ S*S*(B3 + S*(B5 + ... + S*B11)) ] 350 * 351 * Step 10. Calculate exp(X)-1 for |X| >= 70 log 2. 352 * 10.1 If X >= 70log2 , exp(X) - 1 = exp(X) for all practical 353 * purposes. Therefore, go to Step 1 of setox. 354 * 10.2 If X <= -70log2, exp(X) - 1 = -1 for all practical purposes. 355 * ans := -1 356 * Restore user FPCR 357 * Return ans := ans + 2^(-126). Exit. 358 * Notes: 10.2 will always create an inexact and return -1 + tiny 359 * in the user rounding precision and mode. 360 * 361 362 setox IDNT 2,1 Motorola 040 Floating Point Software Package 363 364 section 8 365 366 include fpsp.h 367 368 L2 DC.L $3FDC0000,$82E30865,$4361C4C6,$00000000 369 370 EXPA3 DC.L $3FA55555,$55554431 371 EXPA2 DC.L $3FC55555,$55554018 372 373 HUGE DC.L $7FFE0000,$FFFFFFFF,$FFFFFFFF,$00000000 374 TINY DC.L $00010000,$FFFFFFFF,$FFFFFFFF,$00000000 375 376 EM1A4 DC.L $3F811111,$11174385 377 EM1A3 DC.L $3FA55555,$55554F5A 378 379 EM1A2 DC.L $3FC55555,$55555555,$00000000,$00000000 380 381 EM1B8 DC.L $3EC71DE3,$A5774682 382 EM1B7 DC.L $3EFA01A0,$19D7CB68 383 384 EM1B6 DC.L $3F2A01A0,$1A019DF3 385 EM1B5 DC.L $3F56C16C,$16C170E2 386 387 EM1B4 DC.L $3F811111,$11111111 388 EM1B3 DC.L $3FA55555,$55555555 389 390 EM1B2 DC.L $3FFC0000,$AAAAAAAA,$AAAAAAAB 391 DC.L $00000000 392 393 TWO140 DC.L $48B00000,$00000000 394 TWON140 DC.L $37300000,$00000000 395 396 EXPTBL 397 DC.L $3FFF0000,$80000000,$00000000,$00000000 398 DC.L $3FFF0000,$8164D1F3,$BC030774,$9F841A9B 399 DC.L $3FFF0000,$82CD8698,$AC2BA1D8,$9FC1D5B9 400 DC.L $3FFF0000,$843A28C3,$ACDE4048,$A0728369 401 DC.L $3FFF0000,$85AAC367,$CC487B14,$1FC5C95C 402 DC.L $3FFF0000,$871F6196,$9E8D1010,$1EE85C9F 403 DC.L $3FFF0000,$88980E80,$92DA8528,$9FA20729 404 DC.L $3FFF0000,$8A14D575,$496EFD9C,$A07BF9AF 405 DC.L $3FFF0000,$8B95C1E3,$EA8BD6E8,$A0020DCF 406 DC.L $3FFF0000,$8D1ADF5B,$7E5BA9E4,$205A63DA 407 DC.L $3FFF0000,$8EA4398B,$45CD53C0,$1EB70051 408 DC.L $3FFF0000,$9031DC43,$1466B1DC,$1F6EB029 409 DC.L $3FFF0000,$91C3D373,$AB11C338,$A0781494 410 DC.L $3FFF0000,$935A2B2F,$13E6E92C,$9EB319B0 411 DC.L $3FFF0000,$94F4EFA8,$FEF70960,$2017457D 412 DC.L $3FFF0000,$96942D37,$20185A00,$1F11D537 413 DC.L $3FFF0000,$9837F051,$8DB8A970,$9FB952DD 414 DC.L $3FFF0000,$99E04593,$20B7FA64,$1FE43087 415 DC.L $3FFF0000,$9B8D39B9,$D54E5538,$1FA2A818 416 DC.L $3FFF0000,$9D3ED9A7,$2CFFB750,$1FDE494D 417 DC.L $3FFF0000,$9EF53260,$91A111AC,$20504890 418 DC.L $3FFF0000,$A0B0510F,$B9714FC4,$A073691C 419 DC.L $3FFF0000,$A2704303,$0C496818,$1F9B7A05 420 DC.L $3FFF0000,$A43515AE,$09E680A0,$A0797126 421 DC.L $3FFF0000,$A5FED6A9,$B15138EC,$A071A140 422 DC.L $3FFF0000,$A7CD93B4,$E9653568,$204F62DA 423 DC.L $3FFF0000,$A9A15AB4,$EA7C0EF8,$1F283C4A 424 DC.L $3FFF0000,$AB7A39B5,$A93ED338,$9F9A7FDC 425 DC.L $3FFF0000,$AD583EEA,$42A14AC8,$A05B3FAC 426 DC.L $3FFF0000,$AF3B78AD,$690A4374,$1FDF2610 427 DC.L $3FFF0000,$B123F581,$D2AC2590,$9F705F90 428 DC.L $3FFF0000,$B311C412,$A9112488,$201F678A 429 DC.L $3FFF0000,$B504F333,$F9DE6484,$1F32FB13 430 DC.L $3FFF0000,$B6FD91E3,$28D17790,$20038B30 431 DC.L $3FFF0000,$B8FBAF47,$62FB9EE8,$200DC3CC 432 DC.L $3FFF0000,$BAFF5AB2,$133E45FC,$9F8B2AE6 433 DC.L $3FFF0000,$BD08A39F,$580C36C0,$A02BBF70 434 DC.L $3FFF0000,$BF1799B6,$7A731084,$A00BF518 435 DC.L $3FFF0000,$C12C4CCA,$66709458,$A041DD41 436 DC.L $3FFF0000,$C346CCDA,$24976408,$9FDF137B 437 DC.L $3FFF0000,$C5672A11,$5506DADC,$201F1568 438 DC.L $3FFF0000,$C78D74C8,$ABB9B15C,$1FC13A2E 439 DC.L $3FFF0000,$C9B9BD86,$6E2F27A4,$A03F8F03 440 DC.L $3FFF0000,$CBEC14FE,$F2727C5C,$1FF4907D 441 DC.L $3FFF0000,$CE248C15,$1F8480E4,$9E6E53E4 442 DC.L $3FFF0000,$D06333DA,$EF2B2594,$1FD6D45C 443 DC.L $3FFF0000,$D2A81D91,$F12AE45C,$A076EDB9 444 DC.L $3FFF0000,$D4F35AAB,$CFEDFA20,$9FA6DE21 445 DC.L $3FFF0000,$D744FCCA,$D69D6AF4,$1EE69A2F 446 DC.L $3FFF0000,$D99D15C2,$78AFD7B4,$207F439F 447 DC.L $3FFF0000,$DBFBB797,$DAF23754,$201EC207 448 DC.L $3FFF0000,$DE60F482,$5E0E9124,$9E8BE175 449 DC.L $3FFF0000,$E0CCDEEC,$2A94E110,$20032C4B 450 DC.L $3FFF0000,$E33F8972,$BE8A5A50,$2004DFF5 451 DC.L $3FFF0000,$E5B906E7,$7C8348A8,$1E72F47A 452 DC.L $3FFF0000,$E8396A50,$3C4BDC68,$1F722F22 453 DC.L $3FFF0000,$EAC0C6E7,$DD243930,$A017E945 454 DC.L $3FFF0000,$ED4F301E,$D9942B84,$1F401A5B 455 DC.L $3FFF0000,$EFE4B99B,$DCDAF5CC,$9FB9A9E3 456 DC.L $3FFF0000,$F281773C,$59FFB138,$20744C05 457 DC.L $3FFF0000,$F5257D15,$2486CC2C,$1F773A19 458 DC.L $3FFF0000,$F7D0DF73,$0AD13BB8,$1FFE90D5 459 DC.L $3FFF0000,$FA83B2DB,$722A033C,$A041ED22 460 DC.L $3FFF0000,$FD3E0C0C,$F486C174,$1F853F3A 461 462 ADJFLAG equ L_SCR2 463 SCALE equ FP_SCR1 464 ADJSCALE equ FP_SCR2 465 SC equ FP_SCR3 466 ONEBYSC equ FP_SCR4 467 468 xref t_frcinx 469 xref t_extdnrm 470 xref t_unfl 471 xref t_ovfl 472 473 xdef setoxd 474 setoxd: 475 *--entry point for EXP(X), X is denormalized 476 MOVE.L (a0),d0 477 ANDI.L #$80000000,d0 478 ORI.L #$00800000,d0 ...sign(X)*2^(-126) 479 MOVE.L d0,-(sp) 480 FMOVE.S #:3F800000,fp0 481 fmove.l d1,fpcr 482 FADD.S (sp)+,fp0 483 bra t_frcinx 484 485 xdef setox 486 setox: 487 *--entry point for EXP(X), here X is finite, non-zero, and not NaN's 488 489 *--Step 1. 490 MOVE.L (a0),d0 ...load part of input X 491 ANDI.L #$7FFF0000,d0 ...biased expo. of X 492 CMPI.L #$3FBE0000,d0 ...2^(-65) 493 BGE.B EXPC1 ...normal case 494 BRA.W EXPSM 495 496 EXPC1: 497 *--The case |X| >= 2^(-65) 498 MOVE.W 4(a0),d0 ...expo. and partial sig. of |X| 499 CMPI.L #$400CB167,d0 ...16380 log2 trunc. 16 bits 500 BLT.B EXPMAIN ...normal case 501 BRA.W EXPBIG 502 503 EXPMAIN: 504 *--Step 2. 505 *--This is the normal branch: 2^(-65) <= |X| < 16380 log2. 506 FMOVE.X (a0),fp0 ...load input from (a0) 507 508 FMOVE.X fp0,fp1 509 FMUL.S #:42B8AA3B,fp0 ...64/log2 * X 510 fmovem.x fp2/fp3,-(a7) ...save fp2 511 CLR.L ADJFLAG(a6) 512 FMOVE.L fp0,d0 ...N = int( X * 64/log2 ) 513 LEA EXPTBL,a1 514 FMOVE.L d0,fp0 ...convert to floating-format 515 516 MOVE.L d0,L_SCR1(a6) ...save N temporarily 517 ANDI.L #$3F,d0 ...D0 is J = N mod 64 518 LSL.L #4,d0 519 ADDA.L d0,a1 ...address of 2^(J/64) 520 MOVE.L L_SCR1(a6),d0 521 ASR.L #6,d0 ...D0 is M 522 ADDI.W #$3FFF,d0 ...biased expo. of 2^(M) 523 MOVE.W L2,L_SCR1(a6) ...prefetch L2, no need in CB 524 525 EXPCONT1: 526 *--Step 3. 527 *--fp1,fp2 saved on the stack. fp0 is N, fp1 is X, 528 *--a0 points to 2^(J/64), D0 is biased expo. of 2^(M) 529 FMOVE.X fp0,fp2 530 FMUL.S #:BC317218,fp0 ...N * L1, L1 = lead(-log2/64) 531 FMUL.X L2,fp2 ...N * L2, L1+L2 = -log2/64 532 FADD.X fp1,fp0 ...X + N*L1 533 FADD.X fp2,fp0 ...fp0 is R, reduced arg. 534 * MOVE.W #$3FA5,EXPA3 ...load EXPA3 in cache 535 536 *--Step 4. 537 *--WE NOW COMPUTE EXP(R)-1 BY A POLYNOMIAL 538 *-- R + R*R*(A1 + R*(A2 + R*(A3 + R*(A4 + R*A5)))) 539 *--TO FULLY USE THE PIPELINE, WE COMPUTE S = R*R 540 *--[R+R*S*(A2+S*A4)] + [S*(A1+S*(A3+S*A5))] 541 542 FMOVE.X fp0,fp1 543 FMUL.X fp1,fp1 ...fp1 IS S = R*R 544 545 FMOVE.S #:3AB60B70,fp2 ...fp2 IS A5 546 * CLR.W 2(a1) ...load 2^(J/64) in cache 547 548 FMUL.X fp1,fp2 ...fp2 IS S*A5 549 FMOVE.X fp1,fp3 550 FMUL.S #:3C088895,fp3 ...fp3 IS S*A4 551 552 FADD.D EXPA3,fp2 ...fp2 IS A3+S*A5 553 FADD.D EXPA2,fp3 ...fp3 IS A2+S*A4 554 555 FMUL.X fp1,fp2 ...fp2 IS S*(A3+S*A5) 556 MOVE.W d0,SCALE(a6) ...SCALE is 2^(M) in extended 557 clr.w SCALE+2(a6) 558 move.l #$80000000,SCALE+4(a6) 559 clr.l SCALE+8(a6) 560 561 FMUL.X fp1,fp3 ...fp3 IS S*(A2+S*A4) 562 563 FADD.S #:3F000000,fp2 ...fp2 IS A1+S*(A3+S*A5) 564 FMUL.X fp0,fp3 ...fp3 IS R*S*(A2+S*A4) 565 566 FMUL.X fp1,fp2 ...fp2 IS S*(A1+S*(A3+S*A5)) 567 FADD.X fp3,fp0 ...fp0 IS R+R*S*(A2+S*A4), 568 * ...fp3 released 569 570 FMOVE.X (a1)+,fp1 ...fp1 is lead. pt. of 2^(J/64) 571 FADD.X fp2,fp0 ...fp0 is EXP(R) - 1 572 * ...fp2 released 573 574 *--Step 5 575 *--final reconstruction process 576 *--EXP(X) = 2^M * ( 2^(J/64) + 2^(J/64)*(EXP(R)-1) ) 577 578 FMUL.X fp1,fp0 ...2^(J/64)*(Exp(R)-1) 579 fmovem.x (a7)+,fp2/fp3 ...fp2 restored 580 FADD.S (a1),fp0 ...accurate 2^(J/64) 581 582 FADD.X fp1,fp0 ...2^(J/64) + 2^(J/64)*... 583 MOVE.L ADJFLAG(a6),d0 584 585 *--Step 6 586 TST.L D0 587 BEQ.B NORMAL 588 ADJUST: 589 FMUL.X ADJSCALE(a6),fp0 590 NORMAL: 591 FMOVE.L d1,FPCR ...restore user FPCR 592 FMUL.X SCALE(a6),fp0 ...multiply 2^(M) 593 bra t_frcinx 594 595 EXPSM: 596 *--Step 7 597 FMOVEM.X (a0),fp0 ...in case X is denormalized 598 FMOVE.L d1,FPCR 599 FADD.S #:3F800000,fp0 ...1+X in user mode 600 bra t_frcinx 601 602 EXPBIG: 603 *--Step 8 604 CMPI.L #$400CB27C,d0 ...16480 log2 605 BGT.B EXP2BIG 606 *--Steps 8.2 -- 8.6 607 FMOVE.X (a0),fp0 ...load input from (a0) 608 609 FMOVE.X fp0,fp1 610 FMUL.S #:42B8AA3B,fp0 ...64/log2 * X 611 fmovem.x fp2/fp3,-(a7) ...save fp2 612 MOVE.L #1,ADJFLAG(a6) 613 FMOVE.L fp0,d0 ...N = int( X * 64/log2 ) 614 LEA EXPTBL,a1 615 FMOVE.L d0,fp0 ...convert to floating-format 616 MOVE.L d0,L_SCR1(a6) ...save N temporarily 617 ANDI.L #$3F,d0 ...D0 is J = N mod 64 618 LSL.L #4,d0 619 ADDA.L d0,a1 ...address of 2^(J/64) 620 MOVE.L L_SCR1(a6),d0 621 ASR.L #6,d0 ...D0 is K 622 MOVE.L d0,L_SCR1(a6) ...save K temporarily 623 ASR.L #1,d0 ...D0 is M1 624 SUB.L d0,L_SCR1(a6) ...a1 is M 625 ADDI.W #$3FFF,d0 ...biased expo. of 2^(M1) 626 MOVE.W d0,ADJSCALE(a6) ...ADJSCALE := 2^(M1) 627 clr.w ADJSCALE+2(a6) 628 move.l #$80000000,ADJSCALE+4(a6) 629 clr.l ADJSCALE+8(a6) 630 MOVE.L L_SCR1(a6),d0 ...D0 is M 631 ADDI.W #$3FFF,d0 ...biased expo. of 2^(M) 632 BRA.W EXPCONT1 ...go back to Step 3 633 634 EXP2BIG: 635 *--Step 9 636 FMOVE.L d1,FPCR 637 MOVE.L (a0),d0 638 bclr.b #sign_bit,(a0) ...setox always returns positive 639 TST.L d0 640 BLT t_unfl 641 BRA t_ovfl 642 643 xdef setoxm1d 644 setoxm1d: 645 *--entry point for EXPM1(X), here X is denormalized 646 *--Step 0. 647 bra t_extdnrm 648 649 650 xdef setoxm1 651 setoxm1: 652 *--entry point for EXPM1(X), here X is finite, non-zero, non-NaN 653 654 *--Step 1. 655 *--Step 1.1 656 MOVE.L (a0),d0 ...load part of input X 657 ANDI.L #$7FFF0000,d0 ...biased expo. of X 658 CMPI.L #$3FFD0000,d0 ...1/4 659 BGE.B EM1CON1 ...|X| >= 1/4 660 BRA.W EM1SM 661 662 EM1CON1: 663 *--Step 1.3 664 *--The case |X| >= 1/4 665 MOVE.W 4(a0),d0 ...expo. and partial sig. of |X| 666 CMPI.L #$4004C215,d0 ...70log2 rounded up to 16 bits 667 BLE.B EM1MAIN ...1/4 <= |X| <= 70log2 668 BRA.W EM1BIG 669 670 EM1MAIN: 671 *--Step 2. 672 *--This is the case: 1/4 <= |X| <= 70 log2. 673 FMOVE.X (a0),fp0 ...load input from (a0) 674 675 FMOVE.X fp0,fp1 676 FMUL.S #:42B8AA3B,fp0 ...64/log2 * X 677 fmovem.x fp2/fp3,-(a7) ...save fp2 678 * MOVE.W #$3F81,EM1A4 ...prefetch in CB mode 679 FMOVE.L fp0,d0 ...N = int( X * 64/log2 ) 680 LEA EXPTBL,a1 681 FMOVE.L d0,fp0 ...convert to floating-format 682 683 MOVE.L d0,L_SCR1(a6) ...save N temporarily 684 ANDI.L #$3F,d0 ...D0 is J = N mod 64 685 LSL.L #4,d0 686 ADDA.L d0,a1 ...address of 2^(J/64) 687 MOVE.L L_SCR1(a6),d0 688 ASR.L #6,d0 ...D0 is M 689 MOVE.L d0,L_SCR1(a6) ...save a copy of M 690 * MOVE.W #$3FDC,L2 ...prefetch L2 in CB mode 691 692 *--Step 3. 693 *--fp1,fp2 saved on the stack. fp0 is N, fp1 is X, 694 *--a0 points to 2^(J/64), D0 and a1 both contain M 695 FMOVE.X fp0,fp2 696 FMUL.S #:BC317218,fp0 ...N * L1, L1 = lead(-log2/64) 697 FMUL.X L2,fp2 ...N * L2, L1+L2 = -log2/64 698 FADD.X fp1,fp0 ...X + N*L1 699 FADD.X fp2,fp0 ...fp0 is R, reduced arg. 700 * MOVE.W #$3FC5,EM1A2 ...load EM1A2 in cache 701 ADDI.W #$3FFF,d0 ...D0 is biased expo. of 2^M 702 703 *--Step 4. 704 *--WE NOW COMPUTE EXP(R)-1 BY A POLYNOMIAL 705 *-- R + R*R*(A1 + R*(A2 + R*(A3 + R*(A4 + R*(A5 + R*A6))))) 706 *--TO FULLY USE THE PIPELINE, WE COMPUTE S = R*R 707 *--[R*S*(A2+S*(A4+S*A6))] + [R+S*(A1+S*(A3+S*A5))] 708 709 FMOVE.X fp0,fp1 710 FMUL.X fp1,fp1 ...fp1 IS S = R*R 711 712 FMOVE.S #:3950097B,fp2 ...fp2 IS a6 713 * CLR.W 2(a1) ...load 2^(J/64) in cache 714 715 FMUL.X fp1,fp2 ...fp2 IS S*A6 716 FMOVE.X fp1,fp3 717 FMUL.S #:3AB60B6A,fp3 ...fp3 IS S*A5 718 719 FADD.D EM1A4,fp2 ...fp2 IS A4+S*A6 720 FADD.D EM1A3,fp3 ...fp3 IS A3+S*A5 721 MOVE.W d0,SC(a6) ...SC is 2^(M) in extended 722 clr.w SC+2(a6) 723 move.l #$80000000,SC+4(a6) 724 clr.l SC+8(a6) 725 726 FMUL.X fp1,fp2 ...fp2 IS S*(A4+S*A6) 727 MOVE.L L_SCR1(a6),d0 ...D0 is M 728 NEG.W D0 ...D0 is -M 729 FMUL.X fp1,fp3 ...fp3 IS S*(A3+S*A5) 730 ADDI.W #$3FFF,d0 ...biased expo. of 2^(-M) 731 FADD.D EM1A2,fp2 ...fp2 IS A2+S*(A4+S*A6) 732 FADD.S #:3F000000,fp3 ...fp3 IS A1+S*(A3+S*A5) 733 734 FMUL.X fp1,fp2 ...fp2 IS S*(A2+S*(A4+S*A6)) 735 ORI.W #$8000,d0 ...signed/expo. of -2^(-M) 736 MOVE.W d0,ONEBYSC(a6) ...OnebySc is -2^(-M) 737 clr.w ONEBYSC+2(a6) 738 move.l #$80000000,ONEBYSC+4(a6) 739 clr.l ONEBYSC+8(a6) 740 FMUL.X fp3,fp1 ...fp1 IS S*(A1+S*(A3+S*A5)) 741 * ...fp3 released 742 743 FMUL.X fp0,fp2 ...fp2 IS R*S*(A2+S*(A4+S*A6)) 744 FADD.X fp1,fp0 ...fp0 IS R+S*(A1+S*(A3+S*A5)) 745 * ...fp1 released 746 747 FADD.X fp2,fp0 ...fp0 IS EXP(R)-1 748 * ...fp2 released 749 fmovem.x (a7)+,fp2/fp3 ...fp2 restored 750 751 *--Step 5 752 *--Compute 2^(J/64)*p 753 754 FMUL.X (a1),fp0 ...2^(J/64)*(Exp(R)-1) 755 756 *--Step 6 757 *--Step 6.1 758 MOVE.L L_SCR1(a6),d0 ...retrieve M 759 CMPI.L #63,d0 760 BLE.B MLE63 761 *--Step 6.2 M >= 64 762 FMOVE.S 12(a1),fp1 ...fp1 is t 763 FADD.X ONEBYSC(a6),fp1 ...fp1 is t+OnebySc 764 FADD.X fp1,fp0 ...p+(t+OnebySc), fp1 released 765 FADD.X (a1),fp0 ...T+(p+(t+OnebySc)) 766 BRA.B EM1SCALE 767 MLE63: 768 *--Step 6.3 M <= 63 769 CMPI.L #-3,d0 770 BGE.B MGEN3 771 MLTN3: 772 *--Step 6.4 M <= -4 773 FADD.S 12(a1),fp0 ...p+t 774 FADD.X (a1),fp0 ...T+(p+t) 775 FADD.X ONEBYSC(a6),fp0 ...OnebySc + (T+(p+t)) 776 BRA.B EM1SCALE 777 MGEN3: 778 *--Step 6.5 -3 <= M <= 63 779 FMOVE.X (a1)+,fp1 ...fp1 is T 780 FADD.S (a1),fp0 ...fp0 is p+t 781 FADD.X ONEBYSC(a6),fp1 ...fp1 is T+OnebySc 782 FADD.X fp1,fp0 ...(T+OnebySc)+(p+t) 783 784 EM1SCALE: 785 *--Step 6.6 786 FMOVE.L d1,FPCR 787 FMUL.X SC(a6),fp0 788 789 bra t_frcinx 790 791 EM1SM: 792 *--Step 7 |X| < 1/4. 793 CMPI.L #$3FBE0000,d0 ...2^(-65) 794 BGE.B EM1POLY 795 796 EM1TINY: 797 *--Step 8 |X| < 2^(-65) 798 CMPI.L #$00330000,d0 ...2^(-16312) 799 BLT.B EM12TINY 800 *--Step 8.2 801 MOVE.L #$80010000,SC(a6) ...SC is -2^(-16382) 802 move.l #$80000000,SC+4(a6) 803 clr.l SC+8(a6) 804 FMOVE.X (a0),fp0 805 FMOVE.L d1,FPCR 806 FADD.X SC(a6),fp0 807 808 bra t_frcinx 809 810 EM12TINY: 811 *--Step 8.3 812 FMOVE.X (a0),fp0 813 FMUL.D TWO140,fp0 814 MOVE.L #$80010000,SC(a6) 815 move.l #$80000000,SC+4(a6) 816 clr.l SC+8(a6) 817 FADD.X SC(a6),fp0 818 FMOVE.L d1,FPCR 819 FMUL.D TWON140,fp0 820 821 bra t_frcinx 822 823 EM1POLY: 824 *--Step 9 exp(X)-1 by a simple polynomial 825 FMOVE.X (a0),fp0 ...fp0 is X 826 FMUL.X fp0,fp0 ...fp0 is S := X*X 827 fmovem.x fp2/fp3,-(a7) ...save fp2 828 FMOVE.S #:2F30CAA8,fp1 ...fp1 is B12 829 FMUL.X fp0,fp1 ...fp1 is S*B12 830 FMOVE.S #:310F8290,fp2 ...fp2 is B11 831 FADD.S #:32D73220,fp1 ...fp1 is B10+S*B12 832 833 FMUL.X fp0,fp2 ...fp2 is S*B11 834 FMUL.X fp0,fp1 ...fp1 is S*(B10 + ... 835 836 FADD.S #:3493F281,fp2 ...fp2 is B9+S*... 837 FADD.D EM1B8,fp1 ...fp1 is B8+S*... 838 839 FMUL.X fp0,fp2 ...fp2 is S*(B9+... 840 FMUL.X fp0,fp1 ...fp1 is S*(B8+... 841 842 FADD.D EM1B7,fp2 ...fp2 is B7+S*... 843 FADD.D EM1B6,fp1 ...fp1 is B6+S*... 844 845 FMUL.X fp0,fp2 ...fp2 is S*(B7+... 846 FMUL.X fp0,fp1 ...fp1 is S*(B6+... 847 848 FADD.D EM1B5,fp2 ...fp2 is B5+S*... 849 FADD.D EM1B4,fp1 ...fp1 is B4+S*... 850 851 FMUL.X fp0,fp2 ...fp2 is S*(B5+... 852 FMUL.X fp0,fp1 ...fp1 is S*(B4+... 853 854 FADD.D EM1B3,fp2 ...fp2 is B3+S*... 855 FADD.X EM1B2,fp1 ...fp1 is B2+S*... 856 857 FMUL.X fp0,fp2 ...fp2 is S*(B3+... 858 FMUL.X fp0,fp1 ...fp1 is S*(B2+... 859 860 FMUL.X fp0,fp2 ...fp2 is S*S*(B3+...) 861 FMUL.X (a0),fp1 ...fp1 is X*S*(B2... 862 863 FMUL.S #:3F000000,fp0 ...fp0 is S*B1 864 FADD.X fp2,fp1 ...fp1 is Q 865 * ...fp2 released 866 867 fmovem.x (a7)+,fp2/fp3 ...fp2 restored 868 869 FADD.X fp1,fp0 ...fp0 is S*B1+Q 870 * ...fp1 released 871 872 FMOVE.L d1,FPCR 873 FADD.X (a0),fp0 874 875 bra t_frcinx 876 877 EM1BIG: 878 *--Step 10 |X| > 70 log2 879 MOVE.L (a0),d0 880 TST.L d0 881 BGT.W EXPC1 882 *--Step 10.2 883 FMOVE.S #:BF800000,fp0 ...fp0 is -1 884 FMOVE.L d1,FPCR 885 FADD.S #:00800000,fp0 ...-1 + 2^(-126) 886 887 bra t_frcinx 888 889 end 890
Indexes created Sun Oct 12 05:10:08 GMT 2025