![[BACK]](/icons/back.gif){kind=icon}
Return to ssin.sa CVS log ![[TXT]](/icons/text.gif){kind=icon}
![[DIR]](/icons/dir.gif){kind=icon}
Up to [local] / sys / arch / m68k / fpsp|
Return to ssin.sa CVS log | Up to [local] / sys / arch / m68k / fpsp |
1.1 nbrk 1: * $OpenBSD: ssin.sa,v 1.3 2003/11/07 10:36:10 miod Exp $
2: * $NetBSD: ssin.sa,v 1.3 1994/10/26 07:50:01 cgd Exp $
3:
4: * MOTOROLA MICROPROCESSOR & MEMORY TECHNOLOGY GROUP
5: * M68000 Hi-Performance Microprocessor Division
6: * M68040 Software Package
7: *
8: * M68040 Software Package Copyright (c) 1993, 1994 Motorola Inc.
9: * All rights reserved.
10: *
11: * THE SOFTWARE is provided on an "AS IS" basis and without warranty.
12: * To the maximum extent permitted by applicable law,
13: * MOTOROLA DISCLAIMS ALL WARRANTIES WHETHER EXPRESS OR IMPLIED,
14: * INCLUDING IMPLIED WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A
15: * PARTICULAR PURPOSE and any warranty against infringement with
16: * regard to the SOFTWARE (INCLUDING ANY MODIFIED VERSIONS THEREOF)
17: * and any accompanying written materials.
18: *
19: * To the maximum extent permitted by applicable law,
20: * IN NO EVENT SHALL MOTOROLA BE LIABLE FOR ANY DAMAGES WHATSOEVER
21: * (INCLUDING WITHOUT LIMITATION, DAMAGES FOR LOSS OF BUSINESS
22: * PROFITS, BUSINESS INTERRUPTION, LOSS OF BUSINESS INFORMATION, OR
23: * OTHER PECUNIARY LOSS) ARISING OF THE USE OR INABILITY TO USE THE
24: * SOFTWARE. Motorola assumes no responsibility for the maintenance
25: * and support of the SOFTWARE.
26: *
27: * You are hereby granted a copyright license to use, modify, and
28: * distribute the SOFTWARE so long as this entire notice is retained
29: * without alteration in any modified and/or redistributed versions,
30: * and that such modified versions are clearly identified as such.
31: * No licenses are granted by implication, estoppel or otherwise
32: * under any patents or trademarks of Motorola, Inc.
33:
34: *
35: * ssin.sa 3.3 7/29/91
36: *
37: * The entry point sSIN computes the sine of an input argument
38: * sCOS computes the cosine, and sSINCOS computes both. The
39: * corresponding entry points with a "d" computes the same
40: * corresponding function values for denormalized inputs.
41: *
42: * Input: Double-extended number X in location pointed to
43: * by address register a0.
44: *
45: * Output: The funtion value sin(X) or cos(X) returned in Fp0 if SIN or
46: * COS is requested. Otherwise, for SINCOS, sin(X) is returned
47: * in Fp0, and cos(X) is returned in Fp1.
48: *
49: * Modifies: Fp0 for SIN or COS; both Fp0 and Fp1 for SINCOS.
50: *
51: * Accuracy and Monotonicity: The returned result is within 1 ulp in
52: * 64 significant bit, i.e. within 0.5001 ulp to 53 bits if the
53: * result is subsequently rounded to double precision. The
54: * result is provably monotonic in double precision.
55: *
56: * Speed: The programs sSIN and sCOS take approximately 150 cycles for
57: * input argument X such that |X| < 15Pi, which is the usual
58: * situation. The speed for sSINCOS is approximately 190 cycles.
59: *
60: * Algorithm:
61: *
62: * SIN and COS:
63: * 1. If SIN is invoked, set AdjN := 0; otherwise, set AdjN := 1.
64: *
65: * 2. If |X| >= 15Pi or |X| < 2**(-40), go to 7.
66: *
67: * 3. Decompose X as X = N(Pi/2) + r where |r| <= Pi/4. Let
68: * k = N mod 4, so in particular, k = 0,1,2,or 3. Overwirte
69: * k by k := k + AdjN.
70: *
71: * 4. If k is even, go to 6.
72: *
73: * 5. (k is odd) Set j := (k-1)/2, sgn := (-1)**j. Return sgn*cos(r)
74: * where cos(r) is approximated by an even polynomial in r,
75: * 1 + r*r*(B1+s*(B2+ ... + s*B8)), s = r*r.
76: * Exit.
77: *
78: * 6. (k is even) Set j := k/2, sgn := (-1)**j. Return sgn*sin(r)
79: * where sin(r) is approximated by an odd polynomial in r
80: * r + r*s*(A1+s*(A2+ ... + s*A7)), s = r*r.
81: * Exit.
82: *
83: * 7. If |X| > 1, go to 9.
84: *
85: * 8. (|X|<2**(-40)) If SIN is invoked, return X; otherwise return 1.
86: *
87: * 9. Overwrite X by X := X rem 2Pi. Now that |X| <= Pi, go back to 3.
88: *
89: * SINCOS:
90: * 1. If |X| >= 15Pi or |X| < 2**(-40), go to 6.
91: *
92: * 2. Decompose X as X = N(Pi/2) + r where |r| <= Pi/4. Let
93: * k = N mod 4, so in particular, k = 0,1,2,or 3.
94: *
95: * 3. If k is even, go to 5.
96: *
97: * 4. (k is odd) Set j1 := (k-1)/2, j2 := j1 (EOR) (k mod 2), i.e.
98: * j1 exclusive or with the l.s.b. of k.
99: * sgn1 := (-1)**j1, sgn2 := (-1)**j2.
100: * SIN(X) = sgn1 * cos(r) and COS(X) = sgn2*sin(r) where
101: * sin(r) and cos(r) are computed as odd and even polynomials
102: * in r, respectively. Exit
103: *
104: * 5. (k is even) Set j1 := k/2, sgn1 := (-1)**j1.
105: * SIN(X) = sgn1 * sin(r) and COS(X) = sgn1*cos(r) where
106: * sin(r) and cos(r) are computed as odd and even polynomials
107: * in r, respectively. Exit
108: *
109: * 6. If |X| > 1, go to 8.
110: *
111: * 7. (|X|<2**(-40)) SIN(X) = X and COS(X) = 1. Exit.
112: *
113: * 8. Overwrite X by X := X rem 2Pi. Now that |X| <= Pi, go back to 2.
114: *
115:
116: SSIN IDNT 2,1 Motorola 040 Floating Point Software Package
117:
118: section 8
119:
120: include fpsp.h
121:
122: BOUNDS1 DC.L $3FD78000,$4004BC7E
123: TWOBYPI DC.L $3FE45F30,$6DC9C883
124:
125: SINA7 DC.L $BD6AAA77,$CCC994F5
126: SINA6 DC.L $3DE61209,$7AAE8DA1
127:
128: SINA5 DC.L $BE5AE645,$2A118AE4
129: SINA4 DC.L $3EC71DE3,$A5341531
130:
131: SINA3 DC.L $BF2A01A0,$1A018B59,$00000000,$00000000
132:
133: SINA2 DC.L $3FF80000,$88888888,$888859AF,$00000000
134:
135: SINA1 DC.L $BFFC0000,$AAAAAAAA,$AAAAAA99,$00000000
136:
137: COSB8 DC.L $3D2AC4D0,$D6011EE3
138: COSB7 DC.L $BDA9396F,$9F45AC19
139:
140: COSB6 DC.L $3E21EED9,$0612C972
141: COSB5 DC.L $BE927E4F,$B79D9FCF
142:
143: COSB4 DC.L $3EFA01A0,$1A01D423,$00000000,$00000000
144:
145: COSB3 DC.L $BFF50000,$B60B60B6,$0B61D438,$00000000
146:
147: COSB2 DC.L $3FFA0000,$AAAAAAAA,$AAAAAB5E
148: COSB1 DC.L $BF000000
149:
150: INVTWOPI DC.L $3FFC0000,$A2F9836E,$4E44152A
151:
152: TWOPI1 DC.L $40010000,$C90FDAA2,$00000000,$00000000
153: TWOPI2 DC.L $3FDF0000,$85A308D4,$00000000,$00000000
154:
155: xref PITBL
156:
157: INARG equ FP_SCR4
158:
159: X equ FP_SCR5
160: XDCARE equ X+2
161: XFRAC equ X+4
162:
163: RPRIME equ FP_SCR1
164: SPRIME equ FP_SCR2
165:
166: POSNEG1 equ L_SCR1
167: TWOTO63 equ L_SCR1
168:
169: ENDFLAG equ L_SCR2
170: N equ L_SCR2
171:
172: ADJN equ L_SCR3
173:
174: xref t_frcinx
175: xref t_extdnrm
176: xref sto_cos
177:
178: xdef ssind
179: ssind:
180: *--SIN(X) = X FOR DENORMALIZED X
181: bra t_extdnrm
182:
183: xdef scosd
184: scosd:
185: *--COS(X) = 1 FOR DENORMALIZED X
186:
187: FMOVE.S #:3F800000,FP0
188: *
189: * 9D25B Fix: Sometimes the previous fmove.s sets fpsr bits
190: *
191: fmove.l #0,fpsr
192: *
193: bra t_frcinx
194:
195: xdef ssin
196: ssin:
197: *--SET ADJN TO 0
198: CLR.L ADJN(a6)
199: BRA.B SINBGN
200:
201: xdef scos
202: scos:
203: *--SET ADJN TO 1
204: MOVE.L #1,ADJN(a6)
205:
206: SINBGN:
207: *--SAVE FPCR, FP1. CHECK IF |X| IS TOO SMALL OR LARGE
208:
209: FMOVE.X (a0),FP0 ...LOAD INPUT
210:
211: MOVE.L (A0),D0
212: MOVE.W 4(A0),D0
213: FMOVE.X FP0,X(a6)
214: ANDI.L #$7FFFFFFF,D0 ...COMPACTIFY X
215:
216: CMPI.L #$3FD78000,D0 ...|X| >= 2**(-40)?
217: BGE.B SOK1
218: BRA.W SINSM
219:
220: SOK1:
221: CMPI.L #$4004BC7E,D0 ...|X| < 15 PI?
222: BLT.B SINMAIN
223: BRA.W REDUCEX
224:
225: SINMAIN:
226: *--THIS IS THE USUAL CASE, |X| <= 15 PI.
227: *--THE ARGUMENT REDUCTION IS DONE BY TABLE LOOK UP.
228: FMOVE.X FP0,FP1
229: FMUL.D TWOBYPI,FP1 ...X*2/PI
230:
231: *--HIDE THE NEXT THREE INSTRUCTIONS
232: LEA PITBL+$200,A1 ...TABLE OF N*PI/2, N = -32,...,32
233:
234:
235: *--FP1 IS NOW READY
236: FMOVE.L FP1,N(a6) ...CONVERT TO INTEGER
237:
238: MOVE.L N(a6),D0
239: ASL.L #4,D0
240: ADDA.L D0,A1 ...A1 IS THE ADDRESS OF N*PIBY2
241: * ...WHICH IS IN TWO PIECES Y1 & Y2
242:
243: FSUB.X (A1)+,FP0 ...X-Y1
244: *--HIDE THE NEXT ONE
245: FSUB.S (A1),FP0 ...FP0 IS R = (X-Y1)-Y2
246:
247: SINCONT:
248: *--continuation from REDUCEX
249:
250: *--GET N+ADJN AND SEE IF SIN(R) OR COS(R) IS NEEDED
251: MOVE.L N(a6),D0
252: ADD.L ADJN(a6),D0 ...SEE IF D0 IS ODD OR EVEN
253: ROR.L #1,D0 ...D0 WAS ODD IFF D0 IS NEGATIVE
254: TST.L D0
255: BLT.W COSPOLY
256:
257: SINPOLY:
258: *--LET J BE THE LEAST SIG. BIT OF D0, LET SGN := (-1)**J.
259: *--THEN WE RETURN SGN*SIN(R). SGN*SIN(R) IS COMPUTED BY
260: *--R' + R'*S*(A1 + S(A2 + S(A3 + S(A4 + ... + SA7)))), WHERE
261: *--R' = SGN*R, S=R*R. THIS CAN BE REWRITTEN AS
262: *--R' + R'*S*( [A1+T(A3+T(A5+TA7))] + [S(A2+T(A4+TA6))])
263: *--WHERE T=S*S.
264: *--NOTE THAT A3 THROUGH A7 ARE STORED IN DOUBLE PRECISION
265: *--WHILE A1 AND A2 ARE IN DOUBLE-EXTENDED FORMAT.
266: FMOVE.X FP0,X(a6) ...X IS R
267: FMUL.X FP0,FP0 ...FP0 IS S
268: *---HIDE THE NEXT TWO WHILE WAITING FOR FP0
269: FMOVE.D SINA7,FP3
270: FMOVE.D SINA6,FP2
271: *--FP0 IS NOW READY
272: FMOVE.X FP0,FP1
273: FMUL.X FP1,FP1 ...FP1 IS T
274: *--HIDE THE NEXT TWO WHILE WAITING FOR FP1
275:
276: ROR.L #1,D0
277: ANDI.L #$80000000,D0
278: * ...LEAST SIG. BIT OF D0 IN SIGN POSITION
279: EOR.L D0,X(a6) ...X IS NOW R'= SGN*R
280:
281: FMUL.X FP1,FP3 ...TA7
282: FMUL.X FP1,FP2 ...TA6
283:
284: FADD.D SINA5,FP3 ...A5+TA7
285: FADD.D SINA4,FP2 ...A4+TA6
286:
287: FMUL.X FP1,FP3 ...T(A5+TA7)
288: FMUL.X FP1,FP2 ...T(A4+TA6)
289:
290: FADD.D SINA3,FP3 ...A3+T(A5+TA7)
291: FADD.X SINA2,FP2 ...A2+T(A4+TA6)
292:
293: FMUL.X FP3,FP1 ...T(A3+T(A5+TA7))
294:
295: FMUL.X FP0,FP2 ...S(A2+T(A4+TA6))
296: FADD.X SINA1,FP1 ...A1+T(A3+T(A5+TA7))
297: FMUL.X X(a6),FP0 ...R'*S
298:
299: FADD.X FP2,FP1 ...[A1+T(A3+T(A5+TA7))]+[S(A2+T(A4+TA6))]
300: *--FP3 RELEASED, RESTORE NOW AND TAKE SOME ADVANTAGE OF HIDING
301: *--FP2 RELEASED, RESTORE NOW AND TAKE FULL ADVANTAGE OF HIDING
302:
303:
304: FMUL.X FP1,FP0 ...SIN(R')-R'
305: *--FP1 RELEASED.
306:
307: FMOVE.L d1,FPCR ;restore users exceptions
308: FADD.X X(a6),FP0 ;last inst - possible exception set
309: bra t_frcinx
310:
311:
312: COSPOLY:
313: *--LET J BE THE LEAST SIG. BIT OF D0, LET SGN := (-1)**J.
314: *--THEN WE RETURN SGN*COS(R). SGN*COS(R) IS COMPUTED BY
315: *--SGN + S'*(B1 + S(B2 + S(B3 + S(B4 + ... + SB8)))), WHERE
316: *--S=R*R AND S'=SGN*S. THIS CAN BE REWRITTEN AS
317: *--SGN + S'*([B1+T(B3+T(B5+TB7))] + [S(B2+T(B4+T(B6+TB8)))])
318: *--WHERE T=S*S.
319: *--NOTE THAT B4 THROUGH B8 ARE STORED IN DOUBLE PRECISION
320: *--WHILE B2 AND B3 ARE IN DOUBLE-EXTENDED FORMAT, B1 IS -1/2
321: *--AND IS THEREFORE STORED AS SINGLE PRECISION.
322:
323: FMUL.X FP0,FP0 ...FP0 IS S
324: *---HIDE THE NEXT TWO WHILE WAITING FOR FP0
325: FMOVE.D COSB8,FP2
326: FMOVE.D COSB7,FP3
327: *--FP0 IS NOW READY
328: FMOVE.X FP0,FP1
329: FMUL.X FP1,FP1 ...FP1 IS T
330: *--HIDE THE NEXT TWO WHILE WAITING FOR FP1
331: FMOVE.X FP0,X(a6) ...X IS S
332: ROR.L #1,D0
333: ANDI.L #$80000000,D0
334: * ...LEAST SIG. BIT OF D0 IN SIGN POSITION
335:
336: FMUL.X FP1,FP2 ...TB8
337: *--HIDE THE NEXT TWO WHILE WAITING FOR THE XU
338: EOR.L D0,X(a6) ...X IS NOW S'= SGN*S
339: ANDI.L #$80000000,D0
340:
341: FMUL.X FP1,FP3 ...TB7
342: *--HIDE THE NEXT TWO WHILE WAITING FOR THE XU
343: ORI.L #$3F800000,D0 ...D0 IS SGN IN SINGLE
344: MOVE.L D0,POSNEG1(a6)
345:
346: FADD.D COSB6,FP2 ...B6+TB8
347: FADD.D COSB5,FP3 ...B5+TB7
348:
349: FMUL.X FP1,FP2 ...T(B6+TB8)
350: FMUL.X FP1,FP3 ...T(B5+TB7)
351:
352: FADD.D COSB4,FP2 ...B4+T(B6+TB8)
353: FADD.X COSB3,FP3 ...B3+T(B5+TB7)
354:
355: FMUL.X FP1,FP2 ...T(B4+T(B6+TB8))
356: FMUL.X FP3,FP1 ...T(B3+T(B5+TB7))
357:
358: FADD.X COSB2,FP2 ...B2+T(B4+T(B6+TB8))
359: FADD.S COSB1,FP1 ...B1+T(B3+T(B5+TB7))
360:
361: FMUL.X FP2,FP0 ...S(B2+T(B4+T(B6+TB8)))
362: *--FP3 RELEASED, RESTORE NOW AND TAKE SOME ADVANTAGE OF HIDING
363: *--FP2 RELEASED.
364:
365:
366: FADD.X FP1,FP0
367: *--FP1 RELEASED
368:
369: FMUL.X X(a6),FP0
370:
371: FMOVE.L d1,FPCR ;restore users exceptions
372: FADD.S POSNEG1(a6),FP0 ;last inst - possible exception set
373: bra t_frcinx
374:
375:
376: SINBORS:
377: *--IF |X| > 15PI, WE USE THE GENERAL ARGUMENT REDUCTION.
378: *--IF |X| < 2**(-40), RETURN X OR 1.
379: CMPI.L #$3FFF8000,D0
380: BGT.B REDUCEX
381:
382:
383: SINSM:
384: MOVE.L ADJN(a6),D0
385: TST.L D0
386: BGT.B COSTINY
387:
388: SINTINY:
389: CLR.W XDCARE(a6) ...JUST IN CASE
390: FMOVE.L d1,FPCR ;restore users exceptions
391: FMOVE.X X(a6),FP0 ;last inst - possible exception set
392: bra t_frcinx
393:
394:
395: COSTINY:
396: FMOVE.S #:3F800000,FP0
397:
398: FMOVE.L d1,FPCR ;restore users exceptions
399: FSUB.S #:00800000,FP0 ;last inst - possible exception set
400: bra t_frcinx
401:
402:
403: REDUCEX:
404: *--WHEN REDUCEX IS USED, THE CODE WILL INEVITABLY BE SLOW.
405: *--THIS REDUCTION METHOD, HOWEVER, IS MUCH FASTER THAN USING
406: *--THE REMAINDER INSTRUCTION WHICH IS NOW IN SOFTWARE.
407:
408: FMOVEM.X FP2-FP5,-(A7) ...save FP2 through FP5
409: MOVE.L D2,-(A7)
410: FMOVE.S #:00000000,FP1
411: *--If compact form of abs(arg) in d0=$7ffeffff, argument is so large that
412: *--there is a danger of unwanted overflow in first LOOP iteration. In this
413: *--case, reduce argument by one remainder step to make subsequent reduction
414: *--safe.
415: cmpi.l #$7ffeffff,d0 ;is argument dangerously large?
416: bne.b LOOP
417: move.l #$7ffe0000,FP_SCR2(a6) ;yes
418: * ;create 2**16383*PI/2
419: move.l #$c90fdaa2,FP_SCR2+4(a6)
420: clr.l FP_SCR2+8(a6)
421: ftst.x fp0 ;test sign of argument
422: move.l #$7fdc0000,FP_SCR3(a6) ;create low half of 2**16383*
423: * ;PI/2 at FP_SCR3
424: move.l #$85a308d3,FP_SCR3+4(a6)
425: clr.l FP_SCR3+8(a6)
426: fblt.w red_neg
427: or.w #$8000,FP_SCR2(a6) ;positive arg
428: or.w #$8000,FP_SCR3(a6)
429: red_neg:
430: fadd.x FP_SCR2(a6),fp0 ;high part of reduction is exact
431: fmove.x fp0,fp1 ;save high result in fp1
432: fadd.x FP_SCR3(a6),fp0 ;low part of reduction
433: fsub.x fp0,fp1 ;determine low component of result
434: fadd.x FP_SCR3(a6),fp1 ;fp0/fp1 are reduced argument.
435:
436: *--ON ENTRY, FP0 IS X, ON RETURN, FP0 IS X REM PI/2, |X| <= PI/4.
437: *--integer quotient will be stored in N
438: *--Intermeditate remainder is 66-bit long; (R,r) in (FP0,FP1)
439:
440: LOOP:
441: FMOVE.X FP0,INARG(a6) ...+-2**K * F, 1 <= F < 2
442: MOVE.W INARG(a6),D0
443: MOVE.L D0,A1 ...save a copy of D0
444: ANDI.L #$00007FFF,D0
445: SUBI.L #$00003FFF,D0 ...D0 IS K
446: CMPI.L #28,D0
447: BLE.B LASTLOOP
448: CONTLOOP:
449: SUBI.L #27,D0 ...D0 IS L := K-27
450: CLR.L ENDFLAG(a6)
451: BRA.B WORK
452: LASTLOOP:
453: CLR.L D0 ...D0 IS L := 0
454: MOVE.L #1,ENDFLAG(a6)
455:
456: WORK:
457: *--FIND THE REMAINDER OF (R,r) W.R.T. 2**L * (PI/2). L IS SO CHOSEN
458: *--THAT INT( X * (2/PI) / 2**(L) ) < 2**29.
459:
460: *--CREATE 2**(-L) * (2/PI), SIGN(INARG)*2**(63),
461: *--2**L * (PIby2_1), 2**L * (PIby2_2)
462:
463: MOVE.L #$00003FFE,D2 ...BIASED EXPO OF 2/PI
464: SUB.L D0,D2 ...BIASED EXPO OF 2**(-L)*(2/PI)
465:
466: MOVE.L #$A2F9836E,FP_SCR1+4(a6)
467: MOVE.L #$4E44152A,FP_SCR1+8(a6)
468: MOVE.W D2,FP_SCR1(a6) ...FP_SCR1 is 2**(-L)*(2/PI)
469:
470: FMOVE.X FP0,FP2
471: FMUL.X FP_SCR1(a6),FP2
472: *--WE MUST NOW FIND INT(FP2). SINCE WE NEED THIS VALUE IN
473: *--FLOATING POINT FORMAT, THE TWO FMOVE'S FMOVE.L FP <--> N
474: *--WILL BE TOO INEFFICIENT. THE WAY AROUND IT IS THAT
475: *--(SIGN(INARG)*2**63 + FP2) - SIGN(INARG)*2**63 WILL GIVE
476: *--US THE DESIRED VALUE IN FLOATING POINT.
477:
478: *--HIDE SIX CYCLES OF INSTRUCTION
479: MOVE.L A1,D2
480: SWAP D2
481: ANDI.L #$80000000,D2
482: ORI.L #$5F000000,D2 ...D2 IS SIGN(INARG)*2**63 IN SGL
483: MOVE.L D2,TWOTO63(a6)
484:
485: MOVE.L D0,D2
486: ADDI.L #$00003FFF,D2 ...BIASED EXPO OF 2**L * (PI/2)
487:
488: *--FP2 IS READY
489: FADD.S TWOTO63(a6),FP2 ...THE FRACTIONAL PART OF FP1 IS ROUNDED
490:
491: *--HIDE 4 CYCLES OF INSTRUCTION; creating 2**(L)*Piby2_1 and 2**(L)*Piby2_2
492: MOVE.W D2,FP_SCR2(a6)
493: CLR.W FP_SCR2+2(a6)
494: MOVE.L #$C90FDAA2,FP_SCR2+4(a6)
495: CLR.L FP_SCR2+8(a6) ...FP_SCR2 is 2**(L) * Piby2_1
496:
497: *--FP2 IS READY
498: FSUB.S TWOTO63(a6),FP2 ...FP2 is N
499:
500: ADDI.L #$00003FDD,D0
501: MOVE.W D0,FP_SCR3(a6)
502: CLR.W FP_SCR3+2(a6)
503: MOVE.L #$85A308D3,FP_SCR3+4(a6)
504: CLR.L FP_SCR3+8(a6) ...FP_SCR3 is 2**(L) * Piby2_2
505:
506: MOVE.L ENDFLAG(a6),D0
507:
508: *--We are now ready to perform (R+r) - N*P1 - N*P2, P1 = 2**(L) * Piby2_1 and
509: *--P2 = 2**(L) * Piby2_2
510: FMOVE.X FP2,FP4
511: FMul.X FP_SCR2(a6),FP4 ...W = N*P1
512: FMove.X FP2,FP5
513: FMul.X FP_SCR3(a6),FP5 ...w = N*P2
514: FMove.X FP4,FP3
515: *--we want P+p = W+w but |p| <= half ulp of P
516: *--Then, we need to compute A := R-P and a := r-p
517: FAdd.X FP5,FP3 ...FP3 is P
518: FSub.X FP3,FP4 ...W-P
519:
520: FSub.X FP3,FP0 ...FP0 is A := R - P
521: FAdd.X FP5,FP4 ...FP4 is p = (W-P)+w
522:
523: FMove.X FP0,FP3 ...FP3 A
524: FSub.X FP4,FP1 ...FP1 is a := r - p
525:
526: *--Now we need to normalize (A,a) to "new (R,r)" where R+r = A+a but
527: *--|r| <= half ulp of R.
528: FAdd.X FP1,FP0 ...FP0 is R := A+a
529: *--No need to calculate r if this is the last loop
530: TST.L D0
531: BGT.W RESTORE
532:
533: *--Need to calculate r
534: FSub.X FP0,FP3 ...A-R
535: FAdd.X FP3,FP1 ...FP1 is r := (A-R)+a
536: BRA.W LOOP
537:
538: RESTORE:
539: FMOVE.L FP2,N(a6)
540: MOVE.L (A7)+,D2
541: FMOVEM.X (A7)+,FP2-FP5
542:
543:
544: MOVE.L ADJN(a6),D0
545: CMPI.L #4,D0
546:
547: BLT.W SINCONT
548: BRA.B SCCONT
549:
550: xdef ssincosd
551: ssincosd:
552: *--SIN AND COS OF X FOR DENORMALIZED X
553:
554: FMOVE.S #:3F800000,FP1
555: bsr sto_cos ;store cosine result
556: bra t_extdnrm
557:
558: xdef ssincos
559: ssincos:
560: *--SET ADJN TO 4
561: MOVE.L #4,ADJN(a6)
562:
563: FMOVE.X (a0),FP0 ...LOAD INPUT
564:
565: MOVE.L (A0),D0
566: MOVE.W 4(A0),D0
567: FMOVE.X FP0,X(a6)
568: ANDI.L #$7FFFFFFF,D0 ...COMPACTIFY X
569:
570: CMPI.L #$3FD78000,D0 ...|X| >= 2**(-40)?
571: BGE.B SCOK1
572: BRA.W SCSM
573:
574: SCOK1:
575: CMPI.L #$4004BC7E,D0 ...|X| < 15 PI?
576: BLT.B SCMAIN
577: BRA.W REDUCEX
578:
579:
580: SCMAIN:
581: *--THIS IS THE USUAL CASE, |X| <= 15 PI.
582: *--THE ARGUMENT REDUCTION IS DONE BY TABLE LOOK UP.
583: FMOVE.X FP0,FP1
584: FMUL.D TWOBYPI,FP1 ...X*2/PI
585:
586: *--HIDE THE NEXT THREE INSTRUCTIONS
587: LEA PITBL+$200,A1 ...TABLE OF N*PI/2, N = -32,...,32
588:
589:
590: *--FP1 IS NOW READY
591: FMOVE.L FP1,N(a6) ...CONVERT TO INTEGER
592:
593: MOVE.L N(a6),D0
594: ASL.L #4,D0
595: ADDA.L D0,A1 ...ADDRESS OF N*PIBY2, IN Y1, Y2
596:
597: FSUB.X (A1)+,FP0 ...X-Y1
598: FSUB.S (A1),FP0 ...FP0 IS R = (X-Y1)-Y2
599:
600: SCCONT:
601: *--continuation point from REDUCEX
602:
603: *--HIDE THE NEXT TWO
604: MOVE.L N(a6),D0
605: ROR.L #1,D0
606:
607: TST.L D0 ...D0 < 0 IFF N IS ODD
608: BGE.W NEVEN
609:
610: NODD:
611: *--REGISTERS SAVED SO FAR: D0, A0, FP2.
612:
613: FMOVE.X FP0,RPRIME(a6)
614: FMUL.X FP0,FP0 ...FP0 IS S = R*R
615: FMOVE.D SINA7,FP1 ...A7
616: FMOVE.D COSB8,FP2 ...B8
617: FMUL.X FP0,FP1 ...SA7
618: MOVE.L d2,-(A7)
619: MOVE.L D0,d2
620: FMUL.X FP0,FP2 ...SB8
621: ROR.L #1,d2
622: ANDI.L #$80000000,d2
623:
624: FADD.D SINA6,FP1 ...A6+SA7
625: EOR.L D0,d2
626: ANDI.L #$80000000,d2
627: FADD.D COSB7,FP2 ...B7+SB8
628:
629: FMUL.X FP0,FP1 ...S(A6+SA7)
630: EOR.L d2,RPRIME(a6)
631: MOVE.L (A7)+,d2
632: FMUL.X FP0,FP2 ...S(B7+SB8)
633: ROR.L #1,D0
634: ANDI.L #$80000000,D0
635:
636: FADD.D SINA5,FP1 ...A5+S(A6+SA7)
637: MOVE.L #$3F800000,POSNEG1(a6)
638: EOR.L D0,POSNEG1(a6)
639: FADD.D COSB6,FP2 ...B6+S(B7+SB8)
640:
641: FMUL.X FP0,FP1 ...S(A5+S(A6+SA7))
642: FMUL.X FP0,FP2 ...S(B6+S(B7+SB8))
643: FMOVE.X FP0,SPRIME(a6)
644:
645: FADD.D SINA4,FP1 ...A4+S(A5+S(A6+SA7))
646: EOR.L D0,SPRIME(a6)
647: FADD.D COSB5,FP2 ...B5+S(B6+S(B7+SB8))
648:
649: FMUL.X FP0,FP1 ...S(A4+...)
650: FMUL.X FP0,FP2 ...S(B5+...)
651:
652: FADD.D SINA3,FP1 ...A3+S(A4+...)
653: FADD.D COSB4,FP2 ...B4+S(B5+...)
654:
655: FMUL.X FP0,FP1 ...S(A3+...)
656: FMUL.X FP0,FP2 ...S(B4+...)
657:
658: FADD.X SINA2,FP1 ...A2+S(A3+...)
659: FADD.X COSB3,FP2 ...B3+S(B4+...)
660:
661: FMUL.X FP0,FP1 ...S(A2+...)
662: FMUL.X FP0,FP2 ...S(B3+...)
663:
664: FADD.X SINA1,FP1 ...A1+S(A2+...)
665: FADD.X COSB2,FP2 ...B2+S(B3+...)
666:
667: FMUL.X FP0,FP1 ...S(A1+...)
668: FMUL.X FP2,FP0 ...S(B2+...)
669:
670:
671:
672: FMUL.X RPRIME(a6),FP1 ...R'S(A1+...)
673: FADD.S COSB1,FP0 ...B1+S(B2...)
674: FMUL.X SPRIME(a6),FP0 ...S'(B1+S(B2+...))
675:
676: move.l d1,-(sp) ;restore users mode & precision
677: andi.l #$ff,d1 ;mask off all exceptions
678: fmove.l d1,FPCR
679: FADD.X RPRIME(a6),FP1 ...COS(X)
680: bsr sto_cos ;store cosine result
681: FMOVE.L (sp)+,FPCR ;restore users exceptions
682: FADD.S POSNEG1(a6),FP0 ...SIN(X)
683:
684: bra t_frcinx
685:
686:
687: NEVEN:
688: *--REGISTERS SAVED SO FAR: FP2.
689:
690: FMOVE.X FP0,RPRIME(a6)
691: FMUL.X FP0,FP0 ...FP0 IS S = R*R
692: FMOVE.D COSB8,FP1 ...B8
693: FMOVE.D SINA7,FP2 ...A7
694: FMUL.X FP0,FP1 ...SB8
695: FMOVE.X FP0,SPRIME(a6)
696: FMUL.X FP0,FP2 ...SA7
697: ROR.L #1,D0
698: ANDI.L #$80000000,D0
699: FADD.D COSB7,FP1 ...B7+SB8
700: FADD.D SINA6,FP2 ...A6+SA7
701: EOR.L D0,RPRIME(a6)
702: EOR.L D0,SPRIME(a6)
703: FMUL.X FP0,FP1 ...S(B7+SB8)
704: ORI.L #$3F800000,D0
705: MOVE.L D0,POSNEG1(a6)
706: FMUL.X FP0,FP2 ...S(A6+SA7)
707:
708: FADD.D COSB6,FP1 ...B6+S(B7+SB8)
709: FADD.D SINA5,FP2 ...A5+S(A6+SA7)
710:
711: FMUL.X FP0,FP1 ...S(B6+S(B7+SB8))
712: FMUL.X FP0,FP2 ...S(A5+S(A6+SA7))
713:
714: FADD.D COSB5,FP1 ...B5+S(B6+S(B7+SB8))
715: FADD.D SINA4,FP2 ...A4+S(A5+S(A6+SA7))
716:
717: FMUL.X FP0,FP1 ...S(B5+...)
718: FMUL.X FP0,FP2 ...S(A4+...)
719:
720: FADD.D COSB4,FP1 ...B4+S(B5+...)
721: FADD.D SINA3,FP2 ...A3+S(A4+...)
722:
723: FMUL.X FP0,FP1 ...S(B4+...)
724: FMUL.X FP0,FP2 ...S(A3+...)
725:
726: FADD.X COSB3,FP1 ...B3+S(B4+...)
727: FADD.X SINA2,FP2 ...A2+S(A3+...)
728:
729: FMUL.X FP0,FP1 ...S(B3+...)
730: FMUL.X FP0,FP2 ...S(A2+...)
731:
732: FADD.X COSB2,FP1 ...B2+S(B3+...)
733: FADD.X SINA1,FP2 ...A1+S(A2+...)
734:
735: FMUL.X FP0,FP1 ...S(B2+...)
736: fmul.x fp2,fp0 ...s(a1+...)
737:
738:
739:
740: FADD.S COSB1,FP1 ...B1+S(B2...)
741: FMUL.X RPRIME(a6),FP0 ...R'S(A1+...)
742: FMUL.X SPRIME(a6),FP1 ...S'(B1+S(B2+...))
743:
744: move.l d1,-(sp) ;save users mode & precision
745: andi.l #$ff,d1 ;mask off all exceptions
746: fmove.l d1,FPCR
747: FADD.S POSNEG1(a6),FP1 ...COS(X)
748: bsr sto_cos ;store cosine result
749: FMOVE.L (sp)+,FPCR ;restore users exceptions
750: FADD.X RPRIME(a6),FP0 ...SIN(X)
751:
752: bra t_frcinx
753:
754: SCBORS:
755: CMPI.L #$3FFF8000,D0
756: BGT.W REDUCEX
757:
758:
759: SCSM:
760: CLR.W XDCARE(a6)
761: FMOVE.S #:3F800000,FP1
762:
763: move.l d1,-(sp) ;save users mode & precision
764: andi.l #$ff,d1 ;mask off all exceptions
765: fmove.l d1,FPCR
766: FSUB.S #:00800000,FP1
767: bsr sto_cos ;store cosine result
768: FMOVE.L (sp)+,FPCR ;restore users exceptions
769: FMOVE.X X(a6),FP0
770: bra t_frcinx
771:
772: end