sys7.1-doc-wip/Toolbox/InSANE/HWElemsSinCos.a
2019-07-27 22:37:48 +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