开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 嵌入式/单片机编程 > VxWorks > 查看源码
l_setox.s
资源名称：vxwork_src.rar　[点击查看]
上传用户：nvosite88
上传日期：2007-01-17
资源大小：4983k
文件大小：28k
源码类别：

VxWorks

开发平台：

C/C++

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