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 |       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