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