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