2aaa0725c035d98696e962f2d662359b54bd4ccc
[linux-3.10.git] / 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