home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / vol_100 / 156_01 / float.c < prev    next >
Text File  |  1985-08-21  |  22KB  |  1,200 lines
  1. #include iolib.h
  2. #include float.h
  3. #asm
  4. ;
  5. ;    name...
  6. ;        float
  7. ;
  8. ;    purpose...
  9. ;        floating point routines for C programs
  10. ;
  11. ;    note...
  12. ;        This code uses some of the Z-80's UNDOCUMENTED 
  13. ;    INSTRUCTIONS. The instructions work because the IX and
  14. ;    IY escape bytes are actually a bit more general than
  15. ;    Zilog advertises. Note that in the instructions
  16. ;        E1    POP HL
  17. ;        DD E1    POP IX
  18. ;        FD E1    POP IY
  19. ;    (which Zilog documents), the opcodes for the IX and IY
  20. ;    registers are the same as for the HL register, but are
  21. ;    preceded by the escape bytes DD and FD. Similarly,
  22. ;        6A    LD L,D
  23. ;        DD 6A    LD XL,D        undocumented
  24. ;        FD 6A    LD YL,D        undocumented
  25. ;    where the new instructions deal with the lower half
  26. ;    of the IX and IY registers, respectively. All of the
  27. ;    single-byte transfers, and single byte operations such
  28. ;    as ADD and INC, involving the H or L registers, can be
  29. ;    preceded by escape bytes in this fashion.
  30. ;
  31. ;    history...
  32. ;        The floating point arithmetic routines were
  33. ;        written by Neil Colvin and were included in
  34. ;        Xitan Disk Basic.  He has consented to release
  35. ;        them to the public domain for individual use
  36. ;        but not for sale. Disassembly, comments, ifix(),
  37. ;        float(), and the interface to c programs were by
  38. ;        James R. Van Zandt.
  39. ;
  40. ;        Aug 84        moved atof() and putf() to
  41. ;            PRINT library.
  42. ;        17 Jun 84    including float.h & iolib.h, so
  43. ;            output can be separately assembled.
  44. ;            Calling err() for errors, so walkback
  45. ;            trace can work.
  46. ;        26 Jun 83    Renamed: abs->fabs, mod->fmod.
  47. ;        25 Jun 83    Added 'amax', 'amin',
  48. ;            'putf', 'atof'.
  49. ;        6 Nov 82    Removed 'sqrt' declaration
  50. ;        20 Oct 82    Added 'odd', used it
  51. ;            in 'ceil'
  52. ;        9 Oct 82    Added 'dswap'.
  53. ;        6 Oct 82    Removed my 'floor', renamed
  54. ;            existing 'int' to 'floor'.
  55. ;        5 Oct 82    "int" -> "qint"
  56. ;            declaring int and mod.
  57. ;        29 Sept 82    Added 'floor' code.
  58. ;        23 Sept 82    Removed transcendental
  59. ;            routines.
  60. ;        26 Aug 82    Added comparison routines.
  61. ;        18 Aug 82    Added 'ifix'
  62. ;        15 Aug 82    Added 'float'.
  63. ;        8 Aug 82    AS.COM format source code.
  64. ;        31 Jul 82    Translated to Zilog mnemonics.
  65. ;        30 Jul 82    Disassembled from Xitan Disk
  66. ;            Basic.
  67. ;
  68. ;    double float(),    /* integer to floating point conversion */
  69. ;    mod(),        /* mod(x,y) */
  70. ;    fabs(),        /* absolute value */
  71. ;    floor(),    /* largest integer not greater than */
  72. ;    ceil(),        /* smallest integer not less than */
  73. ;    rand();        /* random number in range 0...1 */
  74. ;    int ifix();    /* floating point to integer
  75. ;            (takes floor first) */
  76. ;
  77. EXTRA:    DEFS    6
  78. FA:    DEFS    6    ;floating point accumulator
  79. FASIGN:    DEFS    1    ;msb indicates sign of FA
  80. ;            0 => negative, 1 => positive
  81.  
  82. ;
  83. L0F2E:    DEFB    0
  84. SEED:    DEFB    80H,80H,0,0,0,0 ;seed for random number generator
  85. ;
  86. DIVZERO: CALL    GRIPE
  87.     DEFB    'can''t /0',0
  88. ILLFCT:    CALL    GRIPE
  89.     DEFB    'Illegal function',0
  90. OFLOW:    CALL    GRIPE
  91.     DEFB    'Arithmetic overflow',0
  92. GRIPE:    CALL    QERR    ;top word on stack points to message
  93.     JP    0    ;error was fatal
  94. ;
  95. ;    push the floating point accumulator
  96. ;    (preserving return address)
  97. ;
  98. DPUSH:    POP    DE
  99.     LD    HL,(FA+4)
  100.     PUSH    HL
  101.     LD    HL,(FA+2)
  102.     PUSH    HL
  103.     LD    HL,(FA)
  104.     PUSH    HL
  105.     EX    DE,HL
  106.     JP    (HL)
  107. ;
  108. ;    push floating point accumulator
  109. ;    (preserve return address and next stacked word)
  110. ;
  111. DPUSH2:    POP    DE    ;save return address
  112.     POP    BC    ;save next word
  113.     LD    HL,(FA+4)
  114.     PUSH    HL
  115.     LD    HL,(FA+2)
  116.     PUSH    HL
  117.     LD    HL,(FA)
  118.     PUSH    HL
  119.     EX    DE,HL
  120.     PUSH    BC    ;restore next word
  121.     JP    (HL)    ;return
  122. ;
  123. ;    convert the integer in hl to
  124. ;    a floating point number in FA
  125. ;
  126. QFLOAT:    LD    A,H    ;fetch MSB
  127.     CPL        ;reverse sign bit
  128.     LD    (FASIGN),A ;save sign (msb)
  129.     RLA        ;move sign into cy
  130.     JR    C,FL4    ;c => nonnegative number
  131.     EX    DE,HL
  132.     SBC    HL,HL    ;clear hl
  133.     SBC    HL,DE    ;get positive number into hl
  134. FL4:    LD    A,L
  135.     DEFB    0DDH
  136.     LD    H,A    ;move LSB to hx
  137.     LD    C,H    ;move MSB to c
  138.     LD    DE,0    ;clear rest of registers
  139.     LD    B,D
  140.     DEFB    0DDH
  141.     LD    L,D    ;clear lx
  142.     LD    A,16+128
  143.     LD    (FA+5),A ;preset exponent
  144.     JP    NORM    ;go normalize c ix de b
  145. ;
  146. ;    convert the floating point number in FA
  147. ;    to an integer in hl  (rounds toward negative numbers)
  148. ;
  149. QIFIX:    CALL    QFLOOR        ;take floor first
  150.     LD    HL,0        ;initialize the result
  151.     LD    A,(FA+5)    ;get the exponent
  152.     LD    B,A        ;  and save it
  153.     OR    A
  154.     RET    Z        ;z => number was zero
  155.     LD    HL,(FA+3)    ;get most significant bits
  156.     LD    C,H        ;save sign bit (msb)
  157.     LD    A,B        ;get exponent again
  158.     CP    80H+16
  159.     JP    M,IFIX5        ;m => fabs(fa) < 32768
  160.     JR    NZ,OFLOW    ;nz => fabs(fa) > 32768
  161. ;                (overflow)
  162.     LD    A,H
  163.     CP    80H
  164.     JR    NZ,OFLOW    ;nz => fa isn't -32768
  165.     LD    A,L
  166.     OR    A
  167.     JR    NZ,OFLOW    ;nz => overflow
  168.     RET            ;return -32768.
  169. ;
  170. IFIX5:    SET    7,H        ;restore hidden bit
  171. IFIX6:    SRL    H        ;shift right (0 fill)
  172.     RR    L        ;shift right (cy fill)
  173.     INC    A
  174.     CP    16+80H
  175.     JR    NZ,IFIX6    ;nz => haven't shifted enough
  176.     RL    C
  177.     RET    NC        ;nc => positive number
  178.     EX    DE,HL
  179.     LD    HL,1        ;compensate for cy bit
  180.     SBC    HL,DE        ;negate result
  181.     RET
  182. ;
  183. ADDHALF: LD    HL,HALF
  184. HLADD:    CALL    LDBCHL
  185.     JR    FADD
  186. HALF:    DEFB    0,0,0,0,0,80H    ;0.5
  187. ;
  188. L247E:    CALL    PUSHFA
  189.     CALL    L27EC
  190.     POP    BC
  191.     POP    IX
  192.     POP    DE
  193.     JR    FADD
  194. ;
  195. HLSUB:    CALL    LDBCHL
  196.     JR    FSUB
  197. ;
  198. ;    fmod(z,x) = z-x*floor(z/x)
  199. ;        if x>0 then  0 <= fmod(z,x) < x
  200. ;        if x<0 then  x < fmod(z,x) <= 0
  201. ;
  202. QFMOD:    POP    HL    ;return addr
  203.     POP    DE    ;discard next number
  204.     POP    DE    ; (already in FA)
  205.     POP    DE
  206.     POP    DE    ;fetch next number
  207.     POP    IX    ; (1st operand, or "z")
  208.     POP    BC
  209.     PUSH    DE    ;restore stack
  210.     PUSH    DE
  211.     PUSH    DE
  212.     PUSH    DE
  213.     PUSH    DE
  214.     PUSH    DE
  215.     PUSH    HL    ;replace return addr
  216.     PUSH    DE    ;save another copy of z
  217.     PUSH    IX
  218.     PUSH    BC
  219.     CALL    PUSHFA    ;save a copy of 2nd operand ("x")
  220.     CALL    FDIV    ;z/x
  221.     CALL    QFLOOR    ;floor(z/x)
  222.     POP    BC
  223.     POP    IX
  224.     POP    DE
  225.     CALL    FMUL    ;x*floor(z/x)
  226.     POP    BC
  227.     POP    IX
  228.     POP    DE
  229. ;        to find mod(z,x)=z-x*floor(z/x), fall into...
  230. FSUB:    CALL    MINUSFA
  231.     JR    FADD
  232. ;
  233. ;    subtract the floating point accumulator from the value
  234. ;    on the stack (under the return address), leave result
  235. ;    in the floating point accumulator.
  236. ;
  237. DSUB:    CALL    MINUSFA
  238. ;
  239. ;    add the value on the stack (under the return address)
  240. ;    to the floating point accumulator
  241. ;
  242. DADD:    POP    HL    ;save return address
  243.     POP    DE
  244.     POP    IX
  245.     POP    BC
  246.     PUSH    HL    ;replace return address
  247. ;
  248. ;    add bc ix de to floating point accumulator
  249. ;
  250. FADD:    LD    A,B
  251.     OR    A
  252.     RET    Z    ;z => number to be added is zero
  253.     LD    A,(FA+5)
  254.     OR    A
  255.     JP    Z,LDFABC ;z => accumulator is zero,
  256. ;                just load number
  257.     SUB    B
  258.     JR    NC,ADD2 ;nc => accumulator has larger number
  259.     NEG        ;reverse accumulator & bc ix de...
  260.     EXX
  261.     PUSH    IX
  262.     CALL    LDBCFA
  263.     EXX
  264.     EX    (SP),IX
  265.     CALL    LDFABC
  266.     EXX
  267.     POP    IX    ;...end of reversing
  268. ADD2:    CP    29H
  269.     RET    NC    ;nc => addition makes no change
  270.     PUSH    AF    ;save difference of exponents
  271.     CALL    UNPACK    ;restore hidden bit & compare signs
  272.     LD    H,A    ;save difference in signs
  273.     POP    AF    ;recall difference of exponents
  274.     CALL    RSHIFT    ;shift  c ix de b  right by (a)
  275.     OR    H
  276.     LD    HL,FA
  277.     JP    P,ADD4    ;p => opposite signs, must subtract
  278.     CALL    FRADD    ;c ix de += FA
  279.     JR    NC,PACK    ;nc => adding caused no carry
  280.     INC    HL
  281.     INC    (HL)    ;increment exponent
  282.     JP    Z,OFLOW
  283.     LD    L,1
  284.     CALL    RSH8    ;shift  c ix de b  right by 1
  285.     JR    PACK    ;round, hide msb, & load into FA
  286. ;
  287. ADD4:    XOR    A    ;negate b...
  288.     SUB    B
  289.     LD    B,A
  290.     LD    A,(HL)    ;c ix de -= FA...
  291.     SBC    A,E
  292.     LD    E,A
  293.     INC    HL
  294.     LD    A,(HL)
  295.     SBC    A,D
  296.     LD    D,A
  297.     INC    HL
  298.     LD    A,(HL)
  299.     DEFB    0DDH
  300.     SBC    A,L
  301.     DEFB    0DDH
  302.     LD    L,A
  303.     INC    HL
  304.     LD    A,(HL)
  305.     DEFB    0DDH
  306.     SBC    A,H
  307.     DEFB    0DDH
  308.     LD    H,A
  309.     INC    HL
  310.     LD    A,(HL)
  311.     SBC    A,C
  312.     LD    C,A    ;...end of subtraction, fall into...
  313. ;
  314. ;    reverse sign if necessary (cy set) and normalize
  315. ;    (sign reversal necessary because we're using
  316. ;    sign-magnitude representation rather than
  317. ;    twos-complement)
  318. ;
  319. NORMA:    CALL    C,MINUSBC
  320. ;
  321. ;    normalize the 48 bit number in c ix de b
  322. ;    current exponent is in FA+5
  323. ;
  324. ;    result loaded into FA
  325. ;
  326. NORM:    LD    L,B
  327.     LD    H,E
  328.     XOR    A
  329. NORM2:    LD    B,A
  330.     LD    A,C
  331.     OR    A
  332.     JR    NZ,NORM12  ;nz => 7 or fewer shifts needed
  333. ;            shift c ix d hl  left by one byte
  334.     DEFB    0DDH
  335.     LD    C,H
  336.     DEFB    0DDH
  337.     LD    A,L
  338.     DEFB    0DDH
  339.     LD    H,A
  340.     DEFB    0DDH
  341.     LD    L,D
  342.     XOR    A
  343.     LD    D,H
  344.     LD    H,L
  345.     LD    L,A    ;...end of shifting
  346. ;
  347.     LD    A,B
  348.     SUB    8    ;adjust exponent
  349.     CP    0D0H
  350.     JR    NZ,NORM2
  351. ;
  352. NORM4:    XOR    A
  353. NORM6:    LD    (FA+5),A
  354.     RET
  355. ;
  356. NORM8:    DEC    B
  357. ;            shift  c ix d hl  left one bit...
  358.     ADD    HL,HL
  359.     RL    D
  360.     EX    AF,AF'
  361.     ADD    IX,IX
  362.     EX    AF,AF'
  363.     JR    NC,NORM10
  364.     INC    IX
  365. NORM10:    EX    AF,AF'
  366.     RL    C    ;...end of shifting
  367. ;
  368. NORM12:    JP    P,NORM8    ;p => high order bit still zero
  369.     LD    A,B
  370. ;            move number to  c ix de b
  371.     LD    E,H
  372.     LD    B,L
  373.     OR    A
  374.     JR    Z,PACK    ;z => exponent unchanged
  375.     LD    HL,FA+5        ;update exponent
  376.     ADD    A,(HL)
  377.     LD    (HL),A
  378.     JR    NC,NORM4    ;nc => underflow (set to 0)
  379.     RET    Z        ;z => underflow (leave as 0)
  380. PACK:    LD    A,B
  381. PACK2:    LD    HL,FA+5    ;round c ix de b to 40 bits
  382.     OR    A
  383.     CALL    M,INCR
  384.     LD    B,(HL)    ;load exponent
  385.     INC    HL
  386.     LD    A,(HL)    ;recover sign
  387.     AND    80H    ;mask out all but sign
  388.     XOR    C    ;add to high
  389.     LD    C,A    ;   order byte
  390.     JP    LDFABC    ;place answer in FA
  391. ;
  392. INCR:    INC    E    ;increment c ix de
  393.     RET    NZ
  394.     INC    D
  395.     RET    NZ
  396.     DEFB    0DDH
  397.     INC    L
  398.     RET    NZ
  399.     DEFB    0DDH
  400.     INC    H
  401.     RET    NZ
  402.     INC    C
  403.     RET    NZ    ;z => carry
  404.     LD    C,80H    ;set high order bit
  405.     INC    (HL)    ;   and increment exponent
  406.     RET    NZ
  407.     JP    OFLOW
  408. ;
  409. ;    fraction add: c ix de += (hl)
  410. ;
  411. FRADD:    LD    A,(HL)
  412.     ADD    A,E
  413.     LD    E,A
  414.     INC    HL
  415.     LD    A,(HL)
  416.     ADC    A,D
  417.     LD    D,A
  418.     INC    HL
  419.     LD    A,(HL)
  420.     DEFB    0DDH
  421.     ADC    A,L
  422.     DEFB    0DDH
  423.     LD    L,A
  424.     INC    HL
  425.     LD    A,(HL)
  426.     DEFB    0DDH
  427.     ADC    A,H
  428.     DEFB    0DDH
  429.     LD    H,A
  430.     INC    HL
  431.     LD    A,(HL)
  432.     ADC    A,C
  433.     LD    C,A
  434.     RET
  435. ;
  436. ;    complement FASIGN and negate the fraction c ix de b
  437. ;
  438. MINUSBC: LD    HL,FASIGN
  439.     LD    A,(HL)
  440.     CPL
  441.     LD    (HL),A
  442.     XOR    A
  443.     LD    L,A
  444.     LD    H,A
  445.     SUB    B
  446.     LD    B,A
  447.     LD    A,L
  448.     SBC    HL,DE
  449.     EX    DE,HL
  450.     LD    L,A
  451.     DEFB    0DDH
  452.     SBC    A,L
  453.     DEFB    0DDH
  454.     LD    L,A
  455.     LD    A,L
  456.     DEFB    0DDH
  457.     SBC    A,H
  458.     DEFB    0DDH
  459.     LD    H,A
  460.     LD    A,L
  461.     SBC    A,C
  462.     LD    C,A
  463.     RET
  464. ;
  465. ;    shift  c ix de b  right by (a)
  466. ;
  467. RSHIFT:    LD    B,0
  468. RSH2:    SUB    8
  469.     JR    C,RSH4    ;c => 7 or fewer shifts remain
  470.     LD    B,E    ;shift  c ix de b  right by 8...
  471.     LD    E,D
  472.     DEFB    0DDH
  473.     LD    D,L
  474.     EX    AF,AF'
  475.     DEFB    0DDH
  476.     LD    A,H
  477.     DEFB    0DDH
  478.     LD    L,A
  479.     EX    AF,AF'
  480.     DEFB    0DDH
  481.     LD    H,C
  482.     LD    C,0    ;...end of shifting
  483.     JR    RSH2
  484. ;
  485. RSH4:    ADD    A,9
  486.     LD    L,A
  487. RSH6:    XOR    A
  488.     DEC    L
  489.     RET    Z    ;z => requested shift is complete
  490.     LD    A,C
  491. RSH8:    RRA        ;shift  c ix de b  right by one...
  492.     LD    C,A
  493.     DEFB    0DDH
  494.     LD    A,H
  495.     RRA
  496.     DEFB    0DDH
  497.     LD    H,A
  498.     DEFB    0DDH
  499.     LD    A,L
  500.     RRA
  501.     DEFB    0DDH
  502.     LD    L,A
  503.     RR    D
  504.     RR    E
  505.     RR    B    ;...end of shifting
  506.     JR    RSH6
  507. ;
  508. ;    multiply the floating point accumulator by the value
  509. ;    on the stack (under the return address), leave result
  510. ;    in the floating point accumulator.
  511. ;
  512. DMUL:    POP    HL    ;return addr
  513.     POP    DE    ;multiplier...
  514.     POP    IX
  515.     POP    BC
  516.     PUSH    HL    ;replace return addr
  517. ;
  518. ;    multiply floating point accumulator by  bc ix de
  519. ;
  520. FMUL:    CALL    SGN
  521.     RET    Z    ; z => accumulator has zero
  522.     LD    L,0    ;"product" flag
  523.     CALL    DIV14    ;find exponent of product
  524.     LD    A,C  ;c' h'l' d'e' (multiplicand) = c ix de...
  525.     PUSH    DE
  526.     EXX
  527.     LD    C,A
  528.     POP    DE
  529.     PUSH    IX
  530.     POP    HL
  531.     EXX        ;...end of multiplicand loading
  532.     LD    BC,0    ; c ix de b (product) = 0...
  533.     LD    D,B
  534.     LD    E,B
  535.     LD    IX,0
  536.     LD    HL,NORM    ; push addr of normalize routine
  537.     PUSH    HL
  538.     LD    HL,MULLOOP    ; push addr of top of loop
  539.     PUSH    HL    ; (5 iterations wanted,
  540.     PUSH    HL    ; once per byte of fraction)
  541.     PUSH    HL
  542.     PUSH    HL
  543.     LD    HL,FA    ;point to LSB
  544. MULLOOP: LD    A,(HL)    ;get next byte of multiplier
  545.     INC    HL
  546.     OR    A
  547.     JR    NZ,MUL2    ; z => next 8 bits of multiplier are 0
  548.     LD    B,E    ;shift  c ix de b  right by 8...
  549.     LD    E,D
  550.     DEFB    0DDH
  551.     LD    D,L
  552.     EX    AF,AF'
  553.     DEFB    0DDH
  554.     LD    A,H
  555.     DEFB    0DDH
  556.     LD    L,A
  557.     EX    AF,AF'
  558.     DEFB    0DDH
  559.     LD    H,C
  560.     LD    C,A    ;...end of shifting
  561.     RET        ;go to top of loop or NORM
  562. ;
  563. MUL2:    PUSH    HL    ;save multiplier pointer
  564.     EX    DE,HL
  565.     LD    E,8    ;8 iterations (once per bit)
  566. MUL4:    RRA        ;rotate a multiplier bit into cy
  567.     LD    D,A
  568.     LD    A,C
  569.     JR    NC,MUL6    ; nc => no addition needed
  570.     PUSH    HL    ; c ix hl (product)  +=
  571.     EXX        ;    c' h'l' d'e' (multiplicand)
  572.     EX    (SP),HL
  573.     ADD    HL,DE
  574.     EX    (SP),HL
  575.     EX    DE,HL
  576.     PUSH    IX
  577.     EX    (SP),HL
  578.     ADC    HL,DE
  579.     EX    (SP),HL
  580.     POP    IX
  581.     EX    DE,HL
  582.     ADC    A,C
  583.     EXX
  584.     POP    HL
  585. ;
  586. MUL6:    RRA       ;shift  c ix hl b (product)  right by 1...
  587.     LD    C,A
  588.     DEFB    0DDH
  589.     LD    A,H
  590.     RRA
  591.     DEFB    0DDH
  592.     LD    H,A
  593.     DEFB    0DDH
  594.     LD    A,L
  595.     RRA
  596.     DEFB    0DDH
  597.     LD    L,A
  598.     RR    H
  599.     RR    L
  600.     RR    B        ;...end of shifting
  601.     DEC    E
  602.     LD    A,D
  603.     JR    NZ,MUL4        ; z => 8 iterations complete
  604.     EX    DE,HL
  605. MUL8:    POP    HL        ;recover multiplier pointer
  606.     RET            ;go to top of loop or NORM
  607. ;
  608. ;    divide floating point accumulator by 10
  609. ;
  610. DIV10:    CALL    PUSHFA
  611.     LD    BC,8420H    ; 10.0
  612.     LD    IX,0
  613.     LD    DE,0
  614.     CALL    LDFABC
  615. DIV1:    POP    BC
  616.     POP    IX
  617.     POP    DE
  618.     JR    FDIV
  619. ;
  620. ;    divide the value on the stack (under the return
  621. ;    address) by the floating point accumulator, leave
  622. ;    result in the floating point accumulator.
  623. ;
  624. DDIV:    POP    HL    ;save return address
  625.     POP    DE
  626.     POP    IX
  627.     POP    BC
  628.     PUSH    HL    ;replace return address
  629. ;
  630. ;    divide  bc ix de  by FA, leave result in FA
  631. ;
  632. FDIV:    CALL    SGN
  633.     JP    Z,DIVZERO ; z => attempting to divide by 0
  634.     LD    L,0FFH    ;"quotient" flag
  635.     CALL    DIV14    ;find quotient exponent
  636.     PUSH    IY
  637.     INC    (HL)
  638.     INC    (HL)
  639.     DEC    HL
  640.     PUSH    HL    ; c' h'l' d'e' (divisor) = FA...
  641.     EXX
  642.     POP    HL
  643.     LD    C,(HL)
  644.     DEC    HL
  645.     LD    D,(HL)
  646.     DEC    HL
  647.     LD    E,(HL)
  648.     DEC    HL
  649.     LD    A,(HL)
  650.     DEC    HL
  651.     LD    L,(HL)
  652.     LD    H,A
  653.     EX    DE,HL
  654.     EXX
  655.     LD    B,C    ; b iy hl (dividend) = c ix de...
  656.     EX    DE,HL
  657.     PUSH    IX
  658.     POP    IY
  659.     XOR    A    ; c ix de (quotient) = 0...
  660.     LD    C,A
  661.     LD    D,A
  662.     LD    E,A
  663.     LD    IX,0
  664.     LD    (EXTRA),A
  665. DIV2:    PUSH    HL    ;save b iy hl in case the subtraction
  666.     PUSH    IY    ; proves to be unnecessary
  667.     PUSH    BC
  668.     PUSH    HL    ; EXTRA b iy hl (dividend)  -=
  669.     LD    A,B    ;    c' h'l' d'e' (divisor)...
  670.     EXX
  671.     EX    (SP),HL
  672.     OR    A
  673.     SBC    HL,DE
  674.     EX    (SP),HL
  675.     EX    DE,HL
  676.     PUSH    IY
  677.     EX    (SP),HL
  678.     SBC    HL,DE
  679.     EX    (SP),HL
  680.     POP    IY
  681.     EX    DE,HL
  682.     SBC    A,C
  683.     EXX
  684.     POP    HL
  685.     LD    B,A
  686.     LD    A,(EXTRA)
  687.     SBC    A,0
  688.     CCF
  689.     JR    NC,DIV4    ; nc => subtraction caused carry
  690.     LD    (EXTRA),A
  691.     POP    AF    ;discard saved value of dividend...
  692.     POP    AF
  693.     POP    AF
  694.     SCF
  695.     JR    DIV6
  696. DIV4:    POP    BC    ;restore dividend...
  697.     POP    IY
  698.     POP    HL
  699. ;
  700. DIV6:    INC    C
  701.     DEC    C
  702.     RRA
  703.     JP    M,DIV12
  704.     RLA      ;shift  c ix de a (quotient)  left by 1...
  705.     RL    E
  706.     RL    D
  707.     EX    AF,AF'    ;(these 6 lines are  adc ix,ix...)
  708.     ADD    IX,IX
  709.     EX    AF,AF'
  710.     JR    NC,DIV8
  711.     INC    IX
  712. DIV8:    EX    AF,AF'
  713.     RL    C    ;...end of  c ix de a  shifting
  714.     ADD    HL,HL    ;shift  EXTRA b iy hl  left by 1...
  715.     EX    AF,AF'
  716.     ADD    IY,IY
  717.     EX    AF,AF'
  718.     JR    NC,DIV9
  719.     INC    IY
  720. DIV9:    EX    AF,AF'
  721.     RL    B
  722.     LD    A,(EXTRA)
  723.     RLA
  724.     LD    (EXTRA),A  ;...end of  EXTRA b iy hl  shifting
  725.     LD    A,C    ;test  c ix de...
  726.     OR    D
  727.     OR    E
  728.     DEFB    0DDH
  729.     OR    H
  730.     DEFB    0DDH
  731.     OR    L    ;...end of  c ix de  testing
  732.     JR    NZ,DIV2    ;nz => dividend nonzero
  733.     PUSH    HL
  734.     LD    HL,FA+5
  735.     DEC    (HL)
  736.     POP    HL
  737.     JR    NZ,DIV2
  738.     JR    OFLOW2
  739. ;
  740. DIV12:    POP    IY
  741.     JP    PACK2
  742. ;
  743. ;    find exponent for product (L=0) or quotient (L=ff)
  744. ;
  745. DIV14:    LD    A,B
  746.     OR    A
  747.     JR    Z,DIV20
  748.     LD    A,L    ;get product/quotient flag
  749.     LD    HL,FA+5
  750.     XOR    (HL)    ;get +-FA exponent
  751.     ADD    A,B    ;find and...
  752.     LD    B,A    ;...load new exponent
  753.     RRA
  754.     XOR    B
  755.     LD    A,B
  756.     JP    P,DIV18
  757.     ADD    A,80H
  758.     LD    (HL),A
  759.     JP    Z,MUL8
  760.     CALL    UNPACK    ;restore hidden bits & compare signs
  761.     LD    (HL),A    ;save difference of signs
  762.     DEC    HL    ;point to MSB of fraction
  763.     RET
  764. ;
  765. DIV17:    CALL    SGN
  766.     CPL
  767.     OR    A
  768.     DEFB    21H    ;"ignore next 2 bytes"
  769. DIV18:    OR    A
  770. DIV20:    POP    HL
  771.     JP    P,NORM4
  772. OFLOW2:    JP    OFLOW
  773. ;
  774. ;    multiply FA by 10
  775. ;
  776. MUL10:    CALL    LDBCFA
  777.     LD    A,B    ;multiply bc ix de by 4...
  778.     OR    A
  779.     RET    Z
  780.     ADD    A,2
  781.     JR    C,OFLOW2
  782.     LD    B,A    ;...end of multiplication
  783.     CALL    FADD    ;add to FA, yields FA*5
  784.     LD    HL,FA+5
  785.     INC    (HL)    ;double again, yielding FA*10
  786.     RET    NZ
  787.     JR    OFLOW2
  788. ;
  789. ; L27DB:    LD    BC,9980H    ; -2**24
  790.     LD    IX,0
  791.     LD    DE,0
  792.     CALL    COMPARE
  793.     RET    Z
  794.     JP    ILLFCT
  795. ;
  796. L27EC:    LD    B,88H    ; 128.
  797.     LD    DE,0
  798. L27F1:    LD    HL,FA+5
  799.     LD    C,A
  800.     PUSH    DE
  801.     POP    IX
  802.     LD    DE,0
  803.     LD    (HL),B    ;store exponent
  804.     LD    B,0
  805.     INC    HL
  806.     LD    (HL),80H    ;store minus sign
  807.     RLA
  808.     JP    NORMA
  809. ;
  810.     EX    DE,HL
  811.     XOR    A
  812.     LD    B,98H
  813.     JR    L27F1
  814.     LD    B,C
  815. L280C:    LD    D,B
  816.     LD    E,0
  817.     LD    HL,L0F2E
  818.     LD    (HL),E
  819.     LD    B,90H
  820.     JR    L27F1
  821.     LD    B,A
  822.     XOR    A
  823.     JR    L280C
  824. ;
  825. L281B:    CALL    SGN
  826.     JP    M,ILLFCT
  827.     LD    A,(FA+5)
  828.     CP    91H
  829.     JP    C,INT2
  830.     LD    BC,9180H    ; -2**16
  831.     LD    IX,0
  832.     LD    DE,0
  833.     CALL    COMPARE
  834.     LD    D,C
  835.     RET    Z
  836. L2838:    JP    ILLFCT
  837.     CALL    L281B
  838.     LD    A,D
  839.     OR    A
  840.     JR    NZ,L2838
  841.     LD    A,E
  842.     RET
  843. ;
  844. ;    set z & s flags per FA
  845. ;
  846. SGN:    LD    A,(FA+5)
  847.     OR    A
  848.     RET    Z
  849.     LD    A,(FA+4)
  850.     DEFB    0FEH    ;"ignore next byte"
  851. SETFLGS: CPL
  852.     RLA
  853.     SBC    A,A
  854.     RET    NZ
  855.     INC    A
  856.     RET
  857. ;
  858. ;    Double precision comparisons
  859. ;
  860. ;    each compares top of stack
  861. ;    (under two return addresses) to FA
  862. ;
  863. ;TOS >= FA?
  864. DGE:    CALL    DCOMPAR
  865.     JR    Z,YES    ;z => equal
  866.     JR    DG    ;remaining tests are shared
  867. ;
  868. ;TOS > FA?
  869. DGT:    CALL    DCOMPAR
  870.     JR    Z,NO    ;z => equal
  871. DG:    JP    P,NO    ;p => not greater than
  872. YES:    LD    HL,1    ;load "true"
  873.     RET
  874. ;
  875. ;TOS <= FA?
  876. DLE:    CALL    DCOMPAR
  877.     JR    Z,YES
  878.     JR    DL
  879. ;
  880. ;TOS < FA?
  881. DLT:    CALL    DCOMPAR
  882.     JR    Z,NO
  883. DL:    JP    P,YES    ;p => less than
  884. NO:    LD    HL,0    ;load "false"
  885.     RET
  886. ;
  887. ;TOS == FA?
  888. DEQ:    CALL    DCOMPAR
  889.     JR    Z,YES
  890.     JR    NO
  891. ;
  892. ;TOS != FA?
  893. DNE:    CALL    DCOMPAR
  894.     JR    NZ,YES
  895.     JR    NO
  896. ;
  897. ;common routine to perform double precision comparisons
  898. DCOMPAR: POP    HL    ;save 1st return addr
  899.     POP    IY    ;save 2nd return addr
  900.     POP    DE    ;get number to compare
  901.     POP    IX
  902.     POP    BC
  903.     PUSH    IY    ;replace 2nd addr
  904.     PUSH    HL    ;replace 1st addr, fall into...
  905. ;
  906. ;    sets flags per FA - (bc ix de)
  907. ;
  908. COMPARE: LD    A,B
  909.     OR    A
  910.     JR    Z,SGN        ;bc ix de = 0, so
  911. ;                sign of FA is result
  912.  
  913.     CALL    SGN
  914.     LD    A,C
  915.     JR    Z,SETFLGS    ;FA = 0, so sign of
  916. ;                -(bc ix de) is result
  917.     LD    HL,FA+4
  918.     XOR    (HL)
  919.     LD    A,C
  920.     JP    M,SETFLGS    ;operands have opposite
  921. ;                signs, so result is sign
  922. ;                of -(bc ix de)
  923.  
  924.     CALL    CPFRAC
  925.     RRA            ;recover cy bit
  926.     XOR    C    ;reverse sign if numbers are negative
  927.     JR    SETFLGS
  928. ;
  929. CPFRAC:    INC    HL    ;compare  bc ix de  to (HL)
  930.     LD    A,B
  931.     CP    (HL)
  932.     RET    NZ
  933.     DEC    HL
  934.     LD    A,C
  935.     CP    (HL)
  936.     RET    NZ
  937.     DEC    HL
  938.     DEFB    0DDH
  939.     LD    A,H
  940.     CP    (HL)
  941.     RET    NZ
  942.     DEC    HL
  943.     DEFB    0DDH
  944.     LD    A,L
  945.     CP    (HL)
  946.     RET    NZ
  947.     DEC    HL
  948.     LD    A,D
  949.     CP    (HL)
  950.     RET    NZ
  951.     DEC    HL
  952.     LD    A,E
  953.     SUB    (HL)
  954.     RET    NZ
  955.     POP    HL    ;return zero to program
  956.     RET        ;that called "COMPARE"
  957. ;
  958. QFABS:    CALL    SGN
  959.     RET    P
  960. MINUSFA: LD    HL,FA+4
  961.     LD    A,(HL)
  962.     XOR    80H
  963.     LD    (HL),A
  964.     RET
  965. ;
  966. PUSHFA:    EX    DE,HL
  967. L2895:    LD    HL,(FA)
  968.     EX    (SP),HL
  969.     PUSH    HL
  970.     LD    HL,(FA+2)
  971.     EX    (SP),HL
  972.     PUSH    HL
  973.     LD    HL,(FA+4)
  974.     EX    (SP),HL
  975.     PUSH    HL
  976.     EX    DE,HL
  977.     RET
  978. ;            This can be reinstated when the compiler
  979. ;            understands "extern". Until then, pi
  980. ;            can't be declared without reserving
  981. ;            storage again.
  982. ; QPI:    DEFW    0A222H    ;3.1415926535 = pi
  983. ;     DEFW    0FDAH
  984. ;     DEFW    8249H
  985. ;
  986. ;    FA = (hl)
  987. ;
  988. DLOAD:    LD    DE,FA
  989.     LD    BC,6
  990.     LDIR
  991.     RET
  992. ;
  993. ;    exchange floating point accumulator with
  994. ;    top of stack (under return address)
  995. ;
  996. DSWAP:    POP    HL    ;return addr
  997.     POP    DE
  998.     POP    IX
  999.     POP    BC
  1000.     EXX        ;protect the values
  1001.     CALL    DPUSH    ;push FA
  1002.     EXX        ;recover the values
  1003.     PUSH    HL    ;replace return addr, fall into...
  1004. ;
  1005. ;    FA = bc ix de
  1006. ;
  1007. LDFABC:    LD    (FA),DE
  1008.     LD    (FA+2),IX
  1009.     LD    (FA+4),BC
  1010.     RET
  1011. ;
  1012. ;    bc ix de = FA
  1013. ;
  1014. LDBCFA:    LD    DE,(FA)
  1015.     LD    IX,(FA+2)
  1016.     LD    BC,(FA+4)
  1017.     RET
  1018. ;
  1019. ;    bc ix de = (hl)
  1020. ;
  1021. LDBCHL:    LD    E,(HL)
  1022.     INC    HL
  1023.     LD    D,(HL)
  1024.     INC    HL
  1025.     LD    C,(HL)
  1026.     DEFB    0DDH
  1027.     LD    L,C
  1028.     INC    HL
  1029.     LD    C,(HL)
  1030.     DEFB    0DDH
  1031.     LD    H,C
  1032.     INC    HL
  1033.     LD    C,(HL)
  1034.     INC    HL
  1035.     LD    B,(HL)
  1036.     INC    HL
  1037.     RET
  1038. ;
  1039. ;    (hl) = FA
  1040. ;
  1041. DSTORE:    LD    DE,FA
  1042.     LD    BC,6
  1043.     EX    DE,HL
  1044.     LDIR
  1045.     EX    DE,HL
  1046.     RET
  1047. ;
  1048. UNPACK:    LD    HL,FA+4
  1049.     LD    A,(HL)    ;get MSB of fraction
  1050.     RLCA        ;rotate sign bit into lsb
  1051.     SCF        ;set carry
  1052.     RRA        ;rotate sign bit into cy, cy into msb
  1053.     LD    (HL),A    ;restore MSB (with hidden bit restored)
  1054.     CCF        ;complement sign bit...
  1055.     RRA
  1056.     INC    HL
  1057.     INC    HL
  1058.     LD    (HL),A    ;...and save in msb of FASIGN
  1059.     LD    A,C    ;similarly, get sign bit of bc ix de...
  1060.     RLCA
  1061.     SCF
  1062.     RRA
  1063.     LD    C,A    ;...restore hidden bit...
  1064.     RRA
  1065.     XOR    (HL)    ;...and compare with sign of FA.
  1066.     RET
  1067. ;
  1068. INT2:    LD    B,A    ;if a==0, return with  bc ix de = 0...
  1069.     LD    C,A
  1070.     LD    D,A
  1071.     LD    E,A
  1072.     DEFB    0DDH
  1073.     LD    H,A
  1074.     DEFB    0DDH
  1075.     LD    L,A
  1076.     OR    A
  1077.     RET    Z
  1078.     PUSH    HL
  1079.     CALL    LDBCFA    ;copy FA into bc ix de,
  1080.     CALL    UNPACK    ; restore hidden bits
  1081.     XOR    (HL)
  1082.     LD    H,A    ;put sign in msb of h
  1083.     JP    P,INT4 ;p => positive number
  1084.     DEC    DE    ;decrement c ix de...
  1085.     LD    A,D
  1086.     AND    E
  1087.     INC    A
  1088.     JR    NZ,INT4
  1089.     DEC    IX
  1090.     DEFB    0DDH
  1091.     LD    A,H
  1092.     DEFB    0DDH
  1093.     AND    L
  1094.     INC    A
  1095.     JR    NZ,INT4
  1096.     DEC    C    ;...end of c ix de decrementing
  1097. ;
  1098. INT4:    LD    A,0A8H    ;shift  c ix de  right so bits to
  1099.     SUB    B    ; the right of the binary point
  1100.     CALL    RSHIFT    ; are discarded
  1101.     LD    A,H
  1102.     RLA
  1103.     CALL    C,INCR    ;c => negative, increment  c ix de
  1104.     LD    B,0
  1105.     CALL    C,MINUSBC ;negate the fraction c ix de
  1106.     POP    HL
  1107.     RET
  1108. ;
  1109.     POP    BC
  1110.     POP    IX
  1111.     POP    DE
  1112. ;
  1113. ;    divide with integer result
  1114. ;    (truncates toward zero)
  1115. ;
  1116. DIVI:    CALL    FDIV
  1117.     CALL    SGN
  1118.     JP    P,QFLOOR
  1119.     CALL    MINUSFA
  1120.     CALL    QFLOOR
  1121.     JP    MINUSFA
  1122. ;
  1123. ;    return -(floor(-x))
  1124. QCEIL:    CALL    ODD
  1125. ;
  1126. ;    return largest integer not greater than
  1127. QFLOOR:    LD    HL,FA+5
  1128.     LD    A,(HL)    ;fetch exponent
  1129.     CP    0A8H
  1130.     LD    A,(FA)
  1131.     RET    NC    ;nc => binary point is right of lsb
  1132.     LD    A,(HL)
  1133.     CALL    INT2
  1134.     LD    (HL),0A8H  ;place binary pt at end of fraction
  1135.     LD    A,E
  1136.     PUSH    AF
  1137.     LD    A,C
  1138.     RLA
  1139.     CALL    NORMA
  1140.     POP    AF
  1141.     RET
  1142. ;
  1143.     LD    HL,FA+5
  1144.     LD    (HL),0A8H
  1145.     INC    HL
  1146.     LD    (HL),80H
  1147.     LD    A,C
  1148.     RLA
  1149.     JP    NORMA
  1150. ;    amax(a,b)    returns the greater of a and b
  1151. QAMAX:    LD    HL,8    ;offset for 1st argument
  1152.     ADD    HL,SP
  1153.     CALL    LDBCHL    ;bcixde := 1st argument
  1154.     CALL    COMPARE
  1155.     JP    M,LDFABC
  1156.     RET
  1157. ;
  1158. ;    amin(a,b)
  1159. QAMIN:    LD    HL,8
  1160.     ADD    HL,SP
  1161.     CALL    LDBCHL
  1162.     CALL    COMPARE
  1163.     JP    P,LDFABC
  1164.     RET
  1165. ;
  1166. ;    negate FA, and push address of MINUSFA
  1167. ;    called to evaluate functions f(x) when the argument is
  1168. ;    negative and f() satisfies f(-x)=-f(x)
  1169. ODD:    CALL    MINUSFA
  1170. L29D1:    LD    HL,MINUSFA
  1171.     EX    (SP),HL
  1172.     JP    (HL)
  1173. ;
  1174. QRAND:    CALL    SGN
  1175.     LD    HL,SEED
  1176.     JP    M,DSTORE
  1177.     CALL    DLOAD
  1178.     RET    Z
  1179.     LD    BC,9835H    ; 11879545.
  1180.     LD    IX,447AH
  1181.     LD    DE,0
  1182.     CALL    FMUL
  1183.     LD    BC,6828H    ; 3.92767775e-8
  1184.     LD    IX,0B146H
  1185.     LD    DE,0
  1186.     CALL    FADD
  1187.     CALL    LDBCFA
  1188.     LD    A,E
  1189.     LD    E,C
  1190.     LD    C,A
  1191.     LD    HL,FASIGN
  1192.     LD    (HL),80H
  1193.     DEC    HL
  1194.     LD    B,(HL)
  1195.     LD    (HL),80H
  1196.     CALL    NORM
  1197.     LD    HL,SEED
  1198.     JP    DSTORE
  1199. #endasm
  1200.