sys/arch/m68k/fpsp/slogn.sa

   1 *       $NetBSD: slogn.sa,v 1.3 1994/10/26 07:49:54 cgd Exp $
   2
   3 *       MOTOROLA MICROPROCESSOR & MEMORY TECHNOLOGY GROUP
   4 *       M68000 Hi-Performance Microprocessor Division
   5 *       M68040 Software Package
   6 *
   7 *       M68040 Software Package Copyright (c) 1993, 1994 Motorola Inc.
   8 *       All rights reserved.
   9 *
  10 *       THE SOFTWARE is provided on an "AS IS" basis and without warranty.
  11 *       To the maximum extent permitted by applicable law,
  12 *       MOTOROLA DISCLAIMS ALL WARRANTIES WHETHER EXPRESS OR IMPLIED,
  13 *       INCLUDING IMPLIED WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A
  14 *       PARTICULAR PURPOSE and any warranty against infringement with
  15 *       regard to the SOFTWARE (INCLUDING ANY MODIFIED VERSIONS THEREOF)
  16 *       and any accompanying written materials.
  17 *
  18 *       To the maximum extent permitted by applicable law,
  19 *       IN NO EVENT SHALL MOTOROLA BE LIABLE FOR ANY DAMAGES WHATSOEVER
  20 *       (INCLUDING WITHOUT LIMITATION, DAMAGES FOR LOSS OF BUSINESS
  21 *       PROFITS, BUSINESS INTERRUPTION, LOSS OF BUSINESS INFORMATION, OR
  22 *       OTHER PECUNIARY LOSS) ARISING OF THE USE OR INABILITY TO USE THE
  23 *       SOFTWARE.  Motorola assumes no responsibility for the maintenance
  24 *       and support of the SOFTWARE.
  25 *
  26 *       You are hereby granted a copyright license to use, modify, and
  27 *       distribute the SOFTWARE so long as this entire notice is retained
  28 *       without alteration in any modified and/or redistributed versions,
  29 *       and that such modified versions are clearly identified as such.
  30 *       No licenses are granted by implication, estoppel or otherwise
  31 *       under any patents or trademarks of Motorola, Inc.
  32
  33 *
  34 *       slogn.sa 3.1 12/10/90
  35 *
  36 *       slogn computes the natural logarithm of an
  37 *       input value. slognd does the same except the input value is a
  38 *       denormalized number. slognp1 computes log(1+X), and slognp1d
  39 *       computes log(1+X) for denormalized X.
  40 *
  41 *       Input: Double-extended value in memory location pointed to by address
  42 *               register a0.
  43 *
  44 *       Output: log(X) or log(1+X) returned in floating-point register Fp0.
  45 *
  46 *       Accuracy and Monotonicity: The returned result is within 2 ulps in
  47 *               64 significant bit, i.e. within 0.5001 ulp to 53 bits if the
  48 *               result is subsequently rounded to double precision. The
  49 *               result is provably monotonic in double precision.
  50 *
  51 *       Speed: The program slogn takes approximately 190 cycles for input
  52 *               argument X such that |X-1| >= 1/16, which is the usual
  53 *               situation. For those arguments, slognp1 takes approximately
  54 *                210 cycles. For the less common arguments, the program will
  55 *                run no worse than 10% slower.
  56 *
  57 *       Algorithm:
  58 *       LOGN:
  59 *       Step 1. If |X-1| < 1/16, approximate log(X) by an odd polynomial in
  60 *               u, where u = 2(X-1)/(X+1). Otherwise, move on to Step 2.
  61 *
  62 *       Step 2. X = 2**k * Y where 1 <= Y < 2. Define F to be the first seven
  63 *               significant bits of Y plus 2**(-7), i.e. F = 1.xxxxxx1 in base
  64 *               2 where the six "x" match those of Y. Note that |Y-F| <= 2**(-7).
  65 *
  66 *       Step 3. Define u = (Y-F)/F. Approximate log(1+u) by a polynomial in u,
  67 *               log(1+u) = poly.
  68 *
  69 *       Step 4. Reconstruct log(X) = log( 2**k * Y ) = k*log(2) + log(F) + log(1+u)
  70 *               by k*log(2) + (log(F) + poly). The values of log(F) are calculated
  71 *               beforehand and stored in the program.
  72 *
  73 *       lognp1:
  74 *       Step 1: If |X| < 1/16, approximate log(1+X) by an odd polynomial in
  75 *               u where u = 2X/(2+X). Otherwise, move on to Step 2.
  76 *
  77 *       Step 2: Let 1+X = 2**k * Y, where 1 <= Y < 2. Define F as done in Step 2
  78 *               of the algorithm for LOGN and compute log(1+X) as
  79 *               k*log(2) + log(F) + poly where poly approximates log(1+u),
  80 *               u = (Y-F)/F.
  81 *
  82 *       Implementation Notes:
  83 *       Note 1. There are 64 different possible values for F, thus 64 log(F)'s
  84 *               need to be tabulated. Moreover, the values of 1/F are also
  85 *               tabulated so that the division in (Y-F)/F can be performed by a
  86 *               multiplication.
  87 *
  88 *       Note 2. In Step 2 of lognp1, in order to preserved accuracy, the value
  89 *               Y-F has to be calculated carefully when 1/2 <= X < 3/2.
  90 *
  91 *       Note 3. To fully exploit the pipeline, polynomials are usually separated
  92 *               into two parts evaluated independently before being added up.
  93 *
  94
  95 slogn   IDNT    2,1 Motorola 040 Floating Point Software Package
  96
  97         section 8
  98
  99         include fpsp.h
 100
 101 BOUNDS1  DC.L $3FFEF07D,$3FFF8841
 102 BOUNDS2  DC.L $3FFE8000,$3FFFC000
 103
 104 LOGOF2  DC.L $3FFE0000,$B17217F7,$D1CF79AC,$00000000
 105
 106 one     DC.L $3F800000
 107 zero    DC.L $00000000
 108 infty   DC.L $7F800000
 109 negone  DC.L $BF800000
 110
 111 LOGA6   DC.L $3FC2499A,$B5E4040B
 112 LOGA5   DC.L $BFC555B5,$848CB7DB
 113
 114 LOGA4   DC.L $3FC99999,$987D8730
 115 LOGA3   DC.L $BFCFFFFF,$FF6F7E97
 116
 117 LOGA2   DC.L $3FD55555,$555555A4
 118 LOGA1   DC.L $BFE00000,$00000008
 119
 120 LOGB5   DC.L $3F175496,$ADD7DAD6
 121 LOGB4   DC.L $3F3C71C2,$FE80C7E0
 122
 123 LOGB3   DC.L $3F624924,$928BCCFF
 124 LOGB2   DC.L $3F899999,$999995EC
 125
 126 LOGB1   DC.L $3FB55555,$55555555
 127 TWO     DC.L $40000000,$00000000
 128
 129 LTHOLD  DC.L $3f990000,$80000000,$00000000,$00000000
 130
 131 LOGTBL:
 132         DC.L  $3FFE0000,$FE03F80F,$E03F80FE,$00000000
 133         DC.L  $3FF70000,$FF015358,$833C47E2,$00000000
 134         DC.L  $3FFE0000,$FA232CF2,$52138AC0,$00000000
 135         DC.L  $3FF90000,$BDC8D83E,$AD88D549,$00000000
 136         DC.L  $3FFE0000,$F6603D98,$0F6603DA,$00000000
 137         DC.L  $3FFA0000,$9CF43DCF,$F5EAFD48,$00000000
 138         DC.L  $3FFE0000,$F2B9D648,$0F2B9D65,$00000000
 139         DC.L  $3FFA0000,$DA16EB88,$CB8DF614,$00000000
 140         DC.L  $3FFE0000,$EF2EB71F,$C4345238,$00000000
 141         DC.L  $3FFB0000,$8B29B775,$1BD70743,$00000000
 142         DC.L  $3FFE0000,$EBBDB2A5,$C1619C8C,$00000000
 143         DC.L  $3FFB0000,$A8D839F8,$30C1FB49,$00000000
 144         DC.L  $3FFE0000,$E865AC7B,$7603A197,$00000000
 145         DC.L  $3FFB0000,$C61A2EB1,$8CD907AD,$00000000
 146         DC.L  $3FFE0000,$E525982A,$F70C880E,$00000000
 147         DC.L  $3FFB0000,$E2F2A47A,$DE3A18AF,$00000000
 148         DC.L  $3FFE0000,$E1FC780E,$1FC780E2,$00000000
 149         DC.L  $3FFB0000,$FF64898E,$DF55D551,$00000000
 150         DC.L  $3FFE0000,$DEE95C4C,$A037BA57,$00000000
 151         DC.L  $3FFC0000,$8DB956A9,$7B3D0148,$00000000
 152         DC.L  $3FFE0000,$DBEB61EE,$D19C5958,$00000000
 153         DC.L  $3FFC0000,$9B8FE100,$F47BA1DE,$00000000
 154         DC.L  $3FFE0000,$D901B203,$6406C80E,$00000000
 155         DC.L  $3FFC0000,$A9372F1D,$0DA1BD17,$00000000
 156         DC.L  $3FFE0000,$D62B80D6,$2B80D62C,$00000000
 157         DC.L  $3FFC0000,$B6B07F38,$CE90E46B,$00000000
 158         DC.L  $3FFE0000,$D3680D36,$80D3680D,$00000000
 159         DC.L  $3FFC0000,$C3FD0329,$06488481,$00000000
 160         DC.L  $3FFE0000,$D0B69FCB,$D2580D0B,$00000000
 161         DC.L  $3FFC0000,$D11DE0FF,$15AB18CA,$00000000
 162         DC.L  $3FFE0000,$CE168A77,$25080CE1,$00000000
 163         DC.L  $3FFC0000,$DE1433A1,$6C66B150,$00000000
 164         DC.L  $3FFE0000,$CB8727C0,$65C393E0,$00000000
 165         DC.L  $3FFC0000,$EAE10B5A,$7DDC8ADD,$00000000
 166         DC.L  $3FFE0000,$C907DA4E,$871146AD,$00000000
 167         DC.L  $3FFC0000,$F7856E5E,$E2C9B291,$00000000
 168         DC.L  $3FFE0000,$C6980C69,$80C6980C,$00000000
 169         DC.L  $3FFD0000,$82012CA5,$A68206D7,$00000000
 170         DC.L  $3FFE0000,$C4372F85,$5D824CA6,$00000000
 171         DC.L  $3FFD0000,$882C5FCD,$7256A8C5,$00000000
 172         DC.L  $3FFE0000,$C1E4BBD5,$95F6E947,$00000000
 173         DC.L  $3FFD0000,$8E44C60B,$4CCFD7DE,$00000000
 174         DC.L  $3FFE0000,$BFA02FE8,$0BFA02FF,$00000000
 175         DC.L  $3FFD0000,$944AD09E,$F4351AF6,$00000000
 176         DC.L  $3FFE0000,$BD691047,$07661AA3,$00000000
 177         DC.L  $3FFD0000,$9A3EECD4,$C3EAA6B2,$00000000
 178         DC.L  $3FFE0000,$BB3EE721,$A54D880C,$00000000
 179         DC.L  $3FFD0000,$A0218434,$353F1DE8,$00000000
 180         DC.L  $3FFE0000,$B92143FA,$36F5E02E,$00000000
 181         DC.L  $3FFD0000,$A5F2FCAB,$BBC506DA,$00000000
 182         DC.L  $3FFE0000,$B70FBB5A,$19BE3659,$00000000
 183         DC.L  $3FFD0000,$ABB3B8BA,$2AD362A5,$00000000
 184         DC.L  $3FFE0000,$B509E68A,$9B94821F,$00000000
 185         DC.L  $3FFD0000,$B1641795,$CE3CA97B,$00000000
 186         DC.L  $3FFE0000,$B30F6352,$8917C80B,$00000000
 187         DC.L  $3FFD0000,$B7047551,$5D0F1C61,$00000000
 188         DC.L  $3FFE0000,$B11FD3B8,$0B11FD3C,$00000000
 189         DC.L  $3FFD0000,$BC952AFE,$EA3D13E1,$00000000
 190         DC.L  $3FFE0000,$AF3ADDC6,$80AF3ADE,$00000000
 191         DC.L  $3FFD0000,$C2168ED0,$F458BA4A,$00000000
 192         DC.L  $3FFE0000,$AD602B58,$0AD602B6,$00000000
 193         DC.L  $3FFD0000,$C788F439,$B3163BF1,$00000000
 194         DC.L  $3FFE0000,$AB8F69E2,$8359CD11,$00000000
 195         DC.L  $3FFD0000,$CCECAC08,$BF04565D,$00000000
 196         DC.L  $3FFE0000,$A9C84A47,$A07F5638,$00000000
 197         DC.L  $3FFD0000,$D2420487,$2DD85160,$00000000
 198         DC.L  $3FFE0000,$A80A80A8,$0A80A80B,$00000000
 199         DC.L  $3FFD0000,$D7894992,$3BC3588A,$00000000
 200         DC.L  $3FFE0000,$A655C439,$2D7B73A8,$00000000
 201         DC.L  $3FFD0000,$DCC2C4B4,$9887DACC,$00000000
 202         DC.L  $3FFE0000,$A4A9CF1D,$96833751,$00000000
 203         DC.L  $3FFD0000,$E1EEBD3E,$6D6A6B9E,$00000000
 204         DC.L  $3FFE0000,$A3065E3F,$AE7CD0E0,$00000000
 205         DC.L  $3FFD0000,$E70D785C,$2F9F5BDC,$00000000
 206         DC.L  $3FFE0000,$A16B312E,$A8FC377D,$00000000
 207         DC.L  $3FFD0000,$EC1F392C,$5179F283,$00000000
 208         DC.L  $3FFE0000,$9FD809FD,$809FD80A,$00000000
 209         DC.L  $3FFD0000,$F12440D3,$E36130E6,$00000000
 210         DC.L  $3FFE0000,$9E4CAD23,$DD5F3A20,$00000000
 211         DC.L  $3FFD0000,$F61CCE92,$346600BB,$00000000
 212         DC.L  $3FFE0000,$9CC8E160,$C3FB19B9,$00000000
 213         DC.L  $3FFD0000,$FB091FD3,$8145630A,$00000000
 214         DC.L  $3FFE0000,$9B4C6F9E,$F03A3CAA,$00000000
 215         DC.L  $3FFD0000,$FFE97042,$BFA4C2AD,$00000000
 216         DC.L  $3FFE0000,$99D722DA,$BDE58F06,$00000000
 217         DC.L  $3FFE0000,$825EFCED,$49369330,$00000000
 218         DC.L  $3FFE0000,$9868C809,$868C8098,$00000000
 219         DC.L  $3FFE0000,$84C37A7A,$B9A905C9,$00000000
 220         DC.L  $3FFE0000,$97012E02,$5C04B809,$00000000
 221         DC.L  $3FFE0000,$87224C2E,$8E645FB7,$00000000
 222         DC.L  $3FFE0000,$95A02568,$095A0257,$00000000
 223         DC.L  $3FFE0000,$897B8CAC,$9F7DE298,$00000000
 224         DC.L  $3FFE0000,$94458094,$45809446,$00000000
 225         DC.L  $3FFE0000,$8BCF55DE,$C4CD05FE,$00000000
 226         DC.L  $3FFE0000,$92F11384,$0497889C,$00000000
 227         DC.L  $3FFE0000,$8E1DC0FB,$89E125E5,$00000000
 228         DC.L  $3FFE0000,$91A2B3C4,$D5E6F809,$00000000
 229         DC.L  $3FFE0000,$9066E68C,$955B6C9B,$00000000
 230         DC.L  $3FFE0000,$905A3863,$3E06C43B,$00000000
 231         DC.L  $3FFE0000,$92AADE74,$C7BE59E0,$00000000
 232         DC.L  $3FFE0000,$8F1779D9,$FDC3A219,$00000000
 233         DC.L  $3FFE0000,$94E9BFF6,$15845643,$00000000
 234         DC.L  $3FFE0000,$8DDA5202,$37694809,$00000000
 235         DC.L  $3FFE0000,$9723A1B7,$20134203,$00000000
 236         DC.L  $3FFE0000,$8CA29C04,$6514E023,$00000000
 237         DC.L  $3FFE0000,$995899C8,$90EB8990,$00000000
 238         DC.L  $3FFE0000,$8B70344A,$139BC75A,$00000000
 239         DC.L  $3FFE0000,$9B88BDAA,$3A3DAE2F,$00000000
 240         DC.L  $3FFE0000,$8A42F870,$5669DB46,$00000000
 241         DC.L  $3FFE0000,$9DB4224F,$FFE1157C,$00000000
 242         DC.L  $3FFE0000,$891AC73A,$E9819B50,$00000000
 243         DC.L  $3FFE0000,$9FDADC26,$8B7A12DA,$00000000
 244         DC.L  $3FFE0000,$87F78087,$F78087F8,$00000000
 245         DC.L  $3FFE0000,$A1FCFF17,$CE733BD4,$00000000
 246         DC.L  $3FFE0000,$86D90544,$7A34ACC6,$00000000
 247         DC.L  $3FFE0000,$A41A9E8F,$5446FB9F,$00000000
 248         DC.L  $3FFE0000,$85BF3761,$2CEE3C9B,$00000000
 249         DC.L  $3FFE0000,$A633CD7E,$6771CD8B,$00000000
 250         DC.L  $3FFE0000,$84A9F9C8,$084A9F9D,$00000000
 251         DC.L  $3FFE0000,$A8489E60,$0B435A5E,$00000000
 252         DC.L  $3FFE0000,$83993052,$3FBE3368,$00000000
 253         DC.L  $3FFE0000,$AA59233C,$CCA4BD49,$00000000
 254         DC.L  $3FFE0000,$828CBFBE,$B9A020A3,$00000000
 255         DC.L  $3FFE0000,$AC656DAE,$6BCC4985,$00000000
 256         DC.L  $3FFE0000,$81848DA8,$FAF0D277,$00000000
 257         DC.L  $3FFE0000,$AE6D8EE3,$60BB2468,$00000000
 258         DC.L  $3FFE0000,$80808080,$80808081,$00000000
 259         DC.L  $3FFE0000,$B07197A2,$3C46C654,$00000000
 260
 261 ADJK    equ     L_SCR1
 262
 263 X       equ     FP_SCR1
 264 XDCARE  equ     X+2
 265 XFRAC   equ     X+4
 266
 267 F       equ     FP_SCR2
 268 FFRAC   equ     F+4
 269
 270 KLOG2   equ     FP_SCR3
 271
 272 SAVEU   equ     FP_SCR4
 273
 274         xref    t_frcinx
 275         xref    t_extdnrm
 276         xref    t_operr
 277         xref    t_dz
 278
 279         xdef    slognd
 280 slognd:
 281 *--ENTRY POINT FOR LOG(X) FOR DENORMALIZED INPUT
 282
 283         MOVE.L          #-100,ADJK(a6)  ...INPUT = 2^(ADJK) * FP0
 284
 285 *----normalize the input value by left shifting k bits (k to be determined
 286 *----below), adjusting exponent and storing -k to  ADJK
 287 *----the value TWOTO100 is no longer needed.
 288 *----Note that this code assumes the denormalized input is NON-ZERO.
 289
 290      MoveM.L    D2-D7,-(A7)             ...save some registers
 291      Clr.L      D3                      ...D3 is exponent of smallest norm. #
 292      Move.L     4(A0),D4
 293      Move.L     8(A0),D5                ...(D4,D5) is (Hi_X,Lo_X)
 294      Clr.L      D2                      ...D2 used for holding K
 295
 296      Tst.L      D4
 297      BNE.B      HiX_not0
 298
 299 HiX_0:
 300      Move.L     D5,D4
 301      Clr.L      D5
 302      Move.L     #32,D2
 303      Clr.L      D6
 304      BFFFO      D4{0:32},D6
 305      LSL.L      D6,D4
 306      Add.L      D6,D2                   ...(D3,D4,D5) is normalized
 307
 308      Move.L     D3,X(a6)
 309      Move.L     D4,XFRAC(a6)
 310      Move.L     D5,XFRAC+4(a6)
 311      Neg.L      D2
 312      Move.L     D2,ADJK(a6)
 313      FMove.X    X(a6),FP0
 314      MoveM.L    (A7)+,D2-D7             ...restore registers
 315      LEA        X(a6),A0
 316      Bra.B      LOGBGN                  ...begin regular log(X)
 317
 318
 319 HiX_not0:
 320      Clr.L      D6
 321      BFFFO      D4{0:32},D6             ...find first 1
 322      Move.L     D6,D2                   ...get k
 323      LSL.L      D6,D4
 324      Move.L     D5,D7                   ...a copy of D5
 325      LSL.L      D6,D5
 326      Neg.L      D6
 327      AddI.L     #32,D6
 328      LSR.L      D6,D7
 329      Or.L       D7,D4                   ...(D3,D4,D5) normalized
 330
 331      Move.L     D3,X(a6)
 332      Move.L     D4,XFRAC(a6)
 333      Move.L     D5,XFRAC+4(a6)
 334      Neg.L      D2
 335      Move.L     D2,ADJK(a6)
 336      FMove.X    X(a6),FP0
 337      MoveM.L    (A7)+,D2-D7             ...restore registers
 338      LEA        X(a6),A0
 339      Bra.B      LOGBGN                  ...begin regular log(X)
 340
 341
 342         xdef    slogn
 343 slogn:
 344 *--ENTRY POINT FOR LOG(X) FOR X FINITE, NON-ZERO, NOT NAN'S
 345
 346         FMOVE.X         (A0),FP0        ...LOAD INPUT
 347         CLR.L           ADJK(a6)
 348
 349 LOGBGN:
 350 *--FPCR SAVED AND CLEARED, INPUT IS 2^(ADJK)*FP0, FP0 CONTAINS
 351 *--A FINITE, NON-ZERO, NORMALIZED NUMBER.
 352
 353         move.l  (a0),d0
 354         move.w  4(a0),d0
 355
 356         move.l  (a0),X(a6)
 357         move.l  4(a0),X+4(a6)
 358         move.l  8(a0),X+8(a6)
 359
 360         TST.L   D0              ...CHECK IF X IS NEGATIVE
 361         BLT.W   LOGNEG          ...LOG OF NEGATIVE ARGUMENT IS INVALID
 362         CMP2.L  BOUNDS1,D0      ...X IS POSITIVE, CHECK IF X IS NEAR 1
 363         BCC.W   LOGNEAR1        ...BOUNDS IS ROUGHLY [15/16, 17/16]
 364
 365 LOGMAIN:
 366 *--THIS SHOULD BE THE USUAL CASE, X NOT VERY CLOSE TO 1
 367
 368 *--X = 2^(K) * Y, 1 <= Y < 2. THUS, Y = 1.XXXXXXXX....XX IN BINARY.
 369 *--WE DEFINE F = 1.XXXXXX1, I.E. FIRST 7 BITS OF Y AND ATTACH A 1.
 370 *--THE IDEA IS THAT LOG(X) = K*LOG2 + LOG(Y)
 371 *--                      = K*LOG2 + LOG(F) + LOG(1 + (Y-F)/F).
 372 *--NOTE THAT U = (Y-F)/F IS VERY SMALL AND THUS APPROXIMATING
 373 *--LOG(1+U) CAN BE VERY EFFICIENT.
 374 *--ALSO NOTE THAT THE VALUE 1/F IS STORED IN A TABLE SO THAT NO
 375 *--DIVISION IS NEEDED TO CALCULATE (Y-F)/F.
 376
 377 *--GET K, Y, F, AND ADDRESS OF 1/F.
 378         ASR.L   #8,D0
 379         ASR.L   #8,D0           ...SHIFTED 16 BITS, BIASED EXPO. OF X
 380         SUBI.L  #$3FFF,D0       ...THIS IS K
 381         ADD.L   ADJK(a6),D0     ...ADJUST K, ORIGINAL INPUT MAY BE  DENORM.
 382         LEA     LOGTBL,A0       ...BASE ADDRESS OF 1/F AND LOG(F)
 383         FMOVE.L D0,FP1          ...CONVERT K TO FLOATING-POINT FORMAT
 384
 385 *--WHILE THE CONVERSION IS GOING ON, WE GET F AND ADDRESS OF 1/F
 386         MOVE.L  #$3FFF0000,X(a6)        ...X IS NOW Y, I.E. 2^(-K)*X
 387         MOVE.L  XFRAC(a6),FFRAC(a6)
 388         ANDI.L  #$FE000000,FFRAC(a6) ...FIRST 7 BITS OF Y
 389         ORI.L   #$01000000,FFRAC(a6) ...GET F: ATTACH A 1 AT THE EIGHTH BIT
 390         MOVE.L  FFRAC(a6),D0    ...READY TO GET ADDRESS OF 1/F
 391         ANDI.L  #$7E000000,D0
 392         ASR.L   #8,D0
 393         ASR.L   #8,D0
 394         ASR.L   #4,D0           ...SHIFTED 20, D0 IS THE DISPLACEMENT
 395         ADDA.L  D0,A0           ...A0 IS THE ADDRESS FOR 1/F
 396
 397         FMOVE.X X(a6),FP0
 398         move.l  #$3fff0000,F(a6)
 399         clr.l   F+8(a6)
 400         FSUB.X  F(a6),FP0               ...Y-F
 401         FMOVEm.X FP2/fp3,-(sp)  ...SAVE FP2 WHILE FP0 IS NOT READY
 402 *--SUMMARY: FP0 IS Y-F, A0 IS ADDRESS OF 1/F, FP1 IS K
 403 *--REGISTERS SAVED: FPCR, FP1, FP2
 404
 405 LP1CONT1:
 406 *--AN RE-ENTRY POINT FOR LOGNP1
 407         FMUL.X  (A0),FP0        ...FP0 IS U = (Y-F)/F
 408         FMUL.X  LOGOF2,FP1      ...GET K*LOG2 WHILE FP0 IS NOT READY
 409         FMOVE.X FP0,FP2
 410         FMUL.X  FP2,FP2         ...FP2 IS V=U*U
 411         FMOVE.X FP1,KLOG2(a6)   ...PUT K*LOG2 IN MEMEORY, FREE FP1
 412
 413 *--LOG(1+U) IS APPROXIMATED BY
 414 *--U + V*(A1+U*(A2+U*(A3+U*(A4+U*(A5+U*A6))))) WHICH IS
 415 *--[U + V*(A1+V*(A3+V*A5))]  +  [U*V*(A2+V*(A4+V*A6))]
 416
 417         FMOVE.X FP2,FP3
 418         FMOVE.X FP2,FP1
 419
 420         FMUL.D  LOGA6,FP1       ...V*A6
 421         FMUL.D  LOGA5,FP2       ...V*A5
 422
 423         FADD.D  LOGA4,FP1       ...A4+V*A6
 424         FADD.D  LOGA3,FP2       ...A3+V*A5
 425
 426         FMUL.X  FP3,FP1         ...V*(A4+V*A6)
 427         FMUL.X  FP3,FP2         ...V*(A3+V*A5)
 428
 429         FADD.D  LOGA2,FP1       ...A2+V*(A4+V*A6)
 430         FADD.D  LOGA1,FP2       ...A1+V*(A3+V*A5)
 431
 432         FMUL.X  FP3,FP1         ...V*(A2+V*(A4+V*A6))
 433         ADDA.L  #16,A0          ...ADDRESS OF LOG(F)
 434         FMUL.X  FP3,FP2         ...V*(A1+V*(A3+V*A5)), FP3 RELEASED
 435
 436         FMUL.X  FP0,FP1         ...U*V*(A2+V*(A4+V*A6))
 437         FADD.X  FP2,FP0         ...U+V*(A1+V*(A3+V*A5)), FP2 RELEASED
 438
 439         FADD.X  (A0),FP1        ...LOG(F)+U*V*(A2+V*(A4+V*A6))
 440         FMOVEm.X  (sp)+,FP2/fp3 ...RESTORE FP2
 441         FADD.X  FP1,FP0         ...FP0 IS LOG(F) + LOG(1+U)
 442
 443         fmove.l d1,fpcr
 444         FADD.X  KLOG2(a6),FP0   ...FINAL ADD
 445         bra     t_frcinx
 446
 447
 448 LOGNEAR1:
 449 *--REGISTERS SAVED: FPCR, FP1. FP0 CONTAINS THE INPUT.
 450         FMOVE.X FP0,FP1
 451         FSUB.S  one,FP1         ...FP1 IS X-1
 452         FADD.S  one,FP0         ...FP0 IS X+1
 453         FADD.X  FP1,FP1         ...FP1 IS 2(X-1)
 454 *--LOG(X) = LOG(1+U/2)-LOG(1-U/2) WHICH IS AN ODD POLYNOMIAL
 455 *--IN U, U = 2(X-1)/(X+1) = FP1/FP0
 456
 457 LP1CONT2:
 458 *--THIS IS AN RE-ENTRY POINT FOR LOGNP1
 459         FDIV.X  FP0,FP1         ...FP1 IS U
 460         FMOVEm.X FP2/fp3,-(sp)   ...SAVE FP2
 461 *--REGISTERS SAVED ARE NOW FPCR,FP1,FP2,FP3
 462 *--LET V=U*U, W=V*V, CALCULATE
 463 *--U + U*V*(B1 + V*(B2 + V*(B3 + V*(B4 + V*B5)))) BY
 464 *--U + U*V*(  [B1 + W*(B3 + W*B5)]  +  [V*(B2 + W*B4)]  )
 465         FMOVE.X FP1,FP0
 466         FMUL.X  FP0,FP0 ...FP0 IS V
 467         FMOVE.X FP1,SAVEU(a6) ...STORE U IN MEMORY, FREE FP1
 468         FMOVE.X FP0,FP1
 469         FMUL.X  FP1,FP1 ...FP1 IS W
 470
 471         FMOVE.D LOGB5,FP3
 472         FMOVE.D LOGB4,FP2
 473
 474         FMUL.X  FP1,FP3 ...W*B5
 475         FMUL.X  FP1,FP2 ...W*B4
 476
 477         FADD.D  LOGB3,FP3 ...B3+W*B5
 478         FADD.D  LOGB2,FP2 ...B2+W*B4
 479
 480         FMUL.X  FP3,FP1 ...W*(B3+W*B5), FP3 RELEASED
 481
 482         FMUL.X  FP0,FP2 ...V*(B2+W*B4)
 483
 484         FADD.D  LOGB1,FP1 ...B1+W*(B3+W*B5)
 485         FMUL.X  SAVEU(a6),FP0 ...FP0 IS U*V
 486
 487         FADD.X  FP2,FP1 ...B1+W*(B3+W*B5) + V*(B2+W*B4), FP2 RELEASED
 488         FMOVEm.X (sp)+,FP2/fp3 ...FP2 RESTORED
 489
 490         FMUL.X  FP1,FP0 ...U*V*( [B1+W*(B3+W*B5)] + [V*(B2+W*B4)] )
 491
 492         fmove.l d1,fpcr
 493         FADD.X  SAVEU(a6),FP0
 494         bra     t_frcinx
 495         rts
 496
 497 LOGNEG:
 498 *--REGISTERS SAVED FPCR. LOG(-VE) IS INVALID
 499         bra     t_operr
 500
 501         xdef    slognp1d
 502 slognp1d:
 503 *--ENTRY POINT FOR LOG(1+Z) FOR DENORMALIZED INPUT
 504 * Simply return the denorm
 505
 506         bra     t_extdnrm
 507
 508         xdef    slognp1
 509 slognp1:
 510 *--ENTRY POINT FOR LOG(1+X) FOR X FINITE, NON-ZERO, NOT NAN'S
 511
 512         FMOVE.X (A0),FP0        ...LOAD INPUT
 513         fabs.x  fp0             ;test magnitude
 514         fcmp.x  LTHOLD,fp0      ;compare with min threshold
 515         fbgt.w  LP1REAL         ;if greater, continue
 516         fmove.l #0,fpsr         ;clr N flag from compare
 517         fmove.l d1,fpcr
 518         fmove.x (a0),fp0        ;return signed argument
 519         bra     t_frcinx
 520
 521 LP1REAL:
 522         FMOVE.X (A0),FP0        ...LOAD INPUT
 523         CLR.L   ADJK(a6)
 524         FMOVE.X FP0,FP1 ...FP1 IS INPUT Z
 525         FADD.S  one,FP0 ...X := ROUND(1+Z)
 526         FMOVE.X FP0,X(a6)
 527         MOVE.W  XFRAC(a6),XDCARE(a6)
 528         MOVE.L  X(a6),D0
 529         TST.L   D0
 530         BLE.W   LP1NEG0 ...LOG OF ZERO OR -VE
 531         CMP2.L  BOUNDS2,D0
 532         BCS.W   LOGMAIN ...BOUNDS2 IS [1/2,3/2]
 533 *--IF 1+Z > 3/2 OR 1+Z < 1/2, THEN X, WHICH IS ROUNDING 1+Z,
 534 *--CONTAINS AT LEAST 63 BITS OF INFORMATION OF Z. IN THAT CASE,
 535 *--SIMPLY INVOKE LOG(X) FOR LOG(1+Z).
 536
 537 LP1NEAR1:
 538 *--NEXT SEE IF EXP(-1/16) < X < EXP(1/16)
 539         CMP2.L  BOUNDS1,D0
 540         BCS.B   LP1CARE
 541
 542 LP1ONE16:
 543 *--EXP(-1/16) < X < EXP(1/16). LOG(1+Z) = LOG(1+U/2) - LOG(1-U/2)
 544 *--WHERE U = 2Z/(2+Z) = 2Z/(1+X).
 545         FADD.X  FP1,FP1 ...FP1 IS 2Z
 546         FADD.S  one,FP0 ...FP0 IS 1+X
 547 *--U = FP1/FP0
 548         BRA.W   LP1CONT2
 549
 550 LP1CARE:
 551 *--HERE WE USE THE USUAL TABLE DRIVEN APPROACH. CARE HAS TO BE
 552 *--TAKEN BECAUSE 1+Z CAN HAVE 67 BITS OF INFORMATION AND WE MUST
 553 *--PRESERVE ALL THE INFORMATION. BECAUSE 1+Z IS IN [1/2,3/2],
 554 *--THERE ARE ONLY TWO CASES.
 555 *--CASE 1: 1+Z < 1, THEN K = -1 AND Y-F = (2-F) + 2Z
 556 *--CASE 2: 1+Z > 1, THEN K = 0  AND Y-F = (1-F) + Z
 557 *--ON RETURNING TO LP1CONT1, WE MUST HAVE K IN FP1, ADDRESS OF
 558 *--(1/F) IN A0, Y-F IN FP0, AND FP2 SAVED.
 559
 560         MOVE.L  XFRAC(a6),FFRAC(a6)
 561         ANDI.L  #$FE000000,FFRAC(a6)
 562         ORI.L   #$01000000,FFRAC(a6)    ...F OBTAINED
 563         CMPI.L  #$3FFF8000,D0   ...SEE IF 1+Z > 1
 564         BGE.B   KISZERO
 565
 566 KISNEG1:
 567         FMOVE.S TWO,FP0
 568         move.l  #$3fff0000,F(a6)
 569         clr.l   F+8(a6)
 570         FSUB.X  F(a6),FP0       ...2-F
 571         MOVE.L  FFRAC(a6),D0
 572         ANDI.L  #$7E000000,D0
 573         ASR.L   #8,D0
 574         ASR.L   #8,D0
 575         ASR.L   #4,D0           ...D0 CONTAINS DISPLACEMENT FOR 1/F
 576         FADD.X  FP1,FP1         ...GET 2Z
 577         FMOVEm.X FP2/fp3,-(sp)  ...SAVE FP2
 578         FADD.X  FP1,FP0         ...FP0 IS Y-F = (2-F)+2Z
 579         LEA     LOGTBL,A0       ...A0 IS ADDRESS OF 1/F
 580         ADDA.L  D0,A0
 581         FMOVE.S negone,FP1      ...FP1 IS K = -1
 582         BRA.W   LP1CONT1
 583
 584 KISZERO:
 585         FMOVE.S one,FP0
 586         move.l  #$3fff0000,F(a6)
 587         clr.l   F+8(a6)
 588         FSUB.X  F(a6),FP0               ...1-F
 589         MOVE.L  FFRAC(a6),D0
 590         ANDI.L  #$7E000000,D0
 591         ASR.L   #8,D0
 592         ASR.L   #8,D0
 593         ASR.L   #4,D0
 594         FADD.X  FP1,FP0         ...FP0 IS Y-F
 595         FMOVEm.X FP2/fp3,-(sp)  ...FP2 SAVED
 596         LEA     LOGTBL,A0
 597         ADDA.L  D0,A0           ...A0 IS ADDRESS OF 1/F
 598         FMOVE.S zero,FP1        ...FP1 IS K = 0
 599         BRA.W   LP1CONT1
 600
 601 LP1NEG0:
 602 *--FPCR SAVED. D0 IS X IN COMPACT FORM.
 603         TST.L   D0
 604         BLT.B   LP1NEG
 605 LP1ZERO:
 606         FMOVE.S negone,FP0
 607
 608         fmove.l d1,fpcr
 609         bra t_dz
 610
 611 LP1NEG:
 612         FMOVE.S zero,FP0
 613
 614         fmove.l d1,fpcr
 615         bra     t_operr
 616
 617         end