Pull bsp-removal into release branch
[linux-2.6] / arch / m68k / fpsp040 / stan.S
1 |
2 |       stan.sa 3.3 7/29/91
3 |
4 |       The entry point stan computes the tangent of
5 |       an input argument;
6 |       stand does the same except for denormalized input.
7 |
8 |       Input: Double-extended number X in location pointed to
9 |               by address register a0.
10 |
11 |       Output: The value tan(X) returned in floating-point register Fp0.
12 |
13 |       Accuracy and Monotonicity: The returned result is within 3 ulp in
14 |               64 significant bit, i.e. within 0.5001 ulp to 53 bits if the
15 |               result is subsequently rounded to double precision. The
16 |               result is provably monotonic in double precision.
17 |
18 |       Speed: The program sTAN takes approximately 170 cycles for
19 |               input argument X such that |X| < 15Pi, which is the usual
20 |               situation.
21 |
22 |       Algorithm:
23 |
24 |       1. If |X| >= 15Pi or |X| < 2**(-40), go to 6.
25 |
26 |       2. Decompose X as X = N(Pi/2) + r where |r| <= Pi/4. Let
27 |               k = N mod 2, so in particular, k = 0 or 1.
28 |
29 |       3. If k is odd, go to 5.
30 |
31 |       4. (k is even) Tan(X) = tan(r) and tan(r) is approximated by a
32 |               rational function U/V where
33 |               U = r + r*s*(P1 + s*(P2 + s*P3)), and
34 |               V = 1 + s*(Q1 + s*(Q2 + s*(Q3 + s*Q4))),  s = r*r.
35 |               Exit.
36 |
37 |       4. (k is odd) Tan(X) = -cot(r). Since tan(r) is approximated by a
38 |               rational function U/V where
39 |               U = r + r*s*(P1 + s*(P2 + s*P3)), and
40 |               V = 1 + s*(Q1 + s*(Q2 + s*(Q3 + s*Q4))), s = r*r,
41 |               -Cot(r) = -V/U. Exit.
42 |
43 |       6. If |X| > 1, go to 8.
44 |
45 |       7. (|X|<2**(-40)) Tan(X) = X. Exit.
46 |
47 |       8. Overwrite X by X := X rem 2Pi. Now that |X| <= Pi, go back to 2.
48 |
49
50 |               Copyright (C) Motorola, Inc. 1990
51 |                       All Rights Reserved
52 |
53 |       For details on the license for this file, please see the
54 |       file, README, in this same directory.
55
56 |STAN   idnt    2,1 | Motorola 040 Floating Point Software Package
57
58         |section        8
59
60 #include "fpsp.h"
61
62 BOUNDS1:        .long 0x3FD78000,0x4004BC7E
63 TWOBYPI:        .long 0x3FE45F30,0x6DC9C883
64
65 TANQ4:  .long 0x3EA0B759,0xF50F8688
66 TANP3:  .long 0xBEF2BAA5,0xA8924F04
67
68 TANQ3:  .long 0xBF346F59,0xB39BA65F,0x00000000,0x00000000
69
70 TANP2:  .long 0x3FF60000,0xE073D3FC,0x199C4A00,0x00000000
71
72 TANQ2:  .long 0x3FF90000,0xD23CD684,0x15D95FA1,0x00000000
73
74 TANP1:  .long 0xBFFC0000,0x8895A6C5,0xFB423BCA,0x00000000
75
76 TANQ1:  .long 0xBFFD0000,0xEEF57E0D,0xA84BC8CE,0x00000000
77
78 INVTWOPI: .long 0x3FFC0000,0xA2F9836E,0x4E44152A,0x00000000
79
80 TWOPI1: .long 0x40010000,0xC90FDAA2,0x00000000,0x00000000
81 TWOPI2: .long 0x3FDF0000,0x85A308D4,0x00000000,0x00000000
82
83 |--N*PI/2, -32 <= N <= 32, IN A LEADING TERM IN EXT. AND TRAILING
84 |--TERM IN SGL. NOTE THAT PI IS 64-BIT LONG, THUS N*PI/2 IS AT
85 |--MOST 69 BITS LONG.
86         .global PITBL
87 PITBL:
88   .long  0xC0040000,0xC90FDAA2,0x2168C235,0x21800000
89   .long  0xC0040000,0xC2C75BCD,0x105D7C23,0xA0D00000
90   .long  0xC0040000,0xBC7EDCF7,0xFF523611,0xA1E80000
91   .long  0xC0040000,0xB6365E22,0xEE46F000,0x21480000
92   .long  0xC0040000,0xAFEDDF4D,0xDD3BA9EE,0xA1200000
93   .long  0xC0040000,0xA9A56078,0xCC3063DD,0x21FC0000
94   .long  0xC0040000,0xA35CE1A3,0xBB251DCB,0x21100000
95   .long  0xC0040000,0x9D1462CE,0xAA19D7B9,0xA1580000
96   .long  0xC0040000,0x96CBE3F9,0x990E91A8,0x21E00000
97   .long  0xC0040000,0x90836524,0x88034B96,0x20B00000
98   .long  0xC0040000,0x8A3AE64F,0x76F80584,0xA1880000
99   .long  0xC0040000,0x83F2677A,0x65ECBF73,0x21C40000
100   .long  0xC0030000,0xFB53D14A,0xA9C2F2C2,0x20000000
101   .long  0xC0030000,0xEEC2D3A0,0x87AC669F,0x21380000
102   .long  0xC0030000,0xE231D5F6,0x6595DA7B,0xA1300000
103   .long  0xC0030000,0xD5A0D84C,0x437F4E58,0x9FC00000
104   .long  0xC0030000,0xC90FDAA2,0x2168C235,0x21000000
105   .long  0xC0030000,0xBC7EDCF7,0xFF523611,0xA1680000
106   .long  0xC0030000,0xAFEDDF4D,0xDD3BA9EE,0xA0A00000
107   .long  0xC0030000,0xA35CE1A3,0xBB251DCB,0x20900000
108   .long  0xC0030000,0x96CBE3F9,0x990E91A8,0x21600000
109   .long  0xC0030000,0x8A3AE64F,0x76F80584,0xA1080000
110   .long  0xC0020000,0xFB53D14A,0xA9C2F2C2,0x1F800000
111   .long  0xC0020000,0xE231D5F6,0x6595DA7B,0xA0B00000
112   .long  0xC0020000,0xC90FDAA2,0x2168C235,0x20800000
113   .long  0xC0020000,0xAFEDDF4D,0xDD3BA9EE,0xA0200000
114   .long  0xC0020000,0x96CBE3F9,0x990E91A8,0x20E00000
115   .long  0xC0010000,0xFB53D14A,0xA9C2F2C2,0x1F000000
116   .long  0xC0010000,0xC90FDAA2,0x2168C235,0x20000000
117   .long  0xC0010000,0x96CBE3F9,0x990E91A8,0x20600000
118   .long  0xC0000000,0xC90FDAA2,0x2168C235,0x1F800000
119   .long  0xBFFF0000,0xC90FDAA2,0x2168C235,0x1F000000
120   .long  0x00000000,0x00000000,0x00000000,0x00000000
121   .long  0x3FFF0000,0xC90FDAA2,0x2168C235,0x9F000000
122   .long  0x40000000,0xC90FDAA2,0x2168C235,0x9F800000
123   .long  0x40010000,0x96CBE3F9,0x990E91A8,0xA0600000
124   .long  0x40010000,0xC90FDAA2,0x2168C235,0xA0000000
125   .long  0x40010000,0xFB53D14A,0xA9C2F2C2,0x9F000000
126   .long  0x40020000,0x96CBE3F9,0x990E91A8,0xA0E00000
127   .long  0x40020000,0xAFEDDF4D,0xDD3BA9EE,0x20200000
128   .long  0x40020000,0xC90FDAA2,0x2168C235,0xA0800000
129   .long  0x40020000,0xE231D5F6,0x6595DA7B,0x20B00000
130   .long  0x40020000,0xFB53D14A,0xA9C2F2C2,0x9F800000
131   .long  0x40030000,0x8A3AE64F,0x76F80584,0x21080000
132   .long  0x40030000,0x96CBE3F9,0x990E91A8,0xA1600000
133   .long  0x40030000,0xA35CE1A3,0xBB251DCB,0xA0900000
134   .long  0x40030000,0xAFEDDF4D,0xDD3BA9EE,0x20A00000
135   .long  0x40030000,0xBC7EDCF7,0xFF523611,0x21680000
136   .long  0x40030000,0xC90FDAA2,0x2168C235,0xA1000000
137   .long  0x40030000,0xD5A0D84C,0x437F4E58,0x1FC00000
138   .long  0x40030000,0xE231D5F6,0x6595DA7B,0x21300000
139   .long  0x40030000,0xEEC2D3A0,0x87AC669F,0xA1380000
140   .long  0x40030000,0xFB53D14A,0xA9C2F2C2,0xA0000000
141   .long  0x40040000,0x83F2677A,0x65ECBF73,0xA1C40000
142   .long  0x40040000,0x8A3AE64F,0x76F80584,0x21880000
143   .long  0x40040000,0x90836524,0x88034B96,0xA0B00000
144   .long  0x40040000,0x96CBE3F9,0x990E91A8,0xA1E00000
145   .long  0x40040000,0x9D1462CE,0xAA19D7B9,0x21580000
146   .long  0x40040000,0xA35CE1A3,0xBB251DCB,0xA1100000
147   .long  0x40040000,0xA9A56078,0xCC3063DD,0xA1FC0000
148   .long  0x40040000,0xAFEDDF4D,0xDD3BA9EE,0x21200000
149   .long  0x40040000,0xB6365E22,0xEE46F000,0xA1480000
150   .long  0x40040000,0xBC7EDCF7,0xFF523611,0x21E80000
151   .long  0x40040000,0xC2C75BCD,0x105D7C23,0x20D00000
152   .long  0x40040000,0xC90FDAA2,0x2168C235,0xA1800000
153
154         .set    INARG,FP_SCR4
155
156         .set    TWOTO63,L_SCR1
157         .set    ENDFLAG,L_SCR2
158         .set    N,L_SCR3
159
160         | xref  t_frcinx
161         |xref   t_extdnrm
162
163         .global stand
164 stand:
165 |--TAN(X) = X FOR DENORMALIZED X
166
167         bra             t_extdnrm
168
169         .global stan
170 stan:
171         fmovex          (%a0),%fp0      | ...LOAD INPUT
172
173         movel           (%a0),%d0
174         movew           4(%a0),%d0
175         andil           #0x7FFFFFFF,%d0
176
177         cmpil           #0x3FD78000,%d0         | ...|X| >= 2**(-40)?
178         bges            TANOK1
179         bra             TANSM
180 TANOK1:
181         cmpil           #0x4004BC7E,%d0         | ...|X| < 15 PI?
182         blts            TANMAIN
183         bra             REDUCEX
184
185
186 TANMAIN:
187 |--THIS IS THE USUAL CASE, |X| <= 15 PI.
188 |--THE ARGUMENT REDUCTION IS DONE BY TABLE LOOK UP.
189         fmovex          %fp0,%fp1
190         fmuld           TWOBYPI,%fp1    | ...X*2/PI
191
192 |--HIDE THE NEXT TWO INSTRUCTIONS
193         leal            PITBL+0x200,%a1 | ...TABLE OF N*PI/2, N = -32,...,32
194
195 |--FP1 IS NOW READY
196         fmovel          %fp1,%d0                | ...CONVERT TO INTEGER
197
198         asll            #4,%d0
199         addal           %d0,%a1         | ...ADDRESS N*PIBY2 IN Y1, Y2
200
201         fsubx           (%a1)+,%fp0     | ...X-Y1
202 |--HIDE THE NEXT ONE
203
204         fsubs           (%a1),%fp0      | ...FP0 IS R = (X-Y1)-Y2
205
206         rorl            #5,%d0
207         andil           #0x80000000,%d0 | ...D0 WAS ODD IFF D0 < 0
208
209 TANCONT:
210
211         cmpil           #0,%d0
212         blt             NODD
213
214         fmovex          %fp0,%fp1
215         fmulx           %fp1,%fp1               | ...S = R*R
216
217         fmoved          TANQ4,%fp3
218         fmoved          TANP3,%fp2
219
220         fmulx           %fp1,%fp3               | ...SQ4
221         fmulx           %fp1,%fp2               | ...SP3
222
223         faddd           TANQ3,%fp3      | ...Q3+SQ4
224         faddx           TANP2,%fp2      | ...P2+SP3
225
226         fmulx           %fp1,%fp3               | ...S(Q3+SQ4)
227         fmulx           %fp1,%fp2               | ...S(P2+SP3)
228
229         faddx           TANQ2,%fp3      | ...Q2+S(Q3+SQ4)
230         faddx           TANP1,%fp2      | ...P1+S(P2+SP3)
231
232         fmulx           %fp1,%fp3               | ...S(Q2+S(Q3+SQ4))
233         fmulx           %fp1,%fp2               | ...S(P1+S(P2+SP3))
234
235         faddx           TANQ1,%fp3      | ...Q1+S(Q2+S(Q3+SQ4))
236         fmulx           %fp0,%fp2               | ...RS(P1+S(P2+SP3))
237
238         fmulx           %fp3,%fp1               | ...S(Q1+S(Q2+S(Q3+SQ4)))
239
240
241         faddx           %fp2,%fp0               | ...R+RS(P1+S(P2+SP3))
242
243
244         fadds           #0x3F800000,%fp1        | ...1+S(Q1+...)
245
246         fmovel          %d1,%fpcr               |restore users exceptions
247         fdivx           %fp1,%fp0               |last inst - possible exception set
248
249         bra             t_frcinx
250
251 NODD:
252         fmovex          %fp0,%fp1
253         fmulx           %fp0,%fp0               | ...S = R*R
254
255         fmoved          TANQ4,%fp3
256         fmoved          TANP3,%fp2
257
258         fmulx           %fp0,%fp3               | ...SQ4
259         fmulx           %fp0,%fp2               | ...SP3
260
261         faddd           TANQ3,%fp3      | ...Q3+SQ4
262         faddx           TANP2,%fp2      | ...P2+SP3
263
264         fmulx           %fp0,%fp3               | ...S(Q3+SQ4)
265         fmulx           %fp0,%fp2               | ...S(P2+SP3)
266
267         faddx           TANQ2,%fp3      | ...Q2+S(Q3+SQ4)
268         faddx           TANP1,%fp2      | ...P1+S(P2+SP3)
269
270         fmulx           %fp0,%fp3               | ...S(Q2+S(Q3+SQ4))
271         fmulx           %fp0,%fp2               | ...S(P1+S(P2+SP3))
272
273         faddx           TANQ1,%fp3      | ...Q1+S(Q2+S(Q3+SQ4))
274         fmulx           %fp1,%fp2               | ...RS(P1+S(P2+SP3))
275
276         fmulx           %fp3,%fp0               | ...S(Q1+S(Q2+S(Q3+SQ4)))
277
278
279         faddx           %fp2,%fp1               | ...R+RS(P1+S(P2+SP3))
280         fadds           #0x3F800000,%fp0        | ...1+S(Q1+...)
281
282
283         fmovex          %fp1,-(%sp)
284         eoril           #0x80000000,(%sp)
285
286         fmovel          %d1,%fpcr               |restore users exceptions
287         fdivx           (%sp)+,%fp0     |last inst - possible exception set
288
289         bra             t_frcinx
290
291 TANBORS:
292 |--IF |X| > 15PI, WE USE THE GENERAL ARGUMENT REDUCTION.
293 |--IF |X| < 2**(-40), RETURN X OR 1.
294         cmpil           #0x3FFF8000,%d0
295         bgts            REDUCEX
296
297 TANSM:
298
299         fmovex          %fp0,-(%sp)
300         fmovel          %d1,%fpcr                |restore users exceptions
301         fmovex          (%sp)+,%fp0     |last inst - possible exception set
302
303         bra             t_frcinx
304
305
306 REDUCEX:
307 |--WHEN REDUCEX IS USED, THE CODE WILL INEVITABLY BE SLOW.
308 |--THIS REDUCTION METHOD, HOWEVER, IS MUCH FASTER THAN USING
309 |--THE REMAINDER INSTRUCTION WHICH IS NOW IN SOFTWARE.
310
311         fmovemx %fp2-%fp5,-(%a7)        | ...save FP2 through FP5
312         movel           %d2,-(%a7)
313         fmoves         #0x00000000,%fp1
314
315 |--If compact form of abs(arg) in d0=$7ffeffff, argument is so large that
316 |--there is a danger of unwanted overflow in first LOOP iteration.  In this
317 |--case, reduce argument by one remainder step to make subsequent reduction
318 |--safe.
319         cmpil   #0x7ffeffff,%d0         |is argument dangerously large?
320         bnes    LOOP
321         movel   #0x7ffe0000,FP_SCR2(%a6)        |yes
322 |                                       ;create 2**16383*PI/2
323         movel   #0xc90fdaa2,FP_SCR2+4(%a6)
324         clrl    FP_SCR2+8(%a6)
325         ftstx   %fp0                    |test sign of argument
326         movel   #0x7fdc0000,FP_SCR3(%a6)        |create low half of 2**16383*
327 |                                       ;PI/2 at FP_SCR3
328         movel   #0x85a308d3,FP_SCR3+4(%a6)
329         clrl   FP_SCR3+8(%a6)
330         fblt    red_neg
331         orw     #0x8000,FP_SCR2(%a6)    |positive arg
332         orw     #0x8000,FP_SCR3(%a6)
333 red_neg:
334         faddx  FP_SCR2(%a6),%fp0                |high part of reduction is exact
335         fmovex  %fp0,%fp1               |save high result in fp1
336         faddx  FP_SCR3(%a6),%fp0                |low part of reduction
337         fsubx  %fp0,%fp1                        |determine low component of result
338         faddx  FP_SCR3(%a6),%fp1                |fp0/fp1 are reduced argument.
339
340 |--ON ENTRY, FP0 IS X, ON RETURN, FP0 IS X REM PI/2, |X| <= PI/4.
341 |--integer quotient will be stored in N
342 |--Intermediate remainder is 66-bit long; (R,r) in (FP0,FP1)
343
344 LOOP:
345         fmovex          %fp0,INARG(%a6) | ...+-2**K * F, 1 <= F < 2
346         movew           INARG(%a6),%d0
347         movel          %d0,%a1          | ...save a copy of D0
348         andil           #0x00007FFF,%d0
349         subil           #0x00003FFF,%d0 | ...D0 IS K
350         cmpil           #28,%d0
351         bles            LASTLOOP
352 CONTLOOP:
353         subil           #27,%d0  | ...D0 IS L := K-27
354         movel           #0,ENDFLAG(%a6)
355         bras            WORK
356 LASTLOOP:
357         clrl            %d0             | ...D0 IS L := 0
358         movel           #1,ENDFLAG(%a6)
359
360 WORK:
361 |--FIND THE REMAINDER OF (R,r) W.R.T.   2**L * (PI/2). L IS SO CHOSEN
362 |--THAT INT( X * (2/PI) / 2**(L) ) < 2**29.
363
364 |--CREATE 2**(-L) * (2/PI), SIGN(INARG)*2**(63),
365 |--2**L * (PIby2_1), 2**L * (PIby2_2)
366
367         movel           #0x00003FFE,%d2 | ...BIASED EXPO OF 2/PI
368         subl            %d0,%d2         | ...BIASED EXPO OF 2**(-L)*(2/PI)
369
370         movel           #0xA2F9836E,FP_SCR1+4(%a6)
371         movel           #0x4E44152A,FP_SCR1+8(%a6)
372         movew           %d2,FP_SCR1(%a6)        | ...FP_SCR1 is 2**(-L)*(2/PI)
373
374         fmovex          %fp0,%fp2
375         fmulx           FP_SCR1(%a6),%fp2
376 |--WE MUST NOW FIND INT(FP2). SINCE WE NEED THIS VALUE IN
377 |--FLOATING POINT FORMAT, THE TWO FMOVE'S       FMOVE.L FP <--> N
378 |--WILL BE TOO INEFFICIENT. THE WAY AROUND IT IS THAT
379 |--(SIGN(INARG)*2**63   +       FP2) - SIGN(INARG)*2**63 WILL GIVE
380 |--US THE DESIRED VALUE IN FLOATING POINT.
381
382 |--HIDE SIX CYCLES OF INSTRUCTION
383         movel           %a1,%d2
384         swap            %d2
385         andil           #0x80000000,%d2
386         oril            #0x5F000000,%d2 | ...D2 IS SIGN(INARG)*2**63 IN SGL
387         movel           %d2,TWOTO63(%a6)
388
389         movel           %d0,%d2
390         addil           #0x00003FFF,%d2 | ...BIASED EXPO OF 2**L * (PI/2)
391
392 |--FP2 IS READY
393         fadds           TWOTO63(%a6),%fp2       | ...THE FRACTIONAL PART OF FP1 IS ROUNDED
394
395 |--HIDE 4 CYCLES OF INSTRUCTION; creating 2**(L)*Piby2_1  and  2**(L)*Piby2_2
396         movew           %d2,FP_SCR2(%a6)
397         clrw           FP_SCR2+2(%a6)
398         movel           #0xC90FDAA2,FP_SCR2+4(%a6)
399         clrl            FP_SCR2+8(%a6)          | ...FP_SCR2 is  2**(L) * Piby2_1
400
401 |--FP2 IS READY
402         fsubs           TWOTO63(%a6),%fp2               | ...FP2 is N
403
404         addil           #0x00003FDD,%d0
405         movew           %d0,FP_SCR3(%a6)
406         clrw           FP_SCR3+2(%a6)
407         movel           #0x85A308D3,FP_SCR3+4(%a6)
408         clrl            FP_SCR3+8(%a6)          | ...FP_SCR3 is 2**(L) * Piby2_2
409
410         movel           ENDFLAG(%a6),%d0
411
412 |--We are now ready to perform (R+r) - N*P1 - N*P2, P1 = 2**(L) * Piby2_1 and
413 |--P2 = 2**(L) * Piby2_2
414         fmovex          %fp2,%fp4
415         fmulx           FP_SCR2(%a6),%fp4               | ...W = N*P1
416         fmovex          %fp2,%fp5
417         fmulx           FP_SCR3(%a6),%fp5               | ...w = N*P2
418         fmovex          %fp4,%fp3
419 |--we want P+p = W+w  but  |p| <= half ulp of P
420 |--Then, we need to compute  A := R-P   and  a := r-p
421         faddx           %fp5,%fp3                       | ...FP3 is P
422         fsubx           %fp3,%fp4                       | ...W-P
423
424         fsubx           %fp3,%fp0                       | ...FP0 is A := R - P
425         faddx           %fp5,%fp4                       | ...FP4 is p = (W-P)+w
426
427         fmovex          %fp0,%fp3                       | ...FP3 A
428         fsubx           %fp4,%fp1                       | ...FP1 is a := r - p
429
430 |--Now we need to normalize (A,a) to  "new (R,r)" where R+r = A+a but
431 |--|r| <= half ulp of R.
432         faddx           %fp1,%fp0                       | ...FP0 is R := A+a
433 |--No need to calculate r if this is the last loop
434         cmpil           #0,%d0
435         bgt             RESTORE
436
437 |--Need to calculate r
438         fsubx           %fp0,%fp3                       | ...A-R
439         faddx           %fp3,%fp1                       | ...FP1 is r := (A-R)+a
440         bra             LOOP
441
442 RESTORE:
443         fmovel          %fp2,N(%a6)
444         movel           (%a7)+,%d2
445         fmovemx (%a7)+,%fp2-%fp5
446
447
448         movel           N(%a6),%d0
449         rorl            #1,%d0
450
451
452         bra             TANCONT
453
454         |end