GNU Linux-libre 4.9.306-gnu1
[releases.git] / arch / m68k / fpsp040 / stwotox.S
1 |
2 |       stwotox.sa 3.1 12/10/90
3 |
4 |       stwotox  --- 2**X
5 |       stwotoxd --- 2**X for denormalized X
6 |       stentox  --- 10**X
7 |       stentoxd --- 10**X for denormalized X
8 |
9 |       Input: Double-extended number X in location pointed to
10 |               by address register a0.
11 |
12 |       Output: The function values are returned in 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 stwotox takes approximately 190 cycles and the
20 |               program stentox takes approximately 200 cycles.
21 |
22 |       Algorithm:
23 |
24 |       twotox
25 |       1. If |X| > 16480, go to ExpBig.
26 |
27 |       2. If |X| < 2**(-70), go to ExpSm.
28 |
29 |       3. Decompose X as X = N/64 + r where |r| <= 1/128. Furthermore
30 |               decompose N as
31 |                N = 64(M + M') + j,  j = 0,1,2,...,63.
32 |
33 |       4. Overwrite r := r * log2. Then
34 |               2**X = 2**(M') * 2**(M) * 2**(j/64) * exp(r).
35 |               Go to expr to compute that expression.
36 |
37 |       tentox
38 |       1. If |X| > 16480*log_10(2) (base 10 log of 2), go to ExpBig.
39 |
40 |       2. If |X| < 2**(-70), go to ExpSm.
41 |
42 |       3. Set y := X*log_2(10)*64 (base 2 log of 10). Set
43 |               N := round-to-int(y). Decompose N as
44 |                N = 64(M + M') + j,  j = 0,1,2,...,63.
45 |
46 |       4. Define r as
47 |               r := ((X - N*L1)-N*L2) * L10
48 |               where L1, L2 are the leading and trailing parts of log_10(2)/64
49 |               and L10 is the natural log of 10. Then
50 |               10**X = 2**(M') * 2**(M) * 2**(j/64) * exp(r).
51 |               Go to expr to compute that expression.
52 |
53 |       expr
54 |       1. Fetch 2**(j/64) from table as Fact1 and Fact2.
55 |
56 |       2. Overwrite Fact1 and Fact2 by
57 |               Fact1 := 2**(M) * Fact1
58 |               Fact2 := 2**(M) * Fact2
59 |               Thus Fact1 + Fact2 = 2**(M) * 2**(j/64).
60 |
61 |       3. Calculate P where 1 + P approximates exp(r):
62 |               P = r + r*r*(A1+r*(A2+...+r*A5)).
63 |
64 |       4. Let AdjFact := 2**(M'). Return
65 |               AdjFact * ( Fact1 + ((Fact1*P) + Fact2) ).
66 |               Exit.
67 |
68 |       ExpBig
69 |       1. Generate overflow by Huge * Huge if X > 0; otherwise, generate
70 |               underflow by Tiny * Tiny.
71 |
72 |       ExpSm
73 |       1. Return 1 + X.
74 |
75
76 |               Copyright (C) Motorola, Inc. 1990
77 |                       All Rights Reserved
78 |
79 |       For details on the license for this file, please see the
80 |       file, README, in this same directory.
81
82 |STWOTOX        idnt    2,1 | Motorola 040 Floating Point Software Package
83
84         |section        8
85
86 #include "fpsp.h"
87
88 BOUNDS1:        .long 0x3FB98000,0x400D80C0 | ... 2^(-70),16480
89 BOUNDS2:        .long 0x3FB98000,0x400B9B07 | ... 2^(-70),16480 LOG2/LOG10
90
91 L2TEN64:        .long 0x406A934F,0x0979A371 | ... 64LOG10/LOG2
92 L10TWO1:        .long 0x3F734413,0x509F8000 | ... LOG2/64LOG10
93
94 L10TWO2:        .long 0xBFCD0000,0xC0219DC1,0xDA994FD2,0x00000000
95
96 LOG10:  .long 0x40000000,0x935D8DDD,0xAAA8AC17,0x00000000
97
98 LOG2:   .long 0x3FFE0000,0xB17217F7,0xD1CF79AC,0x00000000
99
100 EXPA5:  .long 0x3F56C16D,0x6F7BD0B2
101 EXPA4:  .long 0x3F811112,0x302C712C
102 EXPA3:  .long 0x3FA55555,0x55554CC1
103 EXPA2:  .long 0x3FC55555,0x55554A54
104 EXPA1:  .long 0x3FE00000,0x00000000,0x00000000,0x00000000
105
106 HUGE:   .long 0x7FFE0000,0xFFFFFFFF,0xFFFFFFFF,0x00000000
107 TINY:   .long 0x00010000,0xFFFFFFFF,0xFFFFFFFF,0x00000000
108
109 EXPTBL:
110         .long  0x3FFF0000,0x80000000,0x00000000,0x3F738000
111         .long  0x3FFF0000,0x8164D1F3,0xBC030773,0x3FBEF7CA
112         .long  0x3FFF0000,0x82CD8698,0xAC2BA1D7,0x3FBDF8A9
113         .long  0x3FFF0000,0x843A28C3,0xACDE4046,0x3FBCD7C9
114         .long  0x3FFF0000,0x85AAC367,0xCC487B15,0xBFBDE8DA
115         .long  0x3FFF0000,0x871F6196,0x9E8D1010,0x3FBDE85C
116         .long  0x3FFF0000,0x88980E80,0x92DA8527,0x3FBEBBF1
117         .long  0x3FFF0000,0x8A14D575,0x496EFD9A,0x3FBB80CA
118         .long  0x3FFF0000,0x8B95C1E3,0xEA8BD6E7,0xBFBA8373
119         .long  0x3FFF0000,0x8D1ADF5B,0x7E5BA9E6,0xBFBE9670
120         .long  0x3FFF0000,0x8EA4398B,0x45CD53C0,0x3FBDB700
121         .long  0x3FFF0000,0x9031DC43,0x1466B1DC,0x3FBEEEB0
122         .long  0x3FFF0000,0x91C3D373,0xAB11C336,0x3FBBFD6D
123         .long  0x3FFF0000,0x935A2B2F,0x13E6E92C,0xBFBDB319
124         .long  0x3FFF0000,0x94F4EFA8,0xFEF70961,0x3FBDBA2B
125         .long  0x3FFF0000,0x96942D37,0x20185A00,0x3FBE91D5
126         .long  0x3FFF0000,0x9837F051,0x8DB8A96F,0x3FBE8D5A
127         .long  0x3FFF0000,0x99E04593,0x20B7FA65,0xBFBCDE7B
128         .long  0x3FFF0000,0x9B8D39B9,0xD54E5539,0xBFBEBAAF
129         .long  0x3FFF0000,0x9D3ED9A7,0x2CFFB751,0xBFBD86DA
130         .long  0x3FFF0000,0x9EF53260,0x91A111AE,0xBFBEBEDD
131         .long  0x3FFF0000,0xA0B0510F,0xB9714FC2,0x3FBCC96E
132         .long  0x3FFF0000,0xA2704303,0x0C496819,0xBFBEC90B
133         .long  0x3FFF0000,0xA43515AE,0x09E6809E,0x3FBBD1DB
134         .long  0x3FFF0000,0xA5FED6A9,0xB15138EA,0x3FBCE5EB
135         .long  0x3FFF0000,0xA7CD93B4,0xE965356A,0xBFBEC274
136         .long  0x3FFF0000,0xA9A15AB4,0xEA7C0EF8,0x3FBEA83C
137         .long  0x3FFF0000,0xAB7A39B5,0xA93ED337,0x3FBECB00
138         .long  0x3FFF0000,0xAD583EEA,0x42A14AC6,0x3FBE9301
139         .long  0x3FFF0000,0xAF3B78AD,0x690A4375,0xBFBD8367
140         .long  0x3FFF0000,0xB123F581,0xD2AC2590,0xBFBEF05F
141         .long  0x3FFF0000,0xB311C412,0xA9112489,0x3FBDFB3C
142         .long  0x3FFF0000,0xB504F333,0xF9DE6484,0x3FBEB2FB
143         .long  0x3FFF0000,0xB6FD91E3,0x28D17791,0x3FBAE2CB
144         .long  0x3FFF0000,0xB8FBAF47,0x62FB9EE9,0x3FBCDC3C
145         .long  0x3FFF0000,0xBAFF5AB2,0x133E45FB,0x3FBEE9AA
146         .long  0x3FFF0000,0xBD08A39F,0x580C36BF,0xBFBEAEFD
147         .long  0x3FFF0000,0xBF1799B6,0x7A731083,0xBFBCBF51
148         .long  0x3FFF0000,0xC12C4CCA,0x66709456,0x3FBEF88A
149         .long  0x3FFF0000,0xC346CCDA,0x24976407,0x3FBD83B2
150         .long  0x3FFF0000,0xC5672A11,0x5506DADD,0x3FBDF8AB
151         .long  0x3FFF0000,0xC78D74C8,0xABB9B15D,0xBFBDFB17
152         .long  0x3FFF0000,0xC9B9BD86,0x6E2F27A3,0xBFBEFE3C
153         .long  0x3FFF0000,0xCBEC14FE,0xF2727C5D,0xBFBBB6F8
154         .long  0x3FFF0000,0xCE248C15,0x1F8480E4,0xBFBCEE53
155         .long  0x3FFF0000,0xD06333DA,0xEF2B2595,0xBFBDA4AE
156         .long  0x3FFF0000,0xD2A81D91,0xF12AE45A,0x3FBC9124
157         .long  0x3FFF0000,0xD4F35AAB,0xCFEDFA1F,0x3FBEB243
158         .long  0x3FFF0000,0xD744FCCA,0xD69D6AF4,0x3FBDE69A
159         .long  0x3FFF0000,0xD99D15C2,0x78AFD7B6,0xBFB8BC61
160         .long  0x3FFF0000,0xDBFBB797,0xDAF23755,0x3FBDF610
161         .long  0x3FFF0000,0xDE60F482,0x5E0E9124,0xBFBD8BE1
162         .long  0x3FFF0000,0xE0CCDEEC,0x2A94E111,0x3FBACB12
163         .long  0x3FFF0000,0xE33F8972,0xBE8A5A51,0x3FBB9BFE
164         .long  0x3FFF0000,0xE5B906E7,0x7C8348A8,0x3FBCF2F4
165         .long  0x3FFF0000,0xE8396A50,0x3C4BDC68,0x3FBEF22F
166         .long  0x3FFF0000,0xEAC0C6E7,0xDD24392F,0xBFBDBF4A
167         .long  0x3FFF0000,0xED4F301E,0xD9942B84,0x3FBEC01A
168         .long  0x3FFF0000,0xEFE4B99B,0xDCDAF5CB,0x3FBE8CAC
169         .long  0x3FFF0000,0xF281773C,0x59FFB13A,0xBFBCBB3F
170         .long  0x3FFF0000,0xF5257D15,0x2486CC2C,0x3FBEF73A
171         .long  0x3FFF0000,0xF7D0DF73,0x0AD13BB9,0xBFB8B795
172         .long  0x3FFF0000,0xFA83B2DB,0x722A033A,0x3FBEF84B
173         .long  0x3FFF0000,0xFD3E0C0C,0xF486C175,0xBFBEF581
174
175         .set    N,L_SCR1
176
177         .set    X,FP_SCR1
178         .set    XDCARE,X+2
179         .set    XFRAC,X+4
180
181         .set    ADJFACT,FP_SCR2
182
183         .set    FACT1,FP_SCR3
184         .set    FACT1HI,FACT1+4
185         .set    FACT1LOW,FACT1+8
186
187         .set    FACT2,FP_SCR4
188         .set    FACT2HI,FACT2+4
189         .set    FACT2LOW,FACT2+8
190
191         | xref  t_unfl
192         |xref   t_ovfl
193         |xref   t_frcinx
194
195         .global stwotoxd
196 stwotoxd:
197 |--ENTRY POINT FOR 2**(X) FOR DENORMALIZED ARGUMENT
198
199         fmovel          %d1,%fpcr               | ...set user's rounding mode/precision
200         fmoves          #0x3F800000,%fp0  | ...RETURN 1 + X
201         movel           (%a0),%d0
202         orl             #0x00800001,%d0
203         fadds           %d0,%fp0
204         bra             t_frcinx
205
206         .global stwotox
207 stwotox:
208 |--ENTRY POINT FOR 2**(X), HERE X IS FINITE, NON-ZERO, AND NOT NAN'S
209         fmovemx (%a0),%fp0-%fp0 | ...LOAD INPUT, do not set cc's
210
211         movel           (%a0),%d0
212         movew           4(%a0),%d0
213         fmovex          %fp0,X(%a6)
214         andil           #0x7FFFFFFF,%d0
215
216         cmpil           #0x3FB98000,%d0         | ...|X| >= 2**(-70)?
217         bges            TWOOK1
218         bra             EXPBORS
219
220 TWOOK1:
221         cmpil           #0x400D80C0,%d0         | ...|X| > 16480?
222         bles            TWOMAIN
223         bra             EXPBORS
224
225
226 TWOMAIN:
227 |--USUAL CASE, 2^(-70) <= |X| <= 16480
228
229         fmovex          %fp0,%fp1
230         fmuls           #0x42800000,%fp1  | ...64 * X
231
232         fmovel          %fp1,N(%a6)             | ...N = ROUND-TO-INT(64 X)
233         movel           %d2,-(%sp)
234         lea             EXPTBL,%a1      | ...LOAD ADDRESS OF TABLE OF 2^(J/64)
235         fmovel          N(%a6),%fp1             | ...N --> FLOATING FMT
236         movel           N(%a6),%d0
237         movel           %d0,%d2
238         andil           #0x3F,%d0               | ...D0 IS J
239         asll            #4,%d0          | ...DISPLACEMENT FOR 2^(J/64)
240         addal           %d0,%a1         | ...ADDRESS FOR 2^(J/64)
241         asrl            #6,%d2          | ...d2 IS L, N = 64L + J
242         movel           %d2,%d0
243         asrl            #1,%d0          | ...D0 IS M
244         subl            %d0,%d2         | ...d2 IS M', N = 64(M+M') + J
245         addil           #0x3FFF,%d2
246         movew           %d2,ADJFACT(%a6)        | ...ADJFACT IS 2^(M')
247         movel           (%sp)+,%d2
248 |--SUMMARY: a1 IS ADDRESS FOR THE LEADING PORTION OF 2^(J/64),
249 |--D0 IS M WHERE N = 64(M+M') + J. NOTE THAT |M| <= 16140 BY DESIGN.
250 |--ADJFACT = 2^(M').
251 |--REGISTERS SAVED SO FAR ARE (IN ORDER) FPCR, D0, FP1, a1, AND FP2.
252
253         fmuls           #0x3C800000,%fp1  | ...(1/64)*N
254         movel           (%a1)+,FACT1(%a6)
255         movel           (%a1)+,FACT1HI(%a6)
256         movel           (%a1)+,FACT1LOW(%a6)
257         movew           (%a1)+,FACT2(%a6)
258         clrw            FACT2+2(%a6)
259
260         fsubx           %fp1,%fp0               | ...X - (1/64)*INT(64 X)
261
262         movew           (%a1)+,FACT2HI(%a6)
263         clrw            FACT2HI+2(%a6)
264         clrl            FACT2LOW(%a6)
265         addw            %d0,FACT1(%a6)
266
267         fmulx           LOG2,%fp0       | ...FP0 IS R
268         addw            %d0,FACT2(%a6)
269
270         bra             expr
271
272 EXPBORS:
273 |--FPCR, D0 SAVED
274         cmpil           #0x3FFF8000,%d0
275         bgts            EXPBIG
276
277 EXPSM:
278 |--|X| IS SMALL, RETURN 1 + X
279
280         fmovel          %d1,%FPCR               |restore users exceptions
281         fadds           #0x3F800000,%fp0  | ...RETURN 1 + X
282
283         bra             t_frcinx
284
285 EXPBIG:
286 |--|X| IS LARGE, GENERATE OVERFLOW IF X > 0; ELSE GENERATE UNDERFLOW
287 |--REGISTERS SAVE SO FAR ARE FPCR AND  D0
288         movel           X(%a6),%d0
289         cmpil           #0,%d0
290         blts            EXPNEG
291
292         bclrb           #7,(%a0)                |t_ovfl expects positive value
293         bra             t_ovfl
294
295 EXPNEG:
296         bclrb           #7,(%a0)                |t_unfl expects positive value
297         bra             t_unfl
298
299         .global stentoxd
300 stentoxd:
301 |--ENTRY POINT FOR 10**(X) FOR DENORMALIZED ARGUMENT
302
303         fmovel          %d1,%fpcr               | ...set user's rounding mode/precision
304         fmoves          #0x3F800000,%fp0  | ...RETURN 1 + X
305         movel           (%a0),%d0
306         orl             #0x00800001,%d0
307         fadds           %d0,%fp0
308         bra             t_frcinx
309
310         .global stentox
311 stentox:
312 |--ENTRY POINT FOR 10**(X), HERE X IS FINITE, NON-ZERO, AND NOT NAN'S
313         fmovemx (%a0),%fp0-%fp0 | ...LOAD INPUT, do not set cc's
314
315         movel           (%a0),%d0
316         movew           4(%a0),%d0
317         fmovex          %fp0,X(%a6)
318         andil           #0x7FFFFFFF,%d0
319
320         cmpil           #0x3FB98000,%d0         | ...|X| >= 2**(-70)?
321         bges            TENOK1
322         bra             EXPBORS
323
324 TENOK1:
325         cmpil           #0x400B9B07,%d0         | ...|X| <= 16480*log2/log10 ?
326         bles            TENMAIN
327         bra             EXPBORS
328
329 TENMAIN:
330 |--USUAL CASE, 2^(-70) <= |X| <= 16480 LOG 2 / LOG 10
331
332         fmovex          %fp0,%fp1
333         fmuld           L2TEN64,%fp1    | ...X*64*LOG10/LOG2
334
335         fmovel          %fp1,N(%a6)             | ...N=INT(X*64*LOG10/LOG2)
336         movel           %d2,-(%sp)
337         lea             EXPTBL,%a1      | ...LOAD ADDRESS OF TABLE OF 2^(J/64)
338         fmovel          N(%a6),%fp1             | ...N --> FLOATING FMT
339         movel           N(%a6),%d0
340         movel           %d0,%d2
341         andil           #0x3F,%d0               | ...D0 IS J
342         asll            #4,%d0          | ...DISPLACEMENT FOR 2^(J/64)
343         addal           %d0,%a1         | ...ADDRESS FOR 2^(J/64)
344         asrl            #6,%d2          | ...d2 IS L, N = 64L + J
345         movel           %d2,%d0
346         asrl            #1,%d0          | ...D0 IS M
347         subl            %d0,%d2         | ...d2 IS M', N = 64(M+M') + J
348         addil           #0x3FFF,%d2
349         movew           %d2,ADJFACT(%a6)        | ...ADJFACT IS 2^(M')
350         movel           (%sp)+,%d2
351
352 |--SUMMARY: a1 IS ADDRESS FOR THE LEADING PORTION OF 2^(J/64),
353 |--D0 IS M WHERE N = 64(M+M') + J. NOTE THAT |M| <= 16140 BY DESIGN.
354 |--ADJFACT = 2^(M').
355 |--REGISTERS SAVED SO FAR ARE (IN ORDER) FPCR, D0, FP1, a1, AND FP2.
356
357         fmovex          %fp1,%fp2
358
359         fmuld           L10TWO1,%fp1    | ...N*(LOG2/64LOG10)_LEAD
360         movel           (%a1)+,FACT1(%a6)
361
362         fmulx           L10TWO2,%fp2    | ...N*(LOG2/64LOG10)_TRAIL
363
364         movel           (%a1)+,FACT1HI(%a6)
365         movel           (%a1)+,FACT1LOW(%a6)
366         fsubx           %fp1,%fp0               | ...X - N L_LEAD
367         movew           (%a1)+,FACT2(%a6)
368
369         fsubx           %fp2,%fp0               | ...X - N L_TRAIL
370
371         clrw            FACT2+2(%a6)
372         movew           (%a1)+,FACT2HI(%a6)
373         clrw            FACT2HI+2(%a6)
374         clrl            FACT2LOW(%a6)
375
376         fmulx           LOG10,%fp0      | ...FP0 IS R
377
378         addw            %d0,FACT1(%a6)
379         addw            %d0,FACT2(%a6)
380
381 expr:
382 |--FPCR, FP2, FP3 ARE SAVED IN ORDER AS SHOWN.
383 |--ADJFACT CONTAINS 2**(M'), FACT1 + FACT2 = 2**(M) * 2**(J/64).
384 |--FP0 IS R. THE FOLLOWING CODE COMPUTES
385 |--     2**(M'+M) * 2**(J/64) * EXP(R)
386
387         fmovex          %fp0,%fp1
388         fmulx           %fp1,%fp1               | ...FP1 IS S = R*R
389
390         fmoved          EXPA5,%fp2      | ...FP2 IS A5
391         fmoved          EXPA4,%fp3      | ...FP3 IS A4
392
393         fmulx           %fp1,%fp2               | ...FP2 IS S*A5
394         fmulx           %fp1,%fp3               | ...FP3 IS S*A4
395
396         faddd           EXPA3,%fp2      | ...FP2 IS A3+S*A5
397         faddd           EXPA2,%fp3      | ...FP3 IS A2+S*A4
398
399         fmulx           %fp1,%fp2               | ...FP2 IS S*(A3+S*A5)
400         fmulx           %fp1,%fp3               | ...FP3 IS S*(A2+S*A4)
401
402         faddd           EXPA1,%fp2      | ...FP2 IS A1+S*(A3+S*A5)
403         fmulx           %fp0,%fp3               | ...FP3 IS R*S*(A2+S*A4)
404
405         fmulx           %fp1,%fp2               | ...FP2 IS S*(A1+S*(A3+S*A5))
406         faddx           %fp3,%fp0               | ...FP0 IS R+R*S*(A2+S*A4)
407
408         faddx           %fp2,%fp0               | ...FP0 IS EXP(R) - 1
409
410
411 |--FINAL RECONSTRUCTION PROCESS
412 |--EXP(X) = 2^M*2^(J/64) + 2^M*2^(J/64)*(EXP(R)-1)  -  (1 OR 0)
413
414         fmulx           FACT1(%a6),%fp0
415         faddx           FACT2(%a6),%fp0
416         faddx           FACT1(%a6),%fp0
417
418         fmovel          %d1,%FPCR               |restore users exceptions
419         clrw            ADJFACT+2(%a6)
420         movel           #0x80000000,ADJFACT+4(%a6)
421         clrl            ADJFACT+8(%a6)
422         fmulx           ADJFACT(%a6),%fp0       | ...FINAL ADJUSTMENT
423
424         bra             t_frcinx
425
426         |end