home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / vol_100 / 156_01 / transcen.c < prev    next >
Text File  |  1985-08-21  |  9KB  |  605 lines
  1. #include iolib.h
  2. #include float.h
  3. #include transcen.h
  4. #asm
  5. ;
  6. ;    transcendental floating point routines
  7. ;
  8. ;    History...
  9. ;        26 Jun 83    Added 'pow'.
  10. ;        25 Jun83    Added 'sinh', 'cosh',
  11. ;            and 'tanh'.
  12. ;        8 Nov 82    Added 'sqrt'
  13. ;        7 Nov 82    Added 'atan2'
  14. ;        17 Oct 82    removed constant pi.
  15. ;        11 Oct 82    atn => atan
  16. ;            EXP constant fixed, several labels
  17. ;            changed, references to 'int' changed
  18. ;            to 'qfloor', 'qpi' inserted.
  19. ;        23 Sept 82    Added colons after labels,
  20. ;            declaring all coefficients with DBs.
  21. ;        4 Sept 82    Separated from rest of
  22. ;            floating point routines.
  23. ;        8 Aug 82    AS.COM format source code.
  24. ;        31 Jul 82    Translated to Zilog mnemonics.
  25. ;        30 Jul 82    Disassembled from Xitan Disk
  26. ;            Basic.
  27. ;
  28. ONE:    DW    0
  29.     DW    0
  30.     DW    8100H
  31. LOGCOEF: DB    6
  32.     DB    23H,85H,0ACH,0C3H,11H,7FH
  33.     DB    53H,0CBH,9EH,0B7H,23H,7FH
  34.     DB    0CCH,0FEH,0A6H,0DH,53H,7FH
  35.     DB    0CBH,5CH,60H,0BBH,13H,80H
  36.     DB    0DDH,0E3H,4EH,38H,76H,80H
  37.     DB    5CH,29H,3BH,0AAH,38H,82H
  38. ;
  39. QLOG10:    CALL    QLOG    
  40.     LD    BC,7F5EH    ; 1/ln(10)
  41.     LD    IX,5BD8H
  42.     LD    DE,0A938H
  43.     JP    FMUL    
  44. ;
  45. QLOG:    CALL    SGN    
  46.     OR    A
  47.     JP    PE,ILLFCT
  48.     LD    HL,FA+5    
  49.     LD    A,(HL)
  50.     LD    BC,8035H    ; 1/sqrt(2)
  51.     LD    IX,04F3H
  52.     LD    DE,33FAH
  53.     SUB    B
  54.     PUSH    AF
  55.     LD    (HL),B
  56.     PUSH    DE
  57.     PUSH    IX
  58.     PUSH    BC
  59.     CALL    FADD    
  60.     POP    BC
  61.     POP    IX
  62.     POP    DE
  63.     INC    B
  64.     CALL    FDIV    
  65.     LD    HL,ONE    
  66.     CALL    HLSUB    
  67.     LD    HL,LOGCOEF
  68.     CALL    EVENPOL
  69.     LD    BC,8080H    ; -0.5
  70.     LD    IX,0
  71.     LD    DE,0
  72.     CALL    FADD    
  73.     POP    AF
  74.     CALL    L247E    
  75.     LD    BC,8031H    ; ln(2)
  76.     LD    IX,7217H
  77.     LD    DE,0F7D2H
  78.     JP    FMUL    
  79. ;
  80. L265F:    POP    BC
  81.     POP    IX
  82.     POP    DE
  83.     JP    FMUL
  84. ;
  85. ;
  86. QEXP:    LD    BC,8138H    ;1.442695041
  87.     LD    IX,0AA3BH
  88.     LD    DE,295CH    
  89.     CALL    FMUL    
  90.     LD    A,(FA+5)
  91.     CP    88H
  92.     JP    NC,DIV17
  93.     CALL    PUSHFA    
  94.     CALL    QFLOOR    
  95.     POP    BC
  96.     POP    IX
  97.     POP    DE
  98.     PUSH    AF
  99.     CALL    FSUB    
  100.     LD    HL,EXPCOEF
  101.     CALL    POLY    
  102.     LD    HL,FA+5    
  103.     POP    AF
  104.     OR    A
  105.     JP    M,EXP5    
  106.     ADD    A,(HL)
  107.     DB    1    ;"ignore next 2 bytes"
  108. EXP5:    ADD    A,(HL)
  109.     CCF    
  110.     LD    (HL),A
  111.     RET    NC
  112.     JP    DIV17    
  113. ;
  114. EXPCOEF: DB    10
  115.     DB    0CCH,0D5H,45H,56H,15H,6AH
  116.     DB    0CFH,37H,0A0H,92H,27H,6DH
  117.     DB    0F5H,95H,0EEH,93H,00H,71H
  118.     DB    0D0H,0FCH,0A7H,78H,21H,74H
  119.     DB    0B1H,21H,82H,0C4H,2EH,77H
  120.     DB    82H,58H,58H,95H,1DH,7AH
  121.     DB    6DH,0CBH,46H,58H,63H,7CH
  122.     DB    0E9H,0FBH,0EFH,0FDH,75H,7EH
  123.     DB    0D2H,0F7H,17H,72H,31H,80H
  124.     DB    0,0,0,0,0,81H
  125. EVENPOL: CALL    PUSHFA    
  126.     LD    DE,L265F
  127.     PUSH    DE
  128.     PUSH    HL
  129.     CALL    LDBCFA    
  130.     CALL    FMUL    
  131.     POP    HL
  132. POLY:    CALL    PUSHFA    
  133.     LD    A,(HL)
  134.     INC    HL
  135.     CALL    DLOAD    
  136.     DB    0FEH    ;"ignore next byte"
  137. POL3:    POP    AF
  138.     POP    BC
  139.     POP    IX
  140.     POP    DE
  141.     DEC    A
  142.     RET    Z
  143.     PUSH    DE
  144.     PUSH    IX
  145.     PUSH    BC
  146.     PUSH    AF
  147.     PUSH    HL
  148.     CALL    FMUL    
  149.     POP    HL
  150.     CALL    LDBCHL    
  151.     PUSH    HL
  152.     CALL    FADD    
  153.     POP    HL
  154.     JR    POL3    
  155. ;
  156. QCOS:    LD    HL,HALFPI
  157.     CALL    HLADD    
  158. QSIN:    CALL    PUSHFA    
  159.     LD    BC,08349H    ;6.283185308... = 2*pi
  160.     LD    IX,00FDAH
  161.     LD    DE,0A222H
  162.     CALL    LDFABC    
  163.     POP    BC
  164.     POP    IX
  165.     POP    DE
  166.     CALL    FDIV    
  167.     CALL    PUSHFA    
  168.     CALL    QFLOOR    
  169.     POP    BC
  170.     POP    IX
  171.     POP    DE
  172.     CALL    FSUB    
  173.     LD    HL,QUARTER
  174.     CALL    HLSUB    
  175.     CALL    SGN    
  176.     SCF    
  177.     JP    P,SIN5    
  178.     CALL    ADDHALF    
  179.     CALL    SGN    
  180.     OR    A
  181. SIN5:    PUSH    AF
  182.     CALL    P,MINUSFA
  183.     LD    HL,QUARTER
  184.     CALL    HLADD    
  185.     POP    AF
  186.     CALL    NC,MINUSFA
  187.     LD    HL,SINCOEF
  188.     JP    EVENPOL
  189. ;
  190. EE:    DW    058A4H    ;2.718281828... = e
  191.     DW    0F854H
  192.     DW    0822DH
  193. HALFPI:    DB    022H,0A2H,0DAH,00FH,049H,081H ; pi/2
  194. ;    DW    0A222H
  195. ;    DW    00FDAH
  196. ;    DW    08149H
  197. QUARTER: DW    0
  198.     DW    0
  199.     DW    7F00H
  200. SINCOEF: DB    7
  201.     DB    90H,0BAH,34H,76H,6AH,82H
  202.     DB    0E4H,0E9H,0E7H,4BH,0F1H,84H
  203.     DB    0B1H,4FH,7FH,3BH,28H,86H
  204.     DB    31H,0B6H,64H,69H,99H,87H
  205.     DB    0E4H,36H,0E3H,35H,23H,87H
  206.     DB    24H,31H,0E7H,5DH,0A5H,86H
  207.     DB    21H,0A2H,0DAH,0FH,49H,83H
  208. QTAN:    CALL    PUSHFA    
  209.     CALL    QSIN    
  210.     POP    BC
  211.     POP    IX
  212.     POP    DE
  213.     CALL    L2895    
  214.     EX    DE,HL
  215.     CALL    LDFABC    
  216.     CALL    QCOS    
  217.     JP    DIV1    
  218. ;
  219. #endasm
  220. pow(x,y)    /* x to the power y */
  221. double x,y;
  222. {    return exp(log(x)*y);
  223. }
  224. sinh(x) double x;
  225. {    double e;
  226.     e=exp(x);
  227.     return .5*(e-1./e);
  228. }
  229. cosh(x) double x;
  230. {    double e;
  231.     e=exp(x);
  232.     return .5*(e+1./e);
  233. }
  234. tanh(x) double x;
  235. {    double e;
  236.     e=exp(x);
  237.     return (e-1./e)/(e+1./e);
  238. }
  239. #asm
  240. ;
  241. QATAN:    CALL    SGN    
  242.     CALL    M,ODD        ;negate argument & answer
  243.     LD    A,(FA+5)
  244.     CP    81H
  245.     JR    C,ATAN5        ;c => argument less than 1
  246.     LD    BC,8100H    ;1.0
  247.     LD    IX,0
  248.     LD    D,C
  249.     LD    E,C
  250.     CALL    FDIV    
  251.     LD    HL,HLSUB
  252.     PUSH    HL    ;will subtract answer from pi/2
  253. ATAN5:    LD    HL,ATNCOEF
  254.     CALL    EVENPOL
  255.     LD    HL,HALFPI    ;may use for subtraction
  256.     RET    
  257. ;
  258. ATNCOEF: DB    13
  259.     DB    14H,7H,0BAH,0FEH,62H,75H
  260.     DB    51H,16H,0CEH,0D8H,0D6H,78H
  261.     DB    4CH,0BDH,7DH,0D1H,3EH,7AH
  262.     DB    1,0CBH,23H,0C4H,0D7H,7BH
  263.     DB    0DCH,3AH,0AH,17H,34H,7CH
  264.     DB    36H,0C1H,0A3H,81H,0F7H,7CH
  265.     DB    0EBH,16H,61H,0AEH,19H,7DH
  266.     DB    5DH,78H,8FH,60H,0B9H,7DH
  267.     DB    0A2H,44H,12H,72H,63H,7DH
  268.     DB    16H,62H,0FBH,47H,92H,7EH
  269.     DB    0C0H,0F0H,0BFH,0CCH,4CH,7EH
  270.     DB    7EH,8EH,0AAH,0AAH,0AAH,7FH
  271.     DB    0F6H,0FFH,0FFH,0FFH,07FH,80H
  272. ;/*    arc tangent of y/x  */
  273. ;atan2(y,x) double x,y;
  274. QATAN2:
  275. ;{    double a;
  276.     PUSH BC
  277.     PUSH BC
  278.     PUSH BC
  279. ;    if(fabs(x)>=fabs(y))
  280.     LD HL,8
  281.     ADD HL,SP
  282.     CALL DLOAD
  283.     CALL DPUSH
  284.     CALL QFABS
  285.     POP BC
  286.     POP BC
  287.     POP BC
  288.     PUSH HL
  289.     LD HL,16
  290.     ADD HL,SP
  291.     CALL DLOAD
  292.     CALL DPUSH
  293.     CALL QFABS
  294.     POP BC
  295.     POP BC
  296.     POP BC
  297.     CALL CCGE
  298.     LD A,H
  299.     OR L
  300.     JP Z,ATAN23
  301. ;        {a=atan(y/x);
  302.     LD HL,0
  303.     ADD HL,SP
  304.     PUSH HL
  305.     LD HL,16
  306.     ADD HL,SP
  307.     CALL DLOAD
  308.     CALL DPUSH
  309.     LD HL,16
  310.     ADD HL,SP
  311.     CALL DLOAD
  312.     CALL DDIV
  313.     CALL DPUSH
  314.     CALL QATAN
  315.     POP BC
  316.     POP BC
  317.     POP BC
  318.     POP HL
  319.     CALL DSTORE
  320. ;        if(x<0.)
  321.     LD HL,8
  322.     ADD HL,SP
  323.     CALL DLOAD
  324.     CALL DPUSH
  325.     LD HL,ATAN27+0
  326.     CALL DLOAD
  327.     CALL DLT
  328.     LD A,H
  329.     OR L
  330.     JR Z,ATAN22
  331. ;            {if(y>=0.) a=a+3.14159265358979;
  332.     LD HL,14
  333.     ADD HL,SP
  334.     CALL DLOAD
  335.     CALL DPUSH
  336.     LD HL,ATAN27+6
  337.     CALL DLOAD
  338.     CALL DGE
  339.     LD A,H
  340.     OR L
  341.     JR Z,ATAN20
  342.     LD HL,0
  343.     ADD HL,SP
  344.     PUSH HL
  345.     LD HL,2
  346.     ADD HL,SP
  347.     CALL DLOAD
  348.     CALL DPUSH
  349.     LD HL,ATAN27+12
  350.     CALL DLOAD
  351.     CALL DADD
  352.     POP HL
  353.     CALL DSTORE
  354. ;            else a=a-3.14159265358979;
  355.     JR ATAN21
  356. ATAN20:
  357.     LD HL,0
  358.     ADD HL,SP
  359.     PUSH HL
  360.     LD HL,2
  361.     ADD HL,SP
  362.     CALL DLOAD
  363.     CALL DPUSH
  364.     LD HL,ATAN27+18
  365.     CALL DLOAD
  366.     CALL DSUB
  367.     POP HL
  368.     CALL DSTORE
  369. ATAN21:
  370. ;            }
  371. ;        }
  372. ATAN22:
  373. ;    else
  374.     JR ATAN26
  375. ATAN23:
  376. ;        {a=-atan(x/y);
  377.     LD HL,0
  378.     ADD HL,SP
  379.     PUSH HL
  380.     LD HL,10
  381.     ADD HL,SP
  382.     CALL DLOAD
  383.     CALL DPUSH
  384.     LD HL,22
  385.     ADD HL,SP
  386.     CALL DLOAD
  387.     CALL DDIV
  388.     CALL DPUSH
  389.     CALL QATAN
  390.     POP BC
  391.     POP BC
  392.     POP BC
  393.     CALL MINUSFA
  394.     POP HL
  395.     CALL DSTORE
  396. ;        if(y<0.) a=a-1.57079632679489;  /* pi/2 */
  397.     LD HL,14
  398.     ADD HL,SP
  399.     CALL DLOAD
  400.     CALL DPUSH
  401.     LD HL,ATAN27+24
  402.     CALL DLOAD
  403.     CALL DLT
  404.     LD A,H
  405.     OR L
  406.     JR Z,ATAN24
  407.     LD HL,0
  408.     ADD HL,SP
  409.     PUSH HL
  410.     LD HL,2
  411.     ADD HL,SP
  412.     CALL DLOAD
  413.     CALL DPUSH
  414.     LD HL,ATAN27+30
  415.     CALL DLOAD
  416.     CALL DSUB
  417.     POP HL
  418.     CALL DSTORE
  419. ;        else     a=a+1.57079632679489;
  420.     JR ATAN25
  421. ATAN24:
  422.     LD HL,0
  423.     ADD HL,SP
  424.     PUSH HL
  425.     LD HL,2
  426.     ADD HL,SP
  427.     CALL DLOAD
  428.     CALL DPUSH
  429.     LD HL,ATAN27+36
  430.     CALL DLOAD
  431.     CALL DADD
  432.     POP HL
  433.     CALL DSTORE
  434. ATAN25:
  435. ;        }
  436. ATAN26:
  437. ;    return a;
  438.     LD HL,0
  439.     ADD HL,SP
  440.     CALL DLOAD
  441.     POP BC
  442.     POP BC
  443.     POP BC
  444.     RET
  445. ;}
  446. ATAN27:    DB 0,0,0,0,0,0
  447.     DB 0,0,0,0,0,0
  448.     DB 33,162,218,15,73,130
  449.     DB 33,162,218,15,73,130
  450.     DB 0,0,0,0,0,0
  451.     DB 34,162,218,15,73,129
  452.     DB 34,162,218,15,73,129
  453. #endasm
  454. #asm
  455. ;double extra; /* current approximate root */
  456. ;sqrt(x) double x;
  457. QSQRT:
  458. ;{    char *px,    /* points to x */
  459.     PUSH BC
  460. ;    *pextra;        /* points to extra */
  461.     PUSH BC
  462. ;    int i;        /* loop counter */
  463.     PUSH BC
  464. ;    if(x==0.)return 0.;
  465.     LD HL,8
  466.     ADD HL,SP
  467.     CALL DLOAD
  468.     CALL DPUSH
  469.     LD HL,SQRT10
  470.     CALL DLOAD
  471.     CALL DEQ
  472.     LD A,H
  473.     OR L
  474.     JR Z,SQRT2
  475.     LD HL,SQRT10
  476.     JP SQRT9
  477. ;    if(x<0.)
  478. SQRT2:
  479.     LD HL,8
  480.     ADD HL,SP
  481.     CALL DLOAD
  482.     CALL DPUSH
  483.     LD HL,SQRT10
  484.     CALL DLOAD
  485.     CALL DLT
  486.     LD A,H
  487.     OR L
  488.     JR Z,SQRT4
  489. ;        {err("ill sqrt"); return 0.;
  490.     LD HL,SQRT12
  491.     PUSH HL
  492.     CALL QERR
  493.     POP BC
  494.     LD HL,SQRT10
  495.     JP SQRT9
  496. ;        }
  497. ;    px=&x; pextra=&extra;    /* set the pointers */
  498. SQRT4:
  499.     LD HL,4
  500.     ADD HL,SP
  501.     PUSH HL
  502.     LD HL,10
  503.     ADD HL,SP
  504.     CALL CCPINT
  505.     LD HL,2
  506.     ADD HL,SP
  507.     PUSH HL
  508.     LD HL,EXTRA
  509.     CALL CCPINT
  510. ;    extra=.707;    /* initialize fraction at sqrt(.5) */
  511.     LD HL,SQRT14
  512.     CALL DLOAD
  513.     LD HL,EXTRA
  514.     CALL DSTORE
  515. ;    pextra[5]=(px[5]>>1)^64; /* answer exponent is half
  516. ;                of "x" exponent */
  517.     LD HL,2
  518.     ADD HL,SP
  519.     CALL CCGINT
  520.     PUSH HL
  521.     LD HL,5
  522.     POP DE
  523.     ADD HL,DE
  524.     PUSH HL
  525.     LD HL,6
  526.     ADD HL,SP
  527.     CALL CCGINT
  528.     PUSH HL
  529.     LD HL,5
  530.     POP DE
  531.     ADD HL,DE
  532.     CALL CCGCHAR
  533.     PUSH HL
  534.     LD HL,1
  535.     POP DE
  536.     CALL CCASR
  537.     PUSH HL
  538.     LD HL,64
  539.     CALL CCXOR
  540.     POP DE
  541.     LD A,L
  542.     LD (DE),A
  543. ;    i=5;    /* 5 iterations of Newton's algorithm */
  544.     LD HL,0
  545.     ADD HL,SP
  546.     PUSH HL
  547.     LD HL,5
  548.     CALL CCPINT
  549. ;    while(i--)
  550. SQRT6:
  551.     LD HL,0
  552.     ADD HL,SP
  553.     PUSH HL
  554.     CALL CCGINT
  555.     DEC HL
  556.     CALL CCPINT
  557.     INC HL
  558.     LD A,H
  559.     OR L
  560.     JR Z,SQRT8
  561. ;        {extra=(extra+x/extra);
  562.     LD HL,EXTRA
  563.     CALL DLOAD
  564.     CALL DPUSH
  565.     LD HL,14
  566.     ADD HL,SP
  567.     CALL DLOAD
  568.     CALL DPUSH
  569.     LD HL,EXTRA
  570.     CALL DLOAD
  571.     CALL DDIV
  572.     CALL DADD
  573.     LD HL,EXTRA
  574.     CALL DSTORE
  575. ;        pextra[5]--;    /* /2 */
  576.     LD HL,2
  577.     ADD HL,SP
  578.     CALL CCGINT
  579.     PUSH HL
  580.     LD HL,5
  581.     POP DE
  582.     ADD HL,DE
  583.     PUSH HL
  584.     CALL CCGCHAR
  585.     DEC HL
  586.     POP DE
  587.     LD A,L
  588.     LD (DE),A
  589.     INC HL
  590. ;        }
  591.     JR SQRT6
  592. SQRT8:
  593. ;    return extra;
  594.     LD HL,EXTRA
  595. SQRT9:    CALL DLOAD
  596.     POP BC
  597.     POP BC
  598.     POP BC
  599.     RET
  600. ;}
  601. SQRT10:    DB 0,0,0,0,0,0
  602. SQRT12:    DB 'ill sqrt',0
  603. SQRT14:    DB 70,182,243,253,52,128
  604. #endasm
  605.