.TITLE IQSIN, RGC06E72 .IDENT '750214' .GLOBL IQSIN,IQSIN. ;ENTRIES ;DECUS NO: 11-191 ;PURPOSE: COMPUTE 16-BIT SIN(PI/2*I/J) ; FOR ABS(I)<=J. ;SYSTEM: PDP-11/40, EIS, DOS, MACRO ;USAGE: THE FORTRAN-COMPATIBLE FUNCTION CALL IS ; IX=IQSIN(I,J), USING THE R5 CALLING SEQUENCE. ; THE ARGUMENTS ARE 16-BIT INTEGERS, WITH J>0 AND ; ABS(I)<=J. THE INTEGER VALUE RETURNED IN R0, ; AND TO IX ABOVE, IS 2**15*SIN(PI/2*I/J) WITH ; MAGNITUDE NOT EXCEEDING 2**15-1. ;METHOD: COMPUTATION IS FIXED-POINT WITH 16-BIT ; PRECISION. A FOUR-TERM CHEBYSHEV APPROXIMATION ; IS SUFFICIENT FOR AN ABSOLUTE ERROR LESS THAN ; 2**(-15). THOSE FOUR TERMS REPRESENT A SEVENTH- ; ORDER POLYNOMIAL APPROXIMATION: SIN(X*PI/2) = ; (((C3*X**2+C2)*X**2+C1)*X**2+C0)*X. EACH COEF ; IS OF MAGNITUDE LESS THAN UNITY, WITH THE ; EXCEPTION OF C0=1.57079. WE SUBTRACT 1 FROM C0 ; AND ADD X AFTER EVALUATING THE REMAINDER OF THE ; EXPRESSION. WITH X**2<=1, EACH PARTIAL RESULT IN ; THE EVALUATION FROM LEFT TO RIGHT IS LESS THAN ; UNITY IN MAGNITUDE. ;REFERENCE: ABRAMOWITZ, HANDBOOK OF MATHEMATICAL ; FUNCTIONS, NBS-AMS55, 1964, PG.76 (4.3.104) ;SIZE: 65 WORDS. ;ACCURACY: FOR J A POWER OF 2, THE MAX ABSOLUTE ; ERROR IS 1.60 LSB (LEAST SIGNIFICANT BIT, OR PARTS ; IN 32768). FOR J NOT A POWER OF 2, THE MAX ; ABSOLUTE ERROR COULD NOT EXCEED 2.4 LSB (BUT ; PROBABLY DOES NOT EXCEED 1.74 LSB). WITH THIS ; ROUTINE MODIFIED FOR J FIXED AT 32768, AND WITH ; I RUNNING FROM 1 TO 32767, WE FOUND A MAX ABSOLUTE ; ERROR OF 1.60 LSB AND AN RMS ABSOLUTE ERROR LESS ; THAN .45 LSB. WE ROUND X=I/J, SO WE HAVE AN UPPER ; BOUND OF 2.4 LSB (=1.6+PI/4) FOR THE MAX ABSOLUTE ; ERROR FOR ALL I/J. HOWEVER, IN OVER 128,000 TEST ; CASES WE FOUND THE MAX ABSOLUTE ERROR TO BE 1.74 ; LSB AND THE RMS ABSOLUTE ERROR TO BE LESS THAN ; .6 LSB. ALL TEST COMPARISONS WERE BETWEEN ; FLOAT(IQSIN(I,J)) AND 32768.*SIN(PI/2*I/J). ; RELATIVE ERRORS WERE ALSO CALCULATED BUT WERE ; NOT USEFUL SINCE THE ERROR DISTRIBUTION WAS FOUND ; TO BE UNIFORM. SIN(-X) IS CALCULATED AS -SIN(X) ; TO GUARANTEE SYMMETRY TO THE LEAST SIGNIFICANT ; BIT. THE CHEBYSHEV COEFFICIENTS WERE PERTURBED ; DYNAMICALLY (WITH OVER 2 MILLION COMPARISONS) TO ; MINIMIZE THE TRUNCATION ERROR. THE ROUNDING OF ; X AND X**2 (VS. TRUNCATION OR VARIOUS CEIL/FLOOR ; COMBINATIONS) WAS SIMILARLY OPTIMIZED. ;TIMING: 185 MICROSECONDS AVERAGE. ;REMARKS: ALL REGISTERS ARE SAVED, WITH THE EXCEPTION ; OF R0, IN WHICH THE RESULT IS RETURNED. THE GLOBAL ; 'IQSIN.' IS AN ENTRY USED BY HIGHER LEVEL ROUTINES ; WHICH SAVE REGISTERS 1-5, DEFINE I AND J IN R1 ; AND R2, AND JUMP TO IQSIN.. ;PROGRAMMER: NELSON (MAGNET LAB, MIT) IQSIN: MOV R1,-(6) ;SAVE REGS 1-5 MOV R2,-(6) MOV R3,-(6) MOV R4,-(6) MOV R5,-(6) MOV @2(5),R1 ;R1=I MOV @4(5),R2 ;R2=J ;THE FOLLOWING ENTRY IS USED BY HIGHER LEVEL ROUTINES ; THAT SAVE REGISTERS 1-5, DEFINE I&J IN R1&R2, ; AND JUMP TO IQSIN.. IQSIN.: MOV R1,-(6) ;SAVE I ON STACK BNE ABS ;IF I=0, SIN=0 CLR R0 BR POP ABS: BPL .+4 ;R1=ABS(I) NEG R1 CMP R1,R2 ;IF ABS(I)=J, BGE OVR ; ABS(SIN)=2**15-1 ;DEFINE X=ABS(I)/J AND X**2 AS 16-BIT FIXED-POINT ; VALUES WITH BINARY POINT TO THE LEFT OF THE HIGH ; ORDER BIT. NEITHER ROUNDING OF X=I/J NOR ROUNDING ; OF X**2 CAN OVERFLOW. MOV R1,R4 ;X=ABS(I)/J CLR R5 ASHC #-1,R4 DIV R2,R4 ASR R2 ;(ROUNDED) CMP R2,R5 BGE .+4 INC R4 MOV R4,R0 ;R0=X MUL R4,R4 ASHC #1,R4 MOV R4,R2 ;R2=X**2 TST R5 ;(ROUNDED) BPL .+4 INC R2 ;COMPUTE SEVENTH-ORDER POLYNOMIAL APPROXIMATION TO ; SIN(X*PI/2). MOV #COEFS,R1 ;R1 -> COEFS MOV (1)+,R4 ;FIRST COEF MOV #3,R3 ;R3=COUNT MAD: MUL R2,R4 ;MULTIPLY BY X**2 ASHC #1,R4 ADD (1)+,R4 ;ADD NEXT COEF DEC R3 BGT MAD ;REPEAT 3 TIMES MUL R0,R4 ;MULTIPLY BY X ASHC #1,R4 ADD R4,R0 ;ADD X ;CHECK FOR OVERFLOW, AND CORRECT THE SIGN. BVC POP ;IF OVERFLOW, OVR: MOV #77777,R0 ; SET TO 2**15-1 POP: TST (6)+ ;NEGATE SIN IF I<0 BPL RETURN NEG R0 RETURN: MOV (6)+,R5 ;RESTORE REGS 5-1 MOV (6)+,R4 MOV (6)+,R3 MOV (6)+,R2 MOV (6)+,R1 RTS R5 ;RETURN FROM IQSIN ;CHEBYSHEV POLYNOMIAL COEFFICIENTS: COEFS: -143. ;-.00437 0784 = C3 2605. ; .07950 0304 = C2 -21165. ;-.64592 5832 = C1 18705. ; .57079 5134 = C0-1 .END