^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 1) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 2) | setox.sa 3.1 12/10/90
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 3) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 4) | The entry point setox computes the exponential of a value.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 5) | setoxd does the same except the input value is a denormalized
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 6) | number. setoxm1 computes exp(X)-1, and setoxm1d computes
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 7) | exp(X)-1 for denormalized X.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 8) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 9) | INPUT
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 10) | -----
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 11) | Double-extended value in memory location pointed to by address
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 12) | register a0.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 13) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 14) | OUTPUT
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 15) | ------
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 16) | exp(X) or exp(X)-1 returned in floating-point register fp0.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 17) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 18) | ACCURACY and MONOTONICITY
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 19) | -------------------------
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 20) | The returned result is within 0.85 ulps in 64 significant bit, i.e.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 21) | within 0.5001 ulp to 53 bits if the result is subsequently rounded
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 22) | to double precision. The result is provably monotonic in double
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 23) | precision.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 24) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 25) | SPEED
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 26) | -----
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 27) | Two timings are measured, both in the copy-back mode. The
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 28) | first one is measured when the function is invoked the first time
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 29) | (so the instructions and data are not in cache), and the
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 30) | second one is measured when the function is reinvoked at the same
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 31) | input argument.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 32) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 33) | The program setox takes approximately 210/190 cycles for input
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 34) | argument X whose magnitude is less than 16380 log2, which
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 35) | is the usual situation. For the less common arguments,
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 36) | depending on their values, the program may run faster or slower --
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 37) | but no worse than 10% slower even in the extreme cases.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 38) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 39) | The program setoxm1 takes approximately ??? / ??? cycles for input
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 40) | argument X, 0.25 <= |X| < 70log2. For |X| < 0.25, it takes
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 41) | approximately ??? / ??? cycles. For the less common arguments,
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 42) | depending on their values, the program may run faster or slower --
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 43) | but no worse than 10% slower even in the extreme cases.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 44) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 45) | ALGORITHM and IMPLEMENTATION NOTES
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 46) | ----------------------------------
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 47) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 48) | setoxd
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 49) | ------
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 50) | Step 1. Set ans := 1.0
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 51) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 52) | Step 2. Return ans := ans + sign(X)*2^(-126). Exit.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 53) | Notes: This will always generate one exception -- inexact.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 54) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 55) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 56) | setox
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 57) | -----
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 58) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 59) | Step 1. Filter out extreme cases of input argument.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 60) | 1.1 If |X| >= 2^(-65), go to Step 1.3.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 61) | 1.2 Go to Step 7.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 62) | 1.3 If |X| < 16380 log(2), go to Step 2.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 63) | 1.4 Go to Step 8.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 64) | Notes: The usual case should take the branches 1.1 -> 1.3 -> 2.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 65) | To avoid the use of floating-point comparisons, a
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 66) | compact representation of |X| is used. This format is a
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 67) | 32-bit integer, the upper (more significant) 16 bits are
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 68) | the sign and biased exponent field of |X|; the lower 16
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 69) | bits are the 16 most significant fraction (including the
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 70) | explicit bit) bits of |X|. Consequently, the comparisons
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 71) | in Steps 1.1 and 1.3 can be performed by integer comparison.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 72) | Note also that the constant 16380 log(2) used in Step 1.3
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 73) | is also in the compact form. Thus taking the branch
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 74) | to Step 2 guarantees |X| < 16380 log(2). There is no harm
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 75) | to have a small number of cases where |X| is less than,
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 76) | but close to, 16380 log(2) and the branch to Step 9 is
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 77) | taken.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 78) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 79) | Step 2. Calculate N = round-to-nearest-int( X * 64/log2 ).
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 80) | 2.1 Set AdjFlag := 0 (indicates the branch 1.3 -> 2 was taken)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 81) | 2.2 N := round-to-nearest-integer( X * 64/log2 ).
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 82) | 2.3 Calculate J = N mod 64; so J = 0,1,2,..., or 63.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 83) | 2.4 Calculate M = (N - J)/64; so N = 64M + J.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 84) | 2.5 Calculate the address of the stored value of 2^(J/64).
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 85) | 2.6 Create the value Scale = 2^M.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 86) | Notes: The calculation in 2.2 is really performed by
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 87) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 88) | Z := X * constant
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 89) | N := round-to-nearest-integer(Z)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 90) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 91) | where
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 92) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 93) | constant := single-precision( 64/log 2 ).
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 94) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 95) | Using a single-precision constant avoids memory access.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 96) | Another effect of using a single-precision "constant" is
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 97) | that the calculated value Z is
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 98) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 99) | Z = X*(64/log2)*(1+eps), |eps| <= 2^(-24).
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 100) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 101) | This error has to be considered later in Steps 3 and 4.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 102) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 103) | Step 3. Calculate X - N*log2/64.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 104) | 3.1 R := X + N*L1, where L1 := single-precision(-log2/64).
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 105) | 3.2 R := R + N*L2, L2 := extended-precision(-log2/64 - L1).
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 106) | Notes: a) The way L1 and L2 are chosen ensures L1+L2 approximate
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 107) | the value -log2/64 to 88 bits of accuracy.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 108) | b) N*L1 is exact because N is no longer than 22 bits and
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 109) | L1 is no longer than 24 bits.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 110) | c) The calculation X+N*L1 is also exact due to cancellation.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 111) | Thus, R is practically X+N(L1+L2) to full 64 bits.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 112) | d) It is important to estimate how large can |R| be after
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 113) | Step 3.2.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 114) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 115) | N = rnd-to-int( X*64/log2 (1+eps) ), |eps|<=2^(-24)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 116) | X*64/log2 (1+eps) = N + f, |f| <= 0.5
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 117) | X*64/log2 - N = f - eps*X 64/log2
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 118) | X - N*log2/64 = f*log2/64 - eps*X
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 119) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 120) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 121) | Now |X| <= 16446 log2, thus
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 122) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 123) | |X - N*log2/64| <= (0.5 + 16446/2^(18))*log2/64
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 124) | <= 0.57 log2/64.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 125) | This bound will be used in Step 4.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 126) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 127) | Step 4. Approximate exp(R)-1 by a polynomial
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 128) | p = R + R*R*(A1 + R*(A2 + R*(A3 + R*(A4 + R*A5))))
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 129) | Notes: a) In order to reduce memory access, the coefficients are
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 130) | made as "short" as possible: A1 (which is 1/2), A4 and A5
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 131) | are single precision; A2 and A3 are double precision.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 132) | b) Even with the restrictions above,
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 133) | |p - (exp(R)-1)| < 2^(-68.8) for all |R| <= 0.0062.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 134) | Note that 0.0062 is slightly bigger than 0.57 log2/64.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 135) | c) To fully utilize the pipeline, p is separated into
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 136) | two independent pieces of roughly equal complexities
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 137) | p = [ R + R*S*(A2 + S*A4) ] +
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 138) | [ S*(A1 + S*(A3 + S*A5)) ]
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 139) | where S = R*R.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 140) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 141) | Step 5. Compute 2^(J/64)*exp(R) = 2^(J/64)*(1+p) by
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 142) | ans := T + ( T*p + t)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 143) | where T and t are the stored values for 2^(J/64).
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 144) | Notes: 2^(J/64) is stored as T and t where T+t approximates
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 145) | 2^(J/64) to roughly 85 bits; T is in extended precision
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 146) | and t is in single precision. Note also that T is rounded
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 147) | to 62 bits so that the last two bits of T are zero. The
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 148) | reason for such a special form is that T-1, T-2, and T-8
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 149) | will all be exact --- a property that will give much
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 150) | more accurate computation of the function EXPM1.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 151) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 152) | Step 6. Reconstruction of exp(X)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 153) | exp(X) = 2^M * 2^(J/64) * exp(R).
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 154) | 6.1 If AdjFlag = 0, go to 6.3
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 155) | 6.2 ans := ans * AdjScale
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 156) | 6.3 Restore the user FPCR
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 157) | 6.4 Return ans := ans * Scale. Exit.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 158) | Notes: If AdjFlag = 0, we have X = Mlog2 + Jlog2/64 + R,
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 159) | |M| <= 16380, and Scale = 2^M. Moreover, exp(X) will
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 160) | neither overflow nor underflow. If AdjFlag = 1, that
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 161) | means that
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 162) | X = (M1+M)log2 + Jlog2/64 + R, |M1+M| >= 16380.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 163) | Hence, exp(X) may overflow or underflow or neither.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 164) | When that is the case, AdjScale = 2^(M1) where M1 is
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 165) | approximately M. Thus 6.2 will never cause over/underflow.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 166) | Possible exception in 6.4 is overflow or underflow.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 167) | The inexact exception is not generated in 6.4. Although
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 168) | one can argue that the inexact flag should always be
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 169) | raised, to simulate that exception cost to much than the
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 170) | flag is worth in practical uses.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 171) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 172) | Step 7. Return 1 + X.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 173) | 7.1 ans := X
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 174) | 7.2 Restore user FPCR.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 175) | 7.3 Return ans := 1 + ans. Exit
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 176) | Notes: For non-zero X, the inexact exception will always be
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 177) | raised by 7.3. That is the only exception raised by 7.3.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 178) | Note also that we use the FMOVEM instruction to move X
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 179) | in Step 7.1 to avoid unnecessary trapping. (Although
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 180) | the FMOVEM may not seem relevant since X is normalized,
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 181) | the precaution will be useful in the library version of
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 182) | this code where the separate entry for denormalized inputs
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 183) | will be done away with.)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 184) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 185) | Step 8. Handle exp(X) where |X| >= 16380log2.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 186) | 8.1 If |X| > 16480 log2, go to Step 9.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 187) | (mimic 2.2 - 2.6)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 188) | 8.2 N := round-to-integer( X * 64/log2 )
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 189) | 8.3 Calculate J = N mod 64, J = 0,1,...,63
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 190) | 8.4 K := (N-J)/64, M1 := truncate(K/2), M = K-M1, AdjFlag := 1.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 191) | 8.5 Calculate the address of the stored value 2^(J/64).
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 192) | 8.6 Create the values Scale = 2^M, AdjScale = 2^M1.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 193) | 8.7 Go to Step 3.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 194) | Notes: Refer to notes for 2.2 - 2.6.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 195) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 196) | Step 9. Handle exp(X), |X| > 16480 log2.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 197) | 9.1 If X < 0, go to 9.3
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 198) | 9.2 ans := Huge, go to 9.4
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 199) | 9.3 ans := Tiny.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 200) | 9.4 Restore user FPCR.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 201) | 9.5 Return ans := ans * ans. Exit.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 202) | Notes: Exp(X) will surely overflow or underflow, depending on
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 203) | X's sign. "Huge" and "Tiny" are respectively large/tiny
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 204) | extended-precision numbers whose square over/underflow
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 205) | with an inexact result. Thus, 9.5 always raises the
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 206) | inexact together with either overflow or underflow.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 207) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 208) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 209) | setoxm1d
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 210) | --------
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 211) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 212) | Step 1. Set ans := 0
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 213) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 214) | Step 2. Return ans := X + ans. Exit.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 215) | Notes: This will return X with the appropriate rounding
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 216) | precision prescribed by the user FPCR.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 217) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 218) | setoxm1
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 219) | -------
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 220) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 221) | Step 1. Check |X|
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 222) | 1.1 If |X| >= 1/4, go to Step 1.3.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 223) | 1.2 Go to Step 7.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 224) | 1.3 If |X| < 70 log(2), go to Step 2.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 225) | 1.4 Go to Step 10.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 226) | Notes: The usual case should take the branches 1.1 -> 1.3 -> 2.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 227) | However, it is conceivable |X| can be small very often
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 228) | because EXPM1 is intended to evaluate exp(X)-1 accurately
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 229) | when |X| is small. For further details on the comparisons,
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 230) | see the notes on Step 1 of setox.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 231) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 232) | Step 2. Calculate N = round-to-nearest-int( X * 64/log2 ).
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 233) | 2.1 N := round-to-nearest-integer( X * 64/log2 ).
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 234) | 2.2 Calculate J = N mod 64; so J = 0,1,2,..., or 63.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 235) | 2.3 Calculate M = (N - J)/64; so N = 64M + J.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 236) | 2.4 Calculate the address of the stored value of 2^(J/64).
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 237) | 2.5 Create the values Sc = 2^M and OnebySc := -2^(-M).
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 238) | Notes: See the notes on Step 2 of setox.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 239) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 240) | Step 3. Calculate X - N*log2/64.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 241) | 3.1 R := X + N*L1, where L1 := single-precision(-log2/64).
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 242) | 3.2 R := R + N*L2, L2 := extended-precision(-log2/64 - L1).
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 243) | Notes: Applying the analysis of Step 3 of setox in this case
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 244) | shows that |R| <= 0.0055 (note that |X| <= 70 log2 in
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 245) | this case).
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 246) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 247) | Step 4. Approximate exp(R)-1 by a polynomial
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 248) | p = R+R*R*(A1+R*(A2+R*(A3+R*(A4+R*(A5+R*A6)))))
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 249) | Notes: a) In order to reduce memory access, the coefficients are
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 250) | made as "short" as possible: A1 (which is 1/2), A5 and A6
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 251) | are single precision; A2, A3 and A4 are double precision.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 252) | b) Even with the restriction above,
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 253) | |p - (exp(R)-1)| < |R| * 2^(-72.7)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 254) | for all |R| <= 0.0055.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 255) | c) To fully utilize the pipeline, p is separated into
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 256) | two independent pieces of roughly equal complexity
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 257) | p = [ R*S*(A2 + S*(A4 + S*A6)) ] +
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 258) | [ R + S*(A1 + S*(A3 + S*A5)) ]
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 259) | where S = R*R.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 260) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 261) | Step 5. Compute 2^(J/64)*p by
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 262) | p := T*p
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 263) | where T and t are the stored values for 2^(J/64).
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 264) | Notes: 2^(J/64) is stored as T and t where T+t approximates
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 265) | 2^(J/64) to roughly 85 bits; T is in extended precision
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 266) | and t is in single precision. Note also that T is rounded
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 267) | to 62 bits so that the last two bits of T are zero. The
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 268) | reason for such a special form is that T-1, T-2, and T-8
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 269) | will all be exact --- a property that will be exploited
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 270) | in Step 6 below. The total relative error in p is no
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 271) | bigger than 2^(-67.7) compared to the final result.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 272) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 273) | Step 6. Reconstruction of exp(X)-1
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 274) | exp(X)-1 = 2^M * ( 2^(J/64) + p - 2^(-M) ).
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 275) | 6.1 If M <= 63, go to Step 6.3.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 276) | 6.2 ans := T + (p + (t + OnebySc)). Go to 6.6
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 277) | 6.3 If M >= -3, go to 6.5.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 278) | 6.4 ans := (T + (p + t)) + OnebySc. Go to 6.6
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 279) | 6.5 ans := (T + OnebySc) + (p + t).
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 280) | 6.6 Restore user FPCR.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 281) | 6.7 Return ans := Sc * ans. Exit.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 282) | Notes: The various arrangements of the expressions give accurate
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 283) | evaluations.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 284) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 285) | Step 7. exp(X)-1 for |X| < 1/4.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 286) | 7.1 If |X| >= 2^(-65), go to Step 9.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 287) | 7.2 Go to Step 8.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 288) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 289) | Step 8. Calculate exp(X)-1, |X| < 2^(-65).
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 290) | 8.1 If |X| < 2^(-16312), goto 8.3
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 291) | 8.2 Restore FPCR; return ans := X - 2^(-16382). Exit.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 292) | 8.3 X := X * 2^(140).
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 293) | 8.4 Restore FPCR; ans := ans - 2^(-16382).
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 294) | Return ans := ans*2^(140). Exit
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 295) | Notes: The idea is to return "X - tiny" under the user
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 296) | precision and rounding modes. To avoid unnecessary
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 297) | inefficiency, we stay away from denormalized numbers the
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 298) | best we can. For |X| >= 2^(-16312), the straightforward
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 299) | 8.2 generates the inexact exception as the case warrants.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 300) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 301) | Step 9. Calculate exp(X)-1, |X| < 1/4, by a polynomial
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 302) | p = X + X*X*(B1 + X*(B2 + ... + X*B12))
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 303) | Notes: a) In order to reduce memory access, the coefficients are
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 304) | made as "short" as possible: B1 (which is 1/2), B9 to B12
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 305) | are single precision; B3 to B8 are double precision; and
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 306) | B2 is double extended.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 307) | b) Even with the restriction above,
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 308) | |p - (exp(X)-1)| < |X| 2^(-70.6)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 309) | for all |X| <= 0.251.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 310) | Note that 0.251 is slightly bigger than 1/4.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 311) | c) To fully preserve accuracy, the polynomial is computed
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 312) | as X + ( S*B1 + Q ) where S = X*X and
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 313) | Q = X*S*(B2 + X*(B3 + ... + X*B12))
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 314) | d) To fully utilize the pipeline, Q is separated into
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 315) | two independent pieces of roughly equal complexity
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 316) | Q = [ X*S*(B2 + S*(B4 + ... + S*B12)) ] +
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 317) | [ S*S*(B3 + S*(B5 + ... + S*B11)) ]
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 318) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 319) | Step 10. Calculate exp(X)-1 for |X| >= 70 log 2.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 320) | 10.1 If X >= 70log2 , exp(X) - 1 = exp(X) for all practical
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 321) | purposes. Therefore, go to Step 1 of setox.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 322) | 10.2 If X <= -70log2, exp(X) - 1 = -1 for all practical purposes.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 323) | ans := -1
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 324) | Restore user FPCR
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 325) | Return ans := ans + 2^(-126). Exit.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 326) | Notes: 10.2 will always create an inexact and return -1 + tiny
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 327) | in the user rounding precision and mode.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 328) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 329) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 330)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 331) | Copyright (C) Motorola, Inc. 1990
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 332) | All Rights Reserved
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 333) |
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 334) | For details on the license for this file, please see the
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 335) | file, README, in this same directory.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 336)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 337) |setox idnt 2,1 | Motorola 040 Floating Point Software Package
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 338)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 339) |section 8
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 340)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 341) #include "fpsp.h"
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 342)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 343) L2: .long 0x3FDC0000,0x82E30865,0x4361C4C6,0x00000000
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 344)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 345) EXPA3: .long 0x3FA55555,0x55554431
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 346) EXPA2: .long 0x3FC55555,0x55554018
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 347)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 348) HUGE: .long 0x7FFE0000,0xFFFFFFFF,0xFFFFFFFF,0x00000000
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 349) TINY: .long 0x00010000,0xFFFFFFFF,0xFFFFFFFF,0x00000000
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 350)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 351) EM1A4: .long 0x3F811111,0x11174385
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 352) EM1A3: .long 0x3FA55555,0x55554F5A
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 353)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 354) EM1A2: .long 0x3FC55555,0x55555555,0x00000000,0x00000000
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 355)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 356) EM1B8: .long 0x3EC71DE3,0xA5774682
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 357) EM1B7: .long 0x3EFA01A0,0x19D7CB68
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 358)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 359) EM1B6: .long 0x3F2A01A0,0x1A019DF3
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 360) EM1B5: .long 0x3F56C16C,0x16C170E2
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 361)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 362) EM1B4: .long 0x3F811111,0x11111111
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 363) EM1B3: .long 0x3FA55555,0x55555555
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 364)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 365) EM1B2: .long 0x3FFC0000,0xAAAAAAAA,0xAAAAAAAB
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 366) .long 0x00000000
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 367)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 368) TWO140: .long 0x48B00000,0x00000000
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 369) TWON140: .long 0x37300000,0x00000000
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 370)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 371) EXPTBL:
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 372) .long 0x3FFF0000,0x80000000,0x00000000,0x00000000
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 373) .long 0x3FFF0000,0x8164D1F3,0xBC030774,0x9F841A9B
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 374) .long 0x3FFF0000,0x82CD8698,0xAC2BA1D8,0x9FC1D5B9
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 375) .long 0x3FFF0000,0x843A28C3,0xACDE4048,0xA0728369
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 376) .long 0x3FFF0000,0x85AAC367,0xCC487B14,0x1FC5C95C
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 377) .long 0x3FFF0000,0x871F6196,0x9E8D1010,0x1EE85C9F
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 378) .long 0x3FFF0000,0x88980E80,0x92DA8528,0x9FA20729
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 379) .long 0x3FFF0000,0x8A14D575,0x496EFD9C,0xA07BF9AF
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 380) .long 0x3FFF0000,0x8B95C1E3,0xEA8BD6E8,0xA0020DCF
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 381) .long 0x3FFF0000,0x8D1ADF5B,0x7E5BA9E4,0x205A63DA
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 382) .long 0x3FFF0000,0x8EA4398B,0x45CD53C0,0x1EB70051
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 383) .long 0x3FFF0000,0x9031DC43,0x1466B1DC,0x1F6EB029
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 384) .long 0x3FFF0000,0x91C3D373,0xAB11C338,0xA0781494
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 385) .long 0x3FFF0000,0x935A2B2F,0x13E6E92C,0x9EB319B0
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 386) .long 0x3FFF0000,0x94F4EFA8,0xFEF70960,0x2017457D
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 387) .long 0x3FFF0000,0x96942D37,0x20185A00,0x1F11D537
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 388) .long 0x3FFF0000,0x9837F051,0x8DB8A970,0x9FB952DD
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 389) .long 0x3FFF0000,0x99E04593,0x20B7FA64,0x1FE43087
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 390) .long 0x3FFF0000,0x9B8D39B9,0xD54E5538,0x1FA2A818
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 391) .long 0x3FFF0000,0x9D3ED9A7,0x2CFFB750,0x1FDE494D
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 392) .long 0x3FFF0000,0x9EF53260,0x91A111AC,0x20504890
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 393) .long 0x3FFF0000,0xA0B0510F,0xB9714FC4,0xA073691C
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 394) .long 0x3FFF0000,0xA2704303,0x0C496818,0x1F9B7A05
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 395) .long 0x3FFF0000,0xA43515AE,0x09E680A0,0xA0797126
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 396) .long 0x3FFF0000,0xA5FED6A9,0xB15138EC,0xA071A140
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 397) .long 0x3FFF0000,0xA7CD93B4,0xE9653568,0x204F62DA
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 398) .long 0x3FFF0000,0xA9A15AB4,0xEA7C0EF8,0x1F283C4A
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 399) .long 0x3FFF0000,0xAB7A39B5,0xA93ED338,0x9F9A7FDC
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 400) .long 0x3FFF0000,0xAD583EEA,0x42A14AC8,0xA05B3FAC
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 401) .long 0x3FFF0000,0xAF3B78AD,0x690A4374,0x1FDF2610
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 402) .long 0x3FFF0000,0xB123F581,0xD2AC2590,0x9F705F90
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 403) .long 0x3FFF0000,0xB311C412,0xA9112488,0x201F678A
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 404) .long 0x3FFF0000,0xB504F333,0xF9DE6484,0x1F32FB13
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 405) .long 0x3FFF0000,0xB6FD91E3,0x28D17790,0x20038B30
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 406) .long 0x3FFF0000,0xB8FBAF47,0x62FB9EE8,0x200DC3CC
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 407) .long 0x3FFF0000,0xBAFF5AB2,0x133E45FC,0x9F8B2AE6
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 408) .long 0x3FFF0000,0xBD08A39F,0x580C36C0,0xA02BBF70
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 409) .long 0x3FFF0000,0xBF1799B6,0x7A731084,0xA00BF518
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 410) .long 0x3FFF0000,0xC12C4CCA,0x66709458,0xA041DD41
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 411) .long 0x3FFF0000,0xC346CCDA,0x24976408,0x9FDF137B
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 412) .long 0x3FFF0000,0xC5672A11,0x5506DADC,0x201F1568
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 413) .long 0x3FFF0000,0xC78D74C8,0xABB9B15C,0x1FC13A2E
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 414) .long 0x3FFF0000,0xC9B9BD86,0x6E2F27A4,0xA03F8F03
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 415) .long 0x3FFF0000,0xCBEC14FE,0xF2727C5C,0x1FF4907D
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 416) .long 0x3FFF0000,0xCE248C15,0x1F8480E4,0x9E6E53E4
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 417) .long 0x3FFF0000,0xD06333DA,0xEF2B2594,0x1FD6D45C
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 418) .long 0x3FFF0000,0xD2A81D91,0xF12AE45C,0xA076EDB9
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 419) .long 0x3FFF0000,0xD4F35AAB,0xCFEDFA20,0x9FA6DE21
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 420) .long 0x3FFF0000,0xD744FCCA,0xD69D6AF4,0x1EE69A2F
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 421) .long 0x3FFF0000,0xD99D15C2,0x78AFD7B4,0x207F439F
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 422) .long 0x3FFF0000,0xDBFBB797,0xDAF23754,0x201EC207
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 423) .long 0x3FFF0000,0xDE60F482,0x5E0E9124,0x9E8BE175
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 424) .long 0x3FFF0000,0xE0CCDEEC,0x2A94E110,0x20032C4B
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 425) .long 0x3FFF0000,0xE33F8972,0xBE8A5A50,0x2004DFF5
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 426) .long 0x3FFF0000,0xE5B906E7,0x7C8348A8,0x1E72F47A
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 427) .long 0x3FFF0000,0xE8396A50,0x3C4BDC68,0x1F722F22
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 428) .long 0x3FFF0000,0xEAC0C6E7,0xDD243930,0xA017E945
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 429) .long 0x3FFF0000,0xED4F301E,0xD9942B84,0x1F401A5B
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 430) .long 0x3FFF0000,0xEFE4B99B,0xDCDAF5CC,0x9FB9A9E3
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 431) .long 0x3FFF0000,0xF281773C,0x59FFB138,0x20744C05
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 432) .long 0x3FFF0000,0xF5257D15,0x2486CC2C,0x1F773A19
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 433) .long 0x3FFF0000,0xF7D0DF73,0x0AD13BB8,0x1FFE90D5
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 434) .long 0x3FFF0000,0xFA83B2DB,0x722A033C,0xA041ED22
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 435) .long 0x3FFF0000,0xFD3E0C0C,0xF486C174,0x1F853F3A
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 436)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 437) .set ADJFLAG,L_SCR2
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 438) .set SCALE,FP_SCR1
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 439) .set ADJSCALE,FP_SCR2
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 440) .set SC,FP_SCR3
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 441) .set ONEBYSC,FP_SCR4
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 442)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 443) | xref t_frcinx
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 444) |xref t_extdnrm
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 445) |xref t_unfl
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 446) |xref t_ovfl
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 447)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 448) .global setoxd
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 449) setoxd:
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 450) |--entry point for EXP(X), X is denormalized
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 451) movel (%a0),%d0
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 452) andil #0x80000000,%d0
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 453) oril #0x00800000,%d0 | ...sign(X)*2^(-126)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 454) movel %d0,-(%sp)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 455) fmoves #0x3F800000,%fp0
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 456) fmovel %d1,%fpcr
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 457) fadds (%sp)+,%fp0
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 458) bra t_frcinx
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 459)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 460) .global setox
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 461) setox:
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 462) |--entry point for EXP(X), here X is finite, non-zero, and not NaN's
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 463)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 464) |--Step 1.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 465) movel (%a0),%d0 | ...load part of input X
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 466) andil #0x7FFF0000,%d0 | ...biased expo. of X
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 467) cmpil #0x3FBE0000,%d0 | ...2^(-65)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 468) bges EXPC1 | ...normal case
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 469) bra EXPSM
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 470)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 471) EXPC1:
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 472) |--The case |X| >= 2^(-65)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 473) movew 4(%a0),%d0 | ...expo. and partial sig. of |X|
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 474) cmpil #0x400CB167,%d0 | ...16380 log2 trunc. 16 bits
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 475) blts EXPMAIN | ...normal case
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 476) bra EXPBIG
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 477)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 478) EXPMAIN:
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 479) |--Step 2.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 480) |--This is the normal branch: 2^(-65) <= |X| < 16380 log2.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 481) fmovex (%a0),%fp0 | ...load input from (a0)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 482)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 483) fmovex %fp0,%fp1
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 484) fmuls #0x42B8AA3B,%fp0 | ...64/log2 * X
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 485) fmovemx %fp2-%fp2/%fp3,-(%a7) | ...save fp2
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 486) movel #0,ADJFLAG(%a6)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 487) fmovel %fp0,%d0 | ...N = int( X * 64/log2 )
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 488) lea EXPTBL,%a1
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 489) fmovel %d0,%fp0 | ...convert to floating-format
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 490)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 491) movel %d0,L_SCR1(%a6) | ...save N temporarily
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 492) andil #0x3F,%d0 | ...D0 is J = N mod 64
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 493) lsll #4,%d0
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 494) addal %d0,%a1 | ...address of 2^(J/64)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 495) movel L_SCR1(%a6),%d0
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 496) asrl #6,%d0 | ...D0 is M
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 497) addiw #0x3FFF,%d0 | ...biased expo. of 2^(M)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 498) movew L2,L_SCR1(%a6) | ...prefetch L2, no need in CB
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 499)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 500) EXPCONT1:
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 501) |--Step 3.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 502) |--fp1,fp2 saved on the stack. fp0 is N, fp1 is X,
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 503) |--a0 points to 2^(J/64), D0 is biased expo. of 2^(M)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 504) fmovex %fp0,%fp2
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 505) fmuls #0xBC317218,%fp0 | ...N * L1, L1 = lead(-log2/64)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 506) fmulx L2,%fp2 | ...N * L2, L1+L2 = -log2/64
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 507) faddx %fp1,%fp0 | ...X + N*L1
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 508) faddx %fp2,%fp0 | ...fp0 is R, reduced arg.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 509) | MOVE.W #$3FA5,EXPA3 ...load EXPA3 in cache
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 510)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 511) |--Step 4.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 512) |--WE NOW COMPUTE EXP(R)-1 BY A POLYNOMIAL
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 513) |-- R + R*R*(A1 + R*(A2 + R*(A3 + R*(A4 + R*A5))))
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 514) |--TO FULLY UTILIZE THE PIPELINE, WE COMPUTE S = R*R
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 515) |--[R+R*S*(A2+S*A4)] + [S*(A1+S*(A3+S*A5))]
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 516)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 517) fmovex %fp0,%fp1
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 518) fmulx %fp1,%fp1 | ...fp1 IS S = R*R
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 519)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 520) fmoves #0x3AB60B70,%fp2 | ...fp2 IS A5
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 521) | MOVE.W #0,2(%a1) ...load 2^(J/64) in cache
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 522)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 523) fmulx %fp1,%fp2 | ...fp2 IS S*A5
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 524) fmovex %fp1,%fp3
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 525) fmuls #0x3C088895,%fp3 | ...fp3 IS S*A4
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 526)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 527) faddd EXPA3,%fp2 | ...fp2 IS A3+S*A5
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 528) faddd EXPA2,%fp3 | ...fp3 IS A2+S*A4
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 529)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 530) fmulx %fp1,%fp2 | ...fp2 IS S*(A3+S*A5)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 531) movew %d0,SCALE(%a6) | ...SCALE is 2^(M) in extended
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 532) clrw SCALE+2(%a6)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 533) movel #0x80000000,SCALE+4(%a6)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 534) clrl SCALE+8(%a6)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 535)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 536) fmulx %fp1,%fp3 | ...fp3 IS S*(A2+S*A4)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 537)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 538) fadds #0x3F000000,%fp2 | ...fp2 IS A1+S*(A3+S*A5)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 539) fmulx %fp0,%fp3 | ...fp3 IS R*S*(A2+S*A4)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 540)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 541) fmulx %fp1,%fp2 | ...fp2 IS S*(A1+S*(A3+S*A5))
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 542) faddx %fp3,%fp0 | ...fp0 IS R+R*S*(A2+S*A4),
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 543) | ...fp3 released
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 544)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 545) fmovex (%a1)+,%fp1 | ...fp1 is lead. pt. of 2^(J/64)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 546) faddx %fp2,%fp0 | ...fp0 is EXP(R) - 1
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 547) | ...fp2 released
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 548)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 549) |--Step 5
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 550) |--final reconstruction process
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 551) |--EXP(X) = 2^M * ( 2^(J/64) + 2^(J/64)*(EXP(R)-1) )
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 552)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 553) fmulx %fp1,%fp0 | ...2^(J/64)*(Exp(R)-1)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 554) fmovemx (%a7)+,%fp2-%fp2/%fp3 | ...fp2 restored
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 555) fadds (%a1),%fp0 | ...accurate 2^(J/64)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 556)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 557) faddx %fp1,%fp0 | ...2^(J/64) + 2^(J/64)*...
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 558) movel ADJFLAG(%a6),%d0
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 559)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 560) |--Step 6
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 561) tstl %d0
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 562) beqs NORMAL
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 563) ADJUST:
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 564) fmulx ADJSCALE(%a6),%fp0
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 565) NORMAL:
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 566) fmovel %d1,%FPCR | ...restore user FPCR
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 567) fmulx SCALE(%a6),%fp0 | ...multiply 2^(M)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 568) bra t_frcinx
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 569)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 570) EXPSM:
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 571) |--Step 7
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 572) fmovemx (%a0),%fp0-%fp0 | ...in case X is denormalized
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 573) fmovel %d1,%FPCR
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 574) fadds #0x3F800000,%fp0 | ...1+X in user mode
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 575) bra t_frcinx
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 576)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 577) EXPBIG:
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 578) |--Step 8
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 579) cmpil #0x400CB27C,%d0 | ...16480 log2
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 580) bgts EXP2BIG
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 581) |--Steps 8.2 -- 8.6
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 582) fmovex (%a0),%fp0 | ...load input from (a0)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 583)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 584) fmovex %fp0,%fp1
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 585) fmuls #0x42B8AA3B,%fp0 | ...64/log2 * X
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 586) fmovemx %fp2-%fp2/%fp3,-(%a7) | ...save fp2
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 587) movel #1,ADJFLAG(%a6)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 588) fmovel %fp0,%d0 | ...N = int( X * 64/log2 )
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 589) lea EXPTBL,%a1
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 590) fmovel %d0,%fp0 | ...convert to floating-format
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 591) movel %d0,L_SCR1(%a6) | ...save N temporarily
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 592) andil #0x3F,%d0 | ...D0 is J = N mod 64
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 593) lsll #4,%d0
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 594) addal %d0,%a1 | ...address of 2^(J/64)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 595) movel L_SCR1(%a6),%d0
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 596) asrl #6,%d0 | ...D0 is K
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 597) movel %d0,L_SCR1(%a6) | ...save K temporarily
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 598) asrl #1,%d0 | ...D0 is M1
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 599) subl %d0,L_SCR1(%a6) | ...a1 is M
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 600) addiw #0x3FFF,%d0 | ...biased expo. of 2^(M1)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 601) movew %d0,ADJSCALE(%a6) | ...ADJSCALE := 2^(M1)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 602) clrw ADJSCALE+2(%a6)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 603) movel #0x80000000,ADJSCALE+4(%a6)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 604) clrl ADJSCALE+8(%a6)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 605) movel L_SCR1(%a6),%d0 | ...D0 is M
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 606) addiw #0x3FFF,%d0 | ...biased expo. of 2^(M)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 607) bra EXPCONT1 | ...go back to Step 3
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 608)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 609) EXP2BIG:
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 610) |--Step 9
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 611) fmovel %d1,%FPCR
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 612) movel (%a0),%d0
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 613) bclrb #sign_bit,(%a0) | ...setox always returns positive
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 614) cmpil #0,%d0
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 615) blt t_unfl
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 616) bra t_ovfl
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 617)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 618) .global setoxm1d
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 619) setoxm1d:
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 620) |--entry point for EXPM1(X), here X is denormalized
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 621) |--Step 0.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 622) bra t_extdnrm
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 623)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 624)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 625) .global setoxm1
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 626) setoxm1:
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 627) |--entry point for EXPM1(X), here X is finite, non-zero, non-NaN
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 628)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 629) |--Step 1.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 630) |--Step 1.1
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 631) movel (%a0),%d0 | ...load part of input X
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 632) andil #0x7FFF0000,%d0 | ...biased expo. of X
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 633) cmpil #0x3FFD0000,%d0 | ...1/4
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 634) bges EM1CON1 | ...|X| >= 1/4
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 635) bra EM1SM
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 636)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 637) EM1CON1:
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 638) |--Step 1.3
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 639) |--The case |X| >= 1/4
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 640) movew 4(%a0),%d0 | ...expo. and partial sig. of |X|
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 641) cmpil #0x4004C215,%d0 | ...70log2 rounded up to 16 bits
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 642) bles EM1MAIN | ...1/4 <= |X| <= 70log2
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 643) bra EM1BIG
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 644)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 645) EM1MAIN:
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 646) |--Step 2.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 647) |--This is the case: 1/4 <= |X| <= 70 log2.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 648) fmovex (%a0),%fp0 | ...load input from (a0)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 649)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 650) fmovex %fp0,%fp1
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 651) fmuls #0x42B8AA3B,%fp0 | ...64/log2 * X
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 652) fmovemx %fp2-%fp2/%fp3,-(%a7) | ...save fp2
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 653) | MOVE.W #$3F81,EM1A4 ...prefetch in CB mode
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 654) fmovel %fp0,%d0 | ...N = int( X * 64/log2 )
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 655) lea EXPTBL,%a1
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 656) fmovel %d0,%fp0 | ...convert to floating-format
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 657)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 658) movel %d0,L_SCR1(%a6) | ...save N temporarily
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 659) andil #0x3F,%d0 | ...D0 is J = N mod 64
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 660) lsll #4,%d0
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 661) addal %d0,%a1 | ...address of 2^(J/64)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 662) movel L_SCR1(%a6),%d0
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 663) asrl #6,%d0 | ...D0 is M
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 664) movel %d0,L_SCR1(%a6) | ...save a copy of M
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 665) | MOVE.W #$3FDC,L2 ...prefetch L2 in CB mode
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 666)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 667) |--Step 3.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 668) |--fp1,fp2 saved on the stack. fp0 is N, fp1 is X,
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 669) |--a0 points to 2^(J/64), D0 and a1 both contain M
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 670) fmovex %fp0,%fp2
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 671) fmuls #0xBC317218,%fp0 | ...N * L1, L1 = lead(-log2/64)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 672) fmulx L2,%fp2 | ...N * L2, L1+L2 = -log2/64
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 673) faddx %fp1,%fp0 | ...X + N*L1
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 674) faddx %fp2,%fp0 | ...fp0 is R, reduced arg.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 675) | MOVE.W #$3FC5,EM1A2 ...load EM1A2 in cache
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 676) addiw #0x3FFF,%d0 | ...D0 is biased expo. of 2^M
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 677)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 678) |--Step 4.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 679) |--WE NOW COMPUTE EXP(R)-1 BY A POLYNOMIAL
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 680) |-- R + R*R*(A1 + R*(A2 + R*(A3 + R*(A4 + R*(A5 + R*A6)))))
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 681) |--TO FULLY UTILIZE THE PIPELINE, WE COMPUTE S = R*R
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 682) |--[R*S*(A2+S*(A4+S*A6))] + [R+S*(A1+S*(A3+S*A5))]
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 683)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 684) fmovex %fp0,%fp1
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 685) fmulx %fp1,%fp1 | ...fp1 IS S = R*R
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 686)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 687) fmoves #0x3950097B,%fp2 | ...fp2 IS a6
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 688) | MOVE.W #0,2(%a1) ...load 2^(J/64) in cache
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 689)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 690) fmulx %fp1,%fp2 | ...fp2 IS S*A6
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 691) fmovex %fp1,%fp3
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 692) fmuls #0x3AB60B6A,%fp3 | ...fp3 IS S*A5
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 693)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 694) faddd EM1A4,%fp2 | ...fp2 IS A4+S*A6
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 695) faddd EM1A3,%fp3 | ...fp3 IS A3+S*A5
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 696) movew %d0,SC(%a6) | ...SC is 2^(M) in extended
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 697) clrw SC+2(%a6)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 698) movel #0x80000000,SC+4(%a6)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 699) clrl SC+8(%a6)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 700)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 701) fmulx %fp1,%fp2 | ...fp2 IS S*(A4+S*A6)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 702) movel L_SCR1(%a6),%d0 | ...D0 is M
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 703) negw %d0 | ...D0 is -M
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 704) fmulx %fp1,%fp3 | ...fp3 IS S*(A3+S*A5)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 705) addiw #0x3FFF,%d0 | ...biased expo. of 2^(-M)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 706) faddd EM1A2,%fp2 | ...fp2 IS A2+S*(A4+S*A6)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 707) fadds #0x3F000000,%fp3 | ...fp3 IS A1+S*(A3+S*A5)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 708)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 709) fmulx %fp1,%fp2 | ...fp2 IS S*(A2+S*(A4+S*A6))
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 710) oriw #0x8000,%d0 | ...signed/expo. of -2^(-M)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 711) movew %d0,ONEBYSC(%a6) | ...OnebySc is -2^(-M)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 712) clrw ONEBYSC+2(%a6)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 713) movel #0x80000000,ONEBYSC+4(%a6)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 714) clrl ONEBYSC+8(%a6)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 715) fmulx %fp3,%fp1 | ...fp1 IS S*(A1+S*(A3+S*A5))
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 716) | ...fp3 released
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 717)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 718) fmulx %fp0,%fp2 | ...fp2 IS R*S*(A2+S*(A4+S*A6))
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 719) faddx %fp1,%fp0 | ...fp0 IS R+S*(A1+S*(A3+S*A5))
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 720) | ...fp1 released
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 721)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 722) faddx %fp2,%fp0 | ...fp0 IS EXP(R)-1
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 723) | ...fp2 released
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 724) fmovemx (%a7)+,%fp2-%fp2/%fp3 | ...fp2 restored
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 725)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 726) |--Step 5
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 727) |--Compute 2^(J/64)*p
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 728)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 729) fmulx (%a1),%fp0 | ...2^(J/64)*(Exp(R)-1)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 730)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 731) |--Step 6
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 732) |--Step 6.1
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 733) movel L_SCR1(%a6),%d0 | ...retrieve M
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 734) cmpil #63,%d0
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 735) bles MLE63
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 736) |--Step 6.2 M >= 64
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 737) fmoves 12(%a1),%fp1 | ...fp1 is t
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 738) faddx ONEBYSC(%a6),%fp1 | ...fp1 is t+OnebySc
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 739) faddx %fp1,%fp0 | ...p+(t+OnebySc), fp1 released
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 740) faddx (%a1),%fp0 | ...T+(p+(t+OnebySc))
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 741) bras EM1SCALE
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 742) MLE63:
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 743) |--Step 6.3 M <= 63
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 744) cmpil #-3,%d0
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 745) bges MGEN3
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 746) MLTN3:
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 747) |--Step 6.4 M <= -4
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 748) fadds 12(%a1),%fp0 | ...p+t
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 749) faddx (%a1),%fp0 | ...T+(p+t)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 750) faddx ONEBYSC(%a6),%fp0 | ...OnebySc + (T+(p+t))
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 751) bras EM1SCALE
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 752) MGEN3:
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 753) |--Step 6.5 -3 <= M <= 63
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 754) fmovex (%a1)+,%fp1 | ...fp1 is T
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 755) fadds (%a1),%fp0 | ...fp0 is p+t
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 756) faddx ONEBYSC(%a6),%fp1 | ...fp1 is T+OnebySc
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 757) faddx %fp1,%fp0 | ...(T+OnebySc)+(p+t)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 758)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 759) EM1SCALE:
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 760) |--Step 6.6
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 761) fmovel %d1,%FPCR
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 762) fmulx SC(%a6),%fp0
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 763)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 764) bra t_frcinx
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 765)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 766) EM1SM:
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 767) |--Step 7 |X| < 1/4.
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 768) cmpil #0x3FBE0000,%d0 | ...2^(-65)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 769) bges EM1POLY
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 770)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 771) EM1TINY:
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 772) |--Step 8 |X| < 2^(-65)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 773) cmpil #0x00330000,%d0 | ...2^(-16312)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 774) blts EM12TINY
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 775) |--Step 8.2
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 776) movel #0x80010000,SC(%a6) | ...SC is -2^(-16382)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 777) movel #0x80000000,SC+4(%a6)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 778) clrl SC+8(%a6)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 779) fmovex (%a0),%fp0
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 780) fmovel %d1,%FPCR
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 781) faddx SC(%a6),%fp0
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 782)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 783) bra t_frcinx
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 784)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 785) EM12TINY:
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 786) |--Step 8.3
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 787) fmovex (%a0),%fp0
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 788) fmuld TWO140,%fp0
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 789) movel #0x80010000,SC(%a6)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 790) movel #0x80000000,SC+4(%a6)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 791) clrl SC+8(%a6)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 792) faddx SC(%a6),%fp0
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 793) fmovel %d1,%FPCR
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 794) fmuld TWON140,%fp0
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 795)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 796) bra t_frcinx
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 797)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 798) EM1POLY:
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 799) |--Step 9 exp(X)-1 by a simple polynomial
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 800) fmovex (%a0),%fp0 | ...fp0 is X
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 801) fmulx %fp0,%fp0 | ...fp0 is S := X*X
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 802) fmovemx %fp2-%fp2/%fp3,-(%a7) | ...save fp2
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 803) fmoves #0x2F30CAA8,%fp1 | ...fp1 is B12
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 804) fmulx %fp0,%fp1 | ...fp1 is S*B12
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 805) fmoves #0x310F8290,%fp2 | ...fp2 is B11
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 806) fadds #0x32D73220,%fp1 | ...fp1 is B10+S*B12
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 807)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 808) fmulx %fp0,%fp2 | ...fp2 is S*B11
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 809) fmulx %fp0,%fp1 | ...fp1 is S*(B10 + ...
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 810)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 811) fadds #0x3493F281,%fp2 | ...fp2 is B9+S*...
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 812) faddd EM1B8,%fp1 | ...fp1 is B8+S*...
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 813)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 814) fmulx %fp0,%fp2 | ...fp2 is S*(B9+...
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 815) fmulx %fp0,%fp1 | ...fp1 is S*(B8+...
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 816)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 817) faddd EM1B7,%fp2 | ...fp2 is B7+S*...
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 818) faddd EM1B6,%fp1 | ...fp1 is B6+S*...
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 819)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 820) fmulx %fp0,%fp2 | ...fp2 is S*(B7+...
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 821) fmulx %fp0,%fp1 | ...fp1 is S*(B6+...
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 822)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 823) faddd EM1B5,%fp2 | ...fp2 is B5+S*...
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 824) faddd EM1B4,%fp1 | ...fp1 is B4+S*...
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 825)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 826) fmulx %fp0,%fp2 | ...fp2 is S*(B5+...
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 827) fmulx %fp0,%fp1 | ...fp1 is S*(B4+...
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 828)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 829) faddd EM1B3,%fp2 | ...fp2 is B3+S*...
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 830) faddx EM1B2,%fp1 | ...fp1 is B2+S*...
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 831)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 832) fmulx %fp0,%fp2 | ...fp2 is S*(B3+...
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 833) fmulx %fp0,%fp1 | ...fp1 is S*(B2+...
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 834)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 835) fmulx %fp0,%fp2 | ...fp2 is S*S*(B3+...)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 836) fmulx (%a0),%fp1 | ...fp1 is X*S*(B2...
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 837)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 838) fmuls #0x3F000000,%fp0 | ...fp0 is S*B1
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 839) faddx %fp2,%fp1 | ...fp1 is Q
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 840) | ...fp2 released
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 841)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 842) fmovemx (%a7)+,%fp2-%fp2/%fp3 | ...fp2 restored
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 843)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 844) faddx %fp1,%fp0 | ...fp0 is S*B1+Q
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 845) | ...fp1 released
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 846)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 847) fmovel %d1,%FPCR
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 848) faddx (%a0),%fp0
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 849)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 850) bra t_frcinx
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 851)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 852) EM1BIG:
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 853) |--Step 10 |X| > 70 log2
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 854) movel (%a0),%d0
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 855) cmpil #0,%d0
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 856) bgt EXPC1
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 857) |--Step 10.2
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 858) fmoves #0xBF800000,%fp0 | ...fp0 is -1
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 859) fmovel %d1,%FPCR
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 860) fadds #0x00800000,%fp0 | ...-1 + 2^(-126)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 861)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 862) bra t_frcinx
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 863)
^8f3ce5b39 (kx 2023-10-28 12:00:06 +0300 864) |end
Orange Pi5 kernel
Deprecated Linux kernel 5.10.110 for OrangePi 5/5B/5+ boards
3 Commits
0 Branches
0 Tags