stwotox.S 12 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426
  1. |
  2. | stwotox.sa 3.1 12/10/90
  3. |
  4. | stwotox --- 2**X
  5. | stwotoxd --- 2**X for denormalized X
  6. | stentox --- 10**X
  7. | stentoxd --- 10**X for denormalized X
  8. |
  9. | Input: Double-extended number X in location pointed to
  10. | by address register a0.
  11. |
  12. | Output: The function values are returned in Fp0.
  13. |
  14. | Accuracy and Monotonicity: The returned result is within 2 ulps in
  15. | 64 significant bit, i.e. within 0.5001 ulp to 53 bits if the
  16. | result is subsequently rounded to double precision. The
  17. | result is provably monotonic in double precision.
  18. |
  19. | Speed: The program stwotox takes approximately 190 cycles and the
  20. | program stentox takes approximately 200 cycles.
  21. |
  22. | Algorithm:
  23. |
  24. | twotox
  25. | 1. If |X| > 16480, go to ExpBig.
  26. |
  27. | 2. If |X| < 2**(-70), go to ExpSm.
  28. |
  29. | 3. Decompose X as X = N/64 + r where |r| <= 1/128. Furthermore
  30. | decompose N as
  31. | N = 64(M + M') + j, j = 0,1,2,...,63.
  32. |
  33. | 4. Overwrite r := r * log2. Then
  34. | 2**X = 2**(M') * 2**(M) * 2**(j/64) * exp(r).
  35. | Go to expr to compute that expression.
  36. |
  37. | tentox
  38. | 1. If |X| > 16480*log_10(2) (base 10 log of 2), go to ExpBig.
  39. |
  40. | 2. If |X| < 2**(-70), go to ExpSm.
  41. |
  42. | 3. Set y := X*log_2(10)*64 (base 2 log of 10). Set
  43. | N := round-to-int(y). Decompose N as
  44. | N = 64(M + M') + j, j = 0,1,2,...,63.
  45. |
  46. | 4. Define r as
  47. | r := ((X - N*L1)-N*L2) * L10
  48. | where L1, L2 are the leading and trailing parts of log_10(2)/64
  49. | and L10 is the natural log of 10. Then
  50. | 10**X = 2**(M') * 2**(M) * 2**(j/64) * exp(r).
  51. | Go to expr to compute that expression.
  52. |
  53. | expr
  54. | 1. Fetch 2**(j/64) from table as Fact1 and Fact2.
  55. |
  56. | 2. Overwrite Fact1 and Fact2 by
  57. | Fact1 := 2**(M) * Fact1
  58. | Fact2 := 2**(M) * Fact2
  59. | Thus Fact1 + Fact2 = 2**(M) * 2**(j/64).
  60. |
  61. | 3. Calculate P where 1 + P approximates exp(r):
  62. | P = r + r*r*(A1+r*(A2+...+r*A5)).
  63. |
  64. | 4. Let AdjFact := 2**(M'). Return
  65. | AdjFact * ( Fact1 + ((Fact1*P) + Fact2) ).
  66. | Exit.
  67. |
  68. | ExpBig
  69. | 1. Generate overflow by Huge * Huge if X > 0; otherwise, generate
  70. | underflow by Tiny * Tiny.
  71. |
  72. | ExpSm
  73. | 1. Return 1 + X.
  74. |
  75. | Copyright (C) Motorola, Inc. 1990
  76. | All Rights Reserved
  77. |
  78. | For details on the license for this file, please see the
  79. | file, README, in this same directory.
  80. |STWOTOX idnt 2,1 | Motorola 040 Floating Point Software Package
  81. |section 8
  82. #include "fpsp.h"
  83. BOUNDS1: .long 0x3FB98000,0x400D80C0 | ... 2^(-70),16480
  84. BOUNDS2: .long 0x3FB98000,0x400B9B07 | ... 2^(-70),16480 LOG2/LOG10
  85. L2TEN64: .long 0x406A934F,0x0979A371 | ... 64LOG10/LOG2
  86. L10TWO1: .long 0x3F734413,0x509F8000 | ... LOG2/64LOG10
  87. L10TWO2: .long 0xBFCD0000,0xC0219DC1,0xDA994FD2,0x00000000
  88. LOG10: .long 0x40000000,0x935D8DDD,0xAAA8AC17,0x00000000
  89. LOG2: .long 0x3FFE0000,0xB17217F7,0xD1CF79AC,0x00000000
  90. EXPA5: .long 0x3F56C16D,0x6F7BD0B2
  91. EXPA4: .long 0x3F811112,0x302C712C
  92. EXPA3: .long 0x3FA55555,0x55554CC1
  93. EXPA2: .long 0x3FC55555,0x55554A54
  94. EXPA1: .long 0x3FE00000,0x00000000,0x00000000,0x00000000
  95. HUGE: .long 0x7FFE0000,0xFFFFFFFF,0xFFFFFFFF,0x00000000
  96. TINY: .long 0x00010000,0xFFFFFFFF,0xFFFFFFFF,0x00000000
  97. EXPTBL:
  98. .long 0x3FFF0000,0x80000000,0x00000000,0x3F738000
  99. .long 0x3FFF0000,0x8164D1F3,0xBC030773,0x3FBEF7CA
  100. .long 0x3FFF0000,0x82CD8698,0xAC2BA1D7,0x3FBDF8A9
  101. .long 0x3FFF0000,0x843A28C3,0xACDE4046,0x3FBCD7C9
  102. .long 0x3FFF0000,0x85AAC367,0xCC487B15,0xBFBDE8DA
  103. .long 0x3FFF0000,0x871F6196,0x9E8D1010,0x3FBDE85C
  104. .long 0x3FFF0000,0x88980E80,0x92DA8527,0x3FBEBBF1
  105. .long 0x3FFF0000,0x8A14D575,0x496EFD9A,0x3FBB80CA
  106. .long 0x3FFF0000,0x8B95C1E3,0xEA8BD6E7,0xBFBA8373
  107. .long 0x3FFF0000,0x8D1ADF5B,0x7E5BA9E6,0xBFBE9670
  108. .long 0x3FFF0000,0x8EA4398B,0x45CD53C0,0x3FBDB700
  109. .long 0x3FFF0000,0x9031DC43,0x1466B1DC,0x3FBEEEB0
  110. .long 0x3FFF0000,0x91C3D373,0xAB11C336,0x3FBBFD6D
  111. .long 0x3FFF0000,0x935A2B2F,0x13E6E92C,0xBFBDB319
  112. .long 0x3FFF0000,0x94F4EFA8,0xFEF70961,0x3FBDBA2B
  113. .long 0x3FFF0000,0x96942D37,0x20185A00,0x3FBE91D5
  114. .long 0x3FFF0000,0x9837F051,0x8DB8A96F,0x3FBE8D5A
  115. .long 0x3FFF0000,0x99E04593,0x20B7FA65,0xBFBCDE7B
  116. .long 0x3FFF0000,0x9B8D39B9,0xD54E5539,0xBFBEBAAF
  117. .long 0x3FFF0000,0x9D3ED9A7,0x2CFFB751,0xBFBD86DA
  118. .long 0x3FFF0000,0x9EF53260,0x91A111AE,0xBFBEBEDD
  119. .long 0x3FFF0000,0xA0B0510F,0xB9714FC2,0x3FBCC96E
  120. .long 0x3FFF0000,0xA2704303,0x0C496819,0xBFBEC90B
  121. .long 0x3FFF0000,0xA43515AE,0x09E6809E,0x3FBBD1DB
  122. .long 0x3FFF0000,0xA5FED6A9,0xB15138EA,0x3FBCE5EB
  123. .long 0x3FFF0000,0xA7CD93B4,0xE965356A,0xBFBEC274
  124. .long 0x3FFF0000,0xA9A15AB4,0xEA7C0EF8,0x3FBEA83C
  125. .long 0x3FFF0000,0xAB7A39B5,0xA93ED337,0x3FBECB00
  126. .long 0x3FFF0000,0xAD583EEA,0x42A14AC6,0x3FBE9301
  127. .long 0x3FFF0000,0xAF3B78AD,0x690A4375,0xBFBD8367
  128. .long 0x3FFF0000,0xB123F581,0xD2AC2590,0xBFBEF05F
  129. .long 0x3FFF0000,0xB311C412,0xA9112489,0x3FBDFB3C
  130. .long 0x3FFF0000,0xB504F333,0xF9DE6484,0x3FBEB2FB
  131. .long 0x3FFF0000,0xB6FD91E3,0x28D17791,0x3FBAE2CB
  132. .long 0x3FFF0000,0xB8FBAF47,0x62FB9EE9,0x3FBCDC3C
  133. .long 0x3FFF0000,0xBAFF5AB2,0x133E45FB,0x3FBEE9AA
  134. .long 0x3FFF0000,0xBD08A39F,0x580C36BF,0xBFBEAEFD
  135. .long 0x3FFF0000,0xBF1799B6,0x7A731083,0xBFBCBF51
  136. .long 0x3FFF0000,0xC12C4CCA,0x66709456,0x3FBEF88A
  137. .long 0x3FFF0000,0xC346CCDA,0x24976407,0x3FBD83B2
  138. .long 0x3FFF0000,0xC5672A11,0x5506DADD,0x3FBDF8AB
  139. .long 0x3FFF0000,0xC78D74C8,0xABB9B15D,0xBFBDFB17
  140. .long 0x3FFF0000,0xC9B9BD86,0x6E2F27A3,0xBFBEFE3C
  141. .long 0x3FFF0000,0xCBEC14FE,0xF2727C5D,0xBFBBB6F8
  142. .long 0x3FFF0000,0xCE248C15,0x1F8480E4,0xBFBCEE53
  143. .long 0x3FFF0000,0xD06333DA,0xEF2B2595,0xBFBDA4AE
  144. .long 0x3FFF0000,0xD2A81D91,0xF12AE45A,0x3FBC9124
  145. .long 0x3FFF0000,0xD4F35AAB,0xCFEDFA1F,0x3FBEB243
  146. .long 0x3FFF0000,0xD744FCCA,0xD69D6AF4,0x3FBDE69A
  147. .long 0x3FFF0000,0xD99D15C2,0x78AFD7B6,0xBFB8BC61
  148. .long 0x3FFF0000,0xDBFBB797,0xDAF23755,0x3FBDF610
  149. .long 0x3FFF0000,0xDE60F482,0x5E0E9124,0xBFBD8BE1
  150. .long 0x3FFF0000,0xE0CCDEEC,0x2A94E111,0x3FBACB12
  151. .long 0x3FFF0000,0xE33F8972,0xBE8A5A51,0x3FBB9BFE
  152. .long 0x3FFF0000,0xE5B906E7,0x7C8348A8,0x3FBCF2F4
  153. .long 0x3FFF0000,0xE8396A50,0x3C4BDC68,0x3FBEF22F
  154. .long 0x3FFF0000,0xEAC0C6E7,0xDD24392F,0xBFBDBF4A
  155. .long 0x3FFF0000,0xED4F301E,0xD9942B84,0x3FBEC01A
  156. .long 0x3FFF0000,0xEFE4B99B,0xDCDAF5CB,0x3FBE8CAC
  157. .long 0x3FFF0000,0xF281773C,0x59FFB13A,0xBFBCBB3F
  158. .long 0x3FFF0000,0xF5257D15,0x2486CC2C,0x3FBEF73A
  159. .long 0x3FFF0000,0xF7D0DF73,0x0AD13BB9,0xBFB8B795
  160. .long 0x3FFF0000,0xFA83B2DB,0x722A033A,0x3FBEF84B
  161. .long 0x3FFF0000,0xFD3E0C0C,0xF486C175,0xBFBEF581
  162. .set N,L_SCR1
  163. .set X,FP_SCR1
  164. .set XDCARE,X+2
  165. .set XFRAC,X+4
  166. .set ADJFACT,FP_SCR2
  167. .set FACT1,FP_SCR3
  168. .set FACT1HI,FACT1+4
  169. .set FACT1LOW,FACT1+8
  170. .set FACT2,FP_SCR4
  171. .set FACT2HI,FACT2+4
  172. .set FACT2LOW,FACT2+8
  173. | xref t_unfl
  174. |xref t_ovfl
  175. |xref t_frcinx
  176. .global stwotoxd
  177. stwotoxd:
  178. |--ENTRY POINT FOR 2**(X) FOR DENORMALIZED ARGUMENT
  179. fmovel %d1,%fpcr | ...set user's rounding mode/precision
  180. fmoves #0x3F800000,%fp0 | ...RETURN 1 + X
  181. movel (%a0),%d0
  182. orl #0x00800001,%d0
  183. fadds %d0,%fp0
  184. bra t_frcinx
  185. .global stwotox
  186. stwotox:
  187. |--ENTRY POINT FOR 2**(X), HERE X IS FINITE, NON-ZERO, AND NOT NAN'S
  188. fmovemx (%a0),%fp0-%fp0 | ...LOAD INPUT, do not set cc's
  189. movel (%a0),%d0
  190. movew 4(%a0),%d0
  191. fmovex %fp0,X(%a6)
  192. andil #0x7FFFFFFF,%d0
  193. cmpil #0x3FB98000,%d0 | ...|X| >= 2**(-70)?
  194. bges TWOOK1
  195. bra EXPBORS
  196. TWOOK1:
  197. cmpil #0x400D80C0,%d0 | ...|X| > 16480?
  198. bles TWOMAIN
  199. bra EXPBORS
  200. TWOMAIN:
  201. |--USUAL CASE, 2^(-70) <= |X| <= 16480
  202. fmovex %fp0,%fp1
  203. fmuls #0x42800000,%fp1 | ...64 * X
  204. fmovel %fp1,N(%a6) | ...N = ROUND-TO-INT(64 X)
  205. movel %d2,-(%sp)
  206. lea EXPTBL,%a1 | ...LOAD ADDRESS OF TABLE OF 2^(J/64)
  207. fmovel N(%a6),%fp1 | ...N --> FLOATING FMT
  208. movel N(%a6),%d0
  209. movel %d0,%d2
  210. andil #0x3F,%d0 | ...D0 IS J
  211. asll #4,%d0 | ...DISPLACEMENT FOR 2^(J/64)
  212. addal %d0,%a1 | ...ADDRESS FOR 2^(J/64)
  213. asrl #6,%d2 | ...d2 IS L, N = 64L + J
  214. movel %d2,%d0
  215. asrl #1,%d0 | ...D0 IS M
  216. subl %d0,%d2 | ...d2 IS M', N = 64(M+M') + J
  217. addil #0x3FFF,%d2
  218. movew %d2,ADJFACT(%a6) | ...ADJFACT IS 2^(M')
  219. movel (%sp)+,%d2
  220. |--SUMMARY: a1 IS ADDRESS FOR THE LEADING PORTION OF 2^(J/64),
  221. |--D0 IS M WHERE N = 64(M+M') + J. NOTE THAT |M| <= 16140 BY DESIGN.
  222. |--ADJFACT = 2^(M').
  223. |--REGISTERS SAVED SO FAR ARE (IN ORDER) FPCR, D0, FP1, a1, AND FP2.
  224. fmuls #0x3C800000,%fp1 | ...(1/64)*N
  225. movel (%a1)+,FACT1(%a6)
  226. movel (%a1)+,FACT1HI(%a6)
  227. movel (%a1)+,FACT1LOW(%a6)
  228. movew (%a1)+,FACT2(%a6)
  229. clrw FACT2+2(%a6)
  230. fsubx %fp1,%fp0 | ...X - (1/64)*INT(64 X)
  231. movew (%a1)+,FACT2HI(%a6)
  232. clrw FACT2HI+2(%a6)
  233. clrl FACT2LOW(%a6)
  234. addw %d0,FACT1(%a6)
  235. fmulx LOG2,%fp0 | ...FP0 IS R
  236. addw %d0,FACT2(%a6)
  237. bra expr
  238. EXPBORS:
  239. |--FPCR, D0 SAVED
  240. cmpil #0x3FFF8000,%d0
  241. bgts EXPBIG
  242. EXPSM:
  243. |--|X| IS SMALL, RETURN 1 + X
  244. fmovel %d1,%FPCR |restore users exceptions
  245. fadds #0x3F800000,%fp0 | ...RETURN 1 + X
  246. bra t_frcinx
  247. EXPBIG:
  248. |--|X| IS LARGE, GENERATE OVERFLOW IF X > 0; ELSE GENERATE UNDERFLOW
  249. |--REGISTERS SAVE SO FAR ARE FPCR AND D0
  250. movel X(%a6),%d0
  251. cmpil #0,%d0
  252. blts EXPNEG
  253. bclrb #7,(%a0) |t_ovfl expects positive value
  254. bra t_ovfl
  255. EXPNEG:
  256. bclrb #7,(%a0) |t_unfl expects positive value
  257. bra t_unfl
  258. .global stentoxd
  259. stentoxd:
  260. |--ENTRY POINT FOR 10**(X) FOR DENORMALIZED ARGUMENT
  261. fmovel %d1,%fpcr | ...set user's rounding mode/precision
  262. fmoves #0x3F800000,%fp0 | ...RETURN 1 + X
  263. movel (%a0),%d0
  264. orl #0x00800001,%d0
  265. fadds %d0,%fp0
  266. bra t_frcinx
  267. .global stentox
  268. stentox:
  269. |--ENTRY POINT FOR 10**(X), HERE X IS FINITE, NON-ZERO, AND NOT NAN'S
  270. fmovemx (%a0),%fp0-%fp0 | ...LOAD INPUT, do not set cc's
  271. movel (%a0),%d0
  272. movew 4(%a0),%d0
  273. fmovex %fp0,X(%a6)
  274. andil #0x7FFFFFFF,%d0
  275. cmpil #0x3FB98000,%d0 | ...|X| >= 2**(-70)?
  276. bges TENOK1
  277. bra EXPBORS
  278. TENOK1:
  279. cmpil #0x400B9B07,%d0 | ...|X| <= 16480*log2/log10 ?
  280. bles TENMAIN
  281. bra EXPBORS
  282. TENMAIN:
  283. |--USUAL CASE, 2^(-70) <= |X| <= 16480 LOG 2 / LOG 10
  284. fmovex %fp0,%fp1
  285. fmuld L2TEN64,%fp1 | ...X*64*LOG10/LOG2
  286. fmovel %fp1,N(%a6) | ...N=INT(X*64*LOG10/LOG2)
  287. movel %d2,-(%sp)
  288. lea EXPTBL,%a1 | ...LOAD ADDRESS OF TABLE OF 2^(J/64)
  289. fmovel N(%a6),%fp1 | ...N --> FLOATING FMT
  290. movel N(%a6),%d0
  291. movel %d0,%d2
  292. andil #0x3F,%d0 | ...D0 IS J
  293. asll #4,%d0 | ...DISPLACEMENT FOR 2^(J/64)
  294. addal %d0,%a1 | ...ADDRESS FOR 2^(J/64)
  295. asrl #6,%d2 | ...d2 IS L, N = 64L + J
  296. movel %d2,%d0
  297. asrl #1,%d0 | ...D0 IS M
  298. subl %d0,%d2 | ...d2 IS M', N = 64(M+M') + J
  299. addil #0x3FFF,%d2
  300. movew %d2,ADJFACT(%a6) | ...ADJFACT IS 2^(M')
  301. movel (%sp)+,%d2
  302. |--SUMMARY: a1 IS ADDRESS FOR THE LEADING PORTION OF 2^(J/64),
  303. |--D0 IS M WHERE N = 64(M+M') + J. NOTE THAT |M| <= 16140 BY DESIGN.
  304. |--ADJFACT = 2^(M').
  305. |--REGISTERS SAVED SO FAR ARE (IN ORDER) FPCR, D0, FP1, a1, AND FP2.
  306. fmovex %fp1,%fp2
  307. fmuld L10TWO1,%fp1 | ...N*(LOG2/64LOG10)_LEAD
  308. movel (%a1)+,FACT1(%a6)
  309. fmulx L10TWO2,%fp2 | ...N*(LOG2/64LOG10)_TRAIL
  310. movel (%a1)+,FACT1HI(%a6)
  311. movel (%a1)+,FACT1LOW(%a6)
  312. fsubx %fp1,%fp0 | ...X - N L_LEAD
  313. movew (%a1)+,FACT2(%a6)
  314. fsubx %fp2,%fp0 | ...X - N L_TRAIL
  315. clrw FACT2+2(%a6)
  316. movew (%a1)+,FACT2HI(%a6)
  317. clrw FACT2HI+2(%a6)
  318. clrl FACT2LOW(%a6)
  319. fmulx LOG10,%fp0 | ...FP0 IS R
  320. addw %d0,FACT1(%a6)
  321. addw %d0,FACT2(%a6)
  322. expr:
  323. |--FPCR, FP2, FP3 ARE SAVED IN ORDER AS SHOWN.
  324. |--ADJFACT CONTAINS 2**(M'), FACT1 + FACT2 = 2**(M) * 2**(J/64).
  325. |--FP0 IS R. THE FOLLOWING CODE COMPUTES
  326. |-- 2**(M'+M) * 2**(J/64) * EXP(R)
  327. fmovex %fp0,%fp1
  328. fmulx %fp1,%fp1 | ...FP1 IS S = R*R
  329. fmoved EXPA5,%fp2 | ...FP2 IS A5
  330. fmoved EXPA4,%fp3 | ...FP3 IS A4
  331. fmulx %fp1,%fp2 | ...FP2 IS S*A5
  332. fmulx %fp1,%fp3 | ...FP3 IS S*A4
  333. faddd EXPA3,%fp2 | ...FP2 IS A3+S*A5
  334. faddd EXPA2,%fp3 | ...FP3 IS A2+S*A4
  335. fmulx %fp1,%fp2 | ...FP2 IS S*(A3+S*A5)
  336. fmulx %fp1,%fp3 | ...FP3 IS S*(A2+S*A4)
  337. faddd EXPA1,%fp2 | ...FP2 IS A1+S*(A3+S*A5)
  338. fmulx %fp0,%fp3 | ...FP3 IS R*S*(A2+S*A4)
  339. fmulx %fp1,%fp2 | ...FP2 IS S*(A1+S*(A3+S*A5))
  340. faddx %fp3,%fp0 | ...FP0 IS R+R*S*(A2+S*A4)
  341. faddx %fp2,%fp0 | ...FP0 IS EXP(R) - 1
  342. |--FINAL RECONSTRUCTION PROCESS
  343. |--EXP(X) = 2^M*2^(J/64) + 2^M*2^(J/64)*(EXP(R)-1) - (1 OR 0)
  344. fmulx FACT1(%a6),%fp0
  345. faddx FACT2(%a6),%fp0
  346. faddx FACT1(%a6),%fp0
  347. fmovel %d1,%FPCR |restore users exceptions
  348. clrw ADJFACT+2(%a6)
  349. movel #0x80000000,ADJFACT+4(%a6)
  350. clrl ADJFACT+8(%a6)
  351. fmulx ADJFACT(%a6),%fp0 | ...FINAL ADJUSTMENT
  352. bra t_frcinx
  353. |end