home *** CD-ROM | disk | FTP | other *** search
/ Shareware Overload / ShartewareOverload.cdr / windows / winsrc.zip / CALCMAND.ASM < prev    next >
Assembly Source File  |  1990-11-06  |  26KB  |  866 lines

  1. ;    CALCMAND.ASM - Mandelbrot/Julia Set calculation Routines
  2.  
  3. ;    This module runs as part of an overlay with calcfrac.c.
  4. ;    It must not be called from anywhere other than calcfrac.
  5.  
  6. ;    The routines in this code perform Mandelbrot and Julia set
  7. ;    calculations using 32-bit integer math as opposed to the
  8. ;    "traditional" floating-point approach.
  9.  
  10. ;    This code relies on several tricks to run as quickly as it does.
  11.  
  12. ;    One can fake floating point arithmetic by using integer
  13. ;    arithmetic and keeping track of the implied decimal point
  14. ;    if things are reasonable -- and in this case, they are.
  15. ;    I replaced code that looked like: z = x*y with code that
  16. ;    looks like:
  17. ;            ix = x * ifudge         (outside the loops)
  18. ;            iy = y * ifudge
  19. ;            ....
  20. ;            iz = (ix * iy) / ifudge     (inside the loops)
  21. ;    (and keep remembering that all the integers are "ifudged" bigger)
  22.  
  23. ;    The 386 has native 32-bit integer arithmetic, and (briefly) keeps
  24. ;    64-bit values around after 32-bit multiplies.    If the result is
  25. ;    divided down right away, you've got 64-bit arithmetic.   You just
  26. ;    have to ensure that the result after the divide is <= 32 bits long.
  27. ;    CPUs predating the 386 have to emulate 32-bit arithmetic using
  28. ;    16-bit arithmetic, which is significantly slower.
  29.  
  30. ;    Dividing is slow -- but shifting is fast, and we can select our
  31. ;    "fudge factor" to be a power of two, permitting us to use that
  32. ;    method instead.   In addition, the 386 can perform 32-bit wide
  33. ;    shifting -- and even 64-bit shifts with the following logic:
  34. ;            shdr    eax,edx,cl
  35. ;            shr    edx,cl
  36. ;    so we make sure that our "fudge factor" is a power of 2 and shift
  37. ;    it down that way.
  38. ;    Calcmand is hardcoded for a fudge factor of 2**29.
  39.  
  40.  
  41. ;                    Bert Tyler
  42.  
  43.  
  44. ;             required for compatibility if Turbo ASM
  45. IFDEF ??version
  46. MASM51
  47. QUIRKS
  48. ENDIF
  49.  
  50. .MODEL    medium,c
  51. DGROUP          group   _DATA,_DATA2
  52.  
  53. .8086
  54.  
  55.     ; these must NOT be in any segment!!
  56.     ; this get's rid of TURBO-C fixup errors
  57.  
  58.     extrn    keypressed:far        ; this routine is in 'general.asm'
  59.     extrn    getakey:far        ; this routine is in 'general.asm'
  60.     extrn    iplot_orbit:far     ; this routine is in 'calcfrac.c'
  61.     extrn    scrub_orbit:far     ; this routine is in 'calcfrac.c'
  62.  
  63. _DATA2        segment DWORD PUBLIC 'DATA'
  64.  
  65. FUDGEFACTOR    equ    29        ; default (non-potential) fudgefactor
  66.  
  67. ; ************************ External variables *****************************
  68.  
  69.     extrn    fractype:word        ; == 0 if Mandelbrot set, else Julia
  70.     extrn    inside:word        ; "inside" color, normally 1 (blue)
  71.     extrn    creal:dword, cimag:dword ; Julia Set Constant
  72.     extrn    delmin:dword        ; min increment - precision required
  73.     extrn    maxit:word        ; maximum iterations
  74.     extrn    lm:dword        ; magnitude bailout limit
  75.  
  76.     extrn    row:word, col:word    ; current pixel to calc
  77.     extrn    color:word        ; color calculated for the pixel
  78.     extrn    realcolor:word        ; color before inside,etc adjustments
  79.  
  80.     extrn    reset_periodicity:word    ; nonzero if to be reset
  81.     extrn    kbdcount:word        ; keyboard counter
  82.  
  83.     extrn    cpu:word        ; cpu type: 86, 186, 286, or 386
  84.     extrn    dotmode:word
  85.  
  86.     extrn    show_orbit:word     ; "show-orbit" flag
  87.     extrn    orbit_ptr:word        ; "orbit pointer" flag
  88.     extrn    periodicitycheck:word    ; no periodicity if zero
  89.  
  90.     public    linitx,linity        ; caller sets these
  91.     public    savedmask        ; caller sets this
  92.  
  93. ; ************************ Internal variables *****************************
  94.  
  95.         align    4
  96. x        dd    0        ; temp value: x
  97. y        dd    0        ; temp value: y
  98. absx        dd    0        ; temp value: abs(x)
  99. linitx        dd    0        ; initial value, set by calcfrac
  100. linity        dd    0        ; initial value, set by calcfrac
  101. savedmask    dd    0        ; saved values mask
  102. savedx        dd    0        ; saved values of X and Y iterations
  103. savedy        dd    0        ;  (for periodicity checks)
  104. k        dw    0        ; iteration countdown counter
  105. oldcolor    dw    0        ; prior pixel's escape time k value
  106. savedand    dw    0        ; AND value for periodicity checks
  107. savedincr    dw    0        ; flag for incrementing AND value
  108. period        db    0        ; periodicity, if in the lake
  109.  
  110. _DATA2        ends
  111.  
  112. .CODE
  113.  
  114. ; ***************** Function calcmandasm() **********************************
  115.  
  116.     public    calcmandasm
  117.  
  118. FRAME    MACRO regs
  119.     push    bp
  120.     mov    bp, sp
  121.     IRP    reg, <regs>
  122.       push    reg
  123.       ENDM
  124.     ENDM
  125.  
  126. UNFRAME MACRO regs
  127.     IRP    reg, <regs>
  128.       pop reg
  129.       ENDM
  130.     pop bp
  131.     ENDM
  132.  
  133. calcmandasm    proc
  134.     FRAME    <di,si>         ; std frame, for TC++ overlays
  135.     sub    ax,ax            ; clear ax
  136.     cmp    periodicitycheck,ax    ; periodicity checking disabled?
  137.     je    initoldcolor        ;  yup, set oldcolor 0 to disable it
  138.     cmp    reset_periodicity,ax    ; periodicity reset?
  139.     je    short initparms     ; inherit oldcolor from prior invocation
  140.     mov    ax,maxit        ; yup.    reset oldcolor to maxit-250
  141.     sub    ax,250            ; (avoids slowness at high maxits)
  142. initoldcolor:
  143.     mov    oldcolor,ax        ; reset oldcolor
  144.  
  145. initparms:
  146.     mov    ax,word ptr creal    ; initialize x == creal
  147.     mov    dx,word ptr creal+2    ;  ...
  148.     mov    word ptr x,ax        ;  ...
  149.     mov    word ptr x+2,dx     ;  ...
  150.  
  151.     mov    ax,word ptr cimag    ; initialize y == cimag
  152.     mov    dx,word ptr cimag+2    ;  ...
  153.     mov    word ptr y,ax        ;  ...
  154.     mov    word ptr y+2,dx     ;  ...
  155.  
  156.     mov    ax,maxit        ; setup k = maxit
  157.     inc    ax            ; (+ 1)
  158.     mov    k,ax            ;  (decrementing to 0 is faster)
  159.  
  160.     cmp    fractype,1        ; julia or mandelbrot set?
  161.     je    short dojulia        ; julia set - go there
  162.  
  163. ;    (Tim wants this code changed so that, for the Mandelbrot,
  164. ;    Z(1) = (x + iy) + (a + ib).  Affects only "fudged" Mandelbrots.
  165. ;    (for the "normal" case, a = b = 0, and this works, too)
  166. ;    cmp    word ptr x,0        ; Mandelbrot shortcut:
  167. ;    jne    short doeither        ;  if creal = cimag = 0,
  168. ;    cmp    word ptr x+2,0        ; the first iteration can be emulated.
  169. ;    jne    short doeither        ;  ...
  170. ;    cmp    word ptr y,0        ;  ...
  171. ;    jne    short doeither        ;  ...
  172. ;    cmp    word ptr y+2,0        ;  ...
  173. ;    jne    short doeither        ;  ...
  174. ;    dec    k            ; we know the first iteration passed
  175. ;    mov    dx,word ptr linitx+2    ; copy x = linitx
  176. ;    mov    ax,word ptr linitx    ;  ...
  177. ;    mov    word ptr x+2,dx     ;  ...
  178. ;    mov    word ptr x,ax        ;  ...
  179. ;    mov    dx,word ptr linity+2    ; copy y = linity
  180. ;    mov    ax,word ptr linity    ;  ...
  181. ;    mov    word ptr y+2,dx     ;  ...
  182. ;    mov    word ptr y,ax        ;  ...
  183.  
  184.     dec    k            ; we know the first iteration passed
  185.     mov    dx,word ptr linitx+2    ; add x += linitx
  186.     mov    ax,word ptr linitx    ;  ...
  187.     add    word ptr x,ax        ;  ...
  188.     adc    word ptr x+2,dx     ;  ...
  189.     mov    dx,word ptr linity+2    ; add y += linity
  190.     mov    ax,word ptr linity    ;  ...
  191.     add    word ptr y,ax        ;  ...
  192.     adc    word ptr y+2,dx     ;  ...
  193.     jmp    short doeither        ; branch around the julia switch
  194.  
  195. dojulia:                ; Julia Set initialization
  196.                     ; "fudge" Mandelbrot start-up values
  197.     mov    ax,word ptr x        ; switch x with linitx
  198.     mov    dx,word ptr x+2     ;  ...
  199.     mov    bx,word ptr linitx    ;  ...
  200.     mov    cx,word ptr linitx+2    ;  ...
  201.     mov    word ptr x,bx        ;  ...
  202.     mov    word ptr x+2,cx     ;  ...
  203.     mov    word ptr linitx,ax    ;  ...
  204.     mov    word ptr linitx+2,dx    ;  ...
  205.  
  206.     mov    ax,word ptr y        ; switch y with linity
  207.     mov    dx,word ptr y+2     ;  ...
  208.     mov    bx,word ptr linity    ;  ...
  209.     mov    cx,word ptr linity+2    ;  ...
  210.     mov    word ptr y,bx        ;  ...
  211.     mov    word ptr y+2,cx     ;  ...
  212.     mov    word ptr linity,ax    ;  ...
  213.     mov    word ptr linity+2,dx    ;  ...
  214.  
  215. doeither:                ; common Mandelbrot, Julia set code
  216.     mov    period,0        ; claim periodicity of 1
  217.     mov    savedand,1        ; initial periodicity check
  218.     mov    savedincr,1        ;  flag for incrementing periodicity
  219.     mov    word ptr savedx+2,0ffffh; impossible value of "old" x
  220.     mov    word ptr savedy+2,0ffffh; impossible value of "old" y
  221.     mov    orbit_ptr,0        ; clear orbits
  222.  
  223.     dec    kbdcount        ; decrement the keyboard counter
  224.     jns    short nokey        ;  skip keyboard test if still positive
  225.     mov    kbdcount,10        ; stuff in a low kbd count
  226.     cmp    show_orbit,0        ; are we showing orbits?
  227.     jne    quickkbd        ;  yup.  leave it that way.
  228.     mov    kbdcount,5000        ; else, stuff an appropriate count val
  229.     cmp    cpu,386         ; ("appropriate" to the CPU)
  230.     je    short kbddiskadj    ;  ...
  231.     cmp    word ptr delmin+2,1    ; is 16-bit math good enough?
  232.     ja    kbddiskadj        ;  yes. test less often
  233.     mov    kbdcount,500        ;  no.    test more often
  234. kbddiskadj:
  235.     cmp    dotmode,11        ; disk video?
  236.     jne    quickkbd        ;  no, leave as is
  237.     shr    kbdcount,1        ; yes, reduce count
  238.     shr    kbdcount,1        ;  ...
  239. quickkbd:
  240.     call    far ptr keypressed    ; has a key been pressed?
  241.     cmp    ax,0            ;  ...
  242.     je    nokey            ; nope.  proceed
  243.     mov    kbdcount,0        ; make sure it goes negative again
  244.     cmp    ax,'o'                  ; orbit toggle hit?
  245.     je    orbitkey        ;  yup.  show orbits
  246.     cmp    ax,'O'                  ; orbit toggle hit?
  247.     jne    keyhit            ;  nope.  normal key.
  248. orbitkey:
  249.     call    far ptr getakey     ; read the key for real
  250.     mov    ax,1            ; reset orbittoggle = 1 - orbittoggle
  251.     sub    ax,show_orbit        ;  ...
  252.     mov    show_orbit,ax        ;  ...
  253.     jmp    short nokey        ; pretend no key was hit
  254. keyhit: mov    ax,-1            ; return with -1
  255.     mov    color,ax        ; set color to -1
  256.     UNFRAME <si,di>         ; pop stack frame
  257.     ret                ; bail out!
  258.  
  259. nokey:
  260.     cmp    show_orbit,0        ; is orbiting on?
  261.     jne    no16bitcode        ;  yup.  slow down.
  262.     cmp    cpu,386         ; are we on a 386?
  263.     je    short code386bit    ;  YAY!! 386-class speed!
  264.     cmp    word ptr delmin+2,1    ; OK, we're desperate.  16 bits OK?
  265.     ja    yes16bitcode        ;  YAY!  16-bit speed!
  266. no16bitcode:
  267.     call    near ptr code32bit    ; BOO!! nap time.  Full 32 bit math
  268.     jmp    kloopend        ;  bypass the 386-specific code.
  269. yes16bitcode:
  270.     call    near ptr code16bit    ; invoke the 16-bit version
  271.     jmp    kloopend        ;  bypass the 386-specific code.
  272.  
  273. .386                    ; 386-specific code starts here
  274.  
  275. code386bit:
  276.     cmp    word ptr delmin+2,3    ; is 16-bit math good enough?
  277.     jbe    code386_32        ; nope, go do 32 bit stuff
  278. IFDEF ??version
  279.     jmp    code386_32        ; TASM screws up IMUL EBX,EBX!!
  280. ENDIF
  281.  
  282.     ; 16 bit on 386, now we are really gonna move
  283.     movsx    esi,word ptr x+2    ; use SI for X
  284.     movsx    edi,word ptr y+2    ; use DI for Y
  285.     push    ebp
  286.     mov    ebp,-1
  287.     shl    ebp,FUDGEFACTOR-1
  288.     mov    cx,FUDGEFACTOR-16
  289.  
  290. kloop386_16:   ; cx=bitshift-16, ebp=overflow.mask
  291.  
  292.     mov    ebx,esi         ; compute (x * x)
  293.     imul    ebx,ebx         ;  ...
  294.     test    ebx,ebp         ; say, did we overflow? <V20-compat>
  295.     jnz    short end386_16     ;  (oops.  We done.)
  296.     shr    ebx,cl            ; get result down to 16 bits
  297.  
  298.     mov    edx,edi         ; compute (y * y)
  299.     imul    edx,edx         ;  ...
  300.     test    edx,ebp         ; say, did we overflow? <V20-compat>
  301.     jnz    short end386_16     ;  (oops.  We done.)
  302.     shr    edx,cl            ; get result down to 16 bits
  303.  
  304.     mov    ax,bx            ; compute (x*x - y*y) / fudge
  305.     sub    bx,dx            ;  for the next iteration
  306.  
  307.     add    ax,dx            ; compute (x*x + y*y) / fudge
  308.  
  309.     cmp    ax,word ptr lm+2    ; while (xx+yy < lm)
  310.     jae    short end386_16     ;  ...
  311.  
  312.     imul    edi,esi         ; compute (y * x)
  313.     shl    edi,1            ; ( * 2 / fudge)
  314.     sar    edi,cl
  315.     add    di,word ptr linity+2    ; (2*y*x) / fudge + linity
  316.     movsx    edi,di            ; save as y
  317.  
  318.     add    bx,word ptr linitx+2    ; (from above) (x*x - y*y)/fudge + linitx
  319.     movsx    esi,bx            ; save as x
  320.  
  321.     mov    ax,oldcolor        ; recall the old color
  322.     cmp    ax,k            ; check it against this iter
  323.     jge    short chkpd386_16    ;  yup.  do periodicity check.
  324. nonmax386_16:
  325.     dec    k            ; while (k < maxit)
  326.     jnz    short kloop386_16    ; try, try again
  327. end386_16:
  328.     pop    ebp
  329.     jmp    kloopend        ; we done
  330.  
  331. chkpd386_16:
  332.     mov    ax,k            ; set up to test for save-time
  333.     test    ax,savedand        ; save on 0, check on anything else
  334.     jz    short chksv386_16    ;  time to save a new "old" value
  335.     mov    ax,si            ; load up x
  336.     xor    ax,word ptr savedx+2    ; does X match?
  337.     test    ax,word ptr savedmask+2 ;  truncate to appropriate precision
  338.     jne    short nonmax386_16    ;  nope.  forget it.
  339.     mov    ax,di            ; now test y
  340.     xor    ax,word ptr savedy+2    ; does Y match?
  341.     test    ax,word ptr savedmask+2 ;  truncate to appropriate precision
  342.     jne    short nonmax386_16    ;  nope.  forget it.
  343.     mov    period,1        ; note that we have found periodicity
  344.     mov    k,0            ; pretend maxit reached
  345.     jmp    short end386_16
  346. chksv386_16:
  347.     mov    word ptr savedx+2,si    ; save x
  348.     mov    word ptr savedy+2,di    ; save y
  349.     dec    savedincr        ; time to change the periodicity?
  350.     jnz    short nonmax386_16    ;  nope.
  351.     shl    savedand,1        ; well then, let's try this one!
  352.     inc    savedand        ;  (2**n -1)
  353.     mov    savedincr,4        ; and reset the increment flag
  354.     jmp    short nonmax386_16
  355.  
  356.     ; 32bit on 386:
  357. code386_32:
  358.     mov    esi,x            ; use ESI for X
  359.     mov    edi,y            ; use EDI for Y
  360.  
  361. ;    This is the main processing loop.  Here, every T-state counts...
  362.  
  363. kloop:                    ; for (k = 0; k <= maxit; k++)
  364.  
  365.     mov    eax,esi         ; compute (x * x)
  366.     imul    esi            ;  ...
  367.     shrd    eax,edx,FUDGEFACTOR    ; ( / fudge)
  368.     shr    edx,FUDGEFACTOR-1    ; (complete 64-bit shift and check
  369.     jne    short kloopend1     ; bail out if too high
  370.     mov    ebx,eax         ; save this for below
  371.  
  372.     mov    eax,edi         ; compute (y * y)
  373.     imul    edi            ;  ...
  374.     shrd    eax,edx,FUDGEFACTOR    ; ( / fudge)
  375.     shr    edx,FUDGEFACTOR-1    ; (complete 64-bit shift and check
  376.     jne    short kloopend1     ; bail out if too high
  377.  
  378.     mov    ecx,ebx         ; compute (x*x - y*y) / fudge
  379.     sub    ebx,eax         ;  for the next iteration
  380.  
  381.     add    ecx,eax         ; compute (x*x + y*y) / fudge
  382.     cmp    ecx,lm            ; while (lr < lm)
  383.     jae    short kloopend1     ;  ...
  384.  
  385.     mov    eax,edi         ; compute (y * x)
  386.     imul    esi            ;  ...
  387.     shrd    eax,edx,FUDGEFACTOR-1    ;  ( * 2 / fudge)
  388.     add    eax,linity        ;  (above) + linity
  389.     mov    edi,eax         ;  save this as y
  390.  
  391. ;    (from the earlier code)     ; compute (x*x - y*y) / fudge
  392.     add    ebx,linitx        ;    + linitx
  393.     mov    esi,ebx         ; save this as x
  394.  
  395.     mov    ax,oldcolor        ; recall the old color
  396.     cmp    ax,k            ; check it against this iter
  397.     jge    short chkperiod1
  398. nonmax1:
  399.     dec    k            ; while (k < maxit) (dec to 0 is faster)
  400.     jnz    short kloop        ; while (k < maxit) ...
  401. kloopend1:
  402.     jmp    short kloopend32    ; we done.
  403.  
  404. chkperiod1:
  405.     mov    eax,esi
  406.     xor    eax,savedx
  407.     test    eax,savedmask
  408.     jnz    short chksave1
  409.     mov    eax,edi
  410.     xor    eax,savedy
  411.     test    eax,savedmask
  412.     jnz    short chksave1
  413.     mov    period,1        ; note that we have found periodicity
  414.     mov    k,0            ; pretend maxit reached
  415.     jmp    short kloopend1
  416. chksave1:
  417.     mov    ax,k
  418.     test    ax,savedand
  419.     jne    short nonmax1
  420.     mov    savedx,esi
  421.     mov    savedy,edi
  422.     dec    savedincr        ; time to change the periodicity?
  423.     jnz    short nonmax1        ;  nope.
  424.     shl    savedand,1        ; well then, let's try this one!
  425.     inc    savedand        ;  (2**n -1)
  426.     mov    savedincr,4        ; and reset the increment flag
  427.     jmp    short nonmax1
  428.  
  429. kloopend32:
  430.  
  431. .8086                    ; 386-specific code ends here
  432.  
  433. kloopend:
  434.     cmp    orbit_ptr,0        ; any orbits to clear?
  435.     je    noorbit2        ;  nope.
  436.     call    far ptr scrub_orbit    ; clear out any old orbits
  437. noorbit2:
  438.  
  439.     mov    ax,k            ; set old color
  440.     sub    ax,10            ; minus 10, for safety
  441.     mov    oldcolor,ax        ; and save it as the "old" color
  442.     mov    ax,maxit        ; compute color
  443.     sub    ax,k            ;  (first, re-compute "k")
  444.     sub    kbdcount,ax        ; adjust the keyboard count
  445.     cmp    ax,1            ; convert any "outlier" region
  446.     jge    short coloradjust1    ;  (where abs(x) > 2 or abs(y) > 2)
  447.     mov    ax,1            ;   to look like we ran through
  448. coloradjust1:                ;    at least one loop.
  449.     mov    realcolor,ax        ; result before adjustments
  450.     cmp    ax,maxit        ; did we max out on iterations?
  451.     jne    short wedone        ;  nope.
  452.     mov    oldcolor,ax        ; set "oldcolor" to maximum
  453.     cmp    inside,0        ; is "inside" >= 0?
  454.     jl    wedone            ;  nope.  leave it at "maxit"
  455.     mov    ax,inside        ; reset max-out color to default
  456.     cmp    periodicitycheck,0    ; show periodicity matches?
  457.     jge    wedone            ;  nope.
  458.     mov    al,period        ;  reset color to periodicity flag
  459. wedone:                 ;
  460.     mov    color,ax        ; save the color result
  461.     UNFRAME <si,di>         ; pop stack frame
  462.     ret                ; and return with color
  463.  
  464. calcmandasm endp
  465.  
  466.  
  467. ; ******************** Function code16bit() *****************************
  468. ;
  469. ;    Performs "short-cut" 16-bit math where we can get away with it.
  470. ;
  471.  
  472. code16bit    proc    near
  473.  
  474.     mov    si,word ptr x+2     ; use SI for X
  475.     mov    di,word ptr y+2     ; use DI for Y
  476.  
  477. start16bit:
  478.     mov    ax,si            ; compute (x * x)
  479.     imul    si            ;  ...
  480.     cmp    dx,0            ; say, did we overflow? <V20-compat>
  481.     jl    end16bit        ;  (oops.  We done.)
  482.     mov    cx,32-FUDGEFACTOR    ; ( / fudge)
  483. loop16bit1:
  484.     shl    ax,1            ;  ...
  485.     rcl    dx,1            ;  ...
  486.     jo    end16bit        ;  (oops.  overflow)
  487.     loop    loop16bit1        ;  ...
  488.     mov    bx,dx            ; save this for a tad
  489.  
  490.     mov    ax,di            ; compute (y * y)
  491.     imul    di            ;  ...
  492.     cmp    dx,0            ; say, did we overflow? <V20-compat>
  493.     jl    end16bit        ;  (oops.  We done.)
  494.     mov    cx,32-FUDGEFACTOR    ; ( / fudge)
  495. loop16bit2:
  496.     shl    ax,1            ;  ...
  497.     rcl    dx,1            ;  ...
  498.     jo    end16bit        ;  (oops.  overflow)
  499.     loop    loop16bit2        ;  ...
  500.  
  501.     mov    cx,bx            ; compute (x*x - y*y) / fudge
  502.     sub    bx,dx            ;  for the next iteration
  503.  
  504.     add    cx,dx            ; compute (x*x + y*y) / fudge
  505.     jo    end16bit        ; bail out if too high
  506.     js    end16bit        ;  ...
  507.  
  508.     cmp    cx,word ptr lm+2    ; while (xx+yy < lm)
  509.     jae    end16bit        ;  ...
  510.  
  511.     mov    ax,di            ; compute (y * x)
  512.     imul    si            ;  ...
  513.     mov    cx,33-FUDGEFACTOR    ; ( * 2 / fudge)
  514. loop16bit3:
  515.     shl    ax,1            ;  ...
  516.     rcl    dx,1            ;  ...
  517.     loop    loop16bit3        ;  ...
  518.     add    dx,word ptr linity+2    ; (2*y*x) / fudge + linity
  519.     mov    di,dx            ; save as y
  520.  
  521.     add    bx,word ptr linitx+2    ; (from above) (x*x - y*y)/fudge + linitx
  522.     mov    si,bx            ; save as x
  523.  
  524.     mov    ax,oldcolor        ; recall the old color
  525.     cmp    ax,k            ; check it against this iter
  526.     jl    short nonmax3        ;  nope.  bypass periodicity check.
  527.     mov    word ptr x+2,si     ; save x for periodicity check
  528.     mov    word ptr y+2,di     ; save y for periodicity check
  529.     call    checkperiod        ; check for periodicity
  530. nonmax3:
  531.  
  532.     dec    k            ; while (k < maxit)
  533.     jz    end16bit        ;  we done.
  534.  
  535.     jmp    start16bit        ; try, try again.
  536.  
  537. end16bit:                ; we done.
  538.     ret
  539. code16bit    endp
  540.  
  541.  
  542. ;    The following routine checks for periodic loops (a known
  543. ;    method of decay inside "Mandelbrot Lake", and an easy way to
  544. ;    bail out of "lake" points quickly).  For speed, only the
  545. ;    high-order sixteen bits of X and Y are checked for periodicity.
  546. ;    For accuracy, this routine is only fired up if the previous pixel
  547. ;    was in the lake (which means that the FIRST "wet" pixel was
  548. ;    detected by the dull-normal maximum iteration process).
  549.  
  550. checkperiod    proc near        ; periodicity check
  551.     mov    ax,k            ; set up to test for save-time
  552.     test    ax,savedand        ; save on 0, check on anything else
  553.     jz    checksave        ;  time to save a new "old" value
  554.     mov    dx,word ptr x+2     ; load up x
  555.     and    dx,word ptr savedmask+2 ;  truncate to appropriate precision
  556.     cmp    dx,word ptr savedx+2    ; does X match?
  557.     jne    checkdone        ;  nope.  forget it.
  558.     mov    ax,word ptr x        ; load up x
  559.     and    ax,word ptr savedmask    ;  truncate to appropriate precision
  560.     cmp    ax,word ptr savedx    ; does X match?
  561.     jne    checkdone        ;  nope.  forget it.
  562.     mov    dx,word ptr y+2     ; now test y
  563.     and    dx,word ptr savedmask+2 ;  truncate to appropriate precision
  564.     cmp    dx,word ptr savedy+2    ; does Y match?
  565.     jne    checkdone        ;  nope.  forget it.
  566.     mov    ax,word ptr y        ; load up y
  567.     and    ax,word ptr savedmask    ;  truncate to appropriate precision
  568.     cmp    ax,word ptr savedy    ; does Y match?
  569.     jne    checkdone        ;  nope.  forget it.
  570.     mov    period,1        ; note that we have found periodicity
  571.     mov    k,1            ; pretend maxit reached
  572. checksave:
  573.     mov    dx,word ptr x+2     ; load up x
  574.     and    dx,word ptr savedmask+2 ;  truncate to appropriate precision
  575.     mov    word ptr savedx+2,dx    ;  and save it
  576.     mov    ax,word ptr x        ; load up x
  577.     and    ax,word ptr savedmask    ;  truncate to appropriate precision
  578.     mov    word ptr savedx,ax    ;  and save it
  579.     mov    dx,word ptr y+2     ; load up y
  580.     and    dx,word ptr savedmask+2 ;  truncate to appropriate precision
  581.     mov    word ptr savedy+2,dx    ;  and save it
  582.     mov    ax,word ptr y        ; load up y
  583.     and    ax,word ptr savedmask    ;  truncate to appropriate precision
  584.     mov    word ptr savedy,ax    ;  and save it
  585.     dec    savedincr        ; time to change the periodicity?
  586.     jnz    checkdone        ;  nope.
  587.     shl    savedand,1        ; well then, let's try this one!
  588.     inc    savedand        ;  (2**n -1)
  589.     mov    savedincr,4        ; and reset the increment flag
  590. checkdone:
  591.     ret                ; we done.
  592. checkperiod    endp
  593.  
  594.  
  595. ; ******************** Function code32bit() *****************************
  596. ;
  597. ;    Perform the 32-bit logic required using 16-bit logic
  598. ;
  599. ;    New twice as fast logic,
  600. ;       Courtesy of Bill Townsend and Mike Gelvin (CIS:73337,520)
  601.  
  602. square    macro    doneaddr
  603.     local    donejp,skip1y,skip2y,squrey,allzero
  604.     sub    bh,bh
  605.     sub    cx,cx
  606.     mov    bp,cx
  607.     or    ax,ax            ;GET LOW HWORD
  608.     jz    skip1y            ;LOW HWORD IS ZERO, ONLY DO HIGH HWORD
  609.     mov    si,ax
  610.     mul    ax            ;SQUARE LOW HWORD - B*B
  611.     mov    bh,dh
  612.     or    di,di            ;TEST HIGH HWORD
  613.     jz    squrey            ;HIGH HWORD ZERO, SKIP THE FOLLOWING
  614.     mov    ax,di
  615.     mul    si            ;GET MIDDLE PART - A*B
  616.     add    bh,ah            ;add
  617.     adc    cx,dx            ; A*B
  618.     adc    bp,0            ;  twice
  619.     add    bh,ah            ;   -
  620.     adc    cx,dx            ;    -
  621.     adc    bp,0            ;     SI,CX,BH = 2(A*B)+(B*B)
  622.     jmp    short skip2y
  623. donejp: JMP    doneaddr        ;M0=DONEJP
  624. skip1y: or    di,di            ;M1=SKIP1Y
  625.     jz    allzero
  626. skip2y: mov    ax,di            ;M2=SKIP2Y
  627.     mul    ax            ;SQUARE HIGH HWORD - A*A
  628.     add    cx,ax
  629.     adc    bp,dx
  630. squrey: shl    bh,1            ;M3=SQUREY
  631.     rcl    cx,1
  632.     rcl    bp,1
  633. ;    jo    donejp            ; squaring a +ve, top bit known off
  634.     shl    bh,1
  635.     rcl    cx,1
  636.     rcl    bp,1
  637. ;    jo    donejp            ; squaring a +ve, 2nd bit known off
  638.     shl    bh,1
  639.     rcl    cx,1
  640.     rcl    bp,1
  641.     jo    donejp
  642.     add    bp,0
  643.     js    donejp            ; if went neg have overflow also
  644. allzero:
  645.  
  646.     endm    ;#EM
  647.  
  648.  
  649. code32bit    proc near
  650. ;
  651. ; BL IS USED FOR THE FLAGS,TWO LOW ORDER BITS ARE "NEGSWT", THIRD BIT IS XSIGN
  652. ;   XSIGN IS ONE IF THE SIGN OF THE NEW "X" IS NEGATIVE
  653. ;   NEGSWT IS ONE IF EITHER "X" OR "Y" ARE NEGATIVE, BUT NOT BOTH NEGATIVE
  654. ;
  655. xsign    equ    00000100b
  656. negswt    equ    00000001b
  657.  
  658.     push    bp
  659.     sub    si,si            ; make a zero __________
  660.     mov    bx,si            ; BOTH XSIGN & NEGSWT START OFF ZERO
  661.     mov    di,1            ; MAKE A ONE
  662.     mov    cx,word ptr y        ; cx=y
  663.     mov    bp,word ptr y+2     ; bp=y+2
  664.     cmp    bp,si            ; y>0?
  665.     jns    qnotng            ;  yup
  666.     not    cx            ; nope, negate it
  667.     not    bp            ; ...
  668.     add    cx,di            ; ...
  669.     adc    bp,si            ; ...
  670.     inc    bl            ; INCREMENT negswt
  671. qnotng: mov    ax,word ptr x
  672.     mov    dx,word ptr x+2
  673.     cmp    dx,si            ; x>0?
  674.     jns    pnotng            ;  yup
  675.     not    ax            ; nope, negate it
  676.     not    dx            ; ...
  677.     add    ax,di            ; ...
  678.     adc    dx,si            ; ...
  679.     inc    bl            ; INCREMENT negswt
  680. pnotng: mov    word ptr absx,ax    ; save unsigned x
  681.     mov    word ptr absx+2,dx    ; ...
  682.  
  683. ;    iteration loop
  684. nextit: push    bp            ; save y on stack
  685.     push    cx            ; ...
  686.     mov    ax,cx            ; load y for squaring
  687.     mov    di,bp            ;  ...
  688.     square    done2            ; square di,ax & br to jmp_done2 if overflow
  689.     push    cx            ; return results in bp,cx
  690.     push    bp            ; save y*y on stack
  691.     mov    ax,word ptr absx    ; load x for squaring
  692.     mov    di,word ptr absx+2    ;  ...
  693.     square    done            ; square di,ax & br to done if overflow
  694.     mov    di,bp            ; return results in bp,cx
  695.     mov    si,cx            ; save x*x in di,si
  696.     pop    dx            ; unstack y*y into dx,ax
  697.     pop    ax            ;  ...
  698.     ADD    cx,ax            ; calc y*y + x*x
  699.     ADC    bp,dx            ;  ...
  700.     jo    done2            ; overflow?
  701.     js    done2            ; overflow?
  702.     cmp    bp,word ptr lm+2    ; while (lr < lm)
  703.     jb    short nextxy        ;  keep going
  704.     jmp    short done2        ; magnitude is past limit, bailout
  705. done:    add    sp,4            ; discard saved value of "y*y"
  706. done2:    add    sp,4            ; discard saved value of "y"
  707. xydone: pop    bp            ; restore saved bp
  708.     ret
  709.  
  710. ;---------------------------------------------------------------------------
  711. ;y=nz, x=nz, x+2=z
  712. x2is0:    pop    ax            ; get old y+2
  713.     or    ax,ax            ; test y+2
  714.     jz    jgotxy
  715.     mul    si            ; do y+2(ax)*x(si)
  716.     add    bh,ah
  717.     adc    cx,dx
  718.     adc    bp,0
  719. jgotxy: jmp    gotxy
  720.  
  721. ;x=z
  722. xis0:    or    di,di            ; test x+2
  723.     jz    popy_2            ; zero all done, around again please
  724.     mul    di            ; do y(ax)*x+2(di)
  725.     add    bh,ah
  726.     adc    cx,dx
  727.     adc    bp,0
  728.     pop    ax            ; get old y+2
  729.     or    ax,ax
  730.     jz    gotxy
  731. mulxyj: jmp    mulxy2
  732.  
  733. ;x=z, x+2=z
  734. popy_2: pop    ax            ; discard old y+2
  735.     sub    cx,cx            ; x==0, so x*y is zero
  736.     mov    bp,cx            ; ...
  737.     jmp    signok            ; go add b and store new y
  738.  
  739. nextxy: sub    si,ax            ; subtract y*y from x*x
  740.     sbb    di,dx            ;  ...
  741.     add    si,word ptr linitx    ; add "A"
  742.     adc    di,word ptr linitx+2    ;  ...
  743.     mov    word ptr x,si        ; store new x = x*x+y*y+a
  744.     mov    word ptr x+2,di     ;  ...
  745.     jns    swtxxx            ; NO SIGN, NOT NEG
  746.     or    bl,xsign        ; REMEMBER THE NEW "X" IS NEGATIVE
  747.     not    si            ; make it positive in register
  748.     not    di            ;  ...
  749.     add    si,1            ;  ...
  750.     adc    di,0            ;  ...
  751. swtxxx: xchg    di,word ptr absx+2    ; save the new x, load the old
  752.     xchg    si,word ptr absx    ;  ...
  753.  
  754.     ; now calculate x*y
  755.     sub    bh,bh            ; zero some registers
  756.     sub    cx,cx            ;  ...
  757.     mov    bp,cx            ;  ...
  758.     pop    ax            ; get old "y"
  759.     or    ax,ax            ; test "y"
  760.     jz    yis0            ; br if "y" is zero
  761.     or    si,si            ; test old "x"
  762.     jz    xis0            ; br if "x" is zero
  763.     mov    bp,ax            ; save old "y"
  764.     mul    si            ; do y(ax)*x(si)
  765.     mov    ax,bp            ; get old "y" again
  766.     mov    bp,cx            ; RE-ZERO BP
  767.     mov    bh,dh            ; save/set low hword(bx)
  768.     or    di,di            ; test old "x+2"
  769.     jz    x2is0
  770.     mul    di            ; do y(ax)*x+2(di)
  771.     add    bh,ah
  772.     adc    cx,dx
  773.     adc    bp,0
  774.  
  775. yis0:    pop    ax            ; get old "y+2"
  776.     or    ax,ax            ; test y+2
  777.     jz    gotxy            ; br if y+2 is zero
  778.     mov    dx,si            ; dx="x"
  779.     mov    si,ax            ; si=save "y+2"
  780.     mul    dx            ; do y+2(ax)*x(dx)--(si)
  781.     add    bh,ah
  782.     adc    cx,dx
  783.     adc    bp,0
  784.  
  785.     mov    ax,si            ; get old "y+2"
  786. mulxy2: mul    di            ; do y+2(ax)*x+2(di)
  787.     add    cx,ax
  788.     adc    bp,dx
  789.  
  790. gotxy:    shl    bh,1
  791.     rcl    cx,1            ; bp,cx,bx
  792.     rcl    bp,1
  793.     jo    jmpxydone
  794.     shl    bh,1
  795.     rcl    cx,1
  796.     rcl    bp,1
  797.     jo    jmpxydone
  798.     shl    bh,1
  799.     rcl    cx,1
  800.     rcl    bp,1
  801.     jo    jmpxydone
  802.     add    bp,0
  803.     js    jmpxydone        ; if went neg have overflow also
  804.  
  805.     test    bl,negswt        ; ZERO IF NONE OR BOTH X , Y NEG
  806.     jz    signok            ; ONE IF ONLY ONE OF X OR Y IS NEG
  807.     not    cx            ; negate result
  808.     not    bp            ;  ...
  809.     add    cx,1            ;  ...
  810.     adc    bp,0            ;  ...
  811. signok:
  812.     shr    bl,1
  813.     shr    bl,1            ;MOVE "XSIGN"(X=NEG) DOWN TO "NEGSWT"
  814.     add    cx,cx            ;bp,CX = 2(X*Y)
  815.     adc    bp,bp
  816.  
  817.     add    cx,word ptr linity
  818.     adc    bp,word ptr linity+2    ; BP,CX = 2(X*Y)+B
  819.     mov    word ptr y,cx        ; save the new value of y
  820.     mov    word ptr y+2,bp     ;  ...
  821.     jns    jmpnit            ; NO SIGN, NOT NEG
  822.     inc    bl            ; INCREMENT NEGSWT
  823.     not    cx            ; negate
  824.     not    bp            ;  ...
  825.     add    cx,1            ;  ...
  826.     adc    bp,0            ;  ...
  827.  
  828. jmpnit: mov    ax,oldcolor        ; recall the old color
  829.     cmp    ax,k            ; check it against this iter
  830.     jl    short chkmaxit        ;  nope.  bypass periodicity check.
  831.     call    checkperiod        ; check for periodicity
  832.  
  833. chkmaxit:
  834.     dec    k            ; while (k < maxit)
  835.     jz    jmpxydone        ;  we done.
  836.     cmp    show_orbit,0        ; orbiting on?
  837.     jne    horbit            ;  yep.
  838.     jmp    nextit            ;go around again
  839.  
  840. jmpxydone:
  841.     jmp    xydone            ; DOES [(X*X)-(Y*Y)+P] BEFORE THE DEC.
  842.  
  843. horbit: push    bx            ; save my flags
  844.     push    bp            ; and registers
  845.     push    cx            ;  ...
  846.     mov    ax,-1            ; color for plot orbit
  847.     push    ax            ;  ...
  848.     push    word ptr y+2        ; co-ordinates for plot orbit
  849.     push    word ptr y        ;  ...
  850.     push    word ptr x+2        ;  ...
  851.     push    word ptr x        ;  ...
  852.     call    far ptr iplot_orbit    ; display the orbit
  853.     add    sp,5*2            ; clear out the parameters
  854.     pop    cx            ; restore registers
  855.     pop    bp            ;  ...
  856.     pop    bx            ; and flags
  857.     jmp    nextit            ; go around again
  858.  
  859. code32bit    endp
  860.  
  861.  
  862. calcmand_TEXT    ends
  863.  
  864.        end
  865.  
  866.