home *** CD-ROM | disk | FTP | other *** search
/ The C Users' Group Library 1994 August / wc-cdrom-cusersgrouplibrary-1994-08.iso / listings / v_01_01 / 1n01034a < prev    next >
Text File  |  1990-05-17  |  8KB  |  341 lines
  1. ;**********
  2. ;
  3. ; Listing two: COMPLEX.ASM
  4. ; Lee Daniel Crocker
  5. ; 3/15/90
  6. ;
  7.  
  8. .model small, c
  9. .8087
  10.  
  11. @fstat macro
  12.       fstsw [status]          ; Can be replaced with
  13.       mov   ax, [status]      ; fstsw ax on 287
  14.       sahf
  15.       endm
  16.  
  17. qptr  equ   <qword ptr> ; To save typing
  18.  
  19. .data
  20.  
  21. degree      dw    8
  22. degreem1    dw    7
  23. ilimit      dw    1000
  24. status      dw    0
  25.  
  26. floatmin    dq    1.175494351e-38
  27. epsilon     dq    0.000001
  28.  
  29. oldr        dq    ?
  30. oldi        dq    ?
  31. newr        dq    ?
  32. newi        dq    ?
  33.  
  34. extrn       eps2:qword, roots:qword
  35.  
  36. ;**********
  37. ;
  38. ; C-callable functions.  Pretty straightforward stuff:
  39. ; just load the parameters onto the 8087 stack, call
  40. ; the routine that does the work, and pop the result
  41. ; back to memory.
  42. ;
  43.  
  44. .code
  45.  
  46. ;
  47. ; void cmult(COMPLEX *f1, COMPLEX *f2, COMPLEX *result);
  48. ;
  49.  
  50. cmult proc f1:ptr qword, f2:ptr qword, result:ptr qword
  51.       mov   bx, [f1]
  52.       fld   qptr [bx+8]
  53.       fld   qptr [bx]
  54.  
  55.       mov   bx, [f2]
  56.       fld   qptr [bx+8]
  57.       fld   qptr [bx]
  58.  
  59.       call  docmult
  60.  
  61.       mov   bx, [result]
  62.       fstp  qptr [bx]
  63.       fstp  qptr [bx+8]
  64.       fwait             ; synchronize memory write
  65.       ret
  66. cmult endp
  67.  
  68. ;
  69. ; int cdiv(COMPLEX *f1, COMPLEX *f2, COMPLEX *result);
  70. ;
  71.  
  72. cdiv proc f1:ptr qword, f2:ptr qword, result:ptr qword
  73.       mov   bx, [f1]
  74.       fld   qptr [bx+8]
  75.       fld   qptr [bx]
  76.  
  77.       mov   bx, [f2]
  78.       fld   qptr [bx+8]
  79.       fld   qptr [bx]
  80.  
  81.       call  docdiv
  82.  
  83.       mov   bx, [result]
  84.       fstp  qptr [bx]
  85.       fstp  qptr [bx+8]
  86.       fwait             ; synchronize memory write
  87.       ret
  88. cdiv endp
  89.  
  90. ;
  91. ; void cpower(COMPLEX *base, unsigned exp,
  92. ;             COMPLEX *result);
  93. ;
  94.  
  95. cpower proc base:ptr qword, exp:word, result:ptr qword
  96.       mov   bx, [base]
  97.       fld   qptr [bx+8]
  98.       fld   qptr [bx]
  99.  
  100.       mov   ax, [exp]
  101.       call  docpower
  102.  
  103.       mov   bx, [result]
  104.       fstp  qptr [bx]
  105.       fstp  qptr [bx+8]
  106.       fwait             ; synchronize memory write
  107.       ret
  108. cpower endp
  109.  
  110. ;
  111. ;
  112. ;
  113.  
  114. onenewton:
  115.       fld   [oldi]
  116.       fld   [oldr]
  117.       mov   ax, [degreem1]
  118.       call  docpower
  119.  
  120.       fld   st(1)
  121.       fld   st(1)
  122.       fld   [oldi]
  123.       fld   [oldr]
  124.       call  docmult     ; - nr  ni  tr  ti
  125.  
  126.       fild  [degreem1]
  127.       fmul  st(2), st
  128.       fmul
  129.       fld1
  130.       fadd
  131.  
  132.       fild  [degree]
  133.       fmul  st(4), st
  134.       fmulp st(3), st
  135.  
  136.       fincstp
  137.       fxch  st(2)
  138.       fdecstp
  139.       fxch  st(2)
  140.       call  docdiv
  141.  
  142.       sub   dx, dx
  143.       or    ax, ax      ; cdiv's return
  144.       jnz   overf
  145.  
  146.       fld   [oldr]
  147.       fsub  st, st(1)
  148.       fabs
  149.       fcomp [epsilon]
  150.       @fstat
  151.       ja    onout
  152.  
  153.       fld   [oldi]
  154.       fsub  st, st(2)
  155.       fabs
  156.       fcomp [epsilon]
  157.       @fstat
  158.       ja    onout
  159. overf:
  160.       inc   dx
  161. onout:
  162.       ret
  163.  
  164. ;
  165. ; int calcnewton(COMPLEX *point);
  166. ;
  167.  
  168. calcnewton proc point:ptr qword
  169.       mov   bx, [point]
  170.       fld   qptr [bx]
  171.       fstp  [oldr]
  172.       fld   qptr [bx+8]
  173.       fstp  [oldi]
  174.  
  175.       mov   cx, [ilimit]
  176. ltop:
  177.       call  onenewton
  178.  
  179.       or    dx, dx
  180.       jz    lbot
  181.  
  182.       fstp  [newr]
  183.       fstp  [newi]
  184.       jmp   short lbreak
  185. lbot:
  186.       fst   [newr]
  187.       fstp  [oldr]
  188.       fst   [newi]
  189.       fstp  [oldi]
  190.       loop  ltop
  191. lbreak:
  192.       mov   dx, [ilimit]
  193.       sub   dx, cx
  194.  
  195.       mov   bx, [degreem1]
  196.       mov   cl, 4
  197.       shl   bx, cl
  198.       mov   cx, [degreem1]
  199.       add   bx, offset DGROUP:roots
  200.       fld   [eps2]
  201.       fld   [newi]
  202.       fld   [newr]
  203. l2top:
  204.       fld   qptr [bx]
  205.       fsub  st, st(1)
  206.       fabs
  207.       fcomp st(3)
  208.       @fstat
  209.       ja    l2next
  210.  
  211.       fld   qptr [bx+8]
  212.       fsub  st, st(2)
  213.       fabs
  214.       fcomp st(3)
  215.       @fstat
  216.       jna   l2out
  217. l2next:
  218.       sub   bx, 16
  219.       loop  l2top
  220. l2out:
  221.       ffree st(2)
  222.       ffree st(1)
  223.       ffree st    ; Now dx=iterations and
  224.                   ; cx=root.
  225.       mov   ax, cx
  226.       and   ax, 3
  227.       ret
  228. calcnewton endp
  229.  
  230. ;
  231. ; Subroutines used in various functions:
  232. ;
  233. ; Multiply complex [a,b] by [c,d].
  234. ; Uses only two extra stack elements.
  235. ;
  236.                         ; - FORTH-style 8087
  237.                         ; - stack trace:
  238.  
  239. docmult:                ; - a  b  c  d
  240.       fld   st(2)       ; - c  a  b  c  d
  241.       fmul  st, st(1)   ; - ac  a  b  c  d
  242.       fld   st(4)       ; - d  ac  a  b  c  d
  243.       fmul  st, st(3)   ; - bd  ac  a  b  c  d
  244.       fsub              ; - rx=(ac-bd)  a  b  c  d
  245.       fincstp           ; - a  b  c  d  ?  ?  ?  rx
  246.       fmulp st(3), st   ; - b  c  ad  ?  ?  ?  rx
  247.       fmul              ; - bc  ad  ?  ?  ?  rx
  248.       fadd              ; - ry=(ad+bc)  ?  ?  ?  rx
  249.       fld   st(4)       ; - rx  ry  ?  ?  ?  rx
  250.       ffree st(5)       ; - rx  ry
  251.       ret
  252.  
  253. ;
  254. ; Divide complex [c,d] by [a,b].
  255. ; Uses every stack element.
  256. ;
  257.  
  258. docdiv:                 ; - a  b  c  d
  259.       fld   st(1)       ; - b  a  b  c  d
  260.       fmul  st, st(2)   ; - bb  a  b  c  d
  261.       fld   st(1)       ; - a  bb  a  b  c  d
  262.       fmul  st, st(2)   ; - aa  bb  a  b  c  d
  263.       fadd              ; - m  a  b  c  d
  264.  
  265.       fcom  [floatmin]  ; Note: m is always
  266.       @fstat            ; positive here.
  267.       jb    overflow
  268.  
  269.       fld1              ; - 1  m  a  b  c  d
  270.       fdivr             ; - m  a  b  c  d
  271.       fstp  st(5)       ; - a  b  c  d  m
  272.       fincstp           ; - b  c  d  m  ?  ?  ?  a
  273.       fchs              ; - -b  c  d  m  ?  ?  ?  a
  274.       fdecstp           ; - a  -b  c  d  m
  275.       call  docmult     ; - rx  ry  m
  276.       fxch  st(2)       ; - m  ry  rx
  277.       fmul  st(1), st   ; - m  m*ry  rx
  278.       fmul  st, st(2)   ; - m*rx  m*ry  rx
  279.       ffree st(2)       ; - rx  ry
  280.  
  281.       sub   ax, ax
  282.       ret
  283. overflow:
  284.       sbb   ax, ax      ; sets AX=-1
  285.       ret
  286.  
  287. ;
  288. ; Raise complex [a,b] to integer power in AX.
  289. ; Uses every stack element.
  290. ;
  291.  
  292. docpower:               ; - a  b
  293.  
  294.       shr   ax, 1 ; Shift bit 0 of exp into CF.
  295.       jc    bit0
  296.  
  297.       fldz              ; - ry=0  a  b
  298.       fld1              ; - rx=1  ry=0  a  b
  299.       jmp   short whiletop
  300. bit0:
  301.       fld   st(1)       ; - ry=b  a  b
  302.       fld   st(1)       ; - rx=a  ry=b  a  b
  303.  
  304. whiletop:         ; Square [a,b].
  305.  
  306.       fld   st(2)       ; - a  rx  ry  a  b
  307.       fadd  st, st(4)   ; - (a+b)  rx  ry  a  b
  308.       fld   st(3)       ; - a  (a+b)  rx  ry  a  b
  309.       fsub  st, st(5)   ; - (a-b)  (a+b)  rx  ry  a  b
  310.       fmul              ; - t2  rx  ry  a  b
  311.       fxch  st(3)       ; - a  rx  ry  t2  b
  312.       fmul  st, st(4)   ; - ab  rx  ry  t2  b
  313.       fst   st(4)       ; - ab  rx  ry  t2  ab
  314.       faddp st(4), st   ; - rx  ry  t2=a  2ab=b
  315.  
  316.       shr   ax, 1       ; Shift next bit of exp into CF
  317.       jnc   nomult      ; If bit is set, multiply result
  318.                         ; by [a,b].
  319.  
  320.       fld   st(2)       ; - a  rx  ry  a  b
  321.       fmul  st, st(1)   ; - a*rx  rx  ry  a  b
  322.       fld   st(4)       ; - b  a*rx  rx  ry  a  b
  323.       fmul  st, st(3)   ; - b*ry  a*rx  rx  ry  a  b
  324.       fsub              ; - t2  rx  ry  a  b
  325.       fld   st(3)       ; - a  t2  rx  ry  a  b
  326.       fmulp st(3), st   ; - t2  rx  a*ry  a  b
  327.  
  328.       fld   st(4)       ; - b  t2  rx  a*ry  a  b
  329.       fmul  st, st(2)   ; - b*rx  t2  rx  a*ry  a  b
  330.       faddp st(3), st   ; - t2  rx  ry  a  b
  331.  
  332.       fstp  st(1)       ; - rx  ry  a  b
  333. nomult:
  334.       or    ax, ax      ; Set ZF if no more exp.
  335.       jnz   whiletop
  336.  
  337.       ffree st(3)       ; - rx  ry  a
  338.       ffree st(2)       ; - rx  ry
  339.       ret
  340. end
  341.