|  is saved (the argument is forced non-negative).  Computations are|  performed for values quadrant values of -0.5 to 0.5* other values|  are transformed as per the following transformation formulae|  and table:||     Let:     z = shifted quadrant value|     Recall:  sin(-x) = -sin(x)    tan(-x) = -tan(x)|              cos(-x) =  cos(x)||     l========l=================l=================l=================l|     l  QUAD  l       SIN       l       COS       l       TAN       l|     l========l=================l=================l=================l|     l    0   l      sin(z)     l      cos(z)     l      tan(z)     l|     l    1   l      cos(z)     l     -sin(z)     l   -1/tan(z)     l|     l    2   l     -sin(z)     l     -cos(z)     l      tan(z)     l|     l    3   l     -cos(z)     l      sin(z)     l   -1/tan(z)     l|     l========l=================l=================l=================l||  The sign of the argument is, by the equations above, important only|  when computing SIN and TAN.||  Chebyshev Polynomials are used to approximate the trigonometric|  value in the range -0.5 to 0.5 quandrants.  (Note: constants are|  adjusted to reflect calculations done in quadrant mode instead|  of radian mode).||||         C16  =   6.5659 63114 97947 23622E-11*/DCOSCN:        .long   0x3DD20C62,0xC2F2D7F5|         C14  =  -6.3866 03083 79185 22411E-09        .long   0xBE3B6E24,0xF44B128F|         C12  =   4.7108 74778 81817 15037E-07        .long   0x3E9F9D38,0xA3763CC3|         C10  =  -2.5202 04237 30606 05481E-05        .long   0xBEFA6D1F,0x2A204A8C|         C08  =   9.1926 02748 39426 58024E-04        .long   0x3F4E1F50,0x6891BABB|         C06  =  -2.0863 48076 33529 96087E-02        .long   0xBF955D3C,0x7E3CBFFA|         C04  =   2.5366 95079 01048 01364E-01        .long   0x3FD03C1F,0x081B5AC4|         C02  =  -1.2337 00550 13616 98274E+00        .long   0xBFF3BD3C,0xC9BE45DE|         C00  =   1.0000 00000 00000 00000E+00        .long   0x3FF00000,0x00000000||       .set    NDCOSC,9||||         C15  =  -6.6880 35109 81146 72325E-10DSINCN:        .long   0xBE06FADB,0x9F155744|         C13  =   5.6921 72921 96792 68118E-08        .long   0x3E6E8F43,0x4D018D63|         C11  =  -3.5988 43235 21208 53405E-06        .long   0xBECE3074,0xFDE8871F|         C09  =   1.6044 11847 87359 82187E-04        .long   0x3F250783,0x487EE782|         C07 =   -4.6817 54135 31868 81007E-03        .long   0xBF732D2C,0xCE62BD86|         C05  =   7.9692 62624 61670 45121E-02        .long   0x3FB466BC,0x6775AAE2|         C03  =  -6.4596 40975 06246 25366E-01        .long   0xBFE4ABBC,0xE625BE53|         C01  =   1.5707 96326 79489 66192E+00DPIO2:        .long   0x3FF921FB,0x54442D18||       .set    NDSINC,8||||         Q3  =    1.7732 32244 08118 84863E+01DTNQCN:        .long   0x4031BB79,0x7BC569F8|         Q2  =   -8.0485 82645 07638 77427E+02        .long   0xC08926DD,0xB9C83D02|         Q1  =    6.4501 55566 23337 83845E+03        .long   0x40B93227,0xD3304CB9|         Q0  =   -5.6630 29630 56808 22115E+03        .long   0xC0B61F07,0x95DE70E0||       .set    NDTNQC,4|||         P3  =    1.0000 00000 00000 00000E+00DTNPCN:        .long   0x3FF00000,0x00000000|         P2  =   -1.5195 78275 66504 79235E+02        .long   0xC062FEA6,0x85FF2B0D|         P1  =    2.8156 53021 77302 44048E+03        .long   0x40A5FF4E,0x58DEAD6E|         P0  =   -8.8954 66142 22700 41047E+03        .long   0xC0C15FBB,0xAA8C6A22||       .set    NDTNPC,4|||         INV2PI  =  1 / 2*PI  =   1.5915 49430 91895 33577E-01DIN2PI:        .long   0x3FC45F30,0x6DC9C883||| ###   PUBLIC  DPCOS|DPCOS:        bsr     DTRGSV          | Prepare argument for operation|        moveb   d0,d1           | Copy into d1        andib   #3,d0           | Strip sign bit        addqb   #1,d0        rorb    #2,d0           | Move sign bit to B7        moveb   d0,a7@(8)       | Save it        rorb    #1,d1        jcs     DSINOP          | Compute sine||DCOSOP:        pea     DCOSCN        movel   sp@+,a6        moveq   #NDCOSC,d0        bsr     DX2SER|DTRGFN:        moveb   a7@(8),d0        andib   #0x80,d0                | Isolate sign bit        eorb    d0,sp@|        movel   a7@(4),a7@(8)        movel   sp@+,sp@        jmp     a4@|||| ###   PUBLIC  DPSIN|DPSIN:        bsr     DTRGSV          | Prepare argument|        addib   #0x7E,a7@(8)    | Compute sign        rorb    #1,d0        jcs     DCOSOP||DSINOP:        movel   a7@(4),sp@-     | Duplicate argument        movel   a7@(4),sp@-|        pea     DSINCN        movel   sp@+,a6        moveq   #NDSINC,d0        bsr     DX2SER        jsr     DPMUL|        jra     DTRGFN|||| ###   PUBLIC  DPTAN|DPTAN:        bsr     DTRGSV          | Prepare argument|        rorb    #1,d0        andib   #0x80,d0        eorb    d0,a7@(8)|        moveq   #4-1,d1         | Double duplication of argumentDPTN01:        movel   a7@(4),sp@-        dbra    d1,DPTN01|        pea     DTNPCN        movel   sp@+,a6        moveq   #NDTNPC,d0        bsr     DX2SER        jsr     DPMUL|        movel   a7@(8),d0        movel   a7@(12),d1        movel   sp@,a7@(8)        movel   a7@(4),a7@(12)        movel   d0,sp@        movel   d1,a7@(4)|        pea     DTNQCN        movel   sp@+,a6        moveq   #NDTNQC,d0        bsr     DX2SER|        btst    #0,a7@(16)        jne     DPTN20          | J/ reverse division|        jsr     DPDIV        jra     DTRGFN|DPTN20:        jsr     DPRDIV        jra     DTRGFN/*||       page||  Trigonometric Service Routine||  The double precision floating point value to be processed is on the|  stack.  Excess 2 PI's are removed from the argument.  The range of|  the argument is then scaled to within PI/4 of 0.  The shift size|  (in number of quadrants) is returned in D0.  The original sign bit|  is returned in the D0.B sign bit.||  If the argument is NaN, +INF, -INF, or too large (>= 2^16), this|  routine causes a return to the caller with a NaN result.|*/DTRGSV:        moveal  sp@+,a6 | DTRGSV return addr|        bsr     DPFADJ        jcs     DFNANR          | Return NAN if +/- INF        jvs     DFNANR          | Return NaN if NaN|        clrw    a7@(8)          | Create arg type return byte|        roxlw   sp@             | Sign bit into X        roxrw   a7@(8)          | Sign bit into temp byte        roxrw   sp@             | Restore DP value (with sign = +)|        cmpiw   #16*DBIAS+256,sp@        jge     DFNANR          | J/ return NaN if ABS(arg) >= 2^16|        movel   DIN2PI+4,sp@-   | Scale arg by 1/(2*PI)        movel   DIN2PI+0,sp@-        jsr     DPMUL|        cmpiw   #16*DBIAS,sp@        jlt     DTS10           | /* J/ no excess two PI's */|        movel   a7@(4),sp@-     | Double the scaled argument        movel   a7@(4),sp@-        jsr     DINT        negl    d1        negxl   d0        jsr     DFLOAT        jsr     DPADD|DTS10:        tstw    sp@        jeq     DTS20|        addiw   #16*3,sp@       | Multiply value by 8        cmpiw   #16*DBIAS,sp@        jlt     DTS20           | J/ < 1|        movel   a7@(4),sp@-     | Double parameter        movel   a7@(4),sp@-        jsr     DINT        addqb   #1,d1        lsrb    #1,d1        addb    d1,a7@(8)       | Mix quadrant number with sign bit|        lslb    #1,d1        negl    d1        negxl   d0        jsr     DFLOAT        jsr     DPADD|DTS20:        tstw    sp@        jeq     DTS30           | J/ arg = 0.0|        subiw   #16*1,sp@       | Range result to -0.5 to 0.5|DTS30:        andib   #0x83,a7@(8)    | Limit quadrant shift to 0 thru 3        moveb   a7@(8),d0       | Return cntl byte in d0 and at a7@(8)        jmp     a6@/*|| ###   SUBTTL  Polynomial Evaluation Routine|       page||  Polynomial Evaluation Routines.||  A list of constants is pointed to by A6.  X is on the stack.  D0|  contains the number of constants (the polynomial degree plus one).||  A6 enters pointing to C[1].  Upon return, the value on stack has|  been replaced with:|  DXSER   = C[1]*X^(N-1)  + C[2]*X^(N-2)  + ... + C[N-1]*X   + C[N]||  DX2SER  = C[1]*X^(2N-2) + C[2]*X^(2N-4) + ... + C[N-1]*X^2 + C[N]||*/DX2SER:        moveal  a4,a5           | Save a4 in a5        bsr     DPFADJ|        moveb   d0,a7@(8)       | Save count|        movel   a7@(4),sp@-     | Square the argument        movel   a7@(4),sp@-        jsr     DPMUL        jra     DXSR01          | J/ join DXSER routine|DXSER:        moveal  a4,a5        bsr     DPFADJ|        moveb   d0,a7@(8)|DXSR01:        movel   a6@(4),sp@-     | Place C[1] on stack as accum        movel   a6@,sp@-        addql   #8,a6|        subqb   #1,a7@(16)|DXSR02:        movel   a7@(12),sp@-    | Copy X        movel   a7@(12),sp@-        jsr     DPMUL        movel   a6@(4),sp@-     | Add in next coefficient        movel   a6@,sp@-        addql   #8,a6           | Advance to next coefficient        jsr     DPADD|        subqb   #1,a7@(16)        jne     DXSR02          | J/ more coefficients|        movel   a7@(4),a7@(16)  | Adjust return        movel   sp@+,a7@(8)        addql   #8,sp|        moveal  a4,a3        moveal  a5,a4           | Restore return address        jmp     a3@/*|| ###   SUBTTL  Function Completion Segments|       page||  DPFADJ - Adjust the stack for double precision functions.|          Move the return address to the function caller into A4.|          Shift the double precision floating point argument down|          four bytes, leaving a temporary area of four bytes at|          8(A7) thru 11(A7).  Set the condition code bits N, Z,|          C (INF), and V (NaN) to quickly type the argument.|*/DPFADJ:        moveal  sp@+,a3 | Function return addr        moveal  sp@,a4          | Caller return addr        movel   a7@(4),sp@      | Down shift parameter        movel   a7@(8),a7@(4)|        movew   sp@,d2          | Get SEM word        rolw    #1,d2        cmpiw   #32*0x7FF,d2        jcs     DPFA10          | J/ not INF or NaN|        andiw   #0x1E,d2        jeq     DPFA05|        orib    #0x02,ccr       | Set -V- bit -> NaN argument        jmp     a3@|DPFA05:        tstb    sp@             | /* Set N bit as req'd, clear Z bit. */        orib    #0x01,ccr       | Set -C- bit -> INF argument        jmp     a3@|DPFA10:        tstw    sp@             | /* Set N, Z bits as req'd, clear V, C */        jmp     a3@|||DFPINR:        |dsw    0               | Result is positive overflow        moveq   #0,d0           | Sign is positive        jra     DFIN01|DFMINR:        |dsw    0               | Result in negative overflow        moveq   #-1,d0          | Sign is negative|DFIN01:        moveaw  d0,a2           | Sign into a2        addql   #8,sp           | Delete twelve bytes from the stack        addql   #4,sp        moveal  a4,a0           | Return addr in a0        jmp     DINFRS          | (Use DPOPNS exception return)||DFNANR:        |dsw    0               | Result is NaN        addql   #8,sp           | Delete twelve bytes from the stack        addql   #4,sp        moveal  a4,a0        jmp     DNANRS||DFZERR:        |dsw    0               | Result is zero        addql   #8,sp        addql   #4,sp        moveal  a4,a0        jmp     DZERRS||DFUNFR:        |dsw    0               | Underflow result        addql   #8,sp           | Delete twelve bytes from the stack        addql   #4,sp        moveal  a4,a0           | Return address into a0        jmp     DUNFRS||DFONER:        |dsw    0               | Result is 1.0        addql   #4,sp        movel   #0x3FF00000,sp@        clrl    a7@(4)        clrb    FPERR        jmp     a4@||       end