arch/m68k/fpsp040/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 |
  62
  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.
  69
  70 |slogn  idnt    2,1 | Motorola 040 Floating Point Software Package
  71
  72         |section        8
  73
  74 #include "fpsp.h"
  75
  76 BOUNDS1:  .long 0x3FFEF07D,0x3FFF8841
  77 BOUNDS2:  .long 0x3FFE8000,0x3FFFC000
  78
  79 LOGOF2: .long 0x3FFE0000,0xB17217F7,0xD1CF79AC,0x00000000
  80
  81 one:    .long 0x3F800000
  82 zero:   .long 0x00000000
  83 infty:  .long 0x7F800000
  84 negone: .long 0xBF800000
  85
  86 LOGA6:  .long 0x3FC2499A,0xB5E4040B
  87 LOGA5:  .long 0xBFC555B5,0x848CB7DB
  88
  89 LOGA4:  .long 0x3FC99999,0x987D8730
  90 LOGA3:  .long 0xBFCFFFFF,0xFF6F7E97
  91
  92 LOGA2:  .long 0x3FD55555,0x555555a4
  93 LOGA1:  .long 0xBFE00000,0x00000008
  94
  95 LOGB5:  .long 0x3F175496,0xADD7DAD6
  96 LOGB4:  .long 0x3F3C71C2,0xFE80C7E0
  97
  98 LOGB3:  .long 0x3F624924,0x928BCCFF
  99 LOGB2:  .long 0x3F899999,0x999995EC
 100
 101 LOGB1:  .long 0x3FB55555,0x55555555
 102 TWO:    .long 0x40000000,0x00000000
 103
 104 LTHOLD: .long 0x3f990000,0x80000000,0x00000000,0x00000000
 105
 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
 235
 236         .set    ADJK,L_SCR1
 237
 238         .set    X,FP_SCR1
 239         .set    XDCARE,X+2
 240         .set    XFRAC,X+4
 241
 242         .set    F,FP_SCR2
 243         .set    FFRAC,F+4
 244
 245         .set    KLOG2,FP_SCR3
 246
 247         .set    SAVEU,FP_SCR4
 248
 249         | xref  t_frcinx
 250         |xref   t_extdnrm
 251         |xref   t_operr
 252         |xref   t_dz
 253
 254         .global slognd
 255 slognd:
 256 |--ENTRY POINT FOR LOG(X) FOR DENORMALIZED INPUT
 257
 258         movel           #-100,ADJK(%a6) | ...INPUT = 2^(ADJK) * FP0
 259
 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.
 264
 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
 270
 271      tstl       %d4
 272      bnes       HiX_not0
 273
 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
 282
 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)
 292
 293
 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
 305
 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)
 315
 316
 317         .global slogn
 318 slogn:
 319 |--ENTRY POINT FOR LOG(X) FOR X FINITE, NON-ZERO, NOT NAN'S
 320
 321         fmovex          (%a0),%fp0      | ...LOAD INPUT
 322         movel           #0x00000000,ADJK(%a6)
 323
 324 LOGBGN:
 325 |--FPCR SAVED AND CLEARED, INPUT IS 2^(ADJK)*FP0, FP0 CONTAINS
 326 |--A FINITE, NON-ZERO, NORMALIZED NUMBER.
 327
 328         movel   (%a0),%d0
 329         movew   4(%a0),%d0
 330
 331         movel   (%a0),X(%a6)
 332         movel   4(%a0),X+4(%a6)
 333         movel   8(%a0),X+8(%a6)
 334
 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]
 339
 340 LOGMAIN:
 341 |--THIS SHOULD BE THE USUAL CASE, X NOT VERY CLOSE TO 1
 342
 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.
 351
 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
 359
 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
 371
 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
 379
 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
 387
 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))]
 391
 392         fmovex  %fp2,%fp3
 393         fmovex  %fp2,%fp1
 394
 395         fmuld   LOGA6,%fp1      | ...V*A6
 396         fmuld   LOGA5,%fp2      | ...V*A5
 397
 398         faddd   LOGA4,%fp1      | ...A4+V*A6
 399         faddd   LOGA3,%fp2      | ...A3+V*A5
 400
 401         fmulx   %fp3,%fp1               | ...V*(A4+V*A6)
 402         fmulx   %fp3,%fp2               | ...V*(A3+V*A5)
 403
 404         faddd   LOGA2,%fp1      | ...A2+V*(A4+V*A6)
 405         faddd   LOGA1,%fp2      | ...A1+V*(A3+V*A5)
 406
 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
 410
 411         fmulx   %fp0,%fp1               | ...U*V*(A2+V*(A4+V*A6))
 412         faddx   %fp2,%fp0               | ...U+V*(A1+V*(A3+V*A5)), FP2 RELEASED
 413
 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)
 417
 418         fmovel  %d1,%fpcr
 419         faddx   KLOG2(%a6),%fp0 | ...FINAL ADD
 420         bra     t_frcinx
 421
 422
 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
 431
 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
 445
 446         fmoved  LOGB5,%fp3
 447         fmoved  LOGB4,%fp2
 448
 449         fmulx   %fp1,%fp3       | ...W*B5
 450         fmulx   %fp1,%fp2       | ...W*B4
 451
 452         faddd   LOGB3,%fp3 | ...B3+W*B5
 453         faddd   LOGB2,%fp2 | ...B2+W*B4
 454
 455         fmulx   %fp3,%fp1       | ...W*(B3+W*B5), FP3 RELEASED
 456
 457         fmulx   %fp0,%fp2       | ...V*(B2+W*B4)
 458
 459         faddd   LOGB1,%fp1 | ...B1+W*(B3+W*B5)
 460         fmulx   SAVEU(%a6),%fp0 | ...FP0 IS U*V
 461
 462         faddx   %fp2,%fp1       | ...B1+W*(B3+W*B5) + V*(B2+W*B4), FP2 RELEASED
 463         fmovemx (%sp)+,%fp2-%fp2/%fp3 | ...FP2 RESTORED
 464
 465         fmulx   %fp1,%fp0       | ...U*V*( [B1+W*(B3+W*B5)] + [V*(B2+W*B4)] )
 466
 467         fmovel  %d1,%fpcr
 468         faddx   SAVEU(%a6),%fp0
 469         bra     t_frcinx
 470         rts
 471
 472 LOGNEG:
 473 |--REGISTERS SAVED FPCR. LOG(-VE) IS INVALID
 474         bra     t_operr
 475
 476         .global slognp1d
 477 slognp1d:
 478 |--ENTRY POINT FOR LOG(1+Z) FOR DENORMALIZED INPUT
 479 | Simply return the denorm
 480
 481         bra     t_extdnrm
 482
 483         .global slognp1
 484 slognp1:
 485 |--ENTRY POINT FOR LOG(1+X) FOR X FINITE, NON-ZERO, NOT NAN'S
 486
 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
 495
 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).
 511
 512 LP1NEAR1:
 513 |--NEXT SEE IF EXP(-1/16) < X < EXP(1/16)
 514         cmp2l   BOUNDS1,%d0
 515         bcss    LP1CARE
 516
 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
 524
 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.
 534
 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
 540
 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
 558
 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
 575
 576 LP1NEG0:
 577 |--FPCR SAVED. D0 IS X IN COMPACT FORM.
 578         cmpil   #0,%d0
 579         blts    LP1NEG
 580 LP1ZERO:
 581         fmoves  negone,%fp0
 582
 583         fmovel  %d1,%fpcr
 584         bra t_dz
 585
 586 LP1NEG:
 587         fmoves  zero,%fp0
 588
 589         fmovel  %d1,%fpcr
 590         bra     t_operr
 591
 592         |end