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

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

VxWorks

开发平台：

C/C++

setox.s：源码内容
  1. /* setox.s - Motorola 68040 FP exponential routines (EXC) */
  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,10jan92,kdl  added modification history; general cleanup.
  16. 01a,15aug91,kdl  original version, from Motorola FPSP v2.0.
  17. */
  19. /*
  20. DESCRIPTION
  22. setoxsa 3.1 12/10/90
  24. The entry point __x_setox computes the exponential of a value.
  25. __x_setoxd does the same except the input value is a denormalized
  26. number. __x_setoxm1 computes exp(X)-1, and __x_setoxm1d computes
  27. exp(X)-1 for denormalized X.
  29. INPUT
  30. -----
  31. Double-extended value in memory location pointed to by address
  32. register a0.
  34. OUTPUT
  35. ------
  36. exp(X) or exp(X)-1 returned in floating-point register fp0.
  38. ACCURACY and MONOTONICITY
  39. -------------------------
  40. The returned result is within 0.85 ulps in 64 significant bit, i.e.
  41. within 0.5001 ulp to 53 bits if the result is subsequently rounded
  42. to double precision. The result is provably monotonic in double
  43. precision.
  45. SPEED
  46. -----
  47. Two timings are measured, both in the copy-back mode. The
  48. first one is measured when the function is invoked the first time
  49. (so the instructions and data are not in cache), and the
  50. second one is measured when the function is reinvoked at the same
  51. input argument.
  53. The program __x_setox takes approximately 210/190 cycles for input
  54. argument X whose magnitude is less than 16380 log2, which
  55. is the usual situation. For the less common arguments,
  56. depending on their values, the program may run faster or slower --
  57. but no worse than 10 slower even in the extreme cases.
  59. The program __x_setoxm1 takes approximately ???/??? cycles for input
  60. argument X, 0.25 <= |X| < 70log2. For |X| < 0.25, it takes
  61. approximately ???/??? cycles. For the less common arguments,
  62. depending on their values, the program may run faster or slower --
  63. but no worse than 10 slower even in the extreme cases.
  65. ALGORITHM and IMPLEMENTATION NOTES
  66. ----------------------------------
  68. __x_setoxd
  69. ------
  70. Step 1. Set ans := 1.0
  72. Step 2. Return ans := ans + sign(X)*2^(-126). Exit.
  73. Notes: This will always generate one exception -- inexact.
  76. __x_setox
  77. -----
  79. Step 1. Filter out extreme cases of input argument.
  80. 1.1 If |X| >= 2^(-65), go to Step 1.3.
  81. 1.2 Go to Step 7.
  82. 1.3 If |X| < 16380 log(2), go to Step 2.
  83. 1.4 Go to Step 8.
  84. Notes: The usual case should take the branches 1.1 -> 1.3 -> 2.
  85.  To avoid the use of floating-point comparisons, a
  86.  compact representation of |X| is used. This format is a
  87.  32-bit integer, the upper (more significant) 16 bits are
  88.  the sign and biased exponent field of |X||  the lower 16
  89.  bits are the 16 most significant fraction (including the
  90.  explicit bit) bits of |X|. Consequently, the comparisons
  91.  in Steps 1.1 and 1.3 can be performed by integer comparison.
  92.  Note also that the constant 16380 log(2) used in Step 1.3
  93.  is also in the compact form. Thus taking the branch
  94.  to Step 2 guarantees |X| < 16380 log(2). There is no harm
  95.  to have a small number of cases where |X| is less than,
  96.  but close to, 16380 log(2) and the branch to Step 9 is
  97.  taken.
  99. Step 2. Calculate N = round-to-nearest-int( X * 64/log2 ).
  100. 2.1 Set AdjFlag := 0 (indicates the branch 1.3 -> 2
  101. 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,
  212. AdjFlag := 1.
  213. 8.5 Calculate the address of the stored value 2^(J/64).
  214. 8.6 Create the values Scale = 2^M, AdjScale = 2^M1.
  215. 8.7 Go to Step 3.
  216. Notes: Refer to notes for 2.2 - 2.6.
  218. Step 9. Handle exp(X), |X| > 16480 log2.
  219. 9.1 If X < 0, go to 9.3
  220. 9.2 ans := Huge, go to 9.4
  221. 9.3 ans := Tiny.
  222. 9.4 Restore user fpcr.
  223. 9.5 Return ans := ans * ans. Exit.
  224. Notes: Exp(X) will surely overflow or underflow, depending on
  225.  X's sign. "Huge" and "Tiny" are respectively large/tiny
  226.  extended-precision numbers whose square over/underflow
  227.  with an inexact result. Thus, 9.5 always raises the
  228.  inexact together with either overflow or underflow.
  231. __x_setoxm1d
  232. --------
  234. Step 1. Set ans := 0
  236. Step 2. Return ans := X + ans. Exit.
  237. Notes: This will return X with the appropriate rounding
  238.  precision prescribed by the user fpcr.
  240. __x_setoxm1
  241. -------
  243. Step 1. Check |X|
  244. 1.1 If |X| >= 1/4, go to Step 1.3.
  245. 1.2 Go to Step 7.
  246. 1.3 If |X| < 70 log(2), go to Step 2.
  247. 1.4 Go to Step 10.
  248. Notes: The usual case should take the branches 1.1 -> 1.3 -> 2.
  249.  However, it is conceivable |X| can be small very often
  250.  because EXPM1 is intended to evaluate exp(X)-1 accurately
  251.  when |X| is small. For further details on the comparisons,
  252.  see the notes on Step 1 of __x_setox.
  254. Step 2. Calculate N = round-to-nearest-int( X * 64/log2 ).
  255. 2.1 N := round-to-nearest-integer( X * 64/log2 ).
  256. 2.2 Calculate J = N mod 64|  so J = 0,1,2,..., or 63.
  257. 2.3 Calculate M = (N - J)/64|  so N = 64M + J.
  258. 2.4 Calculate the address of the stored value of 2^(J/64).
  259. 2.5 Create the values Sc = 2^M and OnebySc := -2^(-M).
  260. Notes: See the notes on Step 2 of __x_setox.
  262. Step 3. Calculate X - N*log2/64.
  263. 3.1 R := X + N*L1, where L1 := single-precision(-log2/64).
  264. 3.2 R := R + N*L2, L2 := extended-precision(-log2/64 - L1).
  265. Notes: Applying the analysis of Step 3 of __x_setox in this case
  266.  shows that |R| <= 0.0055 (note that |X| <= 70 log2 in
  267.  this case).
  269. Step 4. Approximate exp(R)-1 by a polynomial
  270. p = R+R*R*(A1+R*(A2+R*(A3+R*(A4+R*(A5+R*A6)))))
  271. Notes: a) In order to reduce memory access, the coefficients are
  272.  made as "short" as possible: A1 (which is 1/2), A5 and A6
  273.  are single precision|  A2, A3 and A4 are double precision.
  274.  b) Even with the restriction above,
  275. |p - (exp(R)-1)| < |R| * 2^(-72.7)
  276.  for all |R| <= 0.0055.
  277.  c) To fully utilize the pipeline, p is separated into
  278.  two independent pieces of roughly equal complexity
  279. p = [ R*S*(A2 + S*(A4 + S*A6)) ] +
  280. [ R + S*(A1 + S*(A3 + S*A5)) ]
  281.  where S = R*R.
  283. Step 5. Compute 2^(J/64)*p by
  284. p := T*p
  285.  where T and t are the stored values for 2^(J/64).
  286. Notes: 2^(J/64) is stored as T and t where T+t approximates
  287.  2^(J/64) to roughly 85 bits|  T is in extended precision
  288.  and t is in single precision. Note also that T is rounded
  289.  to 62 bits so that the last two bits of T are zero. The
  290.  reason for such a special form is that T-1, T-2, and T-8
  291.  will all be exact --- a property that will be exploited
  292.  in Step 6 below. The total relative error in p is no
  293.  bigger than 2^(-67.7) compared to the final result.
  295. Step 6. Reconstruction of exp(X)-1
  296. exp(X)-1 = 2^M * ( 2^(J/64) + p - 2^(-M) ).
  297. 6.1 If M <= 63, go to Step 6.3.
  298. 6.2 ans := T + (p + (t + OnebySc)). Go to 6.6
  299. 6.3 If M >= -3, go to 6.5.
  300. 6.4 ans := (T + (p + t)) + OnebySc. Go to 6.6
  301. 6.5 ans := (T + OnebySc) + (p + t).
  302. 6.6 Restore user fpcr.
  303. 6.7 Return ans := Sc * ans. Exit.
  304. Notes: The various arrangements of the expressions give accurate
  305.  evaluations.
  307. Step 7. exp(X)-1 for |X| < 1/4.
  308. 7.1 If |X| >= 2^(-65), go to Step 9.
  309. 7.2 Go to Step 8.
  311. Step 8. Calculate exp(X)-1, |X| < 2^(-65).
  312. 8.1 If |X| < 2^(-16312), goto 8.3
  313. 8.2 Restore fpcr|  return ans := X - 2^(-16382). Exit.
  314. 8.3 X := X * 2^(140).
  315. 8.4 Restore fpcr|  ans := ans - 2^(-16382).
  316.  Return ans := ans*2^(140). Exit
  317. Notes: The idea is to return "X - tiny" under the user
  318.  precision and rounding modes. To avoid unnecessary
  319.  inefficiency, we stay away from denormalized numbers the
  320.  best we can. For |X| >= 2^(-16312), the straightforward
  321.  8.2 generates the inexact exception as the case warrants.
  323. Step 9. Calculate exp(X)-1, |X| < 1/4, by a polynomial
  324. p = X + X*X*(B1 + X*(B2 + |... + X*B12))
  325. Notes: a) In order to reduce memory access, the coefficients are
  326.  made as "short" as possible: B1 (which is 1/2), B9 to B12
  327.  are single precision|  B3 to B8 are double precision|  and
  328.  B2 is double extended.
  329.  b) Even with the restriction above,
  330. |p - (exp(X)-1)| < |X| 2^(-70.6)
  331.  for all |X| <= 0.251.
  332.  Note that 0.251 is slightly bigger than 1/4.
  333.  c) To fully preserve accuracy, the polynomial is computed
  334.  as X + ( S*B1 + Q ) where S = X*X and
  335. Q = X*S*(B2 + X*(B3 + |... + X*B12))
  336.  d) To fully utilize the pipeline, Q is separated into
  337.  two independent pieces of roughly equal complexity
  338. Q = [ X*S*(B2 + S*(B4 + |... + S*B12)) ] +
  339. [ S*S*(B3 + S*(B5 + |... + S*B11)) ]
  341. Step 10. Calculate exp(X)-1 for |X| >= 70 log 2.
  342. 10.1 If X >= 70log2 , exp(X) - 1 = exp(X) for all practical
  343.  purposes. Therefore, go to Step 1 of __x_setox.
  344. 10.2 If X <= -70log2, exp(X) - 1 = -1 for all practical purposes.
  345.  ans := -1
  346.  Restore user fpcr
  347.  Return ans := ans + 2^(-126). Exit.
  348. Notes: 10.2 will always create an inexact and return -1 + tiny
  349.  in the user rounding precision and mode.
  353. Copyright (C) Motorola, Inc. 1990
  354. All Rights Reserved
  356. THIS IS UNPUBLISHED PROPRIETARY SOURCE CODE OF MOTOROLA
  357. The copyright notice above does not evidence any
  358. actual or intended publication of such source code.
  360. __x_setox IDNT 2,1 Motorola 040 Floating Point Software Package
  362. section 8
  364. NOMANUAL
  365. */
  367. #include "fpsp040E.h"
  369. L2: .long 0x3FDC0000,0x82E30865,0x4361C4C6,0x00000000
  371. EXPA3: .long 0x3FA55555,0x55554431
  372. EXPA2: .long 0x3FC55555,0x55554018
  374. HUGE: .long 0x7FFE0000,0xFFFFFFFF,0xFFFFFFFF,0x00000000
  375. TINY: .long 0x00010000,0xFFFFFFFF,0xFFFFFFFF,0x00000000
  377. EM1A4: .long 0x3F811111,0x11174385
  378. EM1A3: .long 0x3FA55555,0x55554F5A
  380. EM1A2: .long 0x3FC55555,0x55555555,0x00000000,0x00000000
  382. EM1B8: .long 0x3EC71DE3,0xA5774682
  383. EM1B7: .long 0x3EFA01A0,0x19D7CB68
  385. EM1B6: .long 0x3F2A01A0,0x1A019DF3
  386. EM1B5: .long 0x3F56C16C,0x16C170E2
  388. EM1B4: .long 0x3F811111,0x11111111
  389. EM1B3: .long 0x3FA55555,0x55555555
  391. EM1B2: .long 0x3FFC0000,0xAAAAAAAA,0xAAAAAAAB
  392. .long 0x00000000
  394. TWO140: .long 0x48B00000,0x00000000
  395. TWON140: .long 0x37300000,0x00000000
  397. EXPTBL:
  398. .long 0x3FFF0000,0x80000000,0x00000000,0x00000000
  399. .long 0x3FFF0000,0x8164D1F3,0xBC030774,0x9F841A9B
  400. .long 0x3FFF0000,0x82CD8698,0xAC2BA1D8,0x9FC1D5B9
  401. .long 0x3FFF0000,0x843A28C3,0xACDE4048,0xA0728369
  402. .long 0x3FFF0000,0x85AAC367,0xCC487B14,0x1FC5C95C
  403. .long 0x3FFF0000,0x871F6196,0x9E8D1010,0x1EE85C9F
  404. .long 0x3FFF0000,0x88980E80,0x92DA8528,0x9FA20729
  405. .long 0x3FFF0000,0x8A14D575,0x496EFD9C,0xA07BF9AF
  406. .long 0x3FFF0000,0x8B95C1E3,0xEA8BD6E8,0xA0020DCF
  407. .long 0x3FFF0000,0x8D1ADF5B,0x7E5BA9E4,0x205A63DA
  408. .long 0x3FFF0000,0x8EA4398B,0x45CD53C0,0x1EB70051
  409. .long 0x3FFF0000,0x9031DC43,0x1466B1DC,0x1F6EB029
  410. .long 0x3FFF0000,0x91C3D373,0xAB11C338,0xA0781494
  411. .long 0x3FFF0000,0x935A2B2F,0x13E6E92C,0x9EB319B0
  412. .long 0x3FFF0000,0x94F4EFA8,0xFEF70960,0x2017457D
  413. .long 0x3FFF0000,0x96942D37,0x20185A00,0x1F11D537
  414. .long 0x3FFF0000,0x9837F051,0x8DB8A970,0x9FB952DD
  415. .long 0x3FFF0000,0x99E04593,0x20B7FA64,0x1FE43087
  416. .long 0x3FFF0000,0x9B8D39B9,0xD54E5538,0x1FA2A818
  417. .long 0x3FFF0000,0x9D3ED9A7,0x2CFFB750,0x1FDE494D
  418. .long 0x3FFF0000,0x9EF53260,0x91A111AC,0x20504890
  419. .long 0x3FFF0000,0xA0B0510F,0xB9714FC4,0xA073691C
  420. .long 0x3FFF0000,0xA2704303,0x0C496818,0x1F9B7A05
  421. .long 0x3FFF0000,0xA43515AE,0x09E680A0,0xA0797126
  422. .long 0x3FFF0000,0xA5FED6A9,0xB15138EC,0xA071A140
  423. .long 0x3FFF0000,0xA7CD93B4,0xE9653568,0x204F62DA
  424. .long 0x3FFF0000,0xA9A15AB4,0xEA7C0EF8,0x1F283C4A
  425. .long 0x3FFF0000,0xAB7A39B5,0xA93ED338,0x9F9A7FDC
  426. .long 0x3FFF0000,0xAD583EEA,0x42A14AC8,0xA05B3FAC
  427. .long 0x3FFF0000,0xAF3B78AD,0x690A4374,0x1FDF2610
  428. .long 0x3FFF0000,0xB123F581,0xD2AC2590,0x9F705F90
  429. .long 0x3FFF0000,0xB311C412,0xA9112488,0x201F678A
  430. .long 0x3FFF0000,0xB504F333,0xF9DE6484,0x1F32FB13
  431. .long 0x3FFF0000,0xB6FD91E3,0x28D17790,0x20038B30
  432. .long 0x3FFF0000,0xB8FBAF47,0x62FB9EE8,0x200DC3CC
  433. .long 0x3FFF0000,0xBAFF5AB2,0x133E45FC,0x9F8B2AE6
  434. .long 0x3FFF0000,0xBD08A39F,0x580C36C0,0xA02BBF70
  435. .long 0x3FFF0000,0xBF1799B6,0x7A731084,0xA00BF518
  436. .long 0x3FFF0000,0xC12C4CCA,0x66709458,0xA041DD41
  437. .long 0x3FFF0000,0xC346CCDA,0x24976408,0x9FDF137B
  438. .long 0x3FFF0000,0xC5672A11,0x5506DADC,0x201F1568
  439. .long 0x3FFF0000,0xC78D74C8,0xABB9B15C,0x1FC13A2E
  440. .long 0x3FFF0000,0xC9B9BD86,0x6E2F27A4,0xA03F8F03
  441. .long 0x3FFF0000,0xCBEC14FE,0xF2727C5C,0x1FF4907D
  442. .long 0x3FFF0000,0xCE248C15,0x1F8480E4,0x9E6E53E4
  443. .long 0x3FFF0000,0xD06333DA,0xEF2B2594,0x1FD6D45C
  444. .long 0x3FFF0000,0xD2A81D91,0xF12AE45C,0xA076EDB9
  445. .long 0x3FFF0000,0xD4F35AAB,0xCFEDFA20,0x9FA6DE21
  446. .long 0x3FFF0000,0xD744FCCA,0xD69D6AF4,0x1EE69A2F
  447. .long 0x3FFF0000,0xD99D15C2,0x78AFD7B4,0x207F439F
  448. .long 0x3FFF0000,0xDBFBB797,0xDAF23754,0x201EC207
  449. .long 0x3FFF0000,0xDE60F482,0x5E0E9124,0x9E8BE175
  450. .long 0x3FFF0000,0xE0CCDEEC,0x2A94E110,0x20032C4B
  451. .long 0x3FFF0000,0xE33F8972,0xBE8A5A50,0x2004DFF5
  452. .long 0x3FFF0000,0xE5B906E7,0x7C8348A8,0x1E72F47A
  453. .long 0x3FFF0000,0xE8396A50,0x3C4BDC68,0x1F722F22
  454. .long 0x3FFF0000,0xEAC0C6E7,0xDD243930,0xA017E945
  455. .long 0x3FFF0000,0xED4F301E,0xD9942B84,0x1F401A5B
  456. .long 0x3FFF0000,0xEFE4B99B,0xDCDAF5CC,0x9FB9A9E3
  457. .long 0x3FFF0000,0xF281773C,0x59FFB138,0x20744C05
  458. .long 0x3FFF0000,0xF5257D15,0x2486CC2C,0x1F773A19
  459. .long 0x3FFF0000,0xF7D0DF73,0x0AD13BB8,0x1FFE90D5
  460. .long 0x3FFF0000,0xFA83B2DB,0x722A033C,0xA041ED22
  461. .long 0x3FFF0000,0xFD3E0C0C,0xF486C174,0x1F853F3A
  463. #define ADJFLAG  L_SCR2
  464. #define SCALE  FP_SCR1
  465. #define ADJSCALE  FP_SCR2
  466. #define SC  FP_SCR3
  467. #define ONEBYSC  FP_SCR4
  469. | xref __x_t_frcinx
  470. | xref __x_t_extdnrm
  471. | xref __x_t_unfl
  472. | xref __x_t_ovfl
  475. .text
  477. .globl __x_setoxd
  478. __x_setoxd:
  479. |--entry point for EXP(X), X is denormalized
  480. movel a0@,d0
  481. andil #0x80000000,d0
  482. oril #0x00800000,d0 |...sign(X)*2^(-126)
  483. movel d0,a7@-
  484. /* fmoves &0x3F800000,fp0 */  .long 0xf23c4400,0x3f800000
  485. fmovel d1,fpcr
  486. fadds a7@+,fp0
  487. jra  __x_t_frcinx
  489. .globl __x_setox
  490. __x_setox:
  491. /* |--entry point for EXP(X), here X is finite, non-zero, and not NaN's */
  493. |--Step 1.
  494. movel a0@,d0  |...load part of input X
  495. andil #0x7FFF0000,d0 |...biased expo. of X
  496. cmpil #0x3FBE0000,d0 |...2^(-65)
  497. jge  EXPC1 |...normal case
  498. jra  EXPSM
  500. EXPC1:
  501. |--The case |X| >= 2^(-65)
  502. movew a0@(4),d0 |...expo. and partial sig. of |X|
  503. cmpil #0x400CB167,d0 |...16380 log2 trunc. 16 bits
  504. jlt  EXPMAIN  |...normal case
  505. jra  EXPBIG
  507. EXPMAIN:
  508. |--Step 2.
  509. |--This is the __x_normal branch: 2^(-65) <= |X| < 16380 log2.
  510. fmovex a0@,fp0 |...load input from a0@
  512. fmovex fp0,fp1
  513. /* fmuls &0x42B8AA3B,fp0 */  .long 0xf23c4423,0x42b8aa3b
  514. fmovemx fp2/fp3,a7@- |...save fp2
  515. movel #0,a6@(ADJFLAG)
  516. fmovel fp0,d0 |...N = int( X * 64/log2 )
  517. lea EXPTBL,a1
  518. fmovel d0,fp0 |...convert to floating-format
  520. movel d0,a6@(L_SCR1) |...save N temporarily
  521. andil #0x3F,d0 |...D0 is J = N mod 64
  522. lsll #4,d0
  523. addal d0,a1 |...address of 2^(J/64)
  524. movel a6@(L_SCR1),d0
  525. asrl #6,d0 |...D0 is M
  526. addiw #0x3FFF,d0 |...biased expo. of 2^(M)
  527. movew L2,a6@(L_SCR1) |...prefetch L2, no need in CB
  529. EXPCONT1:
  530. |--Step 3.
  531. |--fp1,fp2 saved on the stack. fp0 is N, fp1 is X,
  532. |--a0 points to 2^(J/64), D0 is biased expo. of 2^(M)
  533. fmovex fp0,fp2
  534. /* fmuls &0xBC317218,fp0 */  .long 0xf23c4423,0xbc317218
  535. fmulx L2,fp2 |...N * L2, L1+L2 = -log2/64
  536. faddx fp1,fp0   |...X + N*L1
  537. faddx fp2,fp0 |...fp0 is R, reduced arg.
  538. | MOVE.w #0x3FA5,EXPA3 |...load EXPA3 in cache
  540. |--Step 4.
  541. |--WE NOW COMPUTE EXP(R)-1 BY A POLYNOMIAL
  542. |-- R + R*R*(A1 + R*(A2 + R*(A3 + R*(A4 + R*A5))))
  543. |--TO FULLY UTILIZE THE PIPELINE, WE COMPUTE S = R*R
  544. |--[R+R*S*(A2+S*A4)] + [S*(A1+S*(A3+S*A5))]
  546. fmovex fp0,fp1
  547. fmulx fp1,fp1   |...fp1 IS S = R*R
  549. /* fmoves &0x3AB60B70,fp2 */  .long 0xf23c4500,0x3ab60b70
  550. | MOVE.w #0,a1@(2) |...load 2^(J/64) in cache
  552. fmulx fp1,fp2   |...fp2 IS S*A5
  553. fmovex fp1,fp3
  554. /* fmuls &0x3C088895,fp3 */  .long 0xf23c45a3,0x3c088895
  556. faddd EXPA3,fp2 |...fp2 IS a3+S*A5
  557. faddd EXPA2,fp3 |...fp3 IS a2+S*A4
  559. fmulx fp1,fp2   |...fp2 IS S*(A3+S*A5)
  560. movew d0,a6@(SCALE) |...SCALE is 2^(M) in extended
  561. clrw a6@(SCALE+2)
  562. movel #0x80000000,a6@(SCALE+4)
  563. clrl a6@(SCALE+8)
  565. fmulx fp1,fp3   |...fp3 IS S*(A2+S*A4)
  567. /* fadds &0x3F000000,fp2 */  .long 0xf23c4522,0x3f000000
  568. fmulx fp0,fp3   |...fp3 IS R*S*(A2+S*A4)
  570. fmulx fp1,fp2   |...fp2 IS S*(A1+S*(A3+S*A5))
  571. faddx fp3,fp0   |...fp0 IS R+R*S*(A2+S*A4),
  572. | |...fp3 released
  574. fmovex a1@+,fp1 |...fp1 is lead. pt. of 2^(J/64)
  575. faddx fp2,fp0   |...fp0 is EXP(R) - 1
  576. | |...fp2 released
  578. |--Step 5
  579. |--final reconstruction process
  580. |--EXP(X) = 2^M * ( 2^(J/64) + 2^(J/64)*(EXP(R)-1) )
  582. fmulx fp1,fp0   |...2^(J/64)*(Exp(R)-1)
  583. fmovemx a7@+,fp2/fp3 |...fp2 restored
  584. fadds a1@,fp0 |...accurate 2^(J/64)
  586. faddx fp1,fp0   |...2^(J/64) + 2^(J/64)*...
  587. movel a6@(ADJFLAG),d0
  589. |--Step 6
  590. tstl d0
  591. jeq  NORMAL
  592. ADJUST:
  593. fmulx a6@(ADJSCALE),fp0
  594. NORMAL:
  595. fmovel d1,fpcr   |...restore user fpcr
  596. fmulx a6@(SCALE),fp0 |...multiply 2^(M)
  597. jra  __x_t_frcinx
  599. EXPSM:
  600. |--Step 7
  601. fmovemx a0@,fp0-fp0 |...in case X is denormalized
  602. fmovel d1,fpcr
  603. .long 0xf23c4422,0x3f800000 /*  fadds  &0x3F800000,fp0 */
  604. jra  __x_t_frcinx
  606. EXPBIG:
  607. |--Step 8
  608. cmpil #0x400CB27C,d0 |...16480 log2
  609. jgt  EXP2BIG
  610. |--Steps 8.2 -- 8.6
  611. fmovex a0@,fp0 |...load input from a0@
  613. fmovex fp0,fp1
  614. .long 0xf23c4423,0x42b8aa3b /*  fmuls  &0x42B8AA3B,fp0 */
  615. fmovemx  fp2/fp3,a7@- |...save fp2
  616. movel #1,a6@(ADJFLAG)
  617. fmovel fp0,d0 |...N = int( X * 64/log2 )
  618. lea EXPTBL,a1
  619. fmovel d0,fp0 |...convert to floating-format
  620. movel d0,a6@(L_SCR1) |...save N temporarily
  621. andil #0x3F,d0 |...D0 is J = N mod 64
  622. lsll #4,d0
  623. addal d0,a1 |...address of 2^(J/64)
  624. movel a6@(L_SCR1),d0
  625. asrl #6,d0 |...D0 is K
  626. movel d0,a6@(L_SCR1) |...save K temporarily
  627. asrl #1,d0 |...D0 is M1
  628. subl d0,a6@(L_SCR1) |...a1 is M
  629. addiw #0x3FFF,d0 |...biased expo. of 2^(M1)
  630. movew d0,a6@(ADJSCALE) |...ADJSCALE := 2^(M1)
  631. clrw a6@(ADJSCALE+2)
  632. movel #0x80000000,a6@(ADJSCALE+4)
  633. clrl a6@(ADJSCALE+8)
  634. movel a6@(L_SCR1),d0 |...D0 is M
  635. addiw #0x3FFF,d0 |...biased expo. of 2^(M)
  636. jra  EXPCONT1 |...go back to Step 3
  638. EXP2BIG:
  639. |--Step 9
  640. fmovel d1,fpcr
  641. movel a0@,d0
  642. bclr #sign_bit,a0@ |...setox always returns positive
  643. cmpil #0,d0
  644. jlt  __x_t_unfl
  645. jra  __x_t_ovfl
  647. .globl __x_setoxm1d
  648. __x_setoxm1d:
  649. |--entry point for EXPM1(X), here X is denormalized
  650. |--Step 0.
  651. jra  __x_t_extdnrm
  654. .globl __x_setoxm1
  655. __x_setoxm1:
  656. |--entry point for EXPM1(X), here X is finite, non-zero, non-NaN
  658. |--Step 1.
  659. |--Step 1.1
  660. movel a0@,d0   |...load part of input X
  661. andil #0x7FFF0000,d0 |...biased expo. of X
  662. cmpil #0x3FFD0000,d0 |...1/4
  663. jge  EM1CON1   |...|X| >= 1/4
  664. jra  EM1SM
  666. EM1CON1:
  667. |--Step 1.3
  668. |--The case |X| >= 1/4
  669. movew a0@(4),d0 |...expo. and partial sig. of |X|
  670. cmpil #0x4004C215,d0 |...70log2 rounded up to 16 bits
  671. jle  EM1MAIN   |...1/4 <= |X| <= 70log2
  672. jra  EM1BIG
  674. EM1MAIN:
  675. |--Step 2.
  676. |--This is the case: 1/4 <= |X| <= 70 log2.
  677. fmovex a0@,fp0 |...load input from a0@
  679. fmovex fp0,fp1
  680. .long 0xf23c4423,0x42b8aa3b /*  fmuls  &0x42B8AA3B,fp0 */
  681. fmovemx fp2/fp3,a7@- |...save fp2
  682. | MOVE.w #0x3F81,EM1A4 |...prefetch in CB mode
  683. fmovel fp0,d0 |...N = int( X * 64/log2 )
  684. lea EXPTBL,a1
  685. fmovel d0,fp0 |...convert to floating-format
  687. movel d0,a6@(L_SCR1) |...save N temporarily
  688. andil #0x3F,d0  |...D0 is J = N mod 64
  689. lsll #4,d0
  690. addal d0,a1 |...address of 2^(J/64)
  691. movel a6@(L_SCR1),d0
  692. asrl #6,d0 |...D0 is M
  693. movel d0,a6@(L_SCR1) |...save a copy of M
  694. | MOVE.w #0x3FDC,L2 |...prefetch L2 in CB mode
  696. |--Step 3.
  697. |--fp1,fp2 saved on the stack. fp0 is N, fp1 is X,
  698. |--a0 points to 2^(J/64), D0 and a1 both contain M
  699. fmovex fp0,fp2
  700. .long 0xf23c4423,0xbc317218 /*  fmuls  &0xBC317218,fp0 */
  701. fmulx L2,fp2 |...N * L2, L1+L2 = -log2/64
  702. faddx fp1,fp0   |...X + N*L1
  703. faddx fp2,fp0   |...fp0 is R, reduced arg.
  704. | MOVE.w #0x3FC5,EM1A2 |...load EM1A2 in cache
  705. addiw #0x3FFF,d0 |...D0 is biased expo. of 2^M
  707. |--Step 4.
  708. |--WE NOW COMPUTE EXP(R)-1 BY A POLYNOMIAL
  709. |-- R + R*R*(A1 + R*(A2 + R*(A3 + R*(A4 + R*(A5 + R*A6)))))
  710. |--TO FULLY UTILIZE THE PIPELINE, WE COMPUTE S = R*R
  711. |--[R*S*(A2+S*(A4+S*A6))] + [R+S*(A1+S*(A3+S*A5))]
  713. fmovex fp0,fp1
  714. fmulx fp1,fp1 |...fp1 IS S = R*R
  716. .long 0xf23c4500,0x3950097b /*  fmoves  &0x3950097B,fp2 */
  717. | MOVE.w #0,a1@(2) |...load 2^(J/64) in cache
  719. fmulx fp1,fp2 |...fp2 IS S*A6
  720. fmovex fp1,fp3
  721. .long 0xf23c45a3,0x3ab60b6a /*  fmuls  &0x3AB60B6A,fp3 */
  723. faddd EM1A4,fp2 |...fp2 IS a4+S*A6
  724. faddd EM1A3,fp3 |...fp3 IS a3+S*A5
  725. movew d0,a6@(SC) |...SC is 2^(M) in extended
  726. clrw a6@(SC+2)
  727. movel #0x80000000,a6@(SC+4)
  728. clrl a6@(SC+8)
  730. fmulx fp1,fp2 |...fp2 IS S*(A4+S*A6)
  731. movel a6@(L_SCR1),d0 |...D0 is M
  732. negw d0 |...D0 is -M
  733. fmulx fp1,fp3 |...fp3 IS S*(A3+S*A5)
  734. addiw #0x3FFF,d0 |...biased expo. of 2^(-M)
  735. faddd EM1A2,fp2 |...fp2 IS a2+S*(A4+S*A6)
  736. .long 0xf23c45a2,0x3f000000 /*  fadds  &0x3F000000,fp3 */
  738. fmulx fp1,fp2 |...fp2 IS S*(A2+S*(A4+S*A6))
  739. oriw #0x8000,d0 |...signed/expo. of -2^(-M)
  740. movew d0,a6@(ONEBYSC) |...OnebySc is -2^(-M)
  741. clrw a6@(ONEBYSC+2)
  742. movel #0x80000000,a6@(ONEBYSC+4)
  743. clrl a6@(ONEBYSC+8)
  744. fmulx fp3,fp1 |...fp1 IS S*(A1+S*(A3+S*A5))
  745. | |...fp3 released
  747. fmulx fp0,fp2 |...fp2 IS R*S*(A2+S*(A4+S*A6))
  748. faddx fp1,fp0 |...fp0 IS R+S*(A1+S*(A3+S*A5))
  749. | |...fp1 released
  751. faddx fp2,fp0 |...fp0 IS EXP(R)-1
  752. | |...fp2 released
  753. fmovemx a7@+,fp2/fp3 |...fp2 restored
  755. |--Step 5
  756. |--Compute 2^(J/64)*p
  758. fmulx a1@,fp0 |...2^(J/64)*(Exp(R)-1)
  760. |--Step 6
  761. |--Step 6.1
  762. movel a6@(L_SCR1),d0 |...retrieve M
  763. cmpil #63,d0
  764. jle  MLE63
  765. |--Step 6.2 M >= 64
  766. fmoves a1@(12),fp1 |...fp1 is t
  767. faddx a6@(ONEBYSC),fp1 |...fp1 is t+OnebySc
  768. faddx fp1,fp0 |...p+(t+OnebySc), fp1 released
  769. faddx a1@,fp0 |...T+(p+(t+OnebySc))
  770. jra  EM1SCALE
  771. MLE63:
  772. |--Step 6.3 M <= 63
  773. cmpil #-3,d0
  774. jge  MGEN3
  775. MLTN3:
  776. |--Step 6.4 M <= -4
  777. fadds a1@(12),fp0 |...p+t
  778. faddx a1@,fp0 |...T+(p+t)
  779. faddx a6@(ONEBYSC),fp0 |...OnebySc + (T+(p+t))
  780. jra  EM1SCALE
  781. MGEN3:
  782. |--Step 6.5 -3 <= M <= 63
  783. fmovex a1@+,fp1 |...fp1 is T
  784. fadds a1@,fp0 |...fp0 is p+t
  785. faddx a6@(ONEBYSC),fp1 |...fp1 is T+OnebySc
  786. faddx fp1,fp0 |...(T+OnebySc)+(p+t)
  788. EM1SCALE:
  789. |--Step 6.6
  790. fmovel d1,fpcr
  791. fmulx a6@(SC),fp0
  793. jra  __x_t_frcinx
  795. EM1SM:
  796. |--Step 7 |X| < 1/4.
  797. cmpil #0x3FBE0000,d0 |...2^(-65)
  798. jge  EM1POLY
  800. EM1TINY:
  801. |--Step 8 |X| < 2^(-65)
  802. cmpil #0x00330000,d0 |...2^(-16312)
  803. jlt  EM12TINY
  804. |--Step 8.2
  805. movel #0x80010000,a6@(SC) |...SC is -2^(-16382)
  806. movel #0x80000000,a6@(SC+4)
  807. clrl a6@(SC+8)
  808. fmovex a0@,fp0
  809. fmovel d1,fpcr
  810. faddx a6@(SC),fp0
  812. jra  __x_t_frcinx
  814. EM12TINY:
  815. |--Step 8.3
  816. fmovex a0@,fp0
  817. fmuld TWO140,fp0
  818. movel #0x80010000,a6@(SC)
  819. movel #0x80000000,a6@(SC+4)
  820. clrl a6@(SC+8)
  821. faddx a6@(SC),fp0
  822. fmovel d1,fpcr
  823. fmuld TWON140,fp0
  825. jra  __x_t_frcinx
  827. EM1POLY:
  828. |--Step 9 exp(X)-1 by a simple polynomial
  829. fmovex a0@,fp0 |...fp0 is X
  830. fmulx fp0,fp0 |...fp0 is S := X*X
  831. fmovemx fp2/fp3,a7@- |...save fp2
  832. .long 0xf23c4480,0x2f30caa8 /*  fmoves  &0x2F30CAA8,fp1 */
  833. fmulx fp0,fp1 |...fp1 is S*B12
  834. .long 0xf23c4500,0x310f8290 /*  fmoves  &0x310F8290,fp2 */
  835. .long 0xf23c44a2,0x32d73220 /*  fadds  &0x32D73220,fp1 */
  837. fmulx fp0,fp2 |...fp2 is S*B11
  838. fmulx fp0,fp1 |...fp1 is S*(B10 + ...
  840. .long 0xf23c4522,0x3493f281 /*  fadds  &0x3493F281,fp2 */
  841. faddd EM1B8,fp1 |...fp1 is B8+S*...
  843. fmulx fp0,fp2 |...fp2 is S*(B9+...
  844. fmulx fp0,fp1 |...fp1 is S*(B8+...
  846. faddd EM1B7,fp2 |...fp2 is B7+S*...
  847. faddd EM1B6,fp1 |...fp1 is B6+S*...
  849. fmulx fp0,fp2 |...fp2 is S*(B7+...
  850. fmulx fp0,fp1 |...fp1 is S*(B6+...
  852. faddd EM1B5,fp2 |...fp2 is B5+S*...
  853. faddd EM1B4,fp1 |...fp1 is B4+S*...
  855. fmulx fp0,fp2 |...fp2 is S*(B5+...
  856. fmulx fp0,fp1 |...fp1 is S*(B4+...
  858. faddd EM1B3,fp2 |...fp2 is B3+S*...
  859. faddx EM1B2,fp1 |...fp1 is B2+S*...
  861. fmulx fp0,fp2 |...fp2 is S*(B3+...
  862. fmulx fp0,fp1 |...fp1 is S*(B2+...
  864. fmulx fp0,fp2 |...fp2 is S*S*(B3+...)
  865. fmulx a0@,fp1 |...fp1 is X*S*(B2...
  867. .long 0xf23c4423,0x3f000000 /*  fmuls  &0x3F000000,fp0 */
  868. faddx fp2,fp1 |...fp1 is Q
  869. | |...fp2 released
  871. fmovemx a7@+,fp2/fp3 |...fp2 restored
  873. faddx fp1,fp0 |...fp0 is S*B1+Q
  874. | |...fp1 released
  876. fmovel d1,fpcr
  877. faddx a0@,fp0
  879. jra  __x_t_frcinx
  881. EM1BIG:
  882. |--Step 10 |X| > 70 log2
  883. movel a0@,d0
  884. cmpil #0,d0
  885. jgt  EXPC1
  886. |--Step 10.2
  887. .long 0xf23c4400,0xbf800000 /*  fmoves  &0xBF800000,fp0 */
  888. fmovel d1,fpcr
  889. .long 0xf23c4422,0x00800000 /*  fadds  &0x00800000,fp0 */
  891. jra  __x_t_frcinx
  893. | end