source: rtems/bsps/m68k/shared/fpsp/stwotox.S @ d7d66d7

5
Last change on this file since d7d66d7 was 3cf2bf63, checked in by Sebastian Huber <sebastian.huber@…>, on 03/26/18 at 10:17:06

bsps/m68k: Move fpsp support to bsps

This patch is a part of the BSP source reorganization.

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