1|
2| ssin.sa 3.3 7/29/91
3|
4| The entry point sSIN computes the sine of an input argument
5| sCOS computes the cosine, and sSINCOS computes both. The
6| corresponding entry points with a "d" computes the same
7| corresponding function values for denormalized inputs.
8|
9| Input: Double-extended number X in location pointed to
10| by address register a0.
11|
12| Output: The function value sin(X) or cos(X) returned in Fp0 if SIN or
13| COS is requested. Otherwise, for SINCOS, sin(X) is returned
14| in Fp0, and cos(X) is returned in Fp1.
15|
16| Modifies: Fp0 for SIN or COS; both Fp0 and Fp1 for SINCOS.
17|
18| Accuracy and Monotonicity: The returned result is within 1 ulp in
19| 64 significant bit, i.e. within 0.5001 ulp to 53 bits if the
20| result is subsequently rounded to double precision. The
21| result is provably monotonic in double precision.
22|
23| Speed: The programs sSIN and sCOS take approximately 150 cycles for
24| input argument X such that |X| < 15Pi, which is the usual
25| situation. The speed for sSINCOS is approximately 190 cycles.
26|
27| Algorithm:
28|
29| SIN and COS:
30| 1. If SIN is invoked, set AdjN := 0; otherwise, set AdjN := 1.
31|
32| 2. If |X| >= 15Pi or |X| < 2**(-40), go to 7.
33|
34| 3. Decompose X as X = N(Pi/2) + r where |r| <= Pi/4. Let
35| k = N mod 4, so in particular, k = 0,1,2,or 3. Overwrite
36| k by k := k + AdjN.
37|
38| 4. If k is even, go to 6.
39|
40| 5. (k is odd) Set j := (k-1)/2, sgn := (-1)**j. Return sgn*cos(r)
41| where cos(r) is approximated by an even polynomial in r,
42| 1 + r*r*(B1+s*(B2+ ... + s*B8)), s = r*r.
43| Exit.
44|
45| 6. (k is even) Set j := k/2, sgn := (-1)**j. Return sgn*sin(r)
46| where sin(r) is approximated by an odd polynomial in r
47| r + r*s*(A1+s*(A2+ ... + s*A7)), s = r*r.
48| Exit.
49|
50| 7. If |X| > 1, go to 9.
51|
52| 8. (|X|<2**(-40)) If SIN is invoked, return X; otherwise return 1.
53|
54| 9. Overwrite X by X := X rem 2Pi. Now that |X| <= Pi, go back to 3.
55|
56| SINCOS:
57| 1. If |X| >= 15Pi or |X| < 2**(-40), go to 6.
58|
59| 2. Decompose X as X = N(Pi/2) + r where |r| <= Pi/4. Let
60| k = N mod 4, so in particular, k = 0,1,2,or 3.
61|
62| 3. If k is even, go to 5.
63|
64| 4. (k is odd) Set j1 := (k-1)/2, j2 := j1 (EOR) (k mod 2), i.e.
65| j1 exclusive or with the l.s.b. of k.
66| sgn1 := (-1)**j1, sgn2 := (-1)**j2.
67| SIN(X) = sgn1 * cos(r) and COS(X) = sgn2*sin(r) where
68| sin(r) and cos(r) are computed as odd and even polynomials
69| in r, respectively. Exit
70|
71| 5. (k is even) Set j1 := k/2, sgn1 := (-1)**j1.
72| SIN(X) = sgn1 * sin(r) and COS(X) = sgn1*cos(r) where
73| sin(r) and cos(r) are computed as odd and even polynomials
74| in r, respectively. Exit
75|
76| 6. If |X| > 1, go to 8.
77|
78| 7. (|X|<2**(-40)) SIN(X) = X and COS(X) = 1. Exit.
79|
80| 8. Overwrite X by X := X rem 2Pi. Now that |X| <= Pi, go back to 2.
81|
82
83| Copyright (C) Motorola, Inc. 1990
84| All Rights Reserved
85|
86| For details on the license for this file, please see the
87| file, README, in this same directory.
88
89|SSIN idnt 2,1 | Motorola 040 Floating Point Software Package
90
91 |section 8
92
93#include "fpsp.h"
94
95BOUNDS1: .long 0x3FD78000,0x4004BC7E
96TWOBYPI: .long 0x3FE45F30,0x6DC9C883
97
98SINA7: .long 0xBD6AAA77,0xCCC994F5
99SINA6: .long 0x3DE61209,0x7AAE8DA1
100
101SINA5: .long 0xBE5AE645,0x2A118AE4
102SINA4: .long 0x3EC71DE3,0xA5341531
103
104SINA3: .long 0xBF2A01A0,0x1A018B59,0x00000000,0x00000000
105
106SINA2: .long 0x3FF80000,0x88888888,0x888859AF,0x00000000
107
108SINA1: .long 0xBFFC0000,0xAAAAAAAA,0xAAAAAA99,0x00000000
109
110COSB8: .long 0x3D2AC4D0,0xD6011EE3
111COSB7: .long 0xBDA9396F,0x9F45AC19
112
113COSB6: .long 0x3E21EED9,0x0612C972
114COSB5: .long 0xBE927E4F,0xB79D9FCF
115
116COSB4: .long 0x3EFA01A0,0x1A01D423,0x00000000,0x00000000
117
118COSB3: .long 0xBFF50000,0xB60B60B6,0x0B61D438,0x00000000
119
120COSB2: .long 0x3FFA0000,0xAAAAAAAA,0xAAAAAB5E
121COSB1: .long 0xBF000000
122
123INVTWOPI: .long 0x3FFC0000,0xA2F9836E,0x4E44152A
124
125TWOPI1: .long 0x40010000,0xC90FDAA2,0x00000000,0x00000000
126TWOPI2: .long 0x3FDF0000,0x85A308D4,0x00000000,0x00000000
127
128 |xref PITBL
129
130 .set INARG,FP_SCR4
131
132 .set X,FP_SCR5
133 .set XDCARE,X+2
134 .set XFRAC,X+4
135
136 .set RPRIME,FP_SCR1
137 .set SPRIME,FP_SCR2
138
139 .set POSNEG1,L_SCR1
140 .set TWOTO63,L_SCR1
141
142 .set ENDFLAG,L_SCR2
143 .set N,L_SCR2
144
145 .set ADJN,L_SCR3
146
147 | xref t_frcinx
148 |xref t_extdnrm
149 |xref sto_cos
150
151 .global ssind
152ssind:
153|--SIN(X) = X FOR DENORMALIZED X
154 bra t_extdnrm
155
156 .global scosd
157scosd:
158|--COS(X) = 1 FOR DENORMALIZED X
159
160 fmoves
161|
162| 9D25B Fix: Sometimes the previous fmove.s sets fpsr bits
163|
164 fmovel
165|
166 bra t_frcinx
167
168 .global ssin
169ssin:
170|--SET ADJN TO 0
171 movel
172 bras SINBGN
173
174 .global scos
175scos:
176|--SET ADJN TO 1
177 movel
178
179SINBGN:
180|--SAVE FPCR, FP1. CHECK IF |X| IS TOO SMALL OR LARGE
181
182 fmovex (%a0),%fp0 | ...LOAD INPUT
183
184 movel (%a0),%d0
185 movew 4(%a0),%d0
186 fmovex %fp0,X(%a6)
187 andil
188
189 cmpil
190 bges SOK1
191 bra SINSM
192
193SOK1:
194 cmpil
195 blts SINMAIN
196 bra REDUCEX
197
198SINMAIN:
199|--THIS IS THE USUAL CASE, |X| <= 15 PI.
200|--THE ARGUMENT REDUCTION IS DONE BY TABLE LOOK UP.
201 fmovex %fp0,%fp1
202 fmuld TWOBYPI,%fp1 | ...X*2/PI
203
204|--HIDE THE NEXT THREE INSTRUCTIONS
205 lea PITBL+0x200,%a1 | ...TABLE OF N*PI/2, N = -32,...,32
206
207
208|--FP1 IS NOW READY
209 fmovel %fp1,N(%a6) | ...CONVERT TO INTEGER
210
211 movel N(%a6),%d0
212 asll
213 addal %d0,%a1 | ...A1 IS THE ADDRESS OF N*PIBY2
214| ...WHICH IS IN TWO PIECES Y1 & Y2
215
216 fsubx (%a1)+,%fp0 | ...X-Y1
217|--HIDE THE NEXT ONE
218 fsubs (%a1),%fp0 | ...FP0 IS R = (X-Y1)-Y2
219
220SINCONT:
221|--continuation from REDUCEX
222
223|--GET N+ADJN AND SEE IF SIN(R) OR COS(R) IS NEEDED
224 movel N(%a6),%d0
225 addl ADJN(%a6),%d0 | ...SEE IF D0 IS ODD OR EVEN
226 rorl
227 cmpil
228 blt COSPOLY
229
230SINPOLY:
231|--LET J BE THE LEAST SIG. BIT OF D0, LET SGN := (-1)**J.
232|--THEN WE RETURN SGN*SIN(R). SGN*SIN(R) IS COMPUTED BY
233|--R' + R'*S*(A1 + S(A2 + S(A3 + S(A4 + ... + SA7)))), WHERE
234|--R' = SGN*R, S=R*R. THIS CAN BE REWRITTEN AS
235|--R' + R'*S*( [A1+T(A3+T(A5+TA7))] + [S(A2+T(A4+TA6))])
236|--WHERE T=S*S.
237|--NOTE THAT A3 THROUGH A7 ARE STORED IN DOUBLE PRECISION
238|--WHILE A1 AND A2 ARE IN DOUBLE-EXTENDED FORMAT.
239 fmovex %fp0,X(%a6) | ...X IS R
240 fmulx %fp0,%fp0 | ...FP0 IS S
241|---HIDE THE NEXT TWO WHILE WAITING FOR FP0
242 fmoved SINA7,%fp3
243 fmoved SINA6,%fp2
244|--FP0 IS NOW READY
245 fmovex %fp0,%fp1
246 fmulx %fp1,%fp1 | ...FP1 IS T
247|--HIDE THE NEXT TWO WHILE WAITING FOR FP1
248
249 rorl
250 andil
251| ...LEAST SIG. BIT OF D0 IN SIGN POSITION
252 eorl %d0,X(%a6) | ...X IS NOW R'= SGN*R
253
254 fmulx %fp1,%fp3 | ...TA7
255 fmulx %fp1,%fp2 | ...TA6
256
257 faddd SINA5,%fp3 | ...A5+TA7
258 faddd SINA4,%fp2 | ...A4+TA6
259
260 fmulx %fp1,%fp3 | ...T(A5+TA7)
261 fmulx %fp1,%fp2 | ...T(A4+TA6)
262
263 faddd SINA3,%fp3 | ...A3+T(A5+TA7)
264 faddx SINA2,%fp2 | ...A2+T(A4+TA6)
265
266 fmulx %fp3,%fp1 | ...T(A3+T(A5+TA7))
267
268 fmulx %fp0,%fp2 | ...S(A2+T(A4+TA6))
269 faddx SINA1,%fp1 | ...A1+T(A3+T(A5+TA7))
270 fmulx X(%a6),%fp0 | ...R'*S
271
272 faddx %fp2,%fp1 | ...[A1+T(A3+T(A5+TA7))]+[S(A2+T(A4+TA6))]
273|--FP3 RELEASED, RESTORE NOW AND TAKE SOME ADVANTAGE OF HIDING
274|--FP2 RELEASED, RESTORE NOW AND TAKE FULL ADVANTAGE OF HIDING
275
276
277 fmulx %fp1,%fp0 | ...SIN(R')-R'
278|--FP1 RELEASED.
279
280 fmovel %d1,%FPCR |restore users exceptions
281 faddx X(%a6),%fp0 |last inst - possible exception set
282 bra t_frcinx
283
284
285COSPOLY:
286|--LET J BE THE LEAST SIG. BIT OF D0, LET SGN := (-1)**J.
287|--THEN WE RETURN SGN*COS(R). SGN*COS(R) IS COMPUTED BY
288|--SGN + S'*(B1 + S(B2 + S(B3 + S(B4 + ... + SB8)))), WHERE
289|--S=R*R AND S'=SGN*S. THIS CAN BE REWRITTEN AS
290|--SGN + S'*([B1+T(B3+T(B5+TB7))] + [S(B2+T(B4+T(B6+TB8)))])
291|--WHERE T=S*S.
292|--NOTE THAT B4 THROUGH B8 ARE STORED IN DOUBLE PRECISION
293|--WHILE B2 AND B3 ARE IN DOUBLE-EXTENDED FORMAT, B1 IS -1/2
294|--AND IS THEREFORE STORED AS SINGLE PRECISION.
295
296 fmulx %fp0,%fp0 | ...FP0 IS S
297|---HIDE THE NEXT TWO WHILE WAITING FOR FP0
298 fmoved COSB8,%fp2
299 fmoved COSB7,%fp3
300|--FP0 IS NOW READY
301 fmovex %fp0,%fp1
302 fmulx %fp1,%fp1 | ...FP1 IS T
303|--HIDE THE NEXT TWO WHILE WAITING FOR FP1
304 fmovex %fp0,X(%a6) | ...X IS S
305 rorl
306 andil
307| ...LEAST SIG. BIT OF D0 IN SIGN POSITION
308
309 fmulx %fp1,%fp2 | ...TB8
310|--HIDE THE NEXT TWO WHILE WAITING FOR THE XU
311 eorl %d0,X(%a6) | ...X IS NOW S'= SGN*S
312 andil
313
314 fmulx %fp1,%fp3 | ...TB7
315|--HIDE THE NEXT TWO WHILE WAITING FOR THE XU
316 oril
317 movel %d0,POSNEG1(%a6)
318
319 faddd COSB6,%fp2 | ...B6+TB8
320 faddd COSB5,%fp3 | ...B5+TB7
321
322 fmulx %fp1,%fp2 | ...T(B6+TB8)
323 fmulx %fp1,%fp3 | ...T(B5+TB7)
324
325 faddd COSB4,%fp2 | ...B4+T(B6+TB8)
326 faddx COSB3,%fp3 | ...B3+T(B5+TB7)
327
328 fmulx %fp1,%fp2 | ...T(B4+T(B6+TB8))
329 fmulx %fp3,%fp1 | ...T(B3+T(B5+TB7))
330
331 faddx COSB2,%fp2 | ...B2+T(B4+T(B6+TB8))
332 fadds COSB1,%fp1 | ...B1+T(B3+T(B5+TB7))
333
334 fmulx %fp2,%fp0 | ...S(B2+T(B4+T(B6+TB8)))
335|--FP3 RELEASED, RESTORE NOW AND TAKE SOME ADVANTAGE OF HIDING
336|--FP2 RELEASED.
337
338
339 faddx %fp1,%fp0
340|--FP1 RELEASED
341
342 fmulx X(%a6),%fp0
343
344 fmovel %d1,%FPCR |restore users exceptions
345 fadds POSNEG1(%a6),%fp0 |last inst - possible exception set
346 bra t_frcinx
347
348
349SINBORS:
350|--IF |X| > 15PI, WE USE THE GENERAL ARGUMENT REDUCTION.
351|--IF |X| < 2**(-40), RETURN X OR 1.
352 cmpil
353 bgts REDUCEX
354
355
356SINSM:
357 movel ADJN(%a6),%d0
358 cmpil
359 bgts COSTINY
360
361SINTINY:
362 movew
363 fmovel %d1,%FPCR |restore users exceptions
364 fmovex X(%a6),%fp0 |last inst - possible exception set
365 bra t_frcinx
366
367
368COSTINY:
369 fmoves
370
371 fmovel %d1,%FPCR |restore users exceptions
372 fsubs
373 bra t_frcinx
374
375
376REDUCEX:
377|--WHEN REDUCEX IS USED, THE CODE WILL INEVITABLY BE SLOW.
378|--THIS REDUCTION METHOD, HOWEVER, IS MUCH FASTER THAN USING
379|--THE REMAINDER INSTRUCTION WHICH IS NOW IN SOFTWARE.
380
381 fmovemx %fp2-%fp5,-(%a7) | ...save FP2 through FP5
382 movel %d2,-(%a7)
383 fmoves
384|--If compact form of abs(arg) in d0=$7ffeffff, argument is so large that
385|--there is a danger of unwanted overflow in first LOOP iteration. In this
386|--case, reduce argument by one remainder step to make subsequent reduction
387|--safe.
388 cmpil
389 bnes LOOP
390 movel
391| ;create 2**16383*PI/2
392 movel
393 clrl FP_SCR2+8(%a6)
394 ftstx %fp0 |test sign of argument
395 movel
396| ;PI/2 at FP_SCR3
397 movel
398 clrl FP_SCR3+8(%a6)
399 fblt red_neg
400 orw
401 orw
402red_neg:
403 faddx FP_SCR2(%a6),%fp0 |high part of reduction is exact
404 fmovex %fp0,%fp1 |save high result in fp1
405 faddx FP_SCR3(%a6),%fp0 |low part of reduction
406 fsubx %fp0,%fp1 |determine low component of result
407 faddx FP_SCR3(%a6),%fp1 |fp0/fp1 are reduced argument.
408
409|--ON ENTRY, FP0 IS X, ON RETURN, FP0 IS X REM PI/2, |X| <= PI/4.
410|--integer quotient will be stored in N
411|--Intermediate remainder is 66-bit long; (R,r) in (FP0,FP1)
412
413LOOP:
414 fmovex %fp0,INARG(%a6) | ...+-2**K * F, 1 <= F < 2
415 movew INARG(%a6),%d0
416 movel %d0,%a1 | ...save a copy of D0
417 andil
418 subil
419 cmpil
420 bles LASTLOOP
421CONTLOOP:
422 subil
423 movel
424 bras WORK
425LASTLOOP:
426 clrl %d0 | ...D0 IS L := 0
427 movel
428
429WORK:
430|--FIND THE REMAINDER OF (R,r) W.R.T. 2**L * (PI/2). L IS SO CHOSEN
431|--THAT INT( X * (2/PI) / 2**(L) ) < 2**29.
432
433|--CREATE 2**(-L) * (2/PI), SIGN(INARG)*2**(63),
434|--2**L * (PIby2_1), 2**L * (PIby2_2)
435
436 movel
437 subl %d0,%d2 | ...BIASED EXPO OF 2**(-L)*(2/PI)
438
439 movel
440 movel
441 movew %d2,FP_SCR1(%a6) | ...FP_SCR1 is 2**(-L)*(2/PI)
442
443 fmovex %fp0,%fp2
444 fmulx FP_SCR1(%a6),%fp2
445|--WE MUST NOW FIND INT(FP2). SINCE WE NEED THIS VALUE IN
446|--FLOATING POINT FORMAT, THE TWO FMOVE'S FMOVE.L FP <--> N
447|--WILL BE TOO INEFFICIENT. THE WAY AROUND IT IS THAT
448|--(SIGN(INARG)*2**63 + FP2) - SIGN(INARG)*2**63 WILL GIVE
449|--US THE DESIRED VALUE IN FLOATING POINT.
450
451|--HIDE SIX CYCLES OF INSTRUCTION
452 movel %a1,%d2
453 swap %d2
454 andil
455 oril
456 movel %d2,TWOTO63(%a6)
457
458 movel %d0,%d2
459 addil
460
461|--FP2 IS READY
462 fadds TWOTO63(%a6),%fp2 | ...THE FRACTIONAL PART OF FP1 IS ROUNDED
463
464|--HIDE 4 CYCLES OF INSTRUCTION; creating 2**(L)*Piby2_1 and 2**(L)*Piby2_2
465 movew %d2,FP_SCR2(%a6)
466 clrw FP_SCR2+2(%a6)
467 movel
468 clrl FP_SCR2+8(%a6) | ...FP_SCR2 is 2**(L) * Piby2_1
469
470|--FP2 IS READY
471 fsubs TWOTO63(%a6),%fp2 | ...FP2 is N
472
473 addil
474 movew %d0,FP_SCR3(%a6)
475 clrw FP_SCR3+2(%a6)
476 movel
477 clrl FP_SCR3+8(%a6) | ...FP_SCR3 is 2**(L) * Piby2_2
478
479 movel ENDFLAG(%a6),%d0
480
481|--We are now ready to perform (R+r) - N*P1 - N*P2, P1 = 2**(L) * Piby2_1 and
482|--P2 = 2**(L) * Piby2_2
483 fmovex %fp2,%fp4
484 fmulx FP_SCR2(%a6),%fp4 | ...W = N*P1
485 fmovex %fp2,%fp5
486 fmulx FP_SCR3(%a6),%fp5 | ...w = N*P2
487 fmovex %fp4,%fp3
488|--we want P+p = W+w but |p| <= half ulp of P
489|--Then, we need to compute A := R-P and a := r-p
490 faddx %fp5,%fp3 | ...FP3 is P
491 fsubx %fp3,%fp4 | ...W-P
492
493 fsubx %fp3,%fp0 | ...FP0 is A := R - P
494 faddx %fp5,%fp4 | ...FP4 is p = (W-P)+w
495
496 fmovex %fp0,%fp3 | ...FP3 A
497 fsubx %fp4,%fp1 | ...FP1 is a := r - p
498
499|--Now we need to normalize (A,a) to "new (R,r)" where R+r = A+a but
500|--|r| <= half ulp of R.
501 faddx %fp1,%fp0 | ...FP0 is R := A+a
502|--No need to calculate r if this is the last loop
503 cmpil
504 bgt RESTORE
505
506|--Need to calculate r
507 fsubx %fp0,%fp3 | ...A-R
508 faddx %fp3,%fp1 | ...FP1 is r := (A-R)+a
509 bra LOOP
510
511RESTORE:
512 fmovel %fp2,N(%a6)
513 movel (%a7)+,%d2
514 fmovemx (%a7)+,%fp2-%fp5
515
516
517 movel ADJN(%a6),%d0
518 cmpil
519
520 blt SINCONT
521 bras SCCONT
522
523 .global ssincosd
524ssincosd:
525|--SIN AND COS OF X FOR DENORMALIZED X
526
527 fmoves
528 bsr sto_cos |store cosine result
529 bra t_extdnrm
530
531 .global ssincos
532ssincos:
533|--SET ADJN TO 4
534 movel
535
536 fmovex (%a0),%fp0 | ...LOAD INPUT
537
538 movel (%a0),%d0
539 movew 4(%a0),%d0
540 fmovex %fp0,X(%a6)
541 andil
542
543 cmpil
544 bges SCOK1
545 bra SCSM
546
547SCOK1:
548 cmpil
549 blts SCMAIN
550 bra REDUCEX
551
552
553SCMAIN:
554|--THIS IS THE USUAL CASE, |X| <= 15 PI.
555|--THE ARGUMENT REDUCTION IS DONE BY TABLE LOOK UP.
556 fmovex %fp0,%fp1
557 fmuld TWOBYPI,%fp1 | ...X*2/PI
558
559|--HIDE THE NEXT THREE INSTRUCTIONS
560 lea PITBL+0x200,%a1 | ...TABLE OF N*PI/2, N = -32,...,32
561
562
563|--FP1 IS NOW READY
564 fmovel %fp1,N(%a6) | ...CONVERT TO INTEGER
565
566 movel N(%a6),%d0
567 asll
568 addal %d0,%a1 | ...ADDRESS OF N*PIBY2, IN Y1, Y2
569
570 fsubx (%a1)+,%fp0 | ...X-Y1
571 fsubs (%a1),%fp0 | ...FP0 IS R = (X-Y1)-Y2
572
573SCCONT:
574|--continuation point from REDUCEX
575
576|--HIDE THE NEXT TWO
577 movel N(%a6),%d0
578 rorl
579
580 cmpil
581 bge NEVEN
582
583NODD:
584|--REGISTERS SAVED SO FAR: D0, A0, FP2.
585
586 fmovex %fp0,RPRIME(%a6)
587 fmulx %fp0,%fp0 | ...FP0 IS S = R*R
588 fmoved SINA7,%fp1 | ...A7
589 fmoved COSB8,%fp2 | ...B8
590 fmulx %fp0,%fp1 | ...SA7
591 movel %d2,-(%a7)
592 movel %d0,%d2
593 fmulx %fp0,%fp2 | ...SB8
594 rorl
595 andil
596
597 faddd SINA6,%fp1 | ...A6+SA7
598 eorl %d0,%d2
599 andil
600 faddd COSB7,%fp2 | ...B7+SB8
601
602 fmulx %fp0,%fp1 | ...S(A6+SA7)
603 eorl %d2,RPRIME(%a6)
604 movel (%a7)+,%d2
605 fmulx %fp0,%fp2 | ...S(B7+SB8)
606 rorl
607 andil
608
609 faddd SINA5,%fp1 | ...A5+S(A6+SA7)
610 movel
611 eorl %d0,POSNEG1(%a6)
612 faddd COSB6,%fp2 | ...B6+S(B7+SB8)
613
614 fmulx %fp0,%fp1 | ...S(A5+S(A6+SA7))
615 fmulx %fp0,%fp2 | ...S(B6+S(B7+SB8))
616 fmovex %fp0,SPRIME(%a6)
617
618 faddd SINA4,%fp1 | ...A4+S(A5+S(A6+SA7))
619 eorl %d0,SPRIME(%a6)
620 faddd COSB5,%fp2 | ...B5+S(B6+S(B7+SB8))
621
622 fmulx %fp0,%fp1 | ...S(A4+...)
623 fmulx %fp0,%fp2 | ...S(B5+...)
624
625 faddd SINA3,%fp1 | ...A3+S(A4+...)
626 faddd COSB4,%fp2 | ...B4+S(B5+...)
627
628 fmulx %fp0,%fp1 | ...S(A3+...)
629 fmulx %fp0,%fp2 | ...S(B4+...)
630
631 faddx SINA2,%fp1 | ...A2+S(A3+...)
632 faddx COSB3,%fp2 | ...B3+S(B4+...)
633
634 fmulx %fp0,%fp1 | ...S(A2+...)
635 fmulx %fp0,%fp2 | ...S(B3+...)
636
637 faddx SINA1,%fp1 | ...A1+S(A2+...)
638 faddx COSB2,%fp2 | ...B2+S(B3+...)
639
640 fmulx %fp0,%fp1 | ...S(A1+...)
641 fmulx %fp2,%fp0 | ...S(B2+...)
642
643
644
645 fmulx RPRIME(%a6),%fp1 | ...R'S(A1+...)
646 fadds COSB1,%fp0 | ...B1+S(B2...)
647 fmulx SPRIME(%a6),%fp0 | ...S'(B1+S(B2+...))
648
649 movel %d1,-(%sp) |restore users mode & precision
650 andil
651 fmovel %d1,%FPCR
652 faddx RPRIME(%a6),%fp1 | ...COS(X)
653 bsr sto_cos |store cosine result
654 fmovel (%sp)+,%FPCR |restore users exceptions
655 fadds POSNEG1(%a6),%fp0 | ...SIN(X)
656
657 bra t_frcinx
658
659
660NEVEN:
661|--REGISTERS SAVED SO FAR: FP2.
662
663 fmovex %fp0,RPRIME(%a6)
664 fmulx %fp0,%fp0 | ...FP0 IS S = R*R
665 fmoved COSB8,%fp1 | ...B8
666 fmoved SINA7,%fp2 | ...A7
667 fmulx %fp0,%fp1 | ...SB8
668 fmovex %fp0,SPRIME(%a6)
669 fmulx %fp0,%fp2 | ...SA7
670 rorl
671 andil
672 faddd COSB7,%fp1 | ...B7+SB8
673 faddd SINA6,%fp2 | ...A6+SA7
674 eorl %d0,RPRIME(%a6)
675 eorl %d0,SPRIME(%a6)
676 fmulx %fp0,%fp1 | ...S(B7+SB8)
677 oril
678 movel %d0,POSNEG1(%a6)
679 fmulx %fp0,%fp2 | ...S(A6+SA7)
680
681 faddd COSB6,%fp1 | ...B6+S(B7+SB8)
682 faddd SINA5,%fp2 | ...A5+S(A6+SA7)
683
684 fmulx %fp0,%fp1 | ...S(B6+S(B7+SB8))
685 fmulx %fp0,%fp2 | ...S(A5+S(A6+SA7))
686
687 faddd COSB5,%fp1 | ...B5+S(B6+S(B7+SB8))
688 faddd SINA4,%fp2 | ...A4+S(A5+S(A6+SA7))
689
690 fmulx %fp0,%fp1 | ...S(B5+...)
691 fmulx %fp0,%fp2 | ...S(A4+...)
692
693 faddd COSB4,%fp1 | ...B4+S(B5+...)
694 faddd SINA3,%fp2 | ...A3+S(A4+...)
695
696 fmulx %fp0,%fp1 | ...S(B4+...)
697 fmulx %fp0,%fp2 | ...S(A3+...)
698
699 faddx COSB3,%fp1 | ...B3+S(B4+...)
700 faddx SINA2,%fp2 | ...A2+S(A3+...)
701
702 fmulx %fp0,%fp1 | ...S(B3+...)
703 fmulx %fp0,%fp2 | ...S(A2+...)
704
705 faddx COSB2,%fp1 | ...B2+S(B3+...)
706 faddx SINA1,%fp2 | ...A1+S(A2+...)
707
708 fmulx %fp0,%fp1 | ...S(B2+...)
709 fmulx %fp2,%fp0 | ...s(a1+...)
710
711
712
713 fadds COSB1,%fp1 | ...B1+S(B2...)
714 fmulx RPRIME(%a6),%fp0 | ...R'S(A1+...)
715 fmulx SPRIME(%a6),%fp1 | ...S'(B1+S(B2+...))
716
717 movel %d1,-(%sp) |save users mode & precision
718 andil
719 fmovel %d1,%FPCR
720 fadds POSNEG1(%a6),%fp1 | ...COS(X)
721 bsr sto_cos |store cosine result
722 fmovel (%sp)+,%FPCR |restore users exceptions
723 faddx RPRIME(%a6),%fp0 | ...SIN(X)
724
725 bra t_frcinx
726
727SCBORS:
728 cmpil
729 bgt REDUCEX
730
731
732SCSM:
733 movew
734 fmoves
735
736 movel %d1,-(%sp) |save users mode & precision
737 andil
738 fmovel %d1,%FPCR
739 fsubs
740 bsr sto_cos |store cosine result
741 fmovel (%sp)+,%FPCR |restore users exceptions
742 fmovex X(%a6),%fp0
743 bra t_frcinx
744
745 |end
746
