mac-rom/Toolbox/InSANE/HWElemsSinCos.a
Elliot Nunn 4325cdcc78 Bring in CubeE sources
Resource forks are included only for .rsrc files. These are DeRezzed into their data fork. 'ckid' resources, from the Projector VCS, are not included.

The Tools directory, containing mostly junk, is also excluded.
2017-12-26 09:52:23 +08:00

414 lines
20 KiB
Plaintext

;
; File: HWElemsSinCos.a
;
; Contains: Routines to calculate SIN(x) and COS(x)
;
; Written by: Apple Numerics Group, DSG
;
; Copyright: © 1985-1993 by Apple Computer, Inc., all rights reserved.
;
; Change History (most recent first):
;
; <SM2> 2/3/92 CSS Update from Horror:
; <H2> 9/29/92 BG Adding Jon Okada's latest fixes.
; <1> 10/24/91 SAM/KSM Rolled in Regatta file.
;
; Terror Change History:
;
; <1> 01/06/91 BG Added to TERROR/BBS for the time.
;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; ;;
;; File 881ELEMSsincos.a, ;;
;; Functions sin(x) & cos(x), ;;
;; Implementation of sine and cosine for ½SANE. ;;
;; For machines with MC68881/68882 FPU. ;;
;; ;;
;; Copyright Apple Computer, Inc. 1985-7, 1989-90, 1992, 1993. ;;
;; ;;
;; Written by Ali Sazegari, started on September 5 1989, ;;
;; ;;
;; based on Elems881, package code for Macintosh by J. Coonen, ;;
;; who wrote it from pascal programs by D. Hough, C. McMaster and K. Hanson. ;;
;; ;;
;; Modification History: ;;
;; 23 Apr 92 Performance tuning done (better A/D register usage, faster ;;
;; FPU ops used, front-end checks for small operands eliminate ;;
;; spurious underflowing. ... JPO ;;
;; ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
RSin
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Classification of the argument with a neat trick. ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FMOVE.L FPSR,D0 ; what kind of arg? fpsr has the key
MOVEQ #0,D1 ; zero sign/quo flag <5/23/92, JPO>
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Here is the strategy ( commonly known as the scoop ) by C. R. Lewis: ;;
;; Register D0 has N, Z, I, NaN burried in it. Mask off unwanted bits to ;;
;; get them. Then, left shift NZINaN to have Z at MSB position. If D0 = 0 ;;
;; then arg is finite. ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; MOVE.L #ExcMask,D1 ; mask it kid DELETED <5/23/92, JPO>
; AND.L D1,D0 ; we now have the meat DELETED <5/23/92, JPO>
ANDI.L #ExcMask,D0 ; mask FPCC bits <5/23/92, JPO>
LSL.L #ShiftIt,D0 ; N->C, Z->MSB, I->30th & NaN->29th bit
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; If arg is a non zero finite number, great! let's fall through, otherwise ;;
;; take the problem cases. ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
BNE.B BadSinArg ; take care of the nasties (Z bar set)
; MOVEM.L A1,-(SP) ; use a1 for denormal flag use DELETED <5/23/92, JPO>
; MOVEA.L A0,A1 ; keep a copy of the pointer to argument DELETED <5/23/92, JPO>
; this is for exception use at the end
; MOVEM.L D2,-(SP) ; use d2 for sign of argument DELETED <5/23/92, JPO>
; BCC.B @60 ; is the argument negative? DELETED <5/23/92, JPO>
; MOVEQ #1,D2 ; argument was negative, then d2 = 1 DELETED <5/23/92, JPO>
; BRA.B @61 DELETED <5/23/92, JPO>
;@60 MOVEQ #0,D2 ; set to positive at first DELETED <5/23/92, JPO>
;@61 label DELETED <5/23/92, JPO>
BCC.B @smchk ; argument is positive & so is flag <5/23/92, JPO>
BSET.L #31,D1 ; argument is negative & so is flag <5/23/92, JPO>
@smchk: ; check for small magnitude input label ADDED <5/23/92, JPO>
BFEXTU (A0){1:15},D0 ; D0.W <- exponent of input <5/23/92, JPO>
CMPI.W #$3FDF,D0 ; | arg | < 2^(-32)? <5/23/92, JPO>
BGE.B @sinred ; no, reduce argument <5/23/92, JPO>
MOVE.L #SetInexact,D0 ; sin(tiny) = tiny, signal INEXACT <5/23/92, JPO>
BFTST (A0){1:16} ; signal UNDERFLOW if arg was denormal <5/23/92, JPO>
BNE.B @tinydone ; <5/23/92, JPO>
OR.L #SetUFlow,D0 ; <5/23/92, JPO>
@tinydone: ; label ADDED <5/23/92, JPO>
JMP (A4) ; done <5/23/92, JPO>
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; About the argument reduction: t is reduced REM approximate ¹/2, leaving ;;
;; its magnitude no bigger than approximate ¹/4. ;;
;; ;;
;; Recall that: ;;
;; ;;
;; sin(¹/2 + t) = cos(t), q = 1 ;;
;; sin(¹ + t) = - sin(t), q = 2 ;;
;; sin(3¹/2 + t) = - cos(t), q = 3 ;;
;; sin(2¹ + t) = sin(t). q = 0 ;;
;; ;;
;; Then if input t = q * (¹/2) + r, ;;
;; ;;
;; q MOD 2 determines whether to use sine or cosine: ;;
;; ;;
;; if ( q MOD 2 = 0 ) use sine, ;;
;; else ( q MOD 2 = 1 ) then use cosine. ;;
;; ;;
;; q MOD 4 = (2 or 3) determines whether to negate result . ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
@sinred: ; label ADDED <5/23/92, JPO>
; FMOVE.X FPKPI2,FP1 ; get the approximate ¹/2 DELETED <5/23/92, JPO>
; FREM.X FP1,FP0 ; t = arg REM ¹/2, quotient q in fpsr DELETED<5/23/92, JPO>
FREM.X FPKPI2,FP0 ; FP0 <- t = arg REM ¹/2, quo q in FPSR <5/23/92, JPO>
FMOVE.L FPSR,D0 ; q is in 1st byte of higher word
SWAP D0 ; q is in 1st byte of lower word
MOVE.B D0,D1 ; copy q for negation section
LSR.B #1,D0 ; is it 0 or 1? ( odd or even )
BCC.B @62 ; even case uses cosine approx
BSR.W UseCos ; q is odd then use cosine computation
; SUBQ.B #1,D2 ; check for negative argument DELETED <5/23/92, JPO>
; BMI.B ChkNegSin ; GET OUT, result should remain positive DELETED <5/23/92, JPO>
TST.L D1 ; check for negative argument <5/23/92, JPO>
BPL.B ChkNegSin ; positive argument <5/23/92, JPO>
; FMUL.W #-1,FP0 ; d2 = 1, then argument was negative DELETED <5/23/92, JPO>
FNEG.X FP0,FP0 ; negative argument so negate result <5/23/92, JPO>
BRA.B ChkNegSin ; JUMP OUT, do we negate the result?
@62 BSR.W UseSin ; q is even then use sine computation
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Determine if we have to negate the result; i. e., is q MOD 4 = 2 or 3? ;;
;; remember that we saved q ( the REM quotient ) in d1 for this section. ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
ChkNegSin
; MOVEM.L (SP)+,D2 ; restore register used for sign DELETED <5/23/92, JPO>
ANDI.B #3,D1 ; extract q = quo mod 4, mask off unwanted bits
SUBQ.B #2,D1 ; if q = 2 or 3 then D1.B>0, N bit not set
BMI.B @63 ; q = 0 or 1 so sign is OK
; FMUL.W #-1,FP0 ; q = 0 or 1 then negate the result DELETED <5/23/92, JPO>
FNEG.X FP0,FP0 ; q = 2 or 3, so negate result <5/23/92, JPO>
@63:
; BFEXTU (A1) {1:16},D0 ; get the exponent & 1st mantissa bit DELETED <5/23/92, JPO>
; MOVEA.L (SP)+,A1 ; restore the register after use DELETED <5/23/92, JPO>
; BEQ.B @64 ; if D0 = 0 then current s = denorm DELETED <5/23/92, JPO>
MOVE.L #SetInexact,D0 ; set the inexact flag
; MOVE.B #ClrUflow ,D0 ; and clear the underflow bit DELETED <5/23/92, JPO>
JMP (A4)
;@64 MOVE.L #SetInexact + SetUflow,D0 ; signal inexact DELETED <5/23/92, JPO>
; JMP (A4) ; out, out. DELETED <5/23/92, JPO>
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Problem cases: NaN, +° and zero. ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
BadSinArg BMI.B SinOut ; is arg = 0 ? ( N set )
BTST #29,D0 ; is arg a NaN? ( bit position 29 )
BNE.B SinOut ; scream, arg = NaN detected!
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; All right boys & girls, at this time it is apparent that the good old arg ;;
;; is an infinity. ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FMOVE.X #NaNTrig,FP0 ; arg = +° => sin(arg) = cos(arg) = NaN
MOVE.L #SetInvalid,D0 ; signal invalid
JMP (A4) ; return either a NaN or a +°
SinOut MOVEQ #0,D0 ; signal no exceptions
JMP (A4) ; return +-0
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; C O S I N E ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; The cosine section. it follows pretty much in the footsteps of sine with ;;
;; slight twists. Read the comments for argument reduction and compare with ;;
;; the one on sine. ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
RCos
FMOVE.L FPSR,D0 ; what kind of arg? fpsr has the key
MOVEQ #0,D1 ; clear sign/quo flag <5/23/92, JPO>
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Here is the strategy ( commonly known as the scoop ) by C. R. Lewis: ;;
;; Register D0 has N, Z, I, NaN burried in it. Mask off unwanted bits to ;;
;; get them. Then, left shift NZINaN to have Z at MSB position. If D0 = 0 ;;
;; then arg is finite. ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; MOVE.L #ExcMask,D1 ; mask it kid DELETED <5/23/92, JPO>
; AND.L D1,D0 ; we now have the meat DELETED <5/23/92, JPO>
AND.L #ExcMask,D0 ; isolate FPCC bits <5/23/92, JPO>
LSL.L #ShiftIt,D0 ; N->C, Z->MSB, I->30th & NaN->29th bit
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; If arg is a non zero finite number, great! let's fall through, otherwise ;;
;; take the problem cases. ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
BNE.B BadCosArgs ; take care of the nasties (Z bar set)
; MOVEM.L D2,-(SP) ; use d2 for sign of argument DELETED <5/23/92, JPO>
; BCC.B @50 ; is the argument negative? DELETED <5/23/92, JPO>
; MOVEQ #1,D2 ; argument was negative, then d2 = 1 DELETED <5/23/92, JPO>
; BRA.B @51 ; DELETED <5/23/92, JPO>
;@50 MOVEQ #0,D2 ; set to positive at first DELETED<5/23/92, JPO>
;@51 ; label DELETED <5/23/92, JPO>
BCC.B @smchk ; arg is positive & so is flag <5/23/92, JPO>
bset.l #31,D1 ; arg is negative & so is flag <5/23/92, JPO>
@smchk: ; label ADDED <5/23/92, JPO>
BFEXTU (A0){1:15},D0 ; is magnitude of arg < 2^(-32)? <5/23/92, JPO>
CMPI.W #$3FDF,D0 ; <5/23/92, JPO>
BGE.B @cosred ; no, reduce argument <5/23/92, JPO>
FMOVE.W #1,FP0 ; yes, deliver +1.0 for cos <5/23/92, JPO>
MOVE.L #SetInexact,D0 ; signal INEXACT
JMP (A4) ; get out
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; About the argument reduction: just as in the sin case, t is reduced REM ;;
;; approximate ¹/2, leaving its magnitude no bigger than approximate ¹/4. ;;
;; ;;
;; Recall that: ;;
;; ;;
;; cos(¹/2 + t) = - sin(t), q = 1 ;;
;; cos(¹ + t) = - cos(t), q = 2 ;;
;; cos(3¹/2 + t) = sin(t), q = 3 ;;
;; cos(2¹ + t) = cos(t). q = 0 ;;
;; ;;
;; Then if input t = q * (¹/2) + r, ;;
;; ;;
;; q MOD 2 determines whether to use sine or cosine: ;;
;; ;;
;; if ( q MOD 2 = 0 ) use cosine, ;;
;; else ( q MOD 2 = 1 ) then use sine. ;;
;; ;;
;; q MOD 4 = (2 or 3) determines whether to negate result . ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
@cosred: ; label ADDED <5/23/92, JPO>
; FMOVE.X FPKPI2,FP1 ; get the approximate ¹/2 DELETED <5/23/92, JPO>
; FREM.X FP1,FP0 ; t = arg REM ¹/2, quotient q in fpsr DELETED <5/23/92, JPO>
FREM.X FPKPI2,FP0 ; t = arg REM ¹/2, quotient q in fpsr <5/23/92, JPO>
FMOVE.L FPSR,D0 ; q is in 1st byte of higher word
SWAP D0 ; q is in 1st byte of lower word
MOVE.B D0,D1 ; copy q for negation section
LSR.B #1,D0 ; is it 0 or 1? ( odd or even )
BCS.B @52 ; q even then use sine computation
BSR.B UseCos ; q odd then use cosine computation
BRA.B ChkNegCos ; computation done, jump out
@52 BSR.W UseSin
; SUBQ.B #1,D2 ; check for negative argument DELETED <5/23/92, JPO>
; BMI.B ChkNegCos ; result should remain positive, get out DELETED <5/23/92, JPO>
TST.L D1 ; negative argument? <5/23/92, JPO>
BPL.B ChkNegCos ; no, don't negate <5/23/92, JPO>
; FMUL.W #-1,FP0 ; d2 = 1, then argument was negative DELETED <5/23/92, JPO>
FNEG.X FP0,FP0 ; yes, negate result <5/23/92, JPO>
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Determine if we have to negate the result; i. e., is q MOD 4 = 1 or 2? ;;
;; remember that we saved q ( the REM quotient ) in d1 for this section. ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
ChkNegCos
; MOVEM.L (SP)+,D2 ; restore register used for sign DELETED <5/23/92, JPO>
ANDI.B #3,D1 ; extract q, mask off unwanted bits
; SUBQ.B #1,D1 ; if q = 0 then d1<0, N bit not set DELETED <5/23/92, JPO>
; BMI.B @53 ; result is to remain positive DELETED <5/23/92, JPO>
; SUBQ.B #2,D1 ; if q was 3 then d1>0, N bit not set DELETED <5/23/92, JPO>
; BPL.B @53 ; result is to remain positive DELETED <5/23/92, JPO>
; FMUL.W #-1,FP0 ; q = 0 or 1 then negate the result DELETED <5/23/92, JPO>
BEQ.B @53 ; q = 0, don't negate <5/23/92, JPO>
CMPI.B #3,D1 ; <5/23/92, JPO>
BEQ.B @53 ; q = 3, don't negate <5/23/92, JPO>
FNEG.X FP0,FP0 ; q = 1 or 2, negate <5/23/92, JPO>
@53 MOVE.L #SetInexact,D0 ; signal inexact
; MOVE.B #ClrUflow ,D0 ; and always clear the underflow bit DELETED <5/23/92, JPO>
JMP (A4) ; out, out.
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Problem casesfor the cosine: NaN, +° and zero. ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
BadCosArgs BMI.B CosIsOne ; is arg = 0 ? ( N set )
BTST #29,D0 ; is arg a NaN? ( bit position 29 )
BNE.B CosOut ; scream, arg = NaN detected!
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; All right boys & girls, at this time it is apparent that the good old arg ;;
;; is an infinity or a NaN. ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FMOVE.X #NaNTrig,FP0 ; arg = +° => sin(arg) = NaN
MOVE.L #SetInvalid,D0 ; signal invalid
JMP (A4) ; return either a NaN or a +°
;CosOut MOVEQ #0,D0 ; signal no exceptions DELETED <5/23/92, JPO>
; JMP (A4) ; return +-0 DELETED <5/23/92, JPO>
CosIsOne FMOVE.W #1,FP0 ; cos(0) = 1
CosOut: ; label ADDED <5/23/92, JPO>
MOVEQ #0,D0 ; signal no exceptions
JMP (A4) ; return 1
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Subroutine UseCos, called by sine and cosine routines. ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Evaluate cos(t) for reduced |t| <= ¹/4. ;;
;; ;;
;; Use the following approximation ( note that s = t * t ): ;;
;; if ( s < 1/4 ) then: ;;
;; cos(t) = 1 - s/2 + ( s * s ) * ( P(s)/Q(s) ), ;;
;; else: ;;
;; cos(t) = 0.875 - ( z/2 + ( z * (z/2) - ( s * ( s * ( P(s)/Q(s) )))) ;;
;; where z = |t| - 0.5. ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
UseCos ; start of the cosine computation
FMOVE.X FP0,FP3 ; fp0 = t, fp1 = Pi/2, fp3 = t
FMUL.X FP0,FP0 ; fp0 = s = t*t, fp1 = Pi/2, fp3 = t
LEA COSP,A0 ; A0 = coeff table ptr for PolyEval
BSR PolyEval ; fp1 = PolyEval(P, t*t)
FMOVE.X FP1,FP2 ; save P(arg)
LEA COSQ,A0 ; A0 = coeff table ptr for PolyEval
BSR PolyEval ; fp1 = PolyEval(Q, t*t)
FDIV.X FP1,FP2 ; form fp2 = P(s)/Q(s)
FMUL.X FP0,FP2 ; fp2 = ( t * t ) * P(s)/Q(s)
FMUL.X FP0,FP2 ; fp2 = s * ( s * P(s)/Q(s) )
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Check for the formula to be used. ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; FCMP.X #"0.25",FP0 ; is ( s < 1/4 ) ? DELETED <5/23/92, JPO>
FCMP.S #"$3E800000",FP0 ; s < 0.25? <5/23/92, JPO>
FBLT.W SmallCos ; send it off to the small formula
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; The more complicated formula has been chosen ( who chose it? ). ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
FABS.X FP3 ; fp3 = |t|, fp2 = s(s(P/Q)), fp0 = s
; FSUB.X #"0.5",FP3 ; fp3 = |t| - 0.5 = z DELETED <5/23/92, JPO>
FSUB.S #"$3F000000",FP3 ; FP3 <- z = |t| - 0.5 <5/23/92, JPO>
FMOVE.X FP3,FP1 ; duplicate z for the next step
; FDIV.W #2,FP1 ; fp1 = z/2 DELETED <5/23/92, JPO>
FMUL.S #"$3F000000",FP1 ; FP1 <- 0.5 * z <5/23/92, JPO>
FMUL.X FP1,FP3 ; fp3 = z*(z/2)
FSUB.X FP2,FP3 ; fp3 = z*(z/2) - s(s(P/Q))
FADD.X FP1,FP3 ; fp3 = z/2 + ( z*(z/2) - s(s(P/Q)))
; FMUL.W #-1,FP3 ; fp3 = -( z/2 + ( z*(z/2) - s(s(P/Q)))) DELETED <5/23/92, JPO>
FNEG.X FP3,FP3 ; FP3 <- -(z/2 + (z*(z/2)-s*(s*(P/Q)))) <5/23/92, JPO>
; FADD.X #"0.875",FP3 ; the result is now in hand DELETED <5/23/92, JPO>
FADD.S #"$3F600000",FP3 ; FP3 <- 0.875 + (FP3) <5/23/92, JPO>
FMOVE.X FP3,FP0 ; put in the right place
RTS ; return to UseCos caller
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; s < 1/4 then the tiny formula will be used. ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
SmallCos:
; FDIV.W #2,FP0 ; fp0 = s/2 DELETED <5/23/92, JPO>
FMUL.S #"$3F000000",FP0 ; FP0 <- 0.5 * s <5/23/92, JPO>
FSUB.X FP0,FP2 ; fp2 = -s/2 + s * ( s * P(s)/Q(s) )
; FADD.W #1,FP2 ; fp2 = result obtained by formula one DELETED <5/23/92, JPO>
FADD.S #"$3F800000",FP2 ; FP2 <- result by first formula <5/23/92, JPO>
FMOVE.X FP2,FP0 ; in the right place
RTS ; IMPORTANT, return to UseCos caller
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Subroutine UseSin, called by sine and cosine routines. ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Evaluate sin(t) for reduced |t| <= ¹/4. ;;
;; ;;
;; Use the following approximation: ;;
;; sin(t) = t - ( ( t * s ) * ( P(s) / Q(s) ) ), where s = t*t ;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
UseSin ; start of the sine computation
FMOVE.X FP0,FP3 ; fp0 = t, fp1 = Pi/2, fp3 = t
FMUL.X FP0,FP0 ; fp0 = s = t*t, fp1 = Pi/2, fp3 = t
LEA SINP,A0 ; A0 = coeff table ptr for PolyEval
BSR PolyEval ; fp1 = PolyEval(P, t*t)
FMOVE.X FP1,FP2 ; save P(arg)
LEA SINQ,A0 ; A0 = coeff table ptr for PolyEval
BSR PolyEval ; fp1 = PolyEval(Q, t*t)
FDIV.X FP1,FP2 ; form fp2 = P(s)/Q(s)
FMUL.X FP3,FP0 ; fp0 = t^3 = t * s, fp1 = junk, fp3 = t
FMUL.X FP0,FP2 ; form fp2 = ( t * s ) * ( P(s)/Q(s) )
FSUB.X FP2,FP3 ; fp3 = result, transfer it to fp0 & go
FMOVE.X FP3,FP0 ; in the right place
RTS ; return to UseSin Caller