home *** CD-ROM | disk | FTP | other *** search
/ Photo CD Demo 1 / Demo.bin / fractint / fras1611.zip / CALCMAND.ASM < prev    next >
Assembly Source File  |  1991-01-31  |  26KB  |  869 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.     cmp    word ptr delmin+2,8    ; is 16-bit math good enough?
  233.     ja    kbddiskadj        ;  yes. test less often
  234.     mov    kbdcount,500        ;  no.    test more often
  235. kbddiskadj:
  236.     cmp    dotmode,11        ; disk video?
  237.     jne    quickkbd        ;  no, leave as is
  238.     shr    kbdcount,1        ; yes, reduce count
  239.     shr    kbdcount,1        ;  ...
  240. quickkbd:
  241.     call    far ptr keypressed    ; has a key been pressed?
  242.     cmp    ax,0            ;  ...
  243.     je    nokey            ; nope.  proceed
  244.     mov    kbdcount,0        ; make sure it goes negative again
  245.     cmp    ax,'o'                  ; orbit toggle hit?
  246.     je    orbitkey        ;  yup.  show orbits
  247.     cmp    ax,'O'                  ; orbit toggle hit?
  248.     jne    keyhit            ;  nope.  normal key.
  249. orbitkey:
  250.     call    far ptr getakey     ; read the key for real
  251.     mov    ax,1            ; reset orbittoggle = 1 - orbittoggle
  252.     sub    ax,show_orbit        ;  ...
  253.     mov    show_orbit,ax        ;  ...
  254.     jmp    short nokey        ; pretend no key was hit
  255. keyhit: mov    ax,-1            ; return with -1
  256.     mov    color,ax        ; set color to -1
  257.     UNFRAME <si,di>         ; pop stack frame
  258.     ret                ; bail out!
  259.  
  260. nokey:
  261.     cmp    show_orbit,0        ; is orbiting on?
  262.     jne    no16bitcode        ;  yup.  slow down.
  263.     cmp    cpu,386         ; are we on a 386?
  264.     je    short code386bit    ;  YAY!! 386-class speed!
  265. ;;    cmp    word ptr delmin+2,1    ; OK, we're desperate.  16 bits OK?
  266.     cmp    word ptr delmin+2,8    ; OK, we're desperate.  16 bits OK?
  267.     ja    yes16bitcode        ;  YAY!  16-bit speed!
  268. no16bitcode:
  269.     call    near ptr code32bit    ; BOO!! nap time.  Full 32 bit math
  270.     jmp    kloopend        ;  bypass the 386-specific code.
  271. yes16bitcode:
  272.     call    near ptr code16bit    ; invoke the 16-bit version
  273.     jmp    kloopend        ;  bypass the 386-specific code.
  274.  
  275. .386                    ; 386-specific code starts here
  276.  
  277. code386bit:
  278. ;;    cmp    word ptr delmin+2,3    ; is 16-bit math good enough?
  279.     cmp    word ptr delmin+2,8    ; is 16-bit math good enough?
  280.     jbe    code386_32        ; nope, go do 32 bit stuff
  281. IFDEF ??version
  282.     jmp    code386_32        ; TASM screws up IMUL EBX,EBX!!
  283. ENDIF
  284.  
  285.     ; 16 bit on 386, now we are really gonna move
  286.     movsx    esi,word ptr x+2    ; use SI for X
  287.     movsx    edi,word ptr y+2    ; use DI for Y
  288.     push    ebp
  289.     mov    ebp,-1
  290.     shl    ebp,FUDGEFACTOR-1
  291.     mov    cx,FUDGEFACTOR-16
  292.  
  293. kloop386_16:   ; cx=bitshift-16, ebp=overflow.mask
  294.  
  295.     mov    ebx,esi         ; compute (x * x)
  296.     imul    ebx,ebx         ;  ...
  297.     test    ebx,ebp         ; say, did we overflow? <V20-compat>
  298.     jnz    short end386_16     ;  (oops.  We done.)
  299.     shr    ebx,cl            ; get result down to 16 bits
  300.  
  301.     mov    edx,edi         ; compute (y * y)
  302.     imul    edx,edx         ;  ...
  303.     test    edx,ebp         ; say, did we overflow? <V20-compat>
  304.     jnz    short end386_16     ;  (oops.  We done.)
  305.     shr    edx,cl            ; get result down to 16 bits
  306.  
  307.     mov    ax,bx            ; compute (x*x - y*y) / fudge
  308.     sub    bx,dx            ;  for the next iteration
  309.  
  310.     add    ax,dx            ; compute (x*x + y*y) / fudge
  311.  
  312.     cmp    ax,word ptr lm+2    ; while (xx+yy < lm)
  313.     jae    short end386_16     ;  ...
  314.  
  315.     imul    edi,esi         ; compute (y * x)
  316.     shl    edi,1            ; ( * 2 / fudge)
  317.     sar    edi,cl
  318.     add    di,word ptr linity+2    ; (2*y*x) / fudge + linity
  319.     movsx    edi,di            ; save as y
  320.  
  321.     add    bx,word ptr linitx+2    ; (from above) (x*x - y*y)/fudge + linitx
  322.     movsx    esi,bx            ; save as x
  323.  
  324.     mov    ax,oldcolor        ; recall the old color
  325.     cmp    ax,k            ; check it against this iter
  326.     jge    short chkpd386_16    ;  yup.  do periodicity check.
  327. nonmax386_16:
  328.     dec    k            ; while (k < maxit)
  329.     jnz    short kloop386_16    ; try, try again
  330. end386_16:
  331.     pop    ebp
  332.     jmp    kloopend        ; we done
  333.  
  334. chkpd386_16:
  335.     mov    ax,k            ; set up to test for save-time
  336.     test    ax,savedand        ; save on 0, check on anything else
  337.     jz    short chksv386_16    ;  time to save a new "old" value
  338.     mov    ax,si            ; load up x
  339.     xor    ax,word ptr savedx+2    ; does X match?
  340.     test    ax,word ptr savedmask+2 ;  truncate to appropriate precision
  341.     jne    short nonmax386_16    ;  nope.  forget it.
  342.     mov    ax,di            ; now test y
  343.     xor    ax,word ptr savedy+2    ; does Y match?
  344.     test    ax,word ptr savedmask+2 ;  truncate to appropriate precision
  345.     jne    short nonmax386_16    ;  nope.  forget it.
  346.     mov    period,1        ; note that we have found periodicity
  347.     mov    k,0            ; pretend maxit reached
  348.     jmp    short end386_16
  349. chksv386_16:
  350.     mov    word ptr savedx+2,si    ; save x
  351.     mov    word ptr savedy+2,di    ; save y
  352.     dec    savedincr        ; time to change the periodicity?
  353.     jnz    short nonmax386_16    ;  nope.
  354.     shl    savedand,1        ; well then, let's try this one!
  355.     inc    savedand        ;  (2**n -1)
  356.     mov    savedincr,4        ; and reset the increment flag
  357.     jmp    short nonmax386_16
  358.  
  359.     ; 32bit on 386:
  360. code386_32:
  361.     mov    esi,x            ; use ESI for X
  362.     mov    edi,y            ; use EDI for Y
  363.  
  364. ;    This is the main processing loop.  Here, every T-state counts...
  365.  
  366. kloop:                    ; for (k = 0; k <= maxit; k++)
  367.  
  368.     mov    eax,esi         ; compute (x * x)
  369.     imul    esi            ;  ...
  370.     shrd    eax,edx,FUDGEFACTOR    ; ( / fudge)
  371.     shr    edx,FUDGEFACTOR-1    ; (complete 64-bit shift and check
  372.     jne    short kloopend1     ; bail out if too high
  373.     mov    ebx,eax         ; save this for below
  374.  
  375.     mov    eax,edi         ; compute (y * y)
  376.     imul    edi            ;  ...
  377.     shrd    eax,edx,FUDGEFACTOR    ; ( / fudge)
  378.     shr    edx,FUDGEFACTOR-1    ; (complete 64-bit shift and check
  379.     jne    short kloopend1     ; bail out if too high
  380.  
  381.     mov    ecx,ebx         ; compute (x*x - y*y) / fudge
  382.     sub    ebx,eax         ;  for the next iteration
  383.  
  384.     add    ecx,eax         ; compute (x*x + y*y) / fudge
  385.     cmp    ecx,lm            ; while (lr < lm)
  386.     jae    short kloopend1     ;  ...
  387.  
  388.     mov    eax,edi         ; compute (y * x)
  389.     imul    esi            ;  ...
  390.     shrd    eax,edx,FUDGEFACTOR-1    ;  ( * 2 / fudge)
  391.     add    eax,linity        ;  (above) + linity
  392.     mov    edi,eax         ;  save this as y
  393.  
  394. ;    (from the earlier code)     ; compute (x*x - y*y) / fudge
  395.     add    ebx,linitx        ;    + linitx
  396.     mov    esi,ebx         ; save this as x
  397.  
  398.     mov    ax,oldcolor        ; recall the old color
  399.     cmp    ax,k            ; check it against this iter
  400.     jge    short chkperiod1
  401. nonmax1:
  402.     dec    k            ; while (k < maxit) (dec to 0 is faster)
  403.     jnz    short kloop        ; while (k < maxit) ...
  404. kloopend1:
  405.     jmp    short kloopend32    ; we done.
  406.  
  407. chkperiod1:
  408.     mov    eax,esi
  409.     xor    eax,savedx
  410.     test    eax,savedmask
  411.     jnz    short chksave1
  412.     mov    eax,edi
  413.     xor    eax,savedy
  414.     test    eax,savedmask
  415.     jnz    short chksave1
  416.     mov    period,1        ; note that we have found periodicity
  417.     mov    k,0            ; pretend maxit reached
  418.     jmp    short kloopend1
  419. chksave1:
  420.     mov    ax,k
  421.     test    ax,savedand
  422.     jne    short nonmax1
  423.     mov    savedx,esi
  424.     mov    savedy,edi
  425.     dec    savedincr        ; time to change the periodicity?
  426.     jnz    short nonmax1        ;  nope.
  427.     shl    savedand,1        ; well then, let's try this one!
  428.     inc    savedand        ;  (2**n -1)
  429.     mov    savedincr,4        ; and reset the increment flag
  430.     jmp    short nonmax1
  431.  
  432. kloopend32:
  433.  
  434. .8086                    ; 386-specific code ends here
  435.  
  436. kloopend:
  437.     cmp    orbit_ptr,0        ; any orbits to clear?
  438.     je    noorbit2        ;  nope.
  439.     call    far ptr scrub_orbit    ; clear out any old orbits
  440. noorbit2:
  441.  
  442.     mov    ax,k            ; set old color
  443.     sub    ax,10            ; minus 10, for safety
  444.     mov    oldcolor,ax        ; and save it as the "old" color
  445.     mov    ax,maxit        ; compute color
  446.     sub    ax,k            ;  (first, re-compute "k")
  447.     sub    kbdcount,ax        ; adjust the keyboard count
  448.     cmp    ax,1            ; convert any "outlier" region
  449.     jge    short coloradjust1    ;  (where abs(x) > 2 or abs(y) > 2)
  450.     mov    ax,1            ;   to look like we ran through
  451. coloradjust1:                ;    at least one loop.
  452.     mov    realcolor,ax        ; result before adjustments
  453.     cmp    ax,maxit        ; did we max out on iterations?
  454.     jne    short wedone        ;  nope.
  455.     mov    oldcolor,ax        ; set "oldcolor" to maximum
  456.     cmp    inside,0        ; is "inside" >= 0?
  457.     jl    wedone            ;  nope.  leave it at "maxit"
  458.     mov    ax,inside        ; reset max-out color to default
  459.     cmp    periodicitycheck,0    ; show periodicity matches?
  460.     jge    wedone            ;  nope.
  461.     mov    al,period        ;  reset color to periodicity flag
  462. wedone:                 ;
  463.     mov    color,ax        ; save the color result
  464.     UNFRAME <si,di>         ; pop stack frame
  465.     ret                ; and return with color
  466.  
  467. calcmandasm endp
  468.  
  469.  
  470. ; ******************** Function code16bit() *****************************
  471. ;
  472. ;    Performs "short-cut" 16-bit math where we can get away with it.
  473. ;
  474.  
  475. code16bit    proc    near
  476.  
  477.     mov    si,word ptr x+2     ; use SI for X
  478.     mov    di,word ptr y+2     ; use DI for Y
  479.  
  480. start16bit:
  481.     mov    ax,si            ; compute (x * x)
  482.     imul    si            ;  ...
  483.     cmp    dx,0            ; say, did we overflow? <V20-compat>
  484.     jl    end16bit        ;  (oops.  We done.)
  485.     mov    cx,32-FUDGEFACTOR    ; ( / fudge)
  486. loop16bit1:
  487.     shl    ax,1            ;  ...
  488.     rcl    dx,1            ;  ...
  489.     jo    end16bit        ;  (oops.  overflow)
  490.     loop    loop16bit1        ;  ...
  491.     mov    bx,dx            ; save this for a tad
  492.  
  493.     mov    ax,di            ; compute (y * y)
  494.     imul    di            ;  ...
  495.     cmp    dx,0            ; say, did we overflow? <V20-compat>
  496.     jl    end16bit        ;  (oops.  We done.)
  497.     mov    cx,32-FUDGEFACTOR    ; ( / fudge)
  498. loop16bit2:
  499.     shl    ax,1            ;  ...
  500.     rcl    dx,1            ;  ...
  501.     jo    end16bit        ;  (oops.  overflow)
  502.     loop    loop16bit2        ;  ...
  503.  
  504.     mov    cx,bx            ; compute (x*x - y*y) / fudge
  505.     sub    bx,dx            ;  for the next iteration
  506.  
  507.     add    cx,dx            ; compute (x*x + y*y) / fudge
  508.     jo    end16bit        ; bail out if too high
  509.     js    end16bit        ;  ...
  510.  
  511.     cmp    cx,word ptr lm+2    ; while (xx+yy < lm)
  512.     jae    end16bit        ;  ...
  513.  
  514.     mov    ax,di            ; compute (y * x)
  515.     imul    si            ;  ...
  516.     mov    cx,33-FUDGEFACTOR    ; ( * 2 / fudge)
  517. loop16bit3:
  518.     shl    ax,1            ;  ...
  519.     rcl    dx,1            ;  ...
  520.     loop    loop16bit3        ;  ...
  521.     add    dx,word ptr linity+2    ; (2*y*x) / fudge + linity
  522.     mov    di,dx            ; save as y
  523.  
  524.     add    bx,word ptr linitx+2    ; (from above) (x*x - y*y)/fudge + linitx
  525.     mov    si,bx            ; save as x
  526.  
  527.     mov    ax,oldcolor        ; recall the old color
  528.     cmp    ax,k            ; check it against this iter
  529.     jl    short nonmax3        ;  nope.  bypass periodicity check.
  530.     mov    word ptr x+2,si     ; save x for periodicity check
  531.     mov    word ptr y+2,di     ; save y for periodicity check
  532.     call    checkperiod        ; check for periodicity
  533. nonmax3:
  534.  
  535.     dec    k            ; while (k < maxit)
  536.     jz    end16bit        ;  we done.
  537.  
  538.     jmp    start16bit        ; try, try again.
  539.  
  540. end16bit:                ; we done.
  541.     ret
  542. code16bit    endp
  543.  
  544.  
  545. ;    The following routine checks for periodic loops (a known
  546. ;    method of decay inside "Mandelbrot Lake", and an easy way to
  547. ;    bail out of "lake" points quickly).  For speed, only the
  548. ;    high-order sixteen bits of X and Y are checked for periodicity.
  549. ;    For accuracy, this routine is only fired up if the previous pixel
  550. ;    was in the lake (which means that the FIRST "wet" pixel was
  551. ;    detected by the dull-normal maximum iteration process).
  552.  
  553. checkperiod    proc near        ; periodicity check
  554.     mov    ax,k            ; set up to test for save-time
  555.     test    ax,savedand        ; save on 0, check on anything else
  556.     jz    checksave        ;  time to save a new "old" value
  557.     mov    dx,word ptr x+2     ; load up x
  558.     and    dx,word ptr savedmask+2 ;  truncate to appropriate precision
  559.     cmp    dx,word ptr savedx+2    ; does X match?
  560.     jne    checkdone        ;  nope.  forget it.
  561.     mov    ax,word ptr x        ; load up x
  562.     and    ax,word ptr savedmask    ;  truncate to appropriate precision
  563.     cmp    ax,word ptr savedx    ; does X match?
  564.     jne    checkdone        ;  nope.  forget it.
  565.     mov    dx,word ptr y+2     ; now test y
  566.     and    dx,word ptr savedmask+2 ;  truncate to appropriate precision
  567.     cmp    dx,word ptr savedy+2    ; does Y match?
  568.     jne    checkdone        ;  nope.  forget it.
  569.     mov    ax,word ptr y        ; load up y
  570.     and    ax,word ptr savedmask    ;  truncate to appropriate precision
  571.     cmp    ax,word ptr savedy    ; does Y match?
  572.     jne    checkdone        ;  nope.  forget it.
  573.     mov    period,1        ; note that we have found periodicity
  574.     mov    k,1            ; pretend maxit reached
  575. checksave:
  576.     mov    dx,word ptr x+2     ; load up x
  577.     and    dx,word ptr savedmask+2 ;  truncate to appropriate precision
  578.     mov    word ptr savedx+2,dx    ;  and save it
  579.     mov    ax,word ptr x        ; load up x
  580.     and    ax,word ptr savedmask    ;  truncate to appropriate precision
  581.     mov    word ptr savedx,ax    ;  and save it
  582.     mov    dx,word ptr y+2     ; load up y
  583.     and    dx,word ptr savedmask+2 ;  truncate to appropriate precision
  584.     mov    word ptr savedy+2,dx    ;  and save it
  585.     mov    ax,word ptr y        ; load up y
  586.     and    ax,word ptr savedmask    ;  truncate to appropriate precision
  587.     mov    word ptr savedy,ax    ;  and save it
  588.     dec    savedincr        ; time to change the periodicity?
  589.     jnz    checkdone        ;  nope.
  590.     shl    savedand,1        ; well then, let's try this one!
  591.     inc    savedand        ;  (2**n -1)
  592.     mov    savedincr,4        ; and reset the increment flag
  593. checkdone:
  594.     ret                ; we done.
  595. checkperiod    endp
  596.  
  597.  
  598. ; ******************** Function code32bit() *****************************
  599. ;
  600. ;    Perform the 32-bit logic required using 16-bit logic
  601. ;
  602. ;    New twice as fast logic,
  603. ;       Courtesy of Bill Townsend and Mike Gelvin (CIS:73337,520)
  604.  
  605. square    macro    doneaddr
  606.     local    donejp,skip1y,skip2y,squrey,allzero
  607.     sub    bh,bh
  608.     sub    cx,cx
  609.     mov    bp,cx
  610.     or    ax,ax            ;GET LOW HWORD
  611.     jz    skip1y            ;LOW HWORD IS ZERO, ONLY DO HIGH HWORD
  612.     mov    si,ax
  613.     mul    ax            ;SQUARE LOW HWORD - B*B
  614.     mov    bh,dh
  615.     or    di,di            ;TEST HIGH HWORD
  616.     jz    squrey            ;HIGH HWORD ZERO, SKIP THE FOLLOWING
  617.     mov    ax,di
  618.     mul    si            ;GET MIDDLE PART - A*B
  619.     add    bh,ah            ;add
  620.     adc    cx,dx            ; A*B
  621.     adc    bp,0            ;  twice
  622.     add    bh,ah            ;   -
  623.     adc    cx,dx            ;    -
  624.     adc    bp,0            ;     SI,CX,BH = 2(A*B)+(B*B)
  625.     jmp    short skip2y
  626. donejp: JMP    doneaddr        ;M0=DONEJP
  627. skip1y: or    di,di            ;M1=SKIP1Y
  628.     jz    allzero
  629. skip2y: mov    ax,di            ;M2=SKIP2Y
  630.     mul    ax            ;SQUARE HIGH HWORD - A*A
  631.     add    cx,ax
  632.     adc    bp,dx
  633. squrey: shl    bh,1            ;M3=SQUREY
  634.     rcl    cx,1
  635.     rcl    bp,1
  636. ;    jo    donejp            ; squaring a +ve, top bit known off
  637.     shl    bh,1
  638.     rcl    cx,1
  639.     rcl    bp,1
  640. ;    jo    donejp            ; squaring a +ve, 2nd bit known off
  641.     shl    bh,1
  642.     rcl    cx,1
  643.     rcl    bp,1
  644.     jo    donejp
  645.     add    bp,0
  646.     js    donejp            ; if went neg have overflow also
  647. allzero:
  648.  
  649.     endm    ;#EM
  650.  
  651.  
  652. code32bit    proc near
  653. ;
  654. ; BL IS USED FOR THE FLAGS,TWO LOW ORDER BITS ARE "NEGSWT", THIRD BIT IS XSIGN
  655. ;   XSIGN IS ONE IF THE SIGN OF THE NEW "X" IS NEGATIVE
  656. ;   NEGSWT IS ONE IF EITHER "X" OR "Y" ARE NEGATIVE, BUT NOT BOTH NEGATIVE
  657. ;
  658. xsign    equ    00000100b
  659. negswt    equ    00000001b
  660.  
  661.     push    bp
  662.     sub    si,si            ; make a zero __________
  663.     mov    bx,si            ; BOTH XSIGN & NEGSWT START OFF ZERO
  664.     mov    di,1            ; MAKE A ONE
  665.     mov    cx,word ptr y        ; cx=y
  666.     mov    bp,word ptr y+2     ; bp=y+2
  667.     cmp    bp,si            ; y>0?
  668.     jns    qnotng            ;  yup
  669.     not    cx            ; nope, negate it
  670.     not    bp            ; ...
  671.     add    cx,di            ; ...
  672.     adc    bp,si            ; ...
  673.     inc    bl            ; INCREMENT negswt
  674. qnotng: mov    ax,word ptr x
  675.     mov    dx,word ptr x+2
  676.     cmp    dx,si            ; x>0?
  677.     jns    pnotng            ;  yup
  678.     not    ax            ; nope, negate it
  679.     not    dx            ; ...
  680.     add    ax,di            ; ...
  681.     adc    dx,si            ; ...
  682.     inc    bl            ; INCREMENT negswt
  683. pnotng: mov    word ptr absx,ax    ; save unsigned x
  684.     mov    word ptr absx+2,dx    ; ...
  685.  
  686. ;    iteration loop
  687. nextit: push    bp            ; save y on stack
  688.     push    cx            ; ...
  689.     mov    ax,cx            ; load y for squaring
  690.     mov    di,bp            ;  ...
  691.     square    done2            ; square di,ax & br to jmp_done2 if overflow
  692.     push    cx            ; return results in bp,cx
  693.     push    bp            ; save y*y on stack
  694.     mov    ax,word ptr absx    ; load x for squaring
  695.     mov    di,word ptr absx+2    ;  ...
  696.     square    done            ; square di,ax & br to done if overflow
  697.     mov    di,bp            ; return results in bp,cx
  698.     mov    si,cx            ; save x*x in di,si
  699.     pop    dx            ; unstack y*y into dx,ax
  700.     pop    ax            ;  ...
  701.     ADD    cx,ax            ; calc y*y + x*x
  702.     ADC    bp,dx            ;  ...
  703.     jo    done2            ; overflow?
  704.     js    done2            ; overflow?
  705.     cmp    bp,word ptr lm+2    ; while (lr < lm)
  706.     jb    short nextxy        ;  keep going
  707.     jmp    short done2        ; magnitude is past limit, bailout
  708. done:    add    sp,4            ; discard saved value of "y*y"
  709. done2:    add    sp,4            ; discard saved value of "y"
  710. xydone: pop    bp            ; restore saved bp
  711.     ret
  712.  
  713. ;---------------------------------------------------------------------------
  714. ;y=nz, x=nz, x+2=z
  715. x2is0:    pop    ax            ; get old y+2
  716.     or    ax,ax            ; test y+2
  717.     jz    jgotxy
  718.     mul    si            ; do y+2(ax)*x(si)
  719.     add    bh,ah
  720.     adc    cx,dx
  721.     adc    bp,0
  722. jgotxy: jmp    gotxy
  723.  
  724. ;x=z
  725. xis0:    or    di,di            ; test x+2
  726.     jz    popy_2            ; zero all done, around again please
  727.     mul    di            ; do y(ax)*x+2(di)
  728.     add    bh,ah
  729.     adc    cx,dx
  730.     adc    bp,0
  731.     pop    ax            ; get old y+2
  732.     or    ax,ax
  733.     jz    gotxy
  734. mulxyj: jmp    mulxy2
  735.  
  736. ;x=z, x+2=z
  737. popy_2: pop    ax            ; discard old y+2
  738.     sub    cx,cx            ; x==0, so x*y is zero
  739.     mov    bp,cx            ; ...
  740.     jmp    signok            ; go add b and store new y
  741.  
  742. nextxy: sub    si,ax            ; subtract y*y from x*x
  743.     sbb    di,dx            ;  ...
  744.     add    si,word ptr linitx    ; add "A"
  745.     adc    di,word ptr linitx+2    ;  ...
  746.     mov    word ptr x,si        ; store new x = x*x+y*y+a
  747.     mov    word ptr x+2,di     ;  ...
  748.     jns    swtxxx            ; NO SIGN, NOT NEG
  749.     or    bl,xsign        ; REMEMBER THE NEW "X" IS NEGATIVE
  750.     not    si            ; make it positive in register
  751.     not    di            ;  ...
  752.     add    si,1            ;  ...
  753.     adc    di,0            ;  ...
  754. swtxxx: xchg    di,word ptr absx+2    ; save the new x, load the old
  755.     xchg    si,word ptr absx    ;  ...
  756.  
  757.     ; now calculate x*y
  758.     sub    bh,bh            ; zero some registers
  759.     sub    cx,cx            ;  ...
  760.     mov    bp,cx            ;  ...
  761.     pop    ax            ; get old "y"
  762.     or    ax,ax            ; test "y"
  763.     jz    yis0            ; br if "y" is zero
  764.     or    si,si            ; test old "x"
  765.     jz    xis0            ; br if "x" is zero
  766.     mov    bp,ax            ; save old "y"
  767.     mul    si            ; do y(ax)*x(si)
  768.     mov    ax,bp            ; get old "y" again
  769.     mov    bp,cx            ; RE-ZERO BP
  770.     mov    bh,dh            ; save/set low hword(bx)
  771.     or    di,di            ; test old "x+2"
  772.     jz    x2is0
  773.     mul    di            ; do y(ax)*x+2(di)
  774.     add    bh,ah
  775.     adc    cx,dx
  776.     adc    bp,0
  777.  
  778. yis0:    pop    ax            ; get old "y+2"
  779.     or    ax,ax            ; test y+2
  780.     jz    gotxy            ; br if y+2 is zero
  781.     mov    dx,si            ; dx="x"
  782.     mov    si,ax            ; si=save "y+2"
  783.     mul    dx            ; do y+2(ax)*x(dx)--(si)
  784.     add    bh,ah
  785.     adc    cx,dx
  786.     adc    bp,0
  787.  
  788.     mov    ax,si            ; get old "y+2"
  789. mulxy2: mul    di            ; do y+2(ax)*x+2(di)
  790.     add    cx,ax
  791.     adc    bp,dx
  792.  
  793. gotxy:    shl    bh,1
  794.     rcl    cx,1            ; bp,cx,bx
  795.     rcl    bp,1
  796.     jo    jmpxydone
  797.     shl    bh,1
  798.     rcl    cx,1
  799.     rcl    bp,1
  800.     jo    jmpxydone
  801.     shl    bh,1
  802.     rcl    cx,1
  803.     rcl    bp,1
  804.     jo    jmpxydone
  805.     add    bp,0
  806.     js    jmpxydone        ; if went neg have overflow also
  807.  
  808.     test    bl,negswt        ; ZERO IF NONE OR BOTH X , Y NEG
  809.     jz    signok            ; ONE IF ONLY ONE OF X OR Y IS NEG
  810.     not    cx            ; negate result
  811.     not    bp            ;  ...
  812.     add    cx,1            ;  ...
  813.     adc    bp,0            ;  ...
  814. signok:
  815.     shr    bl,1
  816.     shr    bl,1            ;MOVE "XSIGN"(X=NEG) DOWN TO "NEGSWT"
  817.     add    cx,cx            ;bp,CX = 2(X*Y)
  818.     adc    bp,bp
  819.  
  820.     add    cx,word ptr linity
  821.     adc    bp,word ptr linity+2    ; BP,CX = 2(X*Y)+B
  822.     mov    word ptr y,cx        ; save the new value of y
  823.     mov    word ptr y+2,bp     ;  ...
  824.     jns    jmpnit            ; NO SIGN, NOT NEG
  825.     inc    bl            ; INCREMENT NEGSWT
  826.     not    cx            ; negate
  827.     not    bp            ;  ...
  828.     add    cx,1            ;  ...
  829.     adc    bp,0            ;  ...
  830.  
  831. jmpnit: mov    ax,oldcolor        ; recall the old color
  832.     cmp    ax,k            ; check it against this iter
  833.     jl    short chkmaxit        ;  nope.  bypass periodicity check.
  834.     call    checkperiod        ; check for periodicity
  835.  
  836. chkmaxit:
  837.     dec    k            ; while (k < maxit)
  838.     jz    jmpxydone        ;  we done.
  839.     cmp    show_orbit,0        ; orbiting on?
  840.     jne    horbit            ;  yep.
  841.     jmp    nextit            ;go around again
  842.  
  843. jmpxydone:
  844.     jmp    xydone            ; DOES [(X*X)-(Y*Y)+P] BEFORE THE DEC.
  845.  
  846. horbit: push    bx            ; save my flags
  847.     push    bp            ; and registers
  848.     push    cx            ;  ...
  849.     mov    ax,-1            ; color for plot orbit
  850.     push    ax            ;  ...
  851.     push    word ptr y+2        ; co-ordinates for plot orbit
  852.     push    word ptr y        ;  ...
  853.     push    word ptr x+2        ;  ...
  854.     push    word ptr x        ;  ...
  855.     call    far ptr iplot_orbit    ; display the orbit
  856.     add    sp,5*2            ; clear out the parameters
  857.     pop    cx            ; restore registers
  858.     pop    bp            ;  ...
  859.     pop    bx            ; and flags
  860.     jmp    nextit            ; go around again
  861.  
  862. code32bit    endp
  863.  
  864.  
  865. calcmand_TEXT    ends
  866.  
  867.        end
  868.  
  869.