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