home *** CD-ROM | disk | FTP | other *** search
/ Fractal Creations (Second Edition) / FRACTALS_2E.iso / frasrc.exe / FPU087.ASM < prev    next >
Assembly Source File  |  1992-12-19  |  39KB  |  1,538 lines
  1. TITLE fpu087.asm (C) 1989, Mark C. Peterson, CompuServe [70441,3353]
  2. SUBTTL All rights reserved.
  3. ;
  4. ;  Code may be used in any program provided the author is credited
  5. ;    either during program execution or in the documentation.  Source
  6. ;    code may be distributed only in combination with public domain or
  7. ;    shareware source code.  Source code may be modified provided the
  8. ;    copyright notice and this message is left unchanged and all
  9. ;    modifications are clearly documented.
  10. ;
  11. ;    I would appreciate a copy of any work which incorporates this code.
  12. ;
  13. ;    Mark C. Peterson
  14. ;    405-C Queen St., Suite #181
  15. ;    Southington, CT 06489
  16. ;    (203) 276-9721
  17. ;
  18. ;  References:
  19. ;     The VNR Concise Encyclopedia of Mathematics
  20. ;        by W. Gellert, H. Hustner, M. Hellwich, and H. Kastner
  21. ;        Published by Van Nostrand Reinhold Comp, 1975
  22. ;
  23. ;     80386/80286 Assembly Language Programming
  24. ;        by William H. Murray, III and Chris H. Pappas
  25. ;        Published by Osborne McGraw-Hill, 1986
  26. ;
  27. ;  History since Fractint 16.3:
  28. ;     CJLT changed e2x and Log086 algorithms for more speed
  29. ;     CJLT corrected SinCos086 for -ve angles in 2nd and 4th quadrants
  30. ;     CJLT speeded up SinCos086 for angles >45 degrees in any quadrant
  31. ;     (See comments containing the string `CJLT')
  32. ; 14 Aug 91 CJLT removed r16Mul - not called from anywhere
  33. ; 21 Aug 91 CJLT corrected Table[1] from 6 to b
  34. ;                improved bx factors in Log086 for more accuracy
  35. ;                corrected Exp086 overflow detection to 15 from 16 bits.
  36. ; 07 Sep 92 MP   corrected problem in FPUcplxlog
  37. ; 07 Sep 92 MP   added argument test for FPUcplxdiv
  38. ; 06 Nov 92 CAE  made some varibles public for PARSERA.ASM
  39. ; 07 Dec 92 CAE  sped up FPUsinhcosh function
  40. ;
  41. ; CJLT=      Chris J Lusby Taylor
  42. ;            32 Turnpike Road
  43. ;            Newbury, England (where's that?)
  44. ;        Contactable via Compuserve user Stan Chelchowski [100016,351]
  45. ;                     or Tel 011 44 635 33270 (home)
  46. ;
  47. ; CAE=   Chuck Ebbert  CompuServe [76306,1226]
  48. ;
  49. ;
  50. ;PROCs in this module:
  51. ;FPUcplxmul     PROC     x:word, y:word, z:word
  52. ;FPUcplxdiv     PROC     x:word, y:word, z:word
  53. ;FPUcplxlog     PROC     x:word, z:word
  54. ;FPUsinhcosh    PROC     x:word, sinh:word, cosh:word
  55. ;FPUsincos  PROC  x:word, sinx:word, cosx:word
  56. ;r16Mul     PROC    uses si di, x1:word, x2:word, y1:word, y2:word
  57. ;RegFloat2Fg     PROC    x1:word, x2:word, Fudge:word
  58. ;RegFg2Float     PROC   x1:word, x2:word, FudgeFact:byte
  59. ;ExpFudged      PROC    uses si, x_low:word, x_high:word, Fudge:word
  60. ;LogFudged      PROC    uses si di, x_low:word, x_high:word, Fudge:word
  61. ;LogFloat14     PROC     x1:word, x2:word
  62. ;RegSftFloat     PROC   x1:word, x2:word, Shift:byte
  63. ;RegDivFloat     PROC  uses si di, x1:word, x2:word, y1:word, y2:word
  64. ;SinCos086   PROC     uses si di, LoNum:WORD, HiNum:WORD, SinAddr:WORD, \
  65. ;_e2xLT   PROC        ;argument in dx.ax (bitshift=16 is hard coded here)
  66. ;Exp086    PROC     uses si di, LoNum:WORD, HiNum:WORD
  67. ;SinhCosh086    PROC     uses si di, LoNum:WORD, HiNum:WORD, SinhAddr:WORD, \
  68. ;Log086   PROC     uses si di, LoNum:WORD, HiNum:WORD, Fudge:WORD
  69.  
  70.  
  71. IFDEF ??version
  72. MASM51
  73. QUIRKS
  74. ENDIF
  75.  
  76. .model medium, c
  77.  
  78. extrn cos:far
  79. extrn _Loaded387sincos:far
  80. extrn compiled_by_turboc:word
  81.  
  82.  
  83. .data
  84.  
  85. extrn cpu:WORD
  86.  
  87. PUBLIC TrigLimit, TrigOverflow
  88.  
  89. ; CAE 6Nov92 made these public for PARSERA.ASM */
  90. PUBLIC _1_, _2_, PointFive, infinity
  91.  
  92. PiFg13         dw       6487h
  93. InvPiFg33      dd       0a2f9836eh
  94. InvPiFg16      dw       517ch
  95. Ln2Fg16        dw       0b172h ;ln(2) * 2^16 . True value is b172.18
  96. Ln2Fg15        dw       058b9h ;used by e2xLT. True value is 58b9.0C
  97. TrigOverflow   dw       0
  98. TrigLimit      dd       0
  99. ;
  100. ;Table of 2^(n/16) for n=0 to 15. All entries fg15.
  101. ;Used by e2xLT
  102. ;
  103. Table        dw    08000h,085abh,08b96h,091c4h
  104.         dw    09838h,09ef5h,0a5ffh,0ad58h
  105.         dw    0b505h,0bd09h,0c567h,0ce25h
  106.         dw    0d745h,0e0cdh,0eac1h,0f525h
  107. one            dw       ?
  108. expSign        dw       ?
  109. a              dw       ?
  110. SinNeg         dw       ?    ;CJLT - not now needed by SinCos086, but by
  111. CosNeg           dw    ?    ;       ArcTan086
  112. Ans           dq    ?
  113. fake_es        dw       ?         ; <BDT> Windows can't use ES for storage
  114.  
  115. TaylorTerm  MACRO
  116. LOCAL Ratio
  117.    add   Factorial, one
  118.    jnc   SHORT Ratio
  119.  
  120.    rcr   Factorial, 1
  121.    shr   Num, 1
  122.    shr   one, 1
  123.  
  124. Ratio:
  125.    mul   Num
  126.    div   Factorial
  127. ENDM
  128.  
  129.  
  130.  
  131. _4_         dq    4.0
  132. _2_         dq    2.0
  133. _1_         dq    1.0
  134. PointFive   dq    0.5
  135. temp        dq     ?
  136. Sign        dw     ?
  137.  
  138. extrn fpu:word
  139.  
  140. .code
  141.  
  142.  
  143. FPUcplxmul     PROC     x:word, y:word, z:word
  144.    mov   bx, x
  145.    fld   QWORD PTR [bx]       ; x.x
  146.    fld   QWORD PTR [bx+8]     ; x.y, x.x
  147.    mov   bx, y
  148.    fld   QWORD PTR [bx]       ; y.x, x.y, x.x
  149.    fld   QWORD PTR [bx+8]     ; y.y, y.x, x.y, x.x
  150.    mov   bx, z
  151.    fld   st                   ; y.y, y.y, y.x, x.y, x.x
  152.    fmul  st, st(3)            ; y.y*x.y, y.y. y.x, x.y, x.x
  153.    fld   st(2)                ; y.x, y.y*x.y, y.y, y.x, x.y, x.x
  154.    fmul  st, st(5)            ; y.x*x.x, y.y*x.y, y.y, y.x, x.y, x.x
  155.    fsubr                      ; y.x*x.x - y.y*x.y, y.y, y.x, x.y, x.x
  156.    fstp  QWORD PTR [bx]       ; y.y, y.x, x.y, x.x
  157.    fmulp st(3), st            ; y.x, x.y, x.x*y.y
  158.    fmul                       ; y.x*x.y, x.x*y.y
  159.    fadd                       ; y.x*x.y + x.x*y.y
  160.    fstp  QWORD PTR [bx+8]
  161.    ret
  162. FPUcplxmul     ENDP
  163.  
  164.  
  165. .DATA
  166.  
  167. infinity          dq    1.0E+300
  168.  
  169. PUBLIC DivideOverflow
  170. DivideOverflow    dw    0
  171.  
  172. .CODE
  173.  
  174.  
  175. FPUcplxdiv     PROC     x:word, y:word, z:word
  176. LOCAL Status:WORD
  177.    mov   bx, x
  178.    fld   QWORD PTR [bx]       ; x.x
  179.    fld   QWORD PTR [bx+8]     ; x.y, x.x
  180.    mov   bx, y
  181.    fld   QWORD PTR [bx]       ; y.x, x.y, x.x
  182.    fld   QWORD PTR [bx+8]     ; y.y, y.x, x.y, x.x
  183.    fld   st                   ; y.y, y.y, y.x, x.y, x.x
  184.    fmul  st, st               ; y.y*y.y, y.y, y.x, x.y, x.x
  185.    fld   st(2)                ; y.x, y.y*y.y, y.y, y.x, x.y, x.x
  186.    fmul  st, st               ; y.x*y.x, y.y*y.y, y.y, y.x, x.y, x.x
  187.    fadd                       ; mod, y.y, y.x, x.y, x.x
  188.  
  189.    ftst                       ; test whether mod is (0,0)
  190.    fstsw Status
  191.    mov   ax, Status
  192.    and   ah, 01000101b
  193.    cmp   ah, 01000000b
  194.    jne   NotZero
  195.  
  196.    fstp  st
  197.    fstp  st
  198.    fstp  st
  199.    fstp  st
  200.    fstp  st
  201.  
  202.    fld   infinity
  203.    fld   st
  204.    mov   bx, z
  205.    fstp  QWORD PTR [bx]
  206.    fstp  QWORD PTR [bx+8]
  207.    inc   DivideOverflow
  208.    jmp   ExitDiv
  209.  
  210. NotZero:
  211.    fdiv  st(1), st            ; mod, y.y=y.y/mod, y.x, x.y, x.x
  212.    fdivp st(2), st            ; y.y, y.x=y.x/mod, x.y, x.x
  213.    mov   bx, z
  214.    fld   st                   ; y.y, y.y, y.x, x.y, x.x
  215.    fmul  st, st(3)            ; y.y*x.y, y.y. y.x, x.y, x.x
  216.    fld   st(2)                ; y.x, y.y*x.y, y.y, y.x, x.y, x.x
  217.    fmul  st, st(5)            ; y.x*x.x, y.y*x.y, y.y, y.x, x.y, x.x
  218.    fadd                       ; y.x*x.x - y.y*x.y, y.y, y.x, x.y, x.x
  219.    fstp  QWORD PTR [bx]       ; y.y, y.x, x.y, x.x
  220.    fmulp st(3), st            ; y.x, x.y, x.x*y.y
  221.    fmul                       ; y.x*x.y, x.x*y.y
  222.    fsubr                      ; y.x*x.y + x.x*y.y
  223.    fstp  QWORD PTR [bx+8]
  224.  
  225. ExitDiv:
  226.    ret
  227. FPUcplxdiv     ENDP
  228.  
  229.  
  230.  
  231. FPUcplxlog     PROC     x:word, z:word
  232. LOCAL Status:word, ImagZero:WORD
  233.    mov   bx, x
  234.  
  235.    mov   ax, WORD PTR [bx+8+6]
  236.    mov   ImagZero, ax
  237.    or    ax, WORD PTR [bx+6]
  238.    jnz   NotBothZero
  239.  
  240.    fldz
  241.    fldz
  242.    jmp   StoreZX
  243.  
  244. NotBothZero:
  245.    fld   QWORD PTR [bx+8]        ; x.y
  246.    fld   QWORD PTR [bx]          ; x.x, x.y
  247.    mov   bx, z
  248.    fldln2                        ; ln2, x.x, x.y
  249.    fdiv  _2_                     ; ln2/2, x.x, x.y
  250.    fld   st(2)                   ; x.y, ln2/2, x.x, x.y
  251.    fmul  st, st                  ; sqr(x.y), ln2/2, x.x, x.y
  252.    fld   st(2)                   ; x.x, sqr(x.y), ln2/2, x.x, x.y
  253.    fmul  st, st                  ; sqr(x.x), sqr(x.y), ln2/2, x.x, x.y
  254.    fadd                          ; mod, ln2/2, x.x, x.y
  255.    fyl2x                         ; z.x, x.x, x.y
  256.    fxch  st(2)                   ; x.y, x.x, z.x
  257.    fxch                          ; x.x, x.y, z.x
  258.    cmp   fpu, 387
  259.    jne   Restricted
  260.  
  261.    fpatan                        ; z.y, z.x
  262.    jmp   StoreZX
  263.  
  264. Restricted:
  265.    mov   bx, x
  266.    mov   dh, BYTE PTR [bx+7]
  267.    or    dh, dh
  268.    jns   ChkYSign
  269.  
  270.    fchs                          ; |x.x|, x.y, z.x
  271.  
  272. ChkYSign:
  273.    mov   dl, BYTE PTR [bx+8+7]
  274.    or    dl, dl
  275.    jns   ChkMagnitudes
  276.  
  277.    fxch                          ; x.y, |x.x|, z.x
  278.    fchs                          ; |x.y|, |x.x|, z.x
  279.    fxch                          ; |x.x|, |x.y|, z.x
  280.  
  281. ChkMagnitudes:
  282.    fcom  st(1)                   ; x.x, x.y, z.x
  283.    fstsw Status                  ; x.x, x.y, z.x
  284.    test  Status, 4500h
  285.    jz    XisGTY
  286.  
  287.    test  Status, 4000h
  288.    jz    XneY
  289.  
  290.    fstp  st                      ; x.y, z.x
  291.    fstp  st                      ; z.x
  292.    fldpi                         ; Pi, z.x
  293.    fdiv  _4_                     ; Pi/4, z.x
  294.    jmp   ChkSignZ
  295.  
  296. XneY:
  297.    fxch                          ; x.y, x.x, z.x
  298.    fpatan                        ; Pi/2 - Angle, z.x
  299.    fldpi                         ; Pi, Pi/2 - Angle, z.x
  300.    fdiv  _2_                     ; Pi/2, Pi/2 - Angle, z.x
  301.    fsubr                         ; Angle, z.x
  302.    jmp   ChkSignZ
  303.  
  304. XisGTY:
  305.    fpatan                        ; Pi-Angle or Angle+Pi, z.x
  306.  
  307. ChkSignZ:
  308.    or    dh, dh
  309.    js    NegX
  310.  
  311.    or    dl, dl
  312.    jns   StoreZX
  313.  
  314.    fchs
  315.    jmp   StoreZX
  316.  
  317. NegX:
  318.    or    dl, dl
  319.    js    QuadIII
  320.  
  321.    fldpi                        ; Pi, Pi-Angle, z.x
  322.    fsubr                        ; Angle, z.x
  323.    jmp   StoreZX
  324.  
  325. QuadIII:
  326.    fldpi                        ; Pi, Angle+Pi, z.x
  327.    fsub                         ; Angle,  z.x
  328.  
  329. StoreZX:
  330.    mov   bx, z
  331.    fstp  QWORD PTR [bx+8]       ; z.x
  332.    fstp  QWORD PTR [bx]         ; <empty>
  333.    ret
  334. FPUcplxlog     ENDP
  335.  
  336. FPUsinhcosh    PROC     x:word, sinh:word, cosh:word
  337. LOCAL Control:word
  338.    fstcw Control
  339.    push  Control                       ; Save control word on the stack
  340.    or    Control, 0000110000000000b
  341.    fldcw Control                       ; Set control to round towards zero
  342.  
  343.    mov   Sign, 0              ; Assume the sign is positive
  344.    mov   bx, x
  345.  
  346.    fldln2                     ; ln(2)
  347.    fdivr QWORD PTR [bx]       ; x/ln(2)
  348.  
  349.    cmp   BYTE PTR [bx+7], 0
  350.    jns   DuplicateX
  351.  
  352.    fchs                       ; x = |x|
  353.  
  354. DuplicateX:
  355.    fld   st                   ; x/ln(2), x/ln(2)
  356.    frndint                    ; int = integer(|x|/ln(2)), x/ln(2)
  357.    fxch                       ; x/ln(2), int
  358.    fsub  st, st(1)            ; rem < 1.0, int
  359.    fmul  PointFive            ; rem/2 < 0.5, int
  360.       ; CAE 7Dec92 changed above from divide by 2 to multiply by .5
  361.    f2xm1                      ; (2**rem/2)-1, int
  362.    fadd  _1_                  ; 2**rem/2, int
  363.    fmul  st, st               ; 2**rem, int
  364.    fscale                     ; e**|x|, int
  365.    fstp  st(1)                ; e**|x|
  366.  
  367.    cmp   BYTE PTR [bx+7], 0
  368.    jns   ExitFexp
  369.  
  370.    fdivr _1_                  ; e**x      
  371.  
  372. ExitFexp:
  373.    fld   st                   ; e**x, e**x
  374.    fdivr PointFive            ; e**-x/2, e**x
  375.    fld   st                   ; e**-x/2, e**-x/2, e**x
  376.    fxch  st(2)                ; e**x, e**-x/2, e**-x/2
  377.    fmul  PointFive            ; e**x/2,  e**-x/2, e**-x/2
  378.       ; CAE 7Dec92 changed above from divide by 2 to multiply by .5
  379.    fadd  st(2), st            ; e**x/2,  e**-x/2, cosh(x)
  380.    fsubr                      ; sinh(x), cosh(x)
  381.  
  382.    mov   bx, sinh             ; sinh, cosh
  383.    fstp  QWORD PTR [bx]       ; cosh
  384.    mov   bx, cosh
  385.    fstp  QWORD PTR [bx]       ; <empty>
  386.  
  387.    pop   Control
  388.    fldcw Control              ; Restore control word
  389.    ret
  390. FPUsinhcosh    ENDP
  391.  
  392.  
  393. FPUsincos  PROC  x:word, sinx:word, cosx:word
  394. LOCAL Status:word
  395.    mov   bx, x
  396.    fld   QWORD PTR [bx]       ; x
  397.  
  398.    cmp   fpu, 387
  399.    jne   Use387FPUsincos
  400.  
  401.    call  _Loaded387sincos     ; cos(x), sin(x)
  402.    mov   bx, cosx
  403.    fstp  QWORD PTR [bx]       ; sin(x)
  404.    mov   bx, sinx
  405.    fstp  QWORD PTR [bx]       ; <empty>
  406.    ret
  407.  
  408. Use387FPUsincos:
  409.  
  410.    sub   sp, 8                ; save 'x' on the CPU stack
  411.    mov   bx, sp
  412.    fstp  QWORD PTR [bx]       ; FPU stack:  <empty>
  413.  
  414.    call  cos
  415.  
  416.    add   sp, 8                ; take 'cos(x)' off the CPU stack
  417.    mov   bx, ax
  418.    cmp   compiled_by_turboc,0
  419.    jne   turbo_c1
  420.  
  421.    fld   QWORD PTR [bx]       ; FPU stack:  cos(x)
  422.  
  423. turbo_c1:
  424.    fld   st                   ; FPU stack:  cos(x), cos(x)
  425.    fmul  st, st               ; cos(x)**2, cos(x)
  426.    fsubr _1_                  ; sin(x)**2, cos(x)
  427.    fsqrt                      ; +/-sin(x), cos(x)
  428.  
  429.    mov   bx, x
  430.    fld   QWORD PTR [bx]       ; x, +/-sin(x), cos(x)
  431.    fldpi                      ; Pi, x, +/-sin(x), cos(x)
  432.    fadd  st, st               ; 2Pi, x, +/-sin(x), cos(x)
  433.    fxch                       ; |x|, 2Pi, +/-sin(x), cos(x)
  434.    fprem                      ; Angle, 2Pi, +/-sin(x), cos(x)
  435.    fstp  st(1)                ; Angle, +/-sin(x), cos(x)
  436.    fldpi                      ; Pi, Angle, +/-sin(x), cos(x)
  437.  
  438.    cmp   BYTE PTR [bx+7], 0
  439.    jns   SignAlignedPi
  440.  
  441.    fchs                       ; -Pi, Angle, +/-sin(x), cos(x)
  442.  
  443. SignAlignedPi:
  444.    fcompp                     ; +/-sin(x), cos(x)
  445.    fstsw Status               ; +/-sin(x), cos(x)
  446.  
  447.    mov   ax, Status
  448.    and   ah, 1
  449.    jz    StoreSinCos          ; Angle <= Pi
  450.  
  451.    fchs                       ; sin(x), cos(x)
  452.  
  453. StoreSinCos:
  454.    mov   bx, sinx
  455.    fstp  QWORD PTR [bx]       ; cos(x)
  456.    mov   bx, cosx
  457.    fstp  QWORD PTR [bx]       ; <empty>
  458.    ret
  459. FPUsincos   ENDP
  460.  
  461.  
  462. PUBLIC r16Mul
  463. r16Mul     PROC    uses si di, x1:word, x2:word, y1:word, y2:word
  464.       mov   si, x1
  465.       mov   bx, x2
  466.       mov   di, y1
  467.       mov   cx, y2
  468.  
  469.       xor   ax, ax
  470.       shl   bx, 1
  471.       jz    Exitr16Mult          ; Destination is zero
  472.  
  473.       rcl   ah, 1
  474.       shl   cx, 1
  475.       jnz   Chkr16Exp
  476.       xor   bx, bx               ; Source is zero
  477.       xor   si, si
  478.       jmp   Exitr16Mult
  479.  
  480.    Chkr16Exp:
  481.       rcl   al, 1
  482.       xor   ah, al               ; Resulting sign in ah
  483.       stc                        ; Put 'one' bit back into number
  484.       rcr   bl, 1
  485.       stc
  486.       rcr   cl, 1
  487.  
  488.       sub   ch, 127              ; Determine resulting exponent
  489.       add   bh, ch
  490.       mov   al, bh
  491.       mov   fake_es, ax          ; es has the resulting exponent and sign
  492.  
  493.       mov   ax, di
  494.       mov   al, ah
  495.       mov   ah, cl
  496.  
  497.       mov   dx, si
  498.       mov   dl, dh
  499.       mov   dh, bl
  500.  
  501.       mul   dx
  502.       mov   cx, fake_es
  503.  
  504.       shl   ax, 1
  505.       rcl   dx, 1
  506.       jnc   Remr16MulOneBit      ; 'One' bit is the next bit over
  507.  
  508.       inc   cl                   ; 'One' bit removed with previous shift
  509.       jmp   Afterr16MulNorm
  510.  
  511.    Remr16MulOneBit:
  512.       shl   ax, 1
  513.       rcl   dx, 1
  514.  
  515.    Afterr16MulNorm:
  516.       mov   bl, dh               ; Perform remaining 8 bit shift
  517.       mov   dh, dl
  518.       mov   dl, ah
  519.       mov   si, dx
  520.       mov   bh, cl               ; Put in the exponent
  521.       rcr   ch, 1                ; Get the sign
  522.       rcr   bx, 1                ; Normalize the result
  523.       rcr   si, 1
  524.    Exitr16Mult:
  525.       mov   ax, si
  526.       mov   dx, bx
  527.       ret
  528. r16Mul      ENDP
  529.  
  530.  
  531. PUBLIC RegFloat2Fg
  532. RegFloat2Fg     PROC    x1:word, x2:word, Fudge:word
  533.       mov   ax, WORD PTR x1
  534.       mov   dx, WORD PTR x2
  535.       mov   bx, ax
  536.       or    bx, dx
  537.       jz    ExitRegFloat2Fg
  538.  
  539.       xor   bx, bx
  540.       mov   cx, bx
  541.  
  542.       shl   ax, 1
  543.       rcl   dx, 1
  544.       rcl   bx, 1                   ; bx contains the sign
  545.  
  546.       xchg  cl, dh                  ; cx contains the exponent
  547.  
  548.       stc                           ; Put in the One bit
  549.       rcr   dl, 1
  550.       rcr   ax, 1
  551.  
  552.       sub   cx, 127 + 23
  553.       add   cx, Fudge
  554.       jz    ChkFgSign
  555.       jns   ShiftFgLeft
  556.  
  557.       neg   cx
  558.    ShiftFgRight:
  559.       shr   dx, 1
  560.       rcr   ax, 1
  561.       loop  ShiftFgRight
  562.       jmp   ChkFgSign
  563.  
  564.    ShiftFgLeft:
  565.       shl   ax, 1
  566.       rcl   dx, 1
  567.       loop  ShiftFgLeft
  568.  
  569.    ChkFgSign:
  570.       or    bx, bx
  571.       jz    ExitRegFloat2Fg
  572.  
  573.       not   ax
  574.       not   dx
  575.       add   ax, 1
  576.       adc   dx, 0
  577.  
  578.    ExitRegFloat2Fg:
  579.       ret
  580. RegFloat2Fg    ENDP
  581.  
  582. PUBLIC RegFg2Float
  583. RegFg2Float     PROC   x1:word, x2:word, FudgeFact:byte
  584.       mov   ax, x1
  585.       mov   dx, x2
  586.  
  587.       mov   cx, ax
  588.       or    cx, dx
  589.       jz    ExitFudgedToRegFloat
  590.  
  591.       mov   ch, 127 + 32
  592.       sub   ch, FudgeFact
  593.       xor   cl, cl
  594.       shl   ax, 1       ; Get the sign bit
  595.       rcl   dx, 1
  596.       jnc   FindOneBit
  597.  
  598.       inc   cl          ; Fudged < 0, convert to postive
  599.       not   ax
  600.       not   dx
  601.       add   ax, 1
  602.       adc   dx, 0
  603.  
  604.    FindOneBit:
  605.       shl   ax, 1
  606.       rcl   dx, 1
  607.       dec   ch
  608.       jnc   FindOneBit
  609.       dec   ch
  610.  
  611.       mov   al, ah
  612.       mov   ah, dl
  613.       mov   dl, dh
  614.       mov   dh, ch
  615.  
  616.       shr   cl, 1       ; Put sign bit in
  617.       rcr   dx, 1
  618.       rcr   ax, 1
  619.  
  620.    ExitFudgedToRegFloat:
  621.       ret
  622. RegFg2Float      ENDP
  623.  
  624.  
  625. PUBLIC ExpFudged
  626. ExpFudged      PROC    uses si, x_low:word, x_high:word, Fudge:word
  627. LOCAL exp:WORD
  628.       xor   ax, ax
  629.       mov   WORD PTR Ans, ax
  630.       mov   WORD PTR Ans + 2, ax
  631.       mov   ax, WORD PTR x_low
  632.       mov   dx, WORD PTR x_high
  633.       or    dx, dx
  634.       js    NegativeExp
  635.  
  636.       div   Ln2Fg16
  637.       mov   exp, ax
  638.       or    dx, dx
  639.       jz    Raiseexp
  640.  
  641.       mov   ax, dx
  642.       mov   si, dx
  643.       mov   bx, 1
  644.  
  645.    PosExpLoop:
  646.       add   WORD PTR Ans, ax
  647.       adc   WORD PTR Ans + 2, 0
  648.       inc   bx
  649.       mul   si
  650.       mov   ax, dx
  651.       xor   dx, dx
  652.       div   bx
  653.       or    ax, ax
  654.       jnz   PosExpLoop
  655.  
  656.    Raiseexp:
  657.       inc   WORD PTR Ans + 2
  658.       mov   ax, WORD PTR Ans
  659.       mov   dx, WORD PTR Ans + 2
  660.       mov   cx, -16
  661.       add   cx, Fudge
  662.       add   cx, exp
  663.       or    cx, cx
  664.       jz    ExitExpFudged
  665.       jns   LeftShift
  666.       neg   cx
  667.  
  668.    RightShift:
  669.       shr   dx, 1
  670.       rcr   ax, 1
  671.       loop  RightShift
  672.       jmp   ExitExpFudged
  673.  
  674.    NegativeExp:
  675.       not   ax
  676.       not   dx
  677.       add   ax, 1
  678.       adc   dx, 0
  679.       div   Ln2Fg16
  680.       neg   ax
  681.       mov   exp, ax
  682.  
  683.       or    dx, dx
  684.       jz    Raiseexp
  685.  
  686.       mov   ax, dx
  687.       mov   si, dx
  688.       mov   bx, 1
  689.  
  690.    NegExpLoop:
  691.       sub   WORD PTR Ans, ax
  692.       sbb   WORD PTR Ans + 2, 0
  693.       inc   bx
  694.       mul   si
  695.       mov   ax, dx
  696.       xor   dx, dx
  697.       div   bx
  698.       or    ax, ax
  699.       jz    Raiseexp
  700.  
  701.       add   WORD PTR Ans, ax
  702.       adc   WORD PTR Ans + 2, 0
  703.       inc   bx
  704.       mul   si
  705.       mov   ax, dx
  706.       xor   dx, dx
  707.       div   bx
  708.       or    ax, ax
  709.       jnz   NegExpLoop
  710.       jmp   Raiseexp
  711.  
  712.    LeftShift:
  713.       shl   ax, 1
  714.       rcl   dx, 1
  715.       loop  LeftShift
  716.  
  717.    ExitExpFudged:
  718.       ret
  719. ExpFudged      ENDP
  720.  
  721.  
  722.  
  723. PUBLIC   LogFudged
  724. LogFudged      PROC    uses si di, x_low:word, x_high:word, Fudge:word
  725. LOCAL exp:WORD
  726.       xor   bx, bx
  727.       mov   cx, 16
  728.       sub   cx, Fudge
  729.       mov   ax, x_low
  730.       mov   dx, x_high
  731.  
  732.       or    dx, dx
  733.       jz    ChkLowWord
  734.  
  735.    Incexp:
  736.       shr   dx, 1
  737.       jz    DetermineOper
  738.       rcr   ax, 1
  739.       inc   cx
  740.       jmp   Incexp
  741.  
  742.    ChkLowWord:
  743.       or    ax, ax
  744.       jnz   Decexp
  745.       jmp   ExitLogFudged
  746.  
  747.    Decexp:
  748.       dec   cx                      ; Determine power of two
  749.       shl   ax, 1
  750.       jnc   Decexp
  751.  
  752.    DetermineOper:
  753.       mov   exp, cx
  754.       mov   si, ax                  ; si =: x + 1
  755.       shr   si, 1
  756.       stc
  757.       rcr   si, 1
  758.       mov   dx, ax
  759.       xor   ax, ax
  760.       shr   dx, 1
  761.       rcr   ax, 1
  762.       shr   dx, 1
  763.       rcr   ax, 1                   ; dx:ax = x - 1
  764.       div   si
  765.  
  766.       mov   bx, ax                  ; ax, Fudged 16, max of 0.3333333
  767.       shl   ax, 1                   ; x = (x - 1) / (x + 1), Fudged 16
  768.       mul   ax
  769.       shl   ax, 1
  770.       rcl   dx, 1
  771.       mov   ax, dx                  ; dx:ax, Fudged 35, max = 0.1111111
  772.       mov   si, ax                  ; si = (ax * ax), Fudged 19
  773.  
  774.       mov   ax, bx
  775.    ; bx is the accumulator, First term is x
  776.       mul   si                      ; dx:ax, Fudged 35, max of 0.037037
  777.       mov   fake_es, dx             ; Save high word, Fudged (35 - 16) = 19
  778.       mov   di, 0c000h              ; di, 3 Fudged 14
  779.       div   di                      ; ax, Fudged (36 - 14) = 21
  780.       or    ax, ax
  781.       jz    Addexp
  782.  
  783.       mov   cl, 5
  784.       shr   ax, cl
  785.       add   bx, ax                  ; bx, max of 0.345679
  786.    ; x = x + x**3/3
  787.  
  788.       mov   ax, fake_es             ; ax, Fudged 19
  789.       mul   si                      ; dx:ax, Fudged 38, max of 0.004115
  790.       mov   fake_es, dx             ; Save high word, Fudged (38 - 16) = 22
  791.       mov   di, 0a000h              ; di, 5 Fudged 13
  792.       div   di                      ; ax, Fudged (38 - 13) = 25
  793.       or    ax, ax
  794.       jz    Addexp
  795.  
  796.       mov   cl, 9
  797.       shr   ax, cl
  798.       add   bx, ax
  799.    ; x = x + x**3/3 + x**5/5
  800.  
  801.       mov   ax, fake_es             ; ax, Fudged 22
  802.       mul   si                      ; dx:ax, Fudged 41, max of 0.0004572
  803.       mov   di, 0e000h              ; di, 7 Fudged 13
  804.       div   di                      ; ax, Fudged (41 - 13) = 28
  805.       mov   cl, 12
  806.       shr   ax, cl
  807.       add   bx, ax
  808.  
  809.    Addexp:
  810.       shl   bx, 1                   ; bx *= 2, Fudged 16, max of 0.693147
  811.    ; x = 2 * (x + x**3/3 + x**5/5 + x**7/7)
  812.       mov   cx, exp
  813.       mov   ax, Ln2Fg16            ; Answer += exp * Ln2Fg16
  814.       or    cx, cx
  815.       js    SubFromAns
  816.  
  817.       mul   cx
  818.       add   ax, bx
  819.       adc   dx, 0
  820.       jmp   ExitLogFudged
  821.  
  822.    SubFromAns:
  823.       neg   cx
  824.       mul   cx
  825.       xor   cx, cx
  826.       xchg  cx, dx
  827.       xchg  bx, ax
  828.       sub   ax, bx
  829.       sbb   dx, cx
  830.  
  831.    ExitLogFudged:
  832.    ; x = 2 * (x + x**3/3 + x**5/5 + x**7/7) + (exp * Ln2Fg16)
  833.       ret
  834. LogFudged      ENDP
  835.  
  836.  
  837.  
  838.  
  839. PUBLIC LogFloat14
  840. LogFloat14     PROC     x1:word, x2:word
  841.       mov   ax, WORD PTR x1
  842.       mov   dx, WORD PTR x2
  843.       shl   ax, 1
  844.       rcl   dx, 1
  845.       xor   cx, cx
  846.       xchg  cl, dh
  847.  
  848.       stc
  849.         rcr   dl, 1
  850.         rcr   ax, 1
  851.  
  852.       sub   cx, 127 + 23
  853.       neg   cx
  854.       push  cx
  855.       push  dx
  856.       push  ax
  857.       call  LogFudged
  858.       add   sp, 6
  859.       ret
  860. LogFloat14     ENDP
  861.  
  862.  
  863. PUBLIC RegSftFloat
  864. RegSftFloat     PROC   x1:word, x2:word, Shift:byte
  865.       mov   ax, x1
  866.       mov   dx, x2
  867.  
  868.       shl   dx, 1
  869.       rcl   cl, 1
  870.  
  871.       add   dh, Shift
  872.  
  873.       shr   cl, 1
  874.       rcr   dx, 1
  875.  
  876.       ret
  877. RegSftFloat      ENDP
  878.  
  879.  
  880.  
  881.  
  882. PUBLIC RegDivFloat
  883. RegDivFloat     PROC  uses si di, x1:word, x2:word, y1:word, y2:word
  884.       mov   si, x1
  885.       mov   bx, x2
  886.       mov   di, y1
  887.       mov   cx, y2
  888.  
  889.       xor   ax, ax
  890.       shl   bx, 1
  891.       jnz   ChkOtherOp
  892.       jmp   ExitRegDiv           ; Destination is zero
  893.  
  894.    ChkOtherOp:
  895.       rcl   ah, 1
  896.       shl   cx, 1
  897.       jnz   ChkDivExp
  898.       xor   bx, bx               ; Source is zero
  899.       xor   si, si
  900.       jmp   ExitRegDiv
  901.  
  902.    ChkDivExp:
  903.       rcl   al, 1
  904.       xor   ah, al               ; Resulting sign in ah
  905.       stc                        ; Put 'one' bit back into number
  906.       rcr   bl, 1
  907.       stc
  908.       rcr   cl, 1
  909.  
  910.       sub   ch, 127              ; Determine resulting exponent
  911.       sub   bh, ch
  912.       mov   al, bh
  913.       mov   fake_es, ax          ; es has the resulting exponent and sign
  914.  
  915.       mov   ax, si               ; 8 bit shift, bx:si moved to dx:ax
  916.       mov   dh, bl
  917.       mov   dl, ah
  918.       mov   ah, al
  919.       xor   al, al
  920.  
  921.       mov   bh, cl               ; 8 bit shift, cx:di moved to bx:cx
  922.       mov   cx, di
  923.       mov   bl, ch
  924.       mov   ch, cl
  925.       xor   cl, cl
  926.  
  927.       shr   dx, 1
  928.       rcr   ax, 1
  929.  
  930.       div   bx
  931.       mov   si, dx               ; Save (and shift) remainder
  932.       mov   dx, cx               ; Save the quess
  933.       mov   cx, ax
  934.       mul   dx                   ; Mult quess times low word
  935.       xor   di, di
  936.       sub   di, ax               ; Determine remainder
  937.       sbb   si, dx
  938.       mov   ax, di
  939.       mov   dx, si
  940.       jc    RemainderNeg
  941.  
  942.       xor   di, di
  943.       jmp   GetNextDigit
  944.  
  945.    RemainderNeg:
  946.       mov   di, 1                ; Flag digit as negative
  947.       not   ax                   ; Convert remainder to positive
  948.       not   dx
  949.       add   ax, 1
  950.       adc   dx, 0
  951.  
  952.    GetNextDigit:
  953.       shr   dx, 1
  954.       rcr   ax, 1
  955.       div   bx
  956.       xor   bx, bx
  957.       shl   dx, 1
  958.       rcl   ax, 1
  959.       rcl   bl, 1                ; Save high bit
  960.  
  961.       mov   dx, cx               ; Retrieve first digit
  962.       or    di, di
  963.       jz    RemoveDivOneBit
  964.  
  965.       neg   ax                   ; Digit was negative
  966.       neg   bx
  967.       dec   dx
  968.  
  969.    RemoveDivOneBit:
  970.       add   dx, bx
  971.       mov   cx, fake_es
  972.       shl   ax, 1
  973.       rcl   dx, 1
  974.       jc    AfterDivNorm
  975.  
  976.       dec   cl
  977.       shl   ax, 1
  978.       rcl   dx, 1
  979.  
  980.    AfterDivNorm:
  981.       mov   bl, dh               ; Perform remaining 8 bit shift
  982.       mov   dh, dl
  983.       mov   dl, ah
  984.       mov   si, dx
  985.       mov   bh, cl               ; Put in the exponent
  986.       shr   ch, 1                ; Get the sign
  987.       rcr   bx, 1                ; Normalize the result
  988.       rcr   si, 1
  989.  
  990.    ExitRegDiv:
  991.       mov   ax, si
  992.       mov   dx, bx
  993.       ret
  994. RegDivFloat      ENDP
  995.  
  996.  
  997.  
  998. Term        equ      <ax>
  999. Num         equ      <bx>
  1000. Factorial   equ      <cx>
  1001. Sin         equ      <si>
  1002. Cos         equ      <di>
  1003. e           equ      <si>
  1004. Inve        equ      <di>
  1005.          
  1006. SinCos086   PROC     uses si di, LoNum:WORD, HiNum:WORD, SinAddr:WORD, \
  1007.                                 CosAddr:WORD
  1008.    mov   ax, LoNum 
  1009.    mov   dx, HiNum
  1010.    
  1011.    xor   cx, cx
  1012. ;   mov   SinNeg, cx    ;CJLT - not needed now
  1013. ;   mov   CosNeg, cx     ;CJLT - not needed now
  1014.    mov   a, cx        ;CJLT - Not needed by the original code, but it
  1015.                ;       is now!
  1016.    or    dx, dx
  1017.    jns   AnglePositive
  1018.    
  1019.    not   ax
  1020.    not   dx
  1021.    add   ax, 1
  1022.    adc   dx, cx        ;conveniently zero
  1023.    mov   a,8        ;a now has 4 bits: Sign+quadrant+octant
  1024.       
  1025. AnglePositive:
  1026.    mov   si, ax
  1027.    mov   di, dx
  1028.    mul   WORD PTR InvPiFg33
  1029.    mov   bx, dx
  1030.    mov   ax, di
  1031.    mul   WORD PTR InvPiFg33
  1032.    add   bx, ax
  1033.    adc   cx, dx
  1034.    mov   ax, si
  1035.    mul   WORD PTR InvPiFg33 + 2
  1036.    add   bx, ax
  1037.    adc   cx, dx
  1038.    mov   ax, di
  1039.    mul   WORD PTR InvPiFg33 + 2
  1040.    add   ax, cx
  1041.    adc   dx, 0
  1042.  
  1043.    and   dx, 3    ;Remove multiples of 2 pi
  1044.    shl   dx, 1  ;move over to make space for octant number
  1045. ;
  1046. ;CJLT - new code to reduce angle to:  0 <= angle <= 45
  1047. ;
  1048.    or    ax, ax
  1049.    jns   Lessthan45
  1050.    inc   dx    ;octant number
  1051.    not   ax    ;angle=90-angle if it was >45 degrees
  1052. Lessthan45:
  1053.    add   a,  dx    ;remember octant and Sign in a
  1054.    mov   Num, ax ;ax=Term, now
  1055. ;
  1056. ; We do the Taylor series with all numbers scaled by pi/2
  1057. ; so InvPiFg33 represents one. Truly.
  1058. ;
  1059.    mov   Factorial, WORD PTR InvPiFg33 + 2
  1060.    mov   one, Factorial
  1061.    mov   Cos, Factorial          ; Cos = 1
  1062.    mov   Sin, Num                ; Sin = Num
  1063.       
  1064. LoopIntSinCos:
  1065.    TaylorTerm                    ; ax=Term = Num * (Num/2) * (Num/3) * . . .
  1066.    sub   Cos, Term               ; Cos = 1 - Num*(Num/2) + (Num**4)/4! - . . .
  1067.    cmp   Term, WORD PTR TrigLimit
  1068.    jbe   SHORT ExitIntSinCos
  1069.  
  1070.    TaylorTerm
  1071.    sub   Sin, Term               ; Sin = Num - Num*(Num/2)*(Num/3) + (Num**5)/5! - . . .
  1072.    cmp   Term, WORD PTR TrigLimit
  1073.    jbe   SHORT ExitIntSinCos
  1074.  
  1075.    TaylorTerm
  1076.    add   Cos, Term
  1077.    cmp   Term, WORD PTR TrigLimit
  1078.    jbe   SHORT ExitIntSinCos
  1079.  
  1080.    TaylorTerm                    ; Term = Num * (x/2) * (x/3) * . . .
  1081.    add   Sin, Term
  1082.    cmp   Term, WORD PTR TrigLimit
  1083.    jnbe  LoopIntSinCos
  1084.       
  1085. ExitIntSinCos:
  1086.    xor   ax, ax
  1087.    mov   cx, ax
  1088.    cmp   Cos, WORD PTR InvPiFg33 + 2
  1089.    jb    CosDivide               ; Cos < 1.0
  1090.       
  1091.    inc   cx                      ; Cos == 1.0
  1092.    jmp   StoreCos
  1093.  
  1094. CosDivide:
  1095.    mov   dx, Cos
  1096.    div   WORD PTR InvPiFg33 + 2
  1097.       
  1098. StoreCos:
  1099.    mov   Cos, ax                 ; cx:Cos
  1100.  
  1101.    xor   ax, ax
  1102.    mov   bx, ax
  1103.    cmp   Sin, WORD PTR InvPiFg33 + 2
  1104.    jb    SinDivide               ; Sin < 1.0
  1105.    
  1106.    inc   bx                      ; Sin == 1.0
  1107.    jmp   StoreSin
  1108.       
  1109. SinDivide:
  1110.    mov   dx, Sin
  1111.    div   WORD PTR InvPiFg33 + 2
  1112.       
  1113. StoreSin:
  1114.    mov   Sin, ax                 ; bx:Sin
  1115. ; CJLT again. New tests are needed to correct signs and exchange sin/cos values
  1116.    mov   ax, a
  1117.    inc   al    ;forces bit 1 to xor of previous bits 0 and 1
  1118.    test  al, 2
  1119.    jz    ChkNegCos
  1120.  
  1121.    xchg  cx, bx
  1122.    xchg  Sin, Cos
  1123. ;   mov   ax, SinNeg    commented out by CJLT. This was a bug.
  1124. ;   xchg  ax, CosNeg
  1125. ;   mov   CosNeg, ax    and this was meant to be  mov  SinNeg,ax
  1126.  
  1127. ChkNegCos:
  1128.    inc   al    ;forces bit 2 to xor of original bits 1 and 2
  1129.    test  al, 4 
  1130.    jz    ChkNegSin
  1131.  
  1132.    not   Cos    ;negate cos if quadrant 2 or 3
  1133.    not   cx
  1134.    add   Cos, 1
  1135.    adc   cx, 0
  1136.  
  1137. ChkNegSin:
  1138.    inc   al
  1139.    inc   al    ;forces bit 3 to xor of original bits 2 and 3
  1140.    test  al, 8
  1141.    jz    CorrectQuad
  1142.  
  1143.    not   Sin
  1144.    not   bx
  1145.    add   Sin, 1
  1146.    adc   bx, 0
  1147.  
  1148. CorrectQuad:
  1149.  
  1150. CosPolarized:     
  1151.    mov   dx, bx
  1152.    mov   bx, CosAddr
  1153.    mov   WORD PTR [bx], Cos
  1154.    mov   WORD PTR [bx+2], cx
  1155.  
  1156. SinPolarized:
  1157.    mov   bx, SinAddr
  1158.    mov   WORD PTR [bx], Sin
  1159.    mov   WORD PTR [bx+2], dx
  1160.    ret
  1161. SinCos086      ENDP
  1162.       
  1163.  
  1164. PUBLIC ArcTan086
  1165. ;
  1166. ;Used to calculate logs of complex numbers, since log(R,theta)=log(R)+i.theta
  1167. ; in polar coordinates.
  1168. ;
  1169. ;In C we need the prototype:
  1170. ;long ArcTan086(long, long)
  1171.  
  1172. ;The parameters are x and y, the returned value arctan(y/x) in the range 0..2pi.
  1173. ArcTan086     PROC  uses si di, x1:word, x2:word, y1:word, y2:word
  1174.       xor   ax, ax ;ax=0
  1175.       mov   si, x1 ;Lo
  1176.       mov   bx, x2 ;Hi
  1177.       or    bx, bx ;Sign set ?
  1178.       jns   xplus
  1179.       inc   ax
  1180.       not   si
  1181.       not   bx
  1182.       add   bx, 1
  1183.       adc   si, 0    ;We have abs(x) now
  1184.  xplus:           
  1185.       mov   di, y1 ;Lo
  1186.       mov   cx, y2 ;Hi
  1187.       or    cx, cx ;Sign set?
  1188.       jns   yplus
  1189.       inc   ax
  1190.       inc   ax       ;Set bit 1 of ax
  1191.       shl   ax, 1  ; and shift it all left 
  1192.       not  di
  1193.       not   cx
  1194.       add   di, 1
  1195.       adc   cx, 0  ;We have abs(y) now
  1196. yplus:
  1197.       cmp   bx, cx    ;y>x?
  1198.       jl    no_xchg
  1199.       jg    xchg_xy
  1200.       cmp   si, di    ;Hi words zero. What about lo words?
  1201.       jle   no_xchg
  1202. xchg_xy:        ;Exchange them
  1203.       inc   ax        ;Flag the exchange
  1204.       xchg  si, di
  1205.       xchg  bx, cx
  1206. no_xchg:
  1207.       mov   SinNeg, ax    ;Remember signs of x and y, and whether exchanged
  1208.       or    cx, cx    ;y<x now, but is y>=1.0 ?      
  1209.       jz    ynorm    ;no, so continue
  1210. normy:            ;yes, normalise by reducing y to 16 bits max.
  1211.       rcr   cx, 1    ; (We don't really lose any precision)
  1212.       rcr   di, 1
  1213.       clc
  1214.       rcr   bx, 1
  1215.       rcr   si, 1
  1216.       or    cx, cx
  1217.       jnz   normy      
  1218. ynorm:
  1219.  
  1220.       ret
  1221. ArcTan086    ENDP
  1222.  
  1223. ;
  1224. ;There now follows Chris Lusby Taylor's novel (such modesty!) method
  1225. ;of calculating exp(x). Originally, I was going to do it by decomposing x
  1226. ;into a sum of terms of the form:
  1227. ;
  1228. ;  th            -i
  1229. ; i   term=ln(1+2  )
  1230. ;                                     -i
  1231. ;If x>term[i] we subtract term[i] from x and multiply the answer by 1 + 2
  1232. ;
  1233. ;We can multiply by that factor by shifting and adding. Neat eh!
  1234. ;
  1235. ;Unfortunately, however, for 16 bit accuracy, that method turns out to be
  1236. ;slower than what follows, which is also novel. Here, I first divide x not
  1237. ;by ln(2) but ln(2)/16 so that we can look up an initial answer in a table of
  1238. ;2^(n/16). This leaves the remainder so small that the polynomial Taylor series
  1239. ;converges in just 1+x+x^2/2 which is calculated as 1+x(1+x/2).
  1240. ;
  1241. _e2xLT   PROC        ;argument in dx.ax (bitshift=16 is hard coded here)
  1242.    or    dx, dx
  1243.    jns   CalcExpLT
  1244.       
  1245.    not   ax        ; if negative, calc exp(abs(x))
  1246.    not   dx
  1247.    add   ax, 1
  1248.    adc   dx, 0
  1249.    
  1250. CalcExpLT:
  1251.          shl   ax, 1
  1252.       rcl   dx, 1
  1253.       shl   ax, 1
  1254.       rcl   dx, 1
  1255.       shl   ax, 1
  1256.       rcl   dx, 1        ;x is now in fg19 form for accuracy
  1257.                 ; so, relative to fg16, we have:
  1258.    div   Ln2Fg15        ; 8x==(a * Ln(2)/2) + Rem
  1259.                    ;  x=(a * Ln(2)/16) + Rem/8
  1260.                    ;so exp(x)=exp((a * Ln(2)/16) + Rem/8)
  1261.                 ;         =exp(a/16 * Ln(2)) * exp(Rem/8)
  1262.                 ;         =2^(a/16) * exp(Rem)
  1263.                 ;a mod 16 will give us an initial estimate 
  1264.    mov  cl, 4
  1265.    mov  di, ax            ;Remember original
  1266.    shr  ax, cl        
  1267.    mov   a, ax            ;a/16 will give us a bit shift
  1268.    mov  ax, di
  1269.    and  ax, 0000fh        ;a mod 16
  1270.    shl  ax, 1            ;used as an index into 2 byte per element Table
  1271.    mov  di, ax
  1272.    dec  cx            ;3, please
  1273.    add   dx, 4            ;to control rounding up/down
  1274.    shr   dx, cl            ;Num=Rem/8 (convert back to fg16)
  1275.                 ; 
  1276.    mov   ax, dx
  1277.    shr   ax, 1             ;Num/2  fg 16
  1278.    inc   ax            ;rounding control
  1279.    stc
  1280.    rcr   ax, 1            ;1+Num/2   fg15
  1281.    mul   dx            ;dx:ax=Num(1+Num/2) fg31, so dx alone is fg15
  1282.    shl   ax, 1            ;more rounding control
  1283.    adc   dx, 8000H         ;dx=1+Num(1+Num/2) fg15
  1284.    mov   ax, Table[di]        ;Get table entry fg15
  1285.    mul   dx            ;Only two multiplys used!  fg14, now (15+15-16)
  1286.    shl   ax, 1
  1287.    rcl   dx, 1            ;fg15
  1288.    mov   e,  dx 
  1289.    ret                          ; return e=e**x * (2**15), 1 < e**x < 2
  1290. _e2xLT   ENDP            ;        a= bitshift needed 
  1291.  
  1292.  
  1293.       
  1294. Exp086    PROC     uses si di, LoNum:WORD, HiNum:WORD
  1295.    mov   ax, LoNum 
  1296.    mov   dx, HiNum
  1297.    mov   TrigOverflow, 0
  1298.    
  1299.     call  _e2xLT    ;Use Chris Lusby Taylor's e2x
  1300.  
  1301.    cmp   a, 15        ;CJLT - was 16
  1302.    jae   Overflow
  1303.       
  1304. ;   cmp   expSign, 0
  1305. ;   jnz   NegNumber
  1306.    cmp   HiNum, 0    ;CJLT - we don't really need expSign
  1307.    jl   NegNumber
  1308.       
  1309.    mov   ax, e
  1310.    mov   dx, ax
  1311.    inc   a
  1312.    mov   cx, 16
  1313.    sub   cx, a
  1314.    shr   dx, cl
  1315.    mov   cx, a
  1316.    shl   ax, cl
  1317.    jmp   ExitExp086
  1318.       
  1319. Overflow:
  1320.    xor   ax, ax
  1321.    xor   dx, dx
  1322.    mov   TrigOverflow, 1
  1323.    jmp   ExitExp086
  1324.       
  1325. NegNumber:
  1326.    cmp   e, 8000h
  1327.    jne   DivideE
  1328.       
  1329.    mov   ax, e
  1330.    dec   a
  1331.    jmp   ShiftE
  1332.       
  1333. DivideE:
  1334.    xor   ax, ax
  1335.    mov   dx, ax
  1336.    stc
  1337.    rcr   dx, 1
  1338.    div   e
  1339.  
  1340. ShiftE:
  1341.    xor   dx, dx
  1342.    mov   cx, a
  1343.    shr   ax, cl
  1344.       
  1345. ExitExp086:
  1346.    ret
  1347. Exp086    ENDP
  1348.  
  1349.  
  1350.  
  1351. SinhCosh086    PROC     uses si di, LoNum:WORD, HiNum:WORD, SinhAddr:WORD, \
  1352.                                    CoshAddr:WORD
  1353.    mov   ax, LoNum
  1354.    mov   dx, HiNum
  1355.  
  1356.    call  _e2xLT        ;calculate exp(|x|) fg 15
  1357.             ;answer is e*2^a
  1358.    cmp   e, 8000h    ;e==1 ?    
  1359.    jne   InvertE        ; e > 1, so we can invert it.
  1360.  
  1361.    mov   dx, 1
  1362.    xor   ax, ax
  1363.    cmp   a, 0
  1364.    jne   Shiftone
  1365.  
  1366.    mov   e, ax
  1367.    mov   cx, ax
  1368.    jmp   ChkSinhSign
  1369.  
  1370. Shiftone:
  1371.    mov   cx, a
  1372.    shl   dx, cl
  1373.    dec   cx
  1374.    shr   e, cl
  1375.    shr   dx, 1
  1376.    shr   e, 1
  1377.    mov   cx, dx
  1378.    sub   ax, e
  1379.    sbb   dx, 0
  1380.    xchg  ax, e
  1381.    xchg  dx, cx
  1382.    jmp   ChkSinhSign
  1383.  
  1384. InvertE:
  1385.    xor   ax, ax               ; calc 1/e
  1386.    mov   dx, 8000h
  1387.    div   e
  1388.  
  1389.    mov   Inve, ax
  1390.  
  1391. ShiftE:
  1392.    mov   cx, a
  1393.    shr   Inve, cl
  1394.    inc   cl
  1395.    mov   dx, e
  1396.    shl   e, cl
  1397.    neg   cl
  1398.    add   cl, 16
  1399.    shr   dx, cl
  1400.    mov   cx, dx               ; cx:e == e**Exp
  1401.  
  1402.    mov   ax, e                ; dx:e == e**Exp
  1403.    add   ax, Inve
  1404.    adc   dx, 0
  1405.    shr   dx, 1
  1406.    rcr   ax, 1                ; cosh(Num) = (e**Exp + 1/e**Exp) / 2
  1407.  
  1408.    sub   e, Inve
  1409.    sbb   cx, 0
  1410.    sar   cx, 1
  1411.    rcr   e, 1
  1412.  
  1413. ChkSinhSign:
  1414.    or    HiNum, 0
  1415.    jns   StoreHyperbolics
  1416.  
  1417.    not   e
  1418.    not   cx
  1419.    add   e, 1
  1420.    adc   cx, 0
  1421.  
  1422. StoreHyperbolics:
  1423.    mov   bx, CoshAddr
  1424.    mov   WORD PTR [bx], ax
  1425.    mov   WORD PTR [bx+2], dx
  1426.  
  1427.    mov   bx, SinhAddr
  1428.    mov   WORD PTR [bx], e
  1429.    mov   WORD PTR [bx+2], cx
  1430.  
  1431.    ret
  1432. SinhCosh086    ENDP
  1433.  
  1434.  
  1435.  
  1436. Log086   PROC     uses si di, LoNum:WORD, HiNum:WORD, Fudge:WORD
  1437. LOCAL Exp:WORD, Accum:WORD, LoAns:WORD, HiAns:WORD
  1438. ;NOTE: CJLT's method does not need LoAns, HiAns, but he hasn't yet
  1439. ;taken them out
  1440.       xor   bx, bx
  1441.       mov   cx, Fudge
  1442.       mov   ax, LoNum
  1443.       mov   dx, HiNum
  1444.       mov   TrigOverflow, 0
  1445.  
  1446.       or    dx, dx
  1447.       js    Overflow
  1448.       jnz   ShiftUp
  1449.  
  1450.       or    ax, ax
  1451.       jnz   ShiftUp      
  1452.  
  1453.    Overflow:
  1454.       mov   TrigOverflow, 1
  1455.       jmp   ExitLog086
  1456.  
  1457.    ShiftUp:
  1458.       inc   cx                      ; cx = Exp
  1459.       shl   ax, 1
  1460.       rcl   dx, 1
  1461.       or    dx, dx
  1462.       jns   ShiftUp        ;shift until dx in fg15 form
  1463.  
  1464.       neg   cx
  1465.       add   cx, 31
  1466.       mov   Exp, cx
  1467. ;CJLT method starts here. We try to reduce x to make it very close to 1
  1468. ;store LoAns in bx for now (fg 16; initially 0)
  1469.       mov  cl,2        ;shift count
  1470. redoterm2:
  1471.       cmp  dx, 0AAABH    ;x > 4/3 ?
  1472.       jb  doterm3
  1473.       mov  ax, dx
  1474.       shr  ax, cl
  1475.       sub  dx, ax    ;x:=x(1-1/4)
  1476.       add  bx, 49a5H    ;ln(4/3) fg 15
  1477.       jmp  redoterm2
  1478.  doterm3:     
  1479.       inc  cl        ;count=3
  1480. redoterm3:
  1481.       cmp  dx, 9249H    ;x > 8/7 ?
  1482.       jb  doterm4
  1483.       mov  ax, dx
  1484.       shr  ax, cl
  1485.       sub  dx, ax    ;x:=x(1-1/8)
  1486.       add  bx, 222eH    ;ln(8/7)
  1487.       jmp  redoterm3
  1488.  doterm4:
  1489.       inc  cl        ;count=4
  1490.       cmp  dx, 8889H    ;x > 16/15 ?
  1491.       jb  skipterm4
  1492.       mov  ax, dx
  1493.       shr  ax, cl
  1494.       sub  dx, ax    ;x:=x(1-1/16)
  1495.       add  bx, 1085h    ;ln(16/15)
  1496. ;No need to retry term4 as error is acceptably low if we ignore it
  1497. skipterm4:
  1498. ;Now we have reduced x to the range 1 <=x <1.072
  1499. ;
  1500. ;Now we continue with the conventional series, but x is so small we
  1501. ;can ignore all terms except the first!
  1502. ;i.e.:
  1503. ;ln(x)=2(x-1)/(x+1)
  1504.       sub   dx, 8000h        ; dx= x-1, fg 15
  1505.       mov   cx, dx
  1506.       stc    
  1507.       rcr   cx, 1        ; cx = 1 + (x-1)/2   fg 15
  1508.                       ;   = 1+x            fg 14
  1509.       mov   ax,4000H        ;rounding control (Trust me)
  1510.       div   cx            ;ax=ln(x)
  1511.       add   bx, ax        ;so add it to the rest of the Ans. No carry
  1512.    MultExp:
  1513.       mov   cx, Exp
  1514.       mov   ax, Ln2Fg16
  1515.       or    cx, cx
  1516.       js    SubFromAns
  1517.  
  1518.       mul   cx                      ; cx = Exp * Lg2Fg16, fg 16
  1519.       add   ax, bx        ;add bx part of answer
  1520.       adc   dx, 0
  1521.       jmp   ExitLog086
  1522.  
  1523.    SubFromAns:
  1524.       inc   bx        ;Somewhat artificial, but we need to add 1 somewhere
  1525.       neg   cx
  1526.       mul   cx
  1527.    not   ax
  1528.       not   dx
  1529.       add   ax, bx
  1530.       adc   dx, 0
  1531.  
  1532.    ExitLog086:
  1533.       ret
  1534. Log086   ENDP
  1535.  
  1536.  
  1537. END
  1538.