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

源码开发语言/平台
当前位置：首页 > 源码/资料 > Linux/Unix编程 > 查看源码
slogn.S
资源名称：linux-2.4.20.tar.gz　[点击查看]
上传用户：jlfgdled
上传日期：2013-04-10
资源大小：33168k
文件大小：19k
源码类别：

Linux/Unix编程

开发平台：

Unix_Linux

slogn.S：源码内容
  1. |
  2. | slogn.sa 3.1 12/10/90
  3. |
  4. | slogn computes the natural logarithm of an
  5. | input value. slognd does the same except the input value is a
  6. | denormalized number. slognp1 computes log(1+X), and slognp1d
  7. | computes log(1+X) for denormalized X.
  8. |
  9. | Input: Double-extended value in memory location pointed to by address
  10. | register a0.
  11. |
  12. | Output: log(X) or log(1+X) returned in floating-point register 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 slogn takes approximately 190 cycles for input 
  20. | argument X such that |X-1| >= 1/16, which is the usual 
  21. | situation. For those arguments, slognp1 takes approximately
  22. |  210 cycles. For the less common arguments, the program will
  23. |  run no worse than 10% slower.
  24. |
  25. | Algorithm:
  26. | LOGN:
  27. | Step 1. If |X-1| < 1/16, approximate log(X) by an odd polynomial in
  28. | u, where u = 2(X-1)/(X+1). Otherwise, move on to Step 2.
  29. |
  30. | Step 2. X = 2**k * Y where 1 <= Y < 2. Define F to be the first seven
  31. | significant bits of Y plus 2**(-7), i.e. F = 1.xxxxxx1 in base
  32. | 2 where the six "x" match those of Y. Note that |Y-F| <= 2**(-7).
  33. |
  34. | Step 3. Define u = (Y-F)/F. Approximate log(1+u) by a polynomial in u,
  35. | log(1+u) = poly.
  36. |
  37. | Step 4. Reconstruct log(X) = log( 2**k * Y ) = k*log(2) + log(F) + log(1+u)
  38. | by k*log(2) + (log(F) + poly). The values of log(F) are calculated
  39. | beforehand and stored in the program.
  40. |
  41. | lognp1:
  42. | Step 1: If |X| < 1/16, approximate log(1+X) by an odd polynomial in
  43. | u where u = 2X/(2+X). Otherwise, move on to Step 2.
  44. |
  45. | Step 2: Let 1+X = 2**k * Y, where 1 <= Y < 2. Define F as done in Step 2
  46. | of the algorithm for LOGN and compute log(1+X) as
  47. | k*log(2) + log(F) + poly where poly approximates log(1+u),
  48. | u = (Y-F)/F. 
  49. |
  50. | Implementation Notes:
  51. | Note 1. There are 64 different possible values for F, thus 64 log(F)'s
  52. | need to be tabulated. Moreover, the values of 1/F are also 
  53. | tabulated so that the division in (Y-F)/F can be performed by a
  54. | multiplication.
  55. |
  56. | Note 2. In Step 2 of lognp1, in order to preserved accuracy, the value
  57. | Y-F has to be calculated carefully when 1/2 <= X < 3/2. 
  58. |
  59. | Note 3. To fully exploit the pipeline, polynomials are usually separated
  60. | into two parts evaluated independently before being added up.
  61. |
  63. | Copyright (C) Motorola, Inc. 1990
  64. | All Rights Reserved
  65. |
  66. | THIS IS UNPUBLISHED PROPRIETARY SOURCE CODE OF MOTOROLA 
  67. | The copyright notice above does not evidence any  
  68. | actual or intended publication of such source code.
  70. |slogn idnt 2,1 | Motorola 040 Floating Point Software Package
  72. |section 8
  74. .include "fpsp.h"
  76. BOUNDS1:  .long 0x3FFEF07D,0x3FFF8841
  77. BOUNDS2:  .long 0x3FFE8000,0x3FFFC000
  79. LOGOF2: .long 0x3FFE0000,0xB17217F7,0xD1CF79AC,0x00000000
  81. one: .long 0x3F800000
  82. zero: .long 0x00000000
  83. infty: .long 0x7F800000
  84. negone: .long 0xBF800000
  86. LOGA6: .long 0x3FC2499A,0xB5E4040B
  87. LOGA5: .long 0xBFC555B5,0x848CB7DB
  89. LOGA4: .long 0x3FC99999,0x987D8730
  90. LOGA3: .long 0xBFCFFFFF,0xFF6F7E97
  92. LOGA2: .long 0x3FD55555,0x555555a4
  93. LOGA1: .long 0xBFE00000,0x00000008
  95. LOGB5: .long 0x3F175496,0xADD7DAD6
  96. LOGB4: .long 0x3F3C71C2,0xFE80C7E0
  98. LOGB3: .long 0x3F624924,0x928BCCFF
  99. LOGB2: .long 0x3F899999,0x999995EC
  101. LOGB1: .long 0x3FB55555,0x55555555
  102. TWO: .long 0x40000000,0x00000000
  104. LTHOLD: .long 0x3f990000,0x80000000,0x00000000,0x00000000
  106. LOGTBL:
  107. .long  0x3FFE0000,0xFE03F80F,0xE03F80FE,0x00000000
  108. .long  0x3FF70000,0xFF015358,0x833C47E2,0x00000000
  109. .long  0x3FFE0000,0xFA232CF2,0x52138AC0,0x00000000
  110. .long  0x3FF90000,0xBDC8D83E,0xAD88D549,0x00000000
  111. .long  0x3FFE0000,0xF6603D98,0x0F6603DA,0x00000000
  112. .long  0x3FFA0000,0x9CF43DCF,0xF5EAFD48,0x00000000
  113. .long  0x3FFE0000,0xF2B9D648,0x0F2B9D65,0x00000000
  114. .long  0x3FFA0000,0xDA16EB88,0xCB8DF614,0x00000000
  115. .long  0x3FFE0000,0xEF2EB71F,0xC4345238,0x00000000
  116. .long  0x3FFB0000,0x8B29B775,0x1BD70743,0x00000000
  117. .long  0x3FFE0000,0xEBBDB2A5,0xC1619C8C,0x00000000
  118. .long  0x3FFB0000,0xA8D839F8,0x30C1FB49,0x00000000
  119. .long  0x3FFE0000,0xE865AC7B,0x7603A197,0x00000000
  120. .long  0x3FFB0000,0xC61A2EB1,0x8CD907AD,0x00000000
  121. .long  0x3FFE0000,0xE525982A,0xF70C880E,0x00000000
  122. .long  0x3FFB0000,0xE2F2A47A,0xDE3A18AF,0x00000000
  123. .long  0x3FFE0000,0xE1FC780E,0x1FC780E2,0x00000000
  124. .long  0x3FFB0000,0xFF64898E,0xDF55D551,0x00000000
  125. .long  0x3FFE0000,0xDEE95C4C,0xA037BA57,0x00000000
  126. .long  0x3FFC0000,0x8DB956A9,0x7B3D0148,0x00000000
  127. .long  0x3FFE0000,0xDBEB61EE,0xD19C5958,0x00000000
  128. .long  0x3FFC0000,0x9B8FE100,0xF47BA1DE,0x00000000
  129. .long  0x3FFE0000,0xD901B203,0x6406C80E,0x00000000
  130. .long  0x3FFC0000,0xA9372F1D,0x0DA1BD17,0x00000000
  131. .long  0x3FFE0000,0xD62B80D6,0x2B80D62C,0x00000000
  132. .long  0x3FFC0000,0xB6B07F38,0xCE90E46B,0x00000000
  133. .long  0x3FFE0000,0xD3680D36,0x80D3680D,0x00000000
  134. .long  0x3FFC0000,0xC3FD0329,0x06488481,0x00000000
  135. .long  0x3FFE0000,0xD0B69FCB,0xD2580D0B,0x00000000
  136. .long  0x3FFC0000,0xD11DE0FF,0x15AB18CA,0x00000000
  137. .long  0x3FFE0000,0xCE168A77,0x25080CE1,0x00000000
  138. .long  0x3FFC0000,0xDE1433A1,0x6C66B150,0x00000000
  139. .long  0x3FFE0000,0xCB8727C0,0x65C393E0,0x00000000
  140. .long  0x3FFC0000,0xEAE10B5A,0x7DDC8ADD,0x00000000
  141. .long  0x3FFE0000,0xC907DA4E,0x871146AD,0x00000000
  142. .long  0x3FFC0000,0xF7856E5E,0xE2C9B291,0x00000000
  143. .long  0x3FFE0000,0xC6980C69,0x80C6980C,0x00000000
  144. .long  0x3FFD0000,0x82012CA5,0xA68206D7,0x00000000
  145. .long  0x3FFE0000,0xC4372F85,0x5D824CA6,0x00000000
  146. .long  0x3FFD0000,0x882C5FCD,0x7256A8C5,0x00000000
  147. .long  0x3FFE0000,0xC1E4BBD5,0x95F6E947,0x00000000
  148. .long  0x3FFD0000,0x8E44C60B,0x4CCFD7DE,0x00000000
  149. .long  0x3FFE0000,0xBFA02FE8,0x0BFA02FF,0x00000000
  150. .long  0x3FFD0000,0x944AD09E,0xF4351AF6,0x00000000
  151. .long  0x3FFE0000,0xBD691047,0x07661AA3,0x00000000
  152. .long  0x3FFD0000,0x9A3EECD4,0xC3EAA6B2,0x00000000
  153. .long  0x3FFE0000,0xBB3EE721,0xA54D880C,0x00000000
  154. .long  0x3FFD0000,0xA0218434,0x353F1DE8,0x00000000
  155. .long  0x3FFE0000,0xB92143FA,0x36F5E02E,0x00000000
  156. .long  0x3FFD0000,0xA5F2FCAB,0xBBC506DA,0x00000000
  157. .long  0x3FFE0000,0xB70FBB5A,0x19BE3659,0x00000000
  158. .long  0x3FFD0000,0xABB3B8BA,0x2AD362A5,0x00000000
  159. .long  0x3FFE0000,0xB509E68A,0x9B94821F,0x00000000
  160. .long  0x3FFD0000,0xB1641795,0xCE3CA97B,0x00000000
  161. .long  0x3FFE0000,0xB30F6352,0x8917C80B,0x00000000
  162. .long  0x3FFD0000,0xB7047551,0x5D0F1C61,0x00000000
  163. .long  0x3FFE0000,0xB11FD3B8,0x0B11FD3C,0x00000000
  164. .long  0x3FFD0000,0xBC952AFE,0xEA3D13E1,0x00000000
  165. .long  0x3FFE0000,0xAF3ADDC6,0x80AF3ADE,0x00000000
  166. .long  0x3FFD0000,0xC2168ED0,0xF458BA4A,0x00000000
  167. .long  0x3FFE0000,0xAD602B58,0x0AD602B6,0x00000000
  168. .long  0x3FFD0000,0xC788F439,0xB3163BF1,0x00000000
  169. .long  0x3FFE0000,0xAB8F69E2,0x8359CD11,0x00000000
  170. .long  0x3FFD0000,0xCCECAC08,0xBF04565D,0x00000000
  171. .long  0x3FFE0000,0xA9C84A47,0xA07F5638,0x00000000
  172. .long  0x3FFD0000,0xD2420487,0x2DD85160,0x00000000
  173. .long  0x3FFE0000,0xA80A80A8,0x0A80A80B,0x00000000
  174. .long  0x3FFD0000,0xD7894992,0x3BC3588A,0x00000000
  175. .long  0x3FFE0000,0xA655C439,0x2D7B73A8,0x00000000
  176. .long  0x3FFD0000,0xDCC2C4B4,0x9887DACC,0x00000000
  177. .long  0x3FFE0000,0xA4A9CF1D,0x96833751,0x00000000
  178. .long  0x3FFD0000,0xE1EEBD3E,0x6D6A6B9E,0x00000000
  179. .long  0x3FFE0000,0xA3065E3F,0xAE7CD0E0,0x00000000
  180. .long  0x3FFD0000,0xE70D785C,0x2F9F5BDC,0x00000000
  181. .long  0x3FFE0000,0xA16B312E,0xA8FC377D,0x00000000
  182. .long  0x3FFD0000,0xEC1F392C,0x5179F283,0x00000000
  183. .long  0x3FFE0000,0x9FD809FD,0x809FD80A,0x00000000
  184. .long  0x3FFD0000,0xF12440D3,0xE36130E6,0x00000000
  185. .long  0x3FFE0000,0x9E4CAD23,0xDD5F3A20,0x00000000
  186. .long  0x3FFD0000,0xF61CCE92,0x346600BB,0x00000000
  187. .long  0x3FFE0000,0x9CC8E160,0xC3FB19B9,0x00000000
  188. .long  0x3FFD0000,0xFB091FD3,0x8145630A,0x00000000
  189. .long  0x3FFE0000,0x9B4C6F9E,0xF03A3CAA,0x00000000
  190. .long  0x3FFD0000,0xFFE97042,0xBFA4C2AD,0x00000000
  191. .long  0x3FFE0000,0x99D722DA,0xBDE58F06,0x00000000
  192. .long  0x3FFE0000,0x825EFCED,0x49369330,0x00000000
  193. .long  0x3FFE0000,0x9868C809,0x868C8098,0x00000000
  194. .long  0x3FFE0000,0x84C37A7A,0xB9A905C9,0x00000000
  195. .long  0x3FFE0000,0x97012E02,0x5C04B809,0x00000000
  196. .long  0x3FFE0000,0x87224C2E,0x8E645FB7,0x00000000
  197. .long  0x3FFE0000,0x95A02568,0x095A0257,0x00000000
  198. .long  0x3FFE0000,0x897B8CAC,0x9F7DE298,0x00000000
  199. .long  0x3FFE0000,0x94458094,0x45809446,0x00000000
  200. .long  0x3FFE0000,0x8BCF55DE,0xC4CD05FE,0x00000000
  201. .long  0x3FFE0000,0x92F11384,0x0497889C,0x00000000
  202. .long  0x3FFE0000,0x8E1DC0FB,0x89E125E5,0x00000000
  203. .long  0x3FFE0000,0x91A2B3C4,0xD5E6F809,0x00000000
  204. .long  0x3FFE0000,0x9066E68C,0x955B6C9B,0x00000000
  205. .long  0x3FFE0000,0x905A3863,0x3E06C43B,0x00000000
  206. .long  0x3FFE0000,0x92AADE74,0xC7BE59E0,0x00000000
  207. .long  0x3FFE0000,0x8F1779D9,0xFDC3A219,0x00000000
  208. .long  0x3FFE0000,0x94E9BFF6,0x15845643,0x00000000
  209. .long  0x3FFE0000,0x8DDA5202,0x37694809,0x00000000
  210. .long  0x3FFE0000,0x9723A1B7,0x20134203,0x00000000
  211. .long  0x3FFE0000,0x8CA29C04,0x6514E023,0x00000000
  212. .long  0x3FFE0000,0x995899C8,0x90EB8990,0x00000000
  213. .long  0x3FFE0000,0x8B70344A,0x139BC75A,0x00000000
  214. .long  0x3FFE0000,0x9B88BDAA,0x3A3DAE2F,0x00000000
  215. .long  0x3FFE0000,0x8A42F870,0x5669DB46,0x00000000
  216. .long  0x3FFE0000,0x9DB4224F,0xFFE1157C,0x00000000
  217. .long  0x3FFE0000,0x891AC73A,0xE9819B50,0x00000000
  218. .long  0x3FFE0000,0x9FDADC26,0x8B7A12DA,0x00000000
  219. .long  0x3FFE0000,0x87F78087,0xF78087F8,0x00000000
  220. .long  0x3FFE0000,0xA1FCFF17,0xCE733BD4,0x00000000
  221. .long  0x3FFE0000,0x86D90544,0x7A34ACC6,0x00000000
  222. .long  0x3FFE0000,0xA41A9E8F,0x5446FB9F,0x00000000
  223. .long  0x3FFE0000,0x85BF3761,0x2CEE3C9B,0x00000000
  224. .long  0x3FFE0000,0xA633CD7E,0x6771CD8B,0x00000000
  225. .long  0x3FFE0000,0x84A9F9C8,0x084A9F9D,0x00000000
  226. .long  0x3FFE0000,0xA8489E60,0x0B435A5E,0x00000000
  227. .long  0x3FFE0000,0x83993052,0x3FBE3368,0x00000000
  228. .long  0x3FFE0000,0xAA59233C,0xCCA4BD49,0x00000000
  229. .long  0x3FFE0000,0x828CBFBE,0xB9A020A3,0x00000000
  230. .long  0x3FFE0000,0xAC656DAE,0x6BCC4985,0x00000000
  231. .long  0x3FFE0000,0x81848DA8,0xFAF0D277,0x00000000
  232. .long  0x3FFE0000,0xAE6D8EE3,0x60BB2468,0x00000000
  233. .long  0x3FFE0000,0x80808080,0x80808081,0x00000000
  234. .long  0x3FFE0000,0xB07197A2,0x3C46C654,0x00000000
  236. .set ADJK,L_SCR1
  238. .set X,FP_SCR1
  239. .set XDCARE,X+2
  240. .set XFRAC,X+4
  242. .set F,FP_SCR2
  243. .set FFRAC,F+4
  245. .set KLOG2,FP_SCR3
  247. .set SAVEU,FP_SCR4
  249. | xref t_frcinx
  250. |xref t_extdnrm
  251. |xref t_operr
  252. |xref t_dz
  254. .global slognd
  255. slognd:
  256. |--ENTRY POINT FOR LOG(X) FOR DENORMALIZED INPUT
  258. movel #-100,ADJK(%a6) | ...INPUT = 2^(ADJK) * FP0
  260. |----normalize the input value by left shifting k bits (k to be determined
  261. |----below), adjusting exponent and storing -k to  ADJK
  262. |----the value TWOTO100 is no longer needed.
  263. |----Note that this code assumes the denormalized input is NON-ZERO.
  265.      moveml %d2-%d7,-(%a7) | ...save some registers 
  266.      movel #0x00000000,%d3 | ...D3 is exponent of smallest norm. #
  267.      movel 4(%a0),%d4
  268.      movel 8(%a0),%d5 | ...(D4,D5) is (Hi_X,Lo_X)
  269.      clrl %d2 | ...D2 used for holding K
  271.      tstl %d4
  272.      bnes HiX_not0
  274. HiX_0:
  275.      movel %d5,%d4
  276.      clrl %d5
  277.      movel #32,%d2
  278.      clrl %d6
  279.      bfffo      %d4{#0:#32},%d6
  280.      lsll      %d6,%d4
  281.      addl %d6,%d2 | ...(D3,D4,D5) is normalized
  283.      movel %d3,X(%a6)
  284.      movel %d4,XFRAC(%a6)
  285.      movel %d5,XFRAC+4(%a6)
  286.      negl %d2
  287.      movel %d2,ADJK(%a6)
  288.      fmovex X(%a6),%fp0
  289.      moveml (%a7)+,%d2-%d7 | ...restore registers
  290.      lea X(%a6),%a0
  291.      bras LOGBGN | ...begin regular log(X)
  294. HiX_not0:
  295.      clrl %d6
  296.      bfffo %d4{#0:#32},%d6 | ...find first 1
  297.      movel %d6,%d2 | ...get k
  298.      lsll %d6,%d4
  299.      movel %d5,%d7 | ...a copy of D5
  300.      lsll %d6,%d5
  301.      negl %d6
  302.      addil #32,%d6
  303.      lsrl %d6,%d7
  304.      orl %d7,%d4 | ...(D3,D4,D5) normalized
  306.      movel %d3,X(%a6)
  307.      movel %d4,XFRAC(%a6)
  308.      movel %d5,XFRAC+4(%a6)
  309.      negl %d2
  310.      movel %d2,ADJK(%a6)
  311.      fmovex X(%a6),%fp0
  312.      moveml (%a7)+,%d2-%d7 | ...restore registers
  313.      lea X(%a6),%a0
  314.      bras LOGBGN | ...begin regular log(X)
  317. .global slogn
  318. slogn:
  319. |--ENTRY POINT FOR LOG(X) FOR X FINITE, NON-ZERO, NOT NAN'S
  321. fmovex (%a0),%fp0 | ...LOAD INPUT
  322. movel #0x00000000,ADJK(%a6)
  324. LOGBGN:
  325. |--FPCR SAVED AND CLEARED, INPUT IS 2^(ADJK)*FP0, FP0 CONTAINS
  326. |--A FINITE, NON-ZERO, NORMALIZED NUMBER.
  328. movel (%a0),%d0
  329. movew 4(%a0),%d0
  331. movel (%a0),X(%a6)
  332. movel 4(%a0),X+4(%a6)
  333. movel 8(%a0),X+8(%a6)
  335. cmpil #0,%d0 | ...CHECK IF X IS NEGATIVE
  336. blt LOGNEG | ...LOG OF NEGATIVE ARGUMENT IS INVALID
  337. cmp2l BOUNDS1,%d0 | ...X IS POSITIVE, CHECK IF X IS NEAR 1
  338. bcc LOGNEAR1 | ...BOUNDS IS ROUGHLY [15/16, 17/16]
  340. LOGMAIN:
  341. |--THIS SHOULD BE THE USUAL CASE, X NOT VERY CLOSE TO 1
  343. |--X = 2^(K) * Y, 1 <= Y < 2. THUS, Y = 1.XXXXXXXX....XX IN BINARY.
  344. |--WE DEFINE F = 1.XXXXXX1, I.E. FIRST 7 BITS OF Y AND ATTACH A 1.
  345. |--THE IDEA IS THAT LOG(X) = K*LOG2 + LOG(Y)
  346. |--  = K*LOG2 + LOG(F) + LOG(1 + (Y-F)/F).
  347. |--NOTE THAT U = (Y-F)/F IS VERY SMALL AND THUS APPROXIMATING
  348. |--LOG(1+U) CAN BE VERY EFFICIENT.
  349. |--ALSO NOTE THAT THE VALUE 1/F IS STORED IN A TABLE SO THAT NO
  350. |--DIVISION IS NEEDED TO CALCULATE (Y-F)/F. 
  352. |--GET K, Y, F, AND ADDRESS OF 1/F.
  353. asrl #8,%d0
  354. asrl #8,%d0 | ...SHIFTED 16 BITS, BIASED EXPO. OF X
  355. subil #0x3FFF,%d0  | ...THIS IS K
  356. addl ADJK(%a6),%d0 | ...ADJUST K, ORIGINAL INPUT MAY BE  DENORM.
  357. lea LOGTBL,%a0  | ...BASE ADDRESS OF 1/F AND LOG(F)
  358. fmovel %d0,%fp1 | ...CONVERT K TO FLOATING-POINT FORMAT
  360. |--WHILE THE CONVERSION IS GOING ON, WE GET F AND ADDRESS OF 1/F
  361. movel #0x3FFF0000,X(%a6) | ...X IS NOW Y, I.E. 2^(-K)*X
  362. movel XFRAC(%a6),FFRAC(%a6)
  363. andil #0xFE000000,FFRAC(%a6) | ...FIRST 7 BITS OF Y
  364. oril #0x01000000,FFRAC(%a6) | ...GET F: ATTACH A 1 AT THE EIGHTH BIT
  365. movel FFRAC(%a6),%d0 | ...READY TO GET ADDRESS OF 1/F
  366. andil #0x7E000000,%d0
  367. asrl #8,%d0
  368. asrl #8,%d0
  369. asrl #4,%d0 | ...SHIFTED 20, D0 IS THE DISPLACEMENT
  370. addal %d0,%a0 | ...A0 IS THE ADDRESS FOR 1/F
  372. fmovex X(%a6),%fp0
  373. movel #0x3fff0000,F(%a6)
  374. clrl F+8(%a6)
  375. fsubx F(%a6),%fp0 | ...Y-F
  376. fmovemx %fp2-%fp2/%fp3,-(%sp) | ...SAVE FP2 WHILE FP0 IS NOT READY
  377. |--SUMMARY: FP0 IS Y-F, A0 IS ADDRESS OF 1/F, FP1 IS K
  378. |--REGISTERS SAVED: FPCR, FP1, FP2
  380. LP1CONT1:
  381. |--AN RE-ENTRY POINT FOR LOGNP1
  382. fmulx (%a0),%fp0 | ...FP0 IS U = (Y-F)/F
  383. fmulx LOGOF2,%fp1 | ...GET K*LOG2 WHILE FP0 IS NOT READY
  384. fmovex %fp0,%fp2
  385. fmulx %fp2,%fp2 | ...FP2 IS V=U*U
  386. fmovex %fp1,KLOG2(%a6) | ...PUT K*LOG2 IN MEMORY, FREE FP1
  388. |--LOG(1+U) IS APPROXIMATED BY
  389. |--U + V*(A1+U*(A2+U*(A3+U*(A4+U*(A5+U*A6))))) WHICH IS
  390. |--[U + V*(A1+V*(A3+V*A5))]  +  [U*V*(A2+V*(A4+V*A6))]
  392. fmovex %fp2,%fp3
  393. fmovex %fp2,%fp1
  395. fmuld LOGA6,%fp1 | ...V*A6
  396. fmuld LOGA5,%fp2 | ...V*A5
  398. faddd LOGA4,%fp1 | ...A4+V*A6
  399. faddd LOGA3,%fp2 | ...A3+V*A5
  401. fmulx %fp3,%fp1 | ...V*(A4+V*A6)
  402. fmulx %fp3,%fp2 | ...V*(A3+V*A5)
  404. faddd LOGA2,%fp1 | ...A2+V*(A4+V*A6)
  405. faddd LOGA1,%fp2 | ...A1+V*(A3+V*A5)
  407. fmulx %fp3,%fp1 | ...V*(A2+V*(A4+V*A6))
  408. addal #16,%a0 | ...ADDRESS OF LOG(F)
  409. fmulx %fp3,%fp2 | ...V*(A1+V*(A3+V*A5)), FP3 RELEASED
  411. fmulx %fp0,%fp1 | ...U*V*(A2+V*(A4+V*A6))
  412. faddx %fp2,%fp0 | ...U+V*(A1+V*(A3+V*A5)), FP2 RELEASED
  414. faddx (%a0),%fp1 | ...LOG(F)+U*V*(A2+V*(A4+V*A6))
  415. fmovemx  (%sp)+,%fp2-%fp2/%fp3 | ...RESTORE FP2
  416. faddx %fp1,%fp0 | ...FP0 IS LOG(F) + LOG(1+U)
  418. fmovel %d1,%fpcr
  419. faddx KLOG2(%a6),%fp0 | ...FINAL ADD
  420. bra t_frcinx
  423. LOGNEAR1:
  424. |--REGISTERS SAVED: FPCR, FP1. FP0 CONTAINS THE INPUT.
  425. fmovex %fp0,%fp1
  426. fsubs one,%fp1 | ...FP1 IS X-1
  427. fadds one,%fp0 | ...FP0 IS X+1
  428. faddx %fp1,%fp1 | ...FP1 IS 2(X-1)
  429. |--LOG(X) = LOG(1+U/2)-LOG(1-U/2) WHICH IS AN ODD POLYNOMIAL
  430. |--IN U, U = 2(X-1)/(X+1) = FP1/FP0
  432. LP1CONT2:
  433. |--THIS IS AN RE-ENTRY POINT FOR LOGNP1
  434. fdivx %fp0,%fp1 | ...FP1 IS U
  435. fmovemx %fp2-%fp2/%fp3,-(%sp)  | ...SAVE FP2
  436. |--REGISTERS SAVED ARE NOW FPCR,FP1,FP2,FP3
  437. |--LET V=U*U, W=V*V, CALCULATE
  438. |--U + U*V*(B1 + V*(B2 + V*(B3 + V*(B4 + V*B5)))) BY
  439. |--U + U*V*(  [B1 + W*(B3 + W*B5)]  +  [V*(B2 + W*B4)]  )
  440. fmovex %fp1,%fp0
  441. fmulx %fp0,%fp0 | ...FP0 IS V
  442. fmovex %fp1,SAVEU(%a6) | ...STORE U IN MEMORY, FREE FP1
  443. fmovex %fp0,%fp1
  444. fmulx %fp1,%fp1 | ...FP1 IS W
  446. fmoved LOGB5,%fp3
  447. fmoved LOGB4,%fp2
  449. fmulx %fp1,%fp3 | ...W*B5
  450. fmulx %fp1,%fp2 | ...W*B4
  452. faddd LOGB3,%fp3 | ...B3+W*B5
  453. faddd LOGB2,%fp2 | ...B2+W*B4
  455. fmulx %fp3,%fp1 | ...W*(B3+W*B5), FP3 RELEASED
  457. fmulx %fp0,%fp2 | ...V*(B2+W*B4)
  459. faddd LOGB1,%fp1 | ...B1+W*(B3+W*B5)
  460. fmulx SAVEU(%a6),%fp0 | ...FP0 IS U*V
  462. faddx %fp2,%fp1 | ...B1+W*(B3+W*B5) + V*(B2+W*B4), FP2 RELEASED
  463. fmovemx (%sp)+,%fp2-%fp2/%fp3 | ...FP2 RESTORED
  465. fmulx %fp1,%fp0 | ...U*V*( [B1+W*(B3+W*B5)] + [V*(B2+W*B4)] )
  467. fmovel %d1,%fpcr
  468. faddx SAVEU(%a6),%fp0
  469. bra t_frcinx
  470. rts
  472. LOGNEG:
  473. |--REGISTERS SAVED FPCR. LOG(-VE) IS INVALID
  474. bra t_operr
  476. .global slognp1d
  477. slognp1d:
  478. |--ENTRY POINT FOR LOG(1+Z) FOR DENORMALIZED INPUT
  479. | Simply return the denorm
  481. bra t_extdnrm
  483. .global slognp1
  484. slognp1:
  485. |--ENTRY POINT FOR LOG(1+X) FOR X FINITE, NON-ZERO, NOT NAN'S
  487. fmovex (%a0),%fp0 | ...LOAD INPUT
  488. fabsx %fp0 |test magnitude
  489. fcmpx LTHOLD,%fp0 |compare with min threshold
  490. fbgt LP1REAL |if greater, continue
  491. fmovel #0,%fpsr |clr N flag from compare
  492. fmovel %d1,%fpcr
  493. fmovex (%a0),%fp0 |return signed argument
  494. bra t_frcinx
  496. LP1REAL:
  497. fmovex (%a0),%fp0 | ...LOAD INPUT
  498. movel #0x00000000,ADJK(%a6)
  499. fmovex %fp0,%fp1 | ...FP1 IS INPUT Z
  500. fadds one,%fp0 | ...X := ROUND(1+Z)
  501. fmovex %fp0,X(%a6)
  502. movew XFRAC(%a6),XDCARE(%a6)
  503. movel X(%a6),%d0
  504. cmpil #0,%d0
  505. ble LP1NEG0 | ...LOG OF ZERO OR -VE
  506. cmp2l BOUNDS2,%d0
  507. bcs LOGMAIN | ...BOUNDS2 IS [1/2,3/2]
  508. |--IF 1+Z > 3/2 OR 1+Z < 1/2, THEN X, WHICH IS ROUNDING 1+Z,
  509. |--CONTAINS AT LEAST 63 BITS OF INFORMATION OF Z. IN THAT CASE,
  510. |--SIMPLY INVOKE LOG(X) FOR LOG(1+Z).
  512. LP1NEAR1:
  513. |--NEXT SEE IF EXP(-1/16) < X < EXP(1/16)
  514. cmp2l BOUNDS1,%d0
  515. bcss LP1CARE
  517. LP1ONE16:
  518. |--EXP(-1/16) < X < EXP(1/16). LOG(1+Z) = LOG(1+U/2) - LOG(1-U/2)
  519. |--WHERE U = 2Z/(2+Z) = 2Z/(1+X).
  520. faddx %fp1,%fp1 | ...FP1 IS 2Z
  521. fadds one,%fp0 | ...FP0 IS 1+X
  522. |--U = FP1/FP0
  523. bra LP1CONT2
  525. LP1CARE:
  526. |--HERE WE USE THE USUAL TABLE DRIVEN APPROACH. CARE HAS TO BE
  527. |--TAKEN BECAUSE 1+Z CAN HAVE 67 BITS OF INFORMATION AND WE MUST
  528. |--PRESERVE ALL THE INFORMATION. BECAUSE 1+Z IS IN [1/2,3/2],
  529. |--THERE ARE ONLY TWO CASES.
  530. |--CASE 1: 1+Z < 1, THEN K = -1 AND Y-F = (2-F) + 2Z
  531. |--CASE 2: 1+Z > 1, THEN K = 0  AND Y-F = (1-F) + Z
  532. |--ON RETURNING TO LP1CONT1, WE MUST HAVE K IN FP1, ADDRESS OF
  533. |--(1/F) IN A0, Y-F IN FP0, AND FP2 SAVED.
  535. movel XFRAC(%a6),FFRAC(%a6)
  536. andil #0xFE000000,FFRAC(%a6)
  537. oril #0x01000000,FFRAC(%a6) | ...F OBTAINED
  538. cmpil #0x3FFF8000,%d0 | ...SEE IF 1+Z > 1
  539. bges KISZERO
  541. KISNEG1:
  542. fmoves TWO,%fp0
  543. movel #0x3fff0000,F(%a6)
  544. clrl F+8(%a6)
  545. fsubx F(%a6),%fp0 | ...2-F
  546. movel FFRAC(%a6),%d0
  547. andil #0x7E000000,%d0
  548. asrl #8,%d0
  549. asrl #8,%d0
  550. asrl #4,%d0 | ...D0 CONTAINS DISPLACEMENT FOR 1/F
  551. faddx %fp1,%fp1 | ...GET 2Z
  552. fmovemx %fp2-%fp2/%fp3,-(%sp) | ...SAVE FP2 
  553. faddx %fp1,%fp0 | ...FP0 IS Y-F = (2-F)+2Z
  554. lea LOGTBL,%a0 | ...A0 IS ADDRESS OF 1/F
  555. addal %d0,%a0
  556. fmoves negone,%fp1 | ...FP1 IS K = -1
  557. bra LP1CONT1
  559. KISZERO:
  560. fmoves one,%fp0
  561. movel #0x3fff0000,F(%a6)
  562. clrl F+8(%a6)
  563. fsubx F(%a6),%fp0 | ...1-F
  564. movel FFRAC(%a6),%d0
  565. andil #0x7E000000,%d0
  566. asrl #8,%d0
  567. asrl #8,%d0
  568. asrl #4,%d0
  569. faddx %fp1,%fp0 | ...FP0 IS Y-F
  570. fmovemx %fp2-%fp2/%fp3,-(%sp) | ...FP2 SAVED
  571. lea LOGTBL,%a0
  572. addal %d0,%a0   | ...A0 IS ADDRESS OF 1/F
  573. fmoves zero,%fp1 | ...FP1 IS K = 0
  574. bra LP1CONT1
  576. LP1NEG0:
  577. |--FPCR SAVED. D0 IS X IN COMPACT FORM.
  578. cmpil #0,%d0
  579. blts LP1NEG
  580. LP1ZERO:
  581. fmoves negone,%fp0
  583. fmovel %d1,%fpcr
  584. bra t_dz
  586. LP1NEG:
  587. fmoves zero,%fp0
  589. fmovel %d1,%fpcr
  590. bra t_operr
  592. |end