home *** CD-ROM | disk | FTP | other *** search
/ The World of Computer Software / World_Of_Computer_Software-02-385-Vol-1of3.iso / w / winsr173.zip / CALCMAND.ASM < prev    next >
Assembly Source File  |  1991-11-10  |  27KB  |  813 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. ; History since Fractint 16.0
  43. ;  (See comments with CJLT in them)
  44. ;  CJLT=Chris Lusby Taylor who has...
  45. ;
  46. ;   1. Speeded up 16 bit on 16 bit CPU
  47. ;    Minor changes, notably prescaling to fg14 before multiplying
  48. ;    instead of scaling the answer.
  49. ;    Also, I added overflow detection after adding linit, since it
  50. ;    seems this could overflow.  
  51. ;    Overall effect is about 10% faster on 386 with debugflag=8088
  52. ;   2. Speeded up 32 bit on 16 bit CPU
  53. ;    The macro `square' is totally rewritten, as is the logic for 2xy,
  54. ;    by prescaling x and y to fg31, not fg29. This allows us to do a
  55. ;    32 bit multiply in 3, not 4, 16 bit chunks while retaining full
  56. ;    fg29 accuracy.
  57. ;       Also, I removed lots of well-meaning but ineffective code handling
  58. ;    special cases of zeros and tidied up the handling of negative numbers,
  59. ;    so the routine is quite a bit shorter now and overall throughput of
  60. ;    Mandel is now over 40% faster on a 386 with debugflag=8088.
  61. ;       By the way, I was tempted to go the whole hog and replace x*x-y*y
  62. ;    by (x+y)*(x-y) to reduce 4 16-bit multiplys to 3, but it makes
  63. ;    escape detection a bit trickier. Another time, maybe.
  64. ;
  65.  
  66. ;             required for compatibility if Turbo ASM
  67. IFDEF ??version
  68. MASM51
  69. QUIRKS
  70. ENDIF
  71.  
  72. .MODEL    medium,c
  73. DGROUP          group   _DATA,_DATA2
  74.  
  75. .8086
  76.  
  77.     ; these must NOT be in any segment!!
  78.     ; this get's rid of TURBO-C fixup errors
  79.  
  80.     extrn    keypressed:far        ; this routine is in 'general.asm'
  81.     extrn    getakey:far        ; this routine is in 'general.asm'
  82.     extrn    iplot_orbit:far     ; this routine is in 'calcfrac.c'
  83.     extrn    scrub_orbit:far     ; this routine is in 'calcfrac.c'
  84.  
  85. _DATA2        segment DWORD PUBLIC 'DATA'
  86.  
  87. FUDGEFACTOR    equ    29        ; default (non-potential) fudgefactor
  88.  
  89. ; ************************ External variables *****************************
  90.  
  91.     extrn    fractype:word        ; == 0 if Mandelbrot set, else Julia
  92.     extrn    inside:word        ; "inside" color, normally 1 (blue)
  93.     extrn    creal:dword, cimag:dword ; Julia Set Constant
  94.     extrn    delmin:dword        ; min increment - precision required
  95.     extrn    maxit:word        ; maximum iterations
  96.     extrn    lm:dword        ; magnitude bailout limit
  97.  
  98.     extrn    row:word, col:word    ; current pixel to calc
  99.     extrn    color:word        ; color calculated for the pixel
  100.     extrn    realcolor:word        ; color before inside,etc adjustments
  101.  
  102.     extrn    reset_periodicity:word    ; nonzero if to be reset
  103.     extrn    kbdcount:word        ; keyboard counter
  104.  
  105.     extrn    cpu:word        ; cpu type: 86, 186, 286, or 386
  106.     extrn    dotmode:word
  107.  
  108.     extrn    show_orbit:word     ; "show-orbit" flag
  109.     extrn    orbit_ptr:word        ; "orbit pointer" flag
  110.     extrn    periodicitycheck:word    ; no periodicity if zero
  111.  
  112.     public    linitx,linity        ; caller sets these
  113.     public    savedmask        ; caller sets this
  114.  
  115. ; ************************ Internal variables *****************************
  116.  
  117.         align    4
  118. x        dd    0        ; temp value: x
  119. y        dd    0        ; temp value: y
  120. absx        dd    0        ; temp value: abs(x)
  121. linitx        dd    0        ; initial value, set by calcfrac
  122. linity        dd    0        ; initial value, set by calcfrac
  123. savedmask    dd    0        ; saved values mask
  124. savedx        dd    0        ; saved values of X and Y iterations
  125. savedy        dd    0        ;  (for periodicity checks)
  126. k        dw    0        ; iteration countdown counter
  127. oldcolor    dw    0        ; prior pixel's escape time k value
  128. savedand    dw    0        ; AND value for periodicity checks
  129. savedincr    dw    0        ; flag for incrementing AND value
  130. period        db    0        ; periodicity, if in the lake
  131.  
  132. _DATA2        ends
  133.  
  134. .CODE
  135.  
  136. ; ***************** Function calcmandasm() **********************************
  137.  
  138.     public    calcmandasm
  139.  
  140. FRAME    MACRO regs
  141.     push    bp
  142.     mov    bp, sp
  143.     IRP    reg, <regs>
  144.       push    reg
  145.       ENDM
  146.     ENDM
  147.  
  148. UNFRAME MACRO regs
  149.     IRP    reg, <regs>
  150.       pop reg
  151.       ENDM
  152.     pop bp
  153.     ENDM
  154.  
  155. calcmandasm    proc
  156.     FRAME    <di,si>         ; std frame, for TC++ overlays
  157.     sub    ax,ax            ; clear ax
  158.     cmp    periodicitycheck,ax    ; periodicity checking disabled?
  159.     je    initoldcolor        ;  yup, set oldcolor 0 to disable it
  160.     cmp    reset_periodicity,ax    ; periodicity reset?
  161.     je    short initparms     ; inherit oldcolor from prior invocation
  162.     mov    ax,maxit        ; yup.    reset oldcolor to maxit-250
  163.     sub    ax,250            ; (avoids slowness at high maxits)
  164. initoldcolor:
  165.     mov    oldcolor,ax        ; reset oldcolor
  166.  
  167. initparms:
  168.     mov    ax,word ptr creal    ; initialize x == creal
  169.     mov    dx,word ptr creal+2    ;  ...
  170.     mov    word ptr x,ax        ;  ...
  171.     mov    word ptr x+2,dx     ;  ...
  172.  
  173.     mov    ax,word ptr cimag    ; initialize y == cimag
  174.     mov    dx,word ptr cimag+2    ;  ...
  175.     mov    word ptr y,ax        ;  ...
  176.     mov    word ptr y+2,dx     ;  ...
  177.  
  178.     mov    ax,maxit        ; setup k = maxit
  179.     inc    ax            ; (+ 1)
  180.     mov    k,ax            ;  (decrementing to 0 is faster)
  181.  
  182.     cmp    fractype,1        ; julia or mandelbrot set?
  183.     je    short dojulia        ; julia set - go there
  184.  
  185. ;    (Tim wants this code changed so that, for the Mandelbrot,
  186. ;    Z(1) = (x + iy) + (a + ib).  Affects only "fudged" Mandelbrots.
  187. ;    (for the "normal" case, a = b = 0, and this works, too)
  188. ;    cmp    word ptr x,0        ; Mandelbrot shortcut:
  189. ;    jne    short doeither        ;  if creal = cimag = 0,
  190. ;    cmp    word ptr x+2,0        ; the first iteration can be emulated.
  191. ;    jne    short doeither        ;  ...
  192. ;    cmp    word ptr y,0        ;  ...
  193. ;    jne    short doeither        ;  ...
  194. ;    cmp    word ptr y+2,0        ;  ...
  195. ;    jne    short doeither        ;  ...
  196. ;    dec    k            ; we know the first iteration passed
  197. ;    mov    dx,word ptr linitx+2    ; copy x = linitx
  198. ;    mov    ax,word ptr linitx    ;  ...
  199. ;    mov    word ptr x+2,dx     ;  ...
  200. ;    mov    word ptr x,ax        ;  ...
  201. ;    mov    dx,word ptr linity+2    ; copy y = linity
  202. ;    mov    ax,word ptr linity    ;  ...
  203. ;    mov    word ptr y+2,dx     ;  ...
  204. ;    mov    word ptr y,ax        ;  ...
  205.  
  206.     dec    k            ; we know the first iteration passed
  207.     mov    dx,word ptr linitx+2    ; add x += linitx
  208.     mov    ax,word ptr linitx    ;  ...
  209.     add    word ptr x,ax        ;  ...
  210.     adc    word ptr x+2,dx     ;  ...
  211.     mov    dx,word ptr linity+2    ; add y += linity
  212.     mov    ax,word ptr linity    ;  ...
  213.     add    word ptr y,ax        ;  ...
  214.     adc    word ptr y+2,dx     ;  ...
  215.     jmp    short doeither        ; branch around the julia switch
  216.  
  217. dojulia:                ; Julia Set initialization
  218.                     ; "fudge" Mandelbrot start-up values
  219.     mov    ax,word ptr x        ; switch x with linitx
  220.     mov    dx,word ptr x+2     ;  ...
  221.     mov    bx,word ptr linitx    ;  ...
  222.     mov    cx,word ptr linitx+2    ;  ...
  223.     mov    word ptr x,bx        ;  ...
  224.     mov    word ptr x+2,cx     ;  ...
  225.     mov    word ptr linitx,ax    ;  ...
  226.     mov    word ptr linitx+2,dx    ;  ...
  227.  
  228.     mov    ax,word ptr y        ; switch y with linity
  229.     mov    dx,word ptr y+2     ;  ...
  230.     mov    bx,word ptr linity    ;  ...
  231.     mov    cx,word ptr linity+2    ;  ...
  232.     mov    word ptr y,bx        ;  ...
  233.     mov    word ptr y+2,cx     ;  ...
  234.     mov    word ptr linity,ax    ;  ...
  235.     mov    word ptr linity+2,dx    ;  ...
  236.  
  237. doeither:                ; common Mandelbrot, Julia set code
  238.     mov    period,0        ; claim periodicity of 1
  239.     mov    savedand,1        ; initial periodicity check
  240.     mov    savedincr,1        ;  flag for incrementing periodicity
  241.     mov    word ptr savedx+2,0ffffh; impossible value of "old" x
  242.     mov    word ptr savedy+2,0ffffh; impossible value of "old" y
  243.     mov    orbit_ptr,0        ; clear orbits
  244.  
  245.     dec    kbdcount        ; decrement the keyboard counter
  246.     jns    short nokey        ;  skip keyboard test if still positive
  247.     mov    kbdcount,10        ; stuff in a low kbd count
  248.     cmp    show_orbit,0        ; are we showing orbits?
  249.     jne    quickkbd        ;  yup.  leave it that way.
  250.     mov    kbdcount,5000        ; else, stuff an appropriate count val
  251.     cmp    cpu,386         ; ("appropriate" to the CPU)
  252.     je    short kbddiskadj    ;  ...
  253. ;;    cmp    word ptr delmin+2,1    ; is 16-bit math good enough?
  254.     cmp    word ptr delmin+2,8    ; is 16-bit math good enough?
  255.     ja    kbddiskadj        ;  yes. test less often
  256.     mov    kbdcount,500        ;  no.    test more often
  257. kbddiskadj:
  258.     cmp    dotmode,11        ; disk video?
  259.     jne    quickkbd        ;  no, leave as is
  260.     shr    kbdcount,1        ; yes, reduce count
  261.     shr    kbdcount,1        ;  ...
  262. quickkbd:
  263.     call    far ptr keypressed    ; has a key been pressed?
  264.     cmp    ax,0            ;  ...
  265.     je    nokey            ; nope.  proceed
  266.     mov    kbdcount,0        ; make sure it goes negative again
  267.     cmp    ax,'o'                  ; orbit toggle hit?
  268.     je    orbitkey        ;  yup.  show orbits
  269.     cmp    ax,'O'                  ; orbit toggle hit?
  270.     jne    keyhit            ;  nope.  normal key.
  271. orbitkey:
  272.     call    far ptr getakey     ; read the key for real
  273.     mov    ax,1            ; reset orbittoggle = 1 - orbittoggle
  274.     sub    ax,show_orbit        ;  ...
  275.     mov    show_orbit,ax        ;  ...
  276.     jmp    short nokey        ; pretend no key was hit
  277. keyhit: mov    ax,-1            ; return with -1
  278.     mov    color,ax        ; set color to -1
  279.     UNFRAME <si,di>         ; pop stack frame
  280.     ret                ; bail out!
  281.  
  282. nokey:
  283.     cmp    show_orbit,0        ; is orbiting on?
  284.     jne    no16bitcode        ;  yup.  slow down.
  285.     cmp    cpu,386         ; are we on a 386?
  286.     je    short code386bit    ;  YAY!! 386-class speed!
  287. ;;    cmp    word ptr delmin+2,1    ; OK, we're desperate.  16 bits OK?
  288.     cmp    word ptr delmin+2,8    ; OK, we're desperate.  16 bits OK?
  289.     ja    yes16bitcode        ;  YAY!  16-bit speed!
  290. no16bitcode:
  291.     call    near ptr code32bit    ; BOO!! nap time.  Full 32 bit math
  292.     jmp    kloopend        ;  bypass the 386-specific code.
  293. yes16bitcode:
  294.     call    near ptr code16bit    ; invoke the 16-bit version
  295.     jmp    kloopend        ;  bypass the 386-specific code.
  296.  
  297. .386                    ; 386-specific code starts here
  298.  
  299. code386bit:
  300. ;;    cmp    word ptr delmin+2,3    ; is 16-bit math good enough?
  301.     cmp    word ptr delmin+2,8    ; is 16-bit math good enough?
  302.     jbe    code386_32        ; nope, go do 32 bit stuff
  303. IFDEF ??version
  304.     jmp    code386_32        ; TASM screws up IMUL EBX,EBX!!
  305. ENDIF
  306.  
  307.     ; 16 bit on 386, now we are really gonna move
  308.     movsx    esi,word ptr x+2    ; use SI for X
  309.     movsx    edi,word ptr y+2    ; use DI for Y
  310.     push    ebp
  311.     mov    ebp,-1
  312.     shl    ebp,FUDGEFACTOR-1
  313.     mov    cx,FUDGEFACTOR-16
  314.  
  315. kloop386_16:   ; cx=bitshift-16, ebp=overflow.mask
  316.  
  317.     mov    ebx,esi         ; compute (x * x)
  318.     imul    ebx,ebx         ;  ...
  319.     test    ebx,ebp         ;
  320.     jnz    short end386_16     ;  (oops.  We done.)
  321.     shr    ebx,cl            ; get result down to 16 bits
  322.  
  323.     mov    edx,edi         ; compute (y * y)
  324.     imul    edx,edx         ;  ...
  325.     test    edx,ebp         ; say, did we overflow? <V20-compat>
  326.     jnz    short end386_16     ;  (oops.  We done.)
  327.     shr    edx,cl            ; get result down to 16 bits
  328.  
  329.     mov    ax,bx            ; compute (x*x - y*y) / fudge
  330.     sub    bx,dx            ;  for the next iteration
  331.  
  332.     add    ax,dx            ; compute (x*x + y*y) / fudge
  333.  
  334.     cmp    ax,word ptr lm+2    ; while (xx+yy < lm)
  335.     jae    short end386_16     ;  ...
  336.  
  337.     imul    edi,esi         ; compute (y * x)
  338.     shl    edi,1            ; ( * 2 / fudge)
  339.     sar    edi,cl
  340.     add    di,word ptr linity+2    ; (2*y*x) / fudge + linity
  341.     movsx    edi,di            ; save as y
  342.  
  343.     add    bx,word ptr linitx+2    ; (from above) (x*x - y*y)/fudge + linitx
  344.     movsx    esi,bx            ; save as x
  345.  
  346.     mov    ax,oldcolor        ; recall the old color
  347.     cmp    ax,k            ; check it against this iter
  348.     jge    short chkpd386_16    ;  yup.  do periodicity check.
  349. nonmax386_16:
  350.     dec    k            ; while (k < maxit)
  351.     jnz    short kloop386_16    ; try, try again
  352. end386_16:
  353.     pop    ebp
  354.     jmp    kloopend        ; we done
  355.  
  356. chkpd386_16:
  357.     mov    ax,k            ; set up to test for save-time
  358.     test    ax,savedand        ; save on 0, check on anything else
  359.     jz    short chksv386_16    ;  time to save a new "old" value
  360.     mov    ax,si            ; load up x
  361.     xor    ax,word ptr savedx+2    ; does X match?
  362.     test    ax,word ptr savedmask+2 ;  truncate to appropriate precision
  363.     jne    short nonmax386_16    ;  nope.  forget it.
  364.     mov    ax,di            ; now test y
  365.     xor    ax,word ptr savedy+2    ; does Y match?
  366.     test    ax,word ptr savedmask+2 ;  truncate to appropriate precision
  367.     jne    short nonmax386_16    ;  nope.  forget it.
  368.     mov    period,1        ; note that we have found periodicity
  369.     mov    k,0            ; pretend maxit reached
  370.     jmp    short end386_16
  371. chksv386_16:
  372.     mov    word ptr savedx+2,si    ; save x
  373.     mov    word ptr savedy+2,di    ; save y
  374.     dec    savedincr        ; time to change the periodicity?
  375.     jnz    short nonmax386_16    ;  nope.
  376.     shl    savedand,1        ; well then, let's try this one!
  377.     inc    savedand        ;  (2**n -1)
  378.     mov    savedincr,4        ; and reset the increment flag
  379.     jmp    short nonmax386_16
  380.  
  381.     ; 32bit on 386:
  382. code386_32:
  383.     mov    esi,x            ; use ESI for X
  384.     mov    edi,y            ; use EDI for Y
  385.  
  386. ;    This is the main processing loop.  Here, every T-state counts...
  387.  
  388. kloop:                    ; for (k = 0; k <= maxit; k++)
  389.  
  390.     mov    eax,esi         ; compute (x * x)
  391.     imul    esi            ;  ...
  392.     shrd    eax,edx,FUDGEFACTOR    ; ( / fudge)
  393.     shr    edx,FUDGEFACTOR-1    ; (complete 64-bit shift and check
  394.     jne    short kloopend1     ; bail out if too high
  395.     mov    ebx,eax         ; save this for below
  396.  
  397.     mov    eax,edi         ; compute (y * y)
  398.     imul    edi            ;  ...
  399.     shrd    eax,edx,FUDGEFACTOR    ; ( / fudge)
  400.     shr    edx,FUDGEFACTOR-1    ; (complete 64-bit shift and check
  401.     jne    short kloopend1     ; bail out if too high
  402.  
  403.     mov    ecx,ebx         ; compute (x*x - y*y) / fudge
  404.     sub    ebx,eax         ;  for the next iteration
  405.  
  406.     add    ecx,eax         ; compute (x*x + y*y) / fudge
  407.     cmp    ecx,lm            ; while (lr < lm)
  408.     jae    short kloopend1     ;  ...
  409.  
  410.     mov    eax,edi         ; compute (y * x)
  411.     imul    esi            ;  ...
  412.     shrd    eax,edx,FUDGEFACTOR-1    ;  ( * 2 / fudge)
  413.     add    eax,linity        ;  (above) + linity
  414.     mov    edi,eax         ;  save this as y
  415.  
  416. ;    (from the earlier code)     ; compute (x*x - y*y) / fudge
  417.     add    ebx,linitx        ;    + linitx
  418.     mov    esi,ebx         ; save this as x
  419.  
  420.     mov    ax,oldcolor        ; recall the old color
  421.     cmp    ax,k            ; check it against this iter
  422.     jge    short chkperiod1
  423. nonmax1:
  424.     dec    k            ; while (k < maxit) (dec to 0 is faster)
  425.     jnz    short kloop        ; while (k < maxit) ...
  426. kloopend1:
  427.     jmp    short kloopend32    ; we done.
  428.  
  429. chkperiod1:
  430.     mov    eax,esi
  431.     xor    eax,savedx
  432.     test    eax,savedmask
  433.     jnz    short chksave1
  434.     mov    eax,edi
  435.     xor    eax,savedy
  436.     test    eax,savedmask
  437.     jnz    short chksave1
  438.     mov    period,1        ; note that we have found periodicity
  439.     mov    k,0            ; pretend maxit reached
  440.     jmp    short kloopend1
  441. chksave1:
  442.     mov    ax,k
  443.     test    ax,savedand
  444.     jne    short nonmax1
  445.     mov    savedx,esi
  446.     mov    savedy,edi
  447.     dec    savedincr        ; time to change the periodicity?
  448.     jnz    short nonmax1        ;  nope.
  449.     shl    savedand,1        ; well then, let's try this one!
  450.     inc    savedand        ;  (2**n -1)
  451.     mov    savedincr,4        ; and reset the increment flag
  452.     jmp    short nonmax1
  453.  
  454. kloopend32:
  455.  
  456. .8086                    ; 386-specific code ends here
  457.  
  458. kloopend:
  459.     cmp    orbit_ptr,0        ; any orbits to clear?
  460.     je    noorbit2        ;  nope.
  461.     call    far ptr scrub_orbit    ; clear out any old orbits
  462. noorbit2:
  463.  
  464.     mov    ax,k            ; set old color
  465.     sub    ax,10            ; minus 10, for safety
  466.     mov    oldcolor,ax        ; and save it as the "old" color
  467.     mov    ax,maxit        ; compute color
  468.     sub    ax,k            ;  (first, re-compute "k")
  469.     sub    kbdcount,ax        ; adjust the keyboard count
  470.     cmp    ax,1            ; convert any "outlier" region
  471.     jge    short coloradjust1    ;  (where abs(x) > 2 or abs(y) > 2)
  472.     mov    ax,1            ;   to look like we ran through
  473. coloradjust1:                ;    at least one loop.
  474.     mov    realcolor,ax        ; result before adjustments
  475.     cmp    ax,maxit        ; did we max out on iterations?
  476.     jne    short wedone        ;  nope.
  477.     mov    oldcolor,ax        ; set "oldcolor" to maximum
  478.     cmp    inside,0        ; is "inside" >= 0?
  479.     jl    wedone            ;  nope.  leave it at "maxit"
  480.     mov    ax,inside        ; reset max-out color to default
  481.     cmp    periodicitycheck,0    ; show periodicity matches?
  482.     jge    wedone            ;  nope.
  483.     mov    al,period        ;  reset color to periodicity flag
  484. wedone:                 ;
  485.     mov    color,ax        ; save the color result
  486.     UNFRAME <si,di>         ; pop stack frame
  487.     ret                ; and return with color
  488.  
  489. calcmandasm endp
  490.  
  491.  
  492. ; ******************** Function code16bit() *****************************
  493. ;
  494. ;    Performs "short-cut" 16-bit math where we can get away with it.
  495. ; CJLT has modified it, mostly by preshifting x and y to fg30 from fg29
  496. ; or, since we ignore the lower 16 bits, fg14 from fg13.
  497. ; If this shift overflows we are outside x*x+y*y=2, so have escaped.
  498. ; Also, he commented out several conditional jumps which he calculated could
  499. ; never be taken (e.g. mov ax,si / imul si ;cannot overflow).
  500.  
  501. code16bit    proc    near
  502.  
  503.     mov    si,word ptr x+2     ; use SI for X fg13
  504.     mov    di,word ptr y+2     ; use DI for Y fg13
  505.  
  506. start16bit:
  507.     add    si,si            ;CJLT-Convert to fg14
  508.     jo    end16bit        ;overflows if <-2 or >2
  509.     mov    ax,si            ; compute (x * x)
  510.     imul    si            ; Answer is fg14+14-16=fg12
  511. ;    cmp    dx,0            ;CJLT commented out-
  512. ;    jl    end16bit        ;CJLT-  imul CANNOT overflow
  513. ;    mov    cx,32-FUDGEFACTOR    ;CJLT. FUDGEFACTOR=29 is hard coded
  514. loop16bit1:
  515.     shl    ax,1            ;  ...
  516.     rcl    dx,1            ;  ...
  517.     jo    end16bit        ;  (oops.  overflow)
  518. ;    loop    loop16bit1        ;CJLT...do it once only. dx now fg13.
  519.     mov    bx,dx            ; save this for a tad
  520.  
  521. ;ditto for y*y...
  522.  
  523.     add    di,di            ;CJLT-Convert to fg14
  524.     jo    end16bit        ;overflows if <-2 or >2
  525.     mov    ax,di            ; compute (y * y)
  526.     imul    di            ;  ...
  527. ;    cmp    dx,0            ; say, did we overflow? <V20-compat>
  528. ;    jl    end16bit        ;  (oops.  We done.)
  529. ;    mov    cx,32-FUDGEFACTOR    ; ( / fudge)
  530. ;loop16bit2:
  531.     shl    ax,1            ;  ...
  532.     rcl    dx,1            ;  ...
  533.     jo    end16bit        ;  (oops.  overflow)
  534. ;    loop    loop16bit2        ;  ...
  535.  
  536.     mov    cx,bx            ; compute (x*x - y*y) / fudge
  537.     sub    bx,dx            ;  for the next iteration
  538.  
  539.     add    cx,dx            ; compute (x*x + y*y) / fudge
  540.     jo    end16bit        ; bail out if too high
  541. ;    js    end16bit        ;  ...
  542.  
  543.     cmp    cx,word ptr lm+2    ; while (xx+yy < lm)
  544.     jae    end16bit        ;  ...
  545.     dec    k            ; while (k < maxit)
  546.     jz    end16bit        ;  we done.
  547.  
  548.     mov    ax,di            ; compute (y * x) fg14+14=fg28
  549.     imul    si            ;  ...
  550. ;    mov    cx,33-FUDGEFACTOR-2    ; ( * 2 / fudge)
  551. ;loop16bit3:
  552.     shl    ax,1            ;  ...
  553.     rcl    dx,1            ;  ...
  554.     shl    ax,1            ;  shift two bits
  555.     rcl    dx,1            ;  cannot overflow as |x|<=2, |y|<=2
  556. ;    loop    loop16bit3        ;  ...
  557.     add    dx,word ptr linity+2    ; (2*y*x) / fudge + linity
  558.     jo    end16bit        ; bail out if too high
  559.     mov    di,dx            ; save as y
  560.  
  561.     add    bx,word ptr linitx+2    ; (from above) (x*x - y*y)/fudge + linitx
  562.     jo    end16bit        ; bail out if too high
  563.     mov    si,bx            ; save as x
  564.  
  565.     mov    ax,oldcolor        ; recall the old color
  566.     cmp    ax,k            ; check it against this iter
  567.     jle    short nonmax3        ;  nope.  bypass periodicity check.
  568.     mov    word ptr x+2,si     ; save x for periodicity check
  569.     mov    word ptr y+2,di     ; save y for periodicity check
  570.     call    checkperiod        ; check for periodicity
  571. nonmax3:
  572.     jmp    start16bit        ; try, try again.
  573.  
  574. end16bit:                ; we done.
  575.     ret
  576. code16bit    endp
  577.  
  578.  
  579. ;    The following routine checks for periodic loops (a known
  580. ;    method of decay inside "Mandelbrot Lake", and an easy way to
  581. ;    bail out of "lake" points quickly).  For speed, only the
  582. ;    high-order sixteen bits of X and Y are checked for periodicity.
  583. ;    For accuracy, this routine is only fired up if the previous pixel
  584. ;    was in the lake (which means that the FIRST "wet" pixel was
  585. ;    detected by the dull-normal maximum iteration process).
  586.  
  587. checkperiod    proc near        ; periodicity check
  588.     mov    ax,k            ; set up to test for save-time
  589.     test    ax,savedand        ; save on 0, check on anything else
  590.     jz    checksave        ;  time to save a new "old" value
  591.     mov    dx,word ptr x+2     ; load up x
  592.     and    dx,word ptr savedmask+2 ;  truncate to appropriate precision
  593.     cmp    dx,word ptr savedx+2    ; does X match?
  594.     jne    checkdone        ;  nope.  forget it.
  595.     mov    ax,word ptr x        ; load up x
  596.     and    ax,word ptr savedmask    ;  truncate to appropriate precision
  597.     cmp    ax,word ptr savedx    ; does X match?
  598.     jne    checkdone        ;  nope.  forget it.
  599.     mov    dx,word ptr y+2     ; now test y
  600.     and    dx,word ptr savedmask+2 ;  truncate to appropriate precision
  601.     cmp    dx,word ptr savedy+2    ; does Y match?
  602.     jne    checkdone        ;  nope.  forget it.
  603.     mov    ax,word ptr y        ; load up y
  604.     and    ax,word ptr savedmask    ;  truncate to appropriate precision
  605.     cmp    ax,word ptr savedy    ; does Y match?
  606.     jne    checkdone        ;  nope.  forget it.
  607.     mov    period,1        ; note that we have found periodicity
  608.     mov    k,1            ; pretend maxit reached
  609. checksave:
  610.     mov    dx,word ptr x+2     ; load up x
  611.     and    dx,word ptr savedmask+2 ;  truncate to appropriate precision
  612.     mov    word ptr savedx+2,dx    ;  and save it
  613.     mov    ax,word ptr x        ; load up x
  614.     and    ax,word ptr savedmask    ;  truncate to appropriate precision
  615.     mov    word ptr savedx,ax    ;  and save it
  616.     mov    dx,word ptr y+2     ; load up y
  617.     and    dx,word ptr savedmask+2 ;  truncate to appropriate precision
  618.     mov    word ptr savedy+2,dx    ;  and save it
  619.     mov    ax,word ptr y        ; load up y
  620.     and    ax,word ptr savedmask    ;  truncate to appropriate precision
  621.     mov    word ptr savedy,ax    ;  and save it
  622.     dec    savedincr        ; time to change the periodicity?
  623.     jnz    checkdone        ;  nope.
  624.     shl    savedand,1        ; well then, let's try this one!
  625.     inc    savedand        ;  (2**n -1)
  626.     mov    savedincr,4        ; and reset the increment flag
  627. checkdone:
  628.     ret                ; we done.
  629. checkperiod    endp
  630.  
  631.  
  632. ; ******************** Function code32bit() *****************************
  633. ;
  634. ;    Perform the 32-bit logic required using 16-bit logic
  635. ;
  636. ;    New twice as fast logic,
  637. ;       Courtesy of Bill Townsend and Mike Gelvin (CIS:73337,520)
  638. ;    Even newer, faster still by Chris Lusby Taylor
  639. ;     who noted that we needn't square the low word if we first multiply
  640. ;     by 4, since we only need 29 places after the point and this will
  641. ;     give 30. (We divide answer by two to give 29 bit shift in answer)
  642. ;    Also, he removed all testing for special cases where a word of operand
  643. ;    happens to be 0, since testing 65536 times costs more than the saving
  644. ;    1 time out of 65536! (I benchmarked it. Just removing the tests speeds
  645. ;    us up by 3%.)
  646. ;
  647. ;Note that square returns DI,AX squared in DX,AX now.
  648. ; DI,AX is first converted to unsigned fg31 form.
  649. ; (For its square to be representable in fg29 (range -4..+3.999)
  650. ; DI:AX must be in the range 0..+1.999 which fits neatly into unsigned fg31.)
  651. ; This allows us to ignore the part of the answer corresponding to AX*AX as it
  652. ; is less than half a least significant bit of the final answer.
  653. ; I thought you'd like that.
  654. ;
  655. ; As we prescaled DI:AX, we need to shift the answer 1 bit to the right to
  656. ; end up in fg29 form since 29=(29+2)+(29+2)-32-1
  657. ; However, the mid term AX*DI is needed twice, so the shifts cancel.
  658. ;
  659. ; Since abs(x) and abs(y) in fg31 form will be needed in calculating 2*X*Y
  660. ; we push them onto the stack for later use.
  661.  
  662. ; Note that square does nor affect bl,si,bp
  663. ; and leaves highword of argument in di
  664. ; but destroys bh,cx
  665. square    MACRO    donepops
  666.     LOCAL    notneg
  667.     shl    ax,1        ;Multiply by 2 to convert to fg30
  668.     rcl    di,1        ;If this overflows DI:AX was negative
  669.     jnc    notneg
  670.     not    ax            ; so negate it
  671.     not    di            ; ...
  672.     add    ax,1            ; ...
  673.     adc    di,0            ; ...
  674.     not    bl            ; change negswt
  675. notneg:    shl    ax,1        ;Multiply by 2 again to give fg31
  676.     rcl    di,1        ;If this gives a carry then DI:AX was >=2.0
  677.                 ;If its high bit is set then DI:AX was >=1.0
  678.                 ;This is OK, but note that this means that
  679.                 ;DI:AX must now be treated as unsigned.
  680.     jc    donepops
  681.     push    di        ; save y or x (in fg31 form) on stack
  682.     push    ax        ; ...
  683.     mul    di        ;GET MIDDLE PART - 2*A*B
  684.     mov    bh,ah        ;Miraculously, it needs no shifting!
  685.     mov    cx,dx
  686.     mov    ax,di
  687.     mul    ax        ;SQUARE HIGH HWORD - A*A
  688.     shl    bh,1        ;See if we round up
  689.     adc    ax,1        ;Anyway, add 1 to round up/down accurately
  690.     adc    dx,0
  691.     shr    dx,1        ;This needs shifting one bit
  692.     rcr    ax,1
  693.     add    ax,cx        ;Add in the 2*A*B term
  694.     adc    dx,0
  695.     ENDM    ;#EM
  696.  
  697.  
  698. code32bit    proc near
  699. ;
  700. ; BL IS USED FOR THE "NEGSWT" FLAG
  701. ;   NEGSWT IS ONE IF EITHER "X" OR "Y" ARE NEGATIVE, BUT NOT BOTH NEGATIVE
  702. ;
  703.  
  704.     push    bp
  705.     xor    bl,bl            ; NEGSWT STARTS OFF ZERO
  706.  
  707. ;    iteration loop
  708.  
  709. nextit:    mov    ax,word ptr y        ; ax=low(y)
  710.     mov    di,word ptr y+2     ; di=high(y)
  711.     square done0    ;square y and quit via done0 if it overflows
  712.     mov    si,ax        ; square returns results in dx,ax
  713.     mov    bp,dx        ; save y*y in bp,si
  714.     mov    ax,word ptr x
  715.     mov    di,word ptr x+2
  716.     square    done2        ; square x and quit via done2 if it overflows
  717.     mov    cx,ax        ; Save low answer in cx.
  718.     ADD    ax,si        ; calc y*y + x*x
  719.     mov    ax,bp
  720.     ADC    ax,dx        ;  ...
  721.     jno    nextxy        ; overflow?
  722.                 ;NOTE: The original code tests against lm
  723.                 ;here, but as lm=4<<29 this is the same
  724.                 ;as testing for signed overflow
  725. done4:    add    sp,4            ; discard saved value of |x| fg 31
  726. done2:    add    sp,4            ; discard saved value of |y| fg 31
  727. done0: pop    bp            ; restore saved bp
  728.     ret
  729.  
  730. ;---------------------------------------------------------------------------
  731.  
  732. nextxy:    dec    k        ; while (k < maxit)
  733.     jz    done4        ;  we done.
  734.     sub    cx,si            ; subtract y*y from x*x
  735.     sbb    dx,bp            ;  ...
  736.     add    cx,word ptr linitx    ; add "A"
  737.     adc    dx,word ptr linitx+2    ;  ...
  738.     jo    done4            ;CJLT-Must detect overflow here
  739.                     ; but increment loop count first
  740.     mov    word ptr x,cx        ; store new x = x*x-y*y+a
  741.     mov    word ptr x+2,dx     ;  ...
  742.  
  743. ; now calculate x*y
  744. ;
  745. ;More CJLT tricks here. We use the pushed abs(x) and abs(y) in fg31 form
  746. ;which, when multiplied, give x*y in fg30, which, luckily, is the same as...
  747. ;2*x*y fg29.
  748. ;As with squaring, we can ignore the product of the low order words, and still
  749. ;be more accurate than the original algorithm.
  750. ;
  751.     pop    bp        ;Xlow
  752.     pop    di        ;Xhigh (already there, actually)
  753.     pop    ax        ;Ylow
  754.     mul    di        ;Xhigh * Ylow
  755.     mov    bh,ah        ;Discard lowest 8 bits of answer
  756.     mov    cx,dx
  757.     pop    ax        ;Yhigh
  758.     mov    si,ax        ; copy it
  759.     mul    bp        ;Xlow * Yhigh
  760.     xor    bp,bp        ;Clear answer
  761.     add    bh,ah
  762.     adc    cx,dx
  763.     adc    bp,0
  764.     mov    ax,si        ;Yhigh
  765.     mul    di        ;Xhigh * Yhigh
  766.     shl    bh,1        ;round up/down
  767.     adc    ax,cx        ;Answer-low
  768.     adc    dx,bp        ;Answer-high
  769.                 ;NOTE: The answer is 0..3.9999 in fg29
  770.     js    done0        ;Overflow if high bit set
  771.     or    bl,bl        ; ZERO IF NONE OR BOTH X , Y NEG
  772.     jz    signok        ; ONE IF ONLY ONE OF X OR Y IS NEG
  773.     not    ax        ; negate result
  774.     not    dx        ;  ...
  775.     add    ax,1        ;  ...
  776.     adc    dx,0        ;  ...
  777.     xor    bl,bl        ;Clear negswt
  778. signok:
  779.     add    ax,word ptr linity
  780.     adc    dx,word ptr linity+2    ; dx,ax = 2(X*Y)+B
  781.     jo    done0
  782.     mov    word ptr y,ax        ; save the new value of y
  783.     mov    word ptr y+2,dx     ;  ...
  784.     mov    ax,oldcolor        ; recall the old color
  785.     cmp    ax,k            ; check it against this iter
  786.     jle    short chkmaxit        ;  nope.  bypass periodicity check.
  787.     call    checkperiod        ; check for periodicity
  788.  
  789. chkmaxit:
  790.     cmp    show_orbit,0        ; orbiting on?
  791.     jne    horbit            ;  yep.
  792.     jmp    nextit            ;go around again
  793.  
  794. jmpdone0:
  795.     jmp    done0            ; DOES [(X*X)-(Y*Y)+P] BEFORE THE DEC.
  796.  
  797. horbit: push    bx            ; save my flags
  798.     mov    ax,-1            ; color for plot orbit
  799.     push    ax            ;  ...
  800.     push    word ptr y+2        ; co-ordinates for plot orbit
  801.     push    word ptr y        ;  ...
  802.     push    word ptr x+2        ;  ...
  803.     push    word ptr x        ;  ...
  804.     call    far ptr iplot_orbit    ; display the orbit
  805.     add    sp,5*2            ; clear out the parameters
  806.     pop    bx            ; restore flags
  807.     jmp    nextit            ; go around again
  808.  
  809. code32bit    endp
  810.  
  811.        end
  812.  
  813.