home *** CD-ROM | disk | FTP | other *** search
-
- ; *******************************************************
- ; * *
- ; * Turbo Pascal Runtime Library Version 7.0 *
- ; * Real Sine/Cosine *
- ; * *
- ; * Copyright (C) 1989-1993 Norbert Juffa *
- ; * *
- ; *******************************************************
-
- TITLE FPSIN
-
-
- CODE SEGMENT BYTE PUBLIC
-
- ASSUME CS:CODE
-
- ; Externals
-
- EXTRN RealAdd:NEAR,RealSub:NEAR,RealPoly:NEAR,RealMulNoChk:NEAR,
- EXTRN HaltError:NEAR,RealTrunc:NEAR,RealMul:NEAR,RealFloat:NEAR
- EXTRN RealReduceMP:NEAR,RFrac:FAR,RealMulFNoChk:NEAR
-
- ; Publics
-
- PUBLIC RSin,RCos
-
- ;-------------------------------------------------------------------------------
- ; RSin and RCos are different entries into the same routine that computes the
- ; sine as well as the cosine, as cosine is the same as the sine with an offset
- ; of pi/2. The argument is reduced to the range -pi/4..+pi/4. Depending on the
- ; octant into which the argument falls and the function to be computed, the
- ; cosine or sine of the reduced argument is computed by polynomial approxi-
- ; mations to the sine/cosine. To prevent excessive loss of precision in the
- ; argument reduction, the mapping of the argument to the reduced argument is
- ; done via a pseudo extended precision calculation. However this extended
- ; precision calculation is only available for arguments inside the range
- ; -1.68e+9..1.68e+9. Outside this range, a simple argument reduction is used.
- ; This causes a sudden drop in accuracy if the absolute value of the argument
- ; gets larger than 1.68e+9. For arguments bigger than 6.9e12, there is a total
- ; loss of precision. Outside this range, Sin always returns 0, Cos returns 1.
- ; The approximating polynomial for the sine on -pi/4..pi/4 is:
- ;
- ; P(x) := ((((2.476053850842e-8 *x² + 2.755533965768e-6)*x² -
- ; 1.984126372865e-4)*x² + 8.333333325083e-3)*x² -
- ; 1.666666666663e-1)*x²*x + x
- ;
- ; This approximation has a theoretical maximum relative error of 1.13e-14.
- ; Maximum observed error when evaluated in REAL arithmetic is 1.04e-13.
- ;
- ; The approximating polynomial for the cosine on -pi/4..pi/4 is:
- ;
- ; P(x) := ((((-2.715314622037e-7 *x² + 2.479862035332e-5)*x² -
- ; 1.388887879411e-3)*x² + 4.166666651281e-2)*x² -
- ; 4.999999999923e-1)*x² + 1.000000000000e+0
- ;
- ; This approximation has a theoretical maximum relative error of 1.41e-13.
- ; Maximum observed error when evaluated in REAL arithmetic is 1.49e-12.
- ;
- ; INPUT: DX:BX:AX argument
- ;
- ; OUTPUT: DX:BX:AX sin of argument
- ;
- ; DESTROYS: AX,BX,CX,DX,SI,DI,Flags
- ;-------------------------------------------------------------------------------
-
- RCos PROC FAR
- PUSH BP ; save caller's frame pointer
- AND DH, 7Fh ; abs(x), since cos (x) = cos (abs(x))
- MOV BP, 2 ; load function code for cosine
- JMP $cos_entry ; goto sine/cosine computation
- RCos ENDP
-
- ALIGN 4
-
- RSin PROC FAR
- PUSH BP ; save caller's frame pointer
- XOR BP, BP ; load function code for sine
- $cos_entry: MOV CH, 80h ; load mask to extract sign bit
- AND CH, DH ; extract sign bit
- XOR DH, CH ; x := abs(x)
- PUSH CX ; save sign
- $value_ok: MOV CX, 04E81h ; load
- MOV SI, 0836Eh ; real constant
- MOV DI, 022F9h ; 4/pi
- CMP AL, 9Fh ; x >= 2^30 ?
- JAE $big_val ; x/(pi/4) too big to be stored in long
- PUSH AX ; save x
- PUSH BX ; on
- PUSH DX ; stack
- CALL RealMulFNoChk ; compute x/(pi/4)
- XOR CH, CH ; signal truncation desired
- CALL RealTrunc ; n = Trunc (x/(pi/4))
- ROR AL, 1 ; if quadrant odd, set CF=1
- ROL AL, 1 ; restore register, CF=1 if odd quadrant
- ADC AX, 0 ; increment quadrant
- ADC DX, 0 ; if quadrant odd
- ADD BP, AX ; q := quadrant + flag
- CALL RealFloat ; convert quadrant to REAL value
- POP DI ; get
- POP SI ; back
- POP CX ; x
- TEST BP, 3 ; determine octants that use cos,sin
- PUSHF ; save sin/cos flag
- PUSH BP ; save q
- MOV BP, OFFSET c1 ; pointer to reduction constant table
- CALL RealReduceMP ; compute x + n*(pi/2) to full precision
- MOV CX, 5 ; use five coefficients
- XOR SI, SI ; signal P(x²)*x+x wanted
- POP BP ; get back q
- POPF ; get back sin/cos flag
- JPE $sine ; octants 0, 3, 4, 7 take sine
- $cosine: MOV DI,OFFSET CS_COEFF; pointer to 1st coeff. of cos approx.
- INC SI ; signal P(x²) wanted
- CALL RealPoly ; P(x²) in DX:BX:AX
- MOV CX, 00081h ; load
- XOR SI, SI ; floating-point
- MOV DI, SI ; constant 1.0
- CALL RealAdd ; 1 + P(x²)
- JMP $make_sign ; put appropriate sign on result
- $sine: MOV DI,OFFSET SN_COEFF; pointer to 1st coeff. of sin approx.
- CALL RealPoly ; P(x²)*x+x in DX:BX:AX
- $make_sign: POP CX ; get sign mask
- INC BP ; determine if
- TEST BP, 4 ; result has to negated
- JZ $ende ; result ok
- XOR CH, 80h ; negate sign flag
- $ende: POP BP ; restore caller's frame pointer
- OR AL, AL ; result zero ?
- JZ $end_sin ; yes, done
- XOR DH, CH ; make sign bit
- $end_sin: RET ; done
-
- $big_val: SUB CL, 3 ; 1/(2*pi) in DI:SI:CX
- CALL RealMulNoChk ; x/(2*pi)
- CALL RFrac ; frac (x/(2*pi))
- MOV CX, 02183h ; load
- MOV SI, 0DAA2h ; constant
- MOV DI, 0490Fh ; 2*pi
- CALL RealMulNoChk ; 2*pi * frac (x/(2*pi)) = remainder
- JMP $value_ok ; reduction to 0..2*pi done
-
- SN_COEFF DB 067h, 0B1h,0D4h ; -2.476053850842e-8
- DB 06Eh,05Ch,08Bh,0B6h,0EBh,038h ; 2.755533965768e-6
- DB 074h,0B5h,09Ch,0FCh,00Ch,0D0h ; -1.984126372865e-4
- DB 07Ah,044h,086h,088h,088h,008h ; 8.333333325083e-3
- DB 07Eh,0A9h,0AAh,0AAh,0AAh,0AAh ; -1.666666666663e-1
- CS_COEFF DB 06Bh, 0C7h,091h ; -2.715314622037e-7
- DB 071h,034h,0B7h,0A1h,006h,050h ; 2.479862035332e-5
- DB 077h,02Eh,00Ah,058h,00Bh,0B6h ; -1.388887879411e-3
- DB 07Ch,018h,0A0h,0AAh,0AAh,02Ah ; 4.166666651281e-2
- DB 07Fh,0EFh,0FFh,0FFh,0FFh,0FFh ; -4.999999999923e-1
-
- C1 DB 080h,000h,000h,000h,00Fh,0C9h ; pi/4-eps =
- C2 DB 070h,0c2h,068h,021h,0a2h,0dah ; eps =
-
- RSIN ENDP
-
- ALIGN 4
-
- CODE ENDS
-
- END
