slogn.s
上传用户:nvosite88
上传日期:2007-01-17
资源大小:4983k
文件大小:19k
源码类别:

VxWorks

开发平台:

C/C++

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