home *** CD-ROM | disk | FTP | other *** search
/ The Fred Fish Collection 1.5 / ffcollection-1-5-1992-11.iso / ff_disks / 200-299 / ff239.lzh / JGoodies / Mandelbrot / Mandelshift2.ASM < prev    next >
Assembly Source File  |  1989-08-21  |  8KB  |  245 lines
  1. \ Copyright 1989 NerveWare
  2. \ No portion of this code may used for commercial purposes,
  3. \ nor may any executable version of this code be disributed for 
  4. \ commercial purposes without the author's express written permission.
  5. \ This code is shareware, all rights reserved.
  6. \ Nick Didkovsky
  7. \ 12/21/88
  8. \ 
  9. \                            MandelShift2.ASM 
  10.  
  11. \ MOD: Uses 16384 as max possible scaler!             12/20/88
  12. \ MOD: 350 as mandelmax seems optimal at this scale   2/11/89
  13. \ MOD: storing mandelmax iterations in register A1    2/11/89
  14. \ MOD: storing registers d2,d5 in scratch registers a0,a2  2/11/89
  15.  
  16. \ 4 := 65536
  17. \ 2 := 32768
  18.  
  19. \ This is the speedy assembler Mandelbrot routine.
  20. \ A little roundoff error is introduced by integer math,
  21. \ visible as a little roughness along the high-contrast boundaries of the 
  22. \ image generated. But it's insanely fast.
  23.  
  24. \ MandelMax iterations of the loop qualifies a given complex 
  25. \ number as belonging to the Mandelbrot Set.
  26.  
  27. \                          THE ALGORITHM
  28. \ The iterative function applied to a given complex number C is: 
  29. \ Z := Z^2+C
  30. \ In terms of the variables below, C = aconst + (bconst)i , where i
  31. \ is the square root of -1.
  32. \ The function is initialized with Z=0+0i 
  33. \ The function's result is fed into Z again, recalculated, fed in, 
  34. \ recalculated, over and over.
  35. \ Each intermediate result's magnitude is tested.  If the magnitude >= 2,
  36. \ break out of the loop, leaving the number of iterations it took to 
  37. \ "blow up" to this value on the stack.
  38. \ If, after 254 iterations of this function, the magnitude is still < 2,
  39. \ tested number C is considered a member of the Mandelbrot Set, and 254 is left
  40. \ on the stack.
  41.  
  42. \ SPEED SHORTCUTS:
  43.  
  44. \ 1) The magnitude of a complex number equals the square root of the sum
  45. \    of the squares of its real and imaginary parts (simple "distance formula"
  46. \    from high school).  
  47. \    Shortcut: Don't bother taking the square root and checking it
  48. \    against 2, rather test the sum of the squares against 4.
  49.  
  50. \ 2) As mentioned before, I use scaled integer math.  This is MUCH faster
  51. \    than floating point.  Additionally, instead of scaling by a nice looking 
  52. \    number like 10000, I use 16384, which is a power of 2.  So, instead of
  53. \    scaling results with a DIVISION by 10000 (slow), I scale it by
  54. \    SHIFTING register contents to the right 14 places (very very fast).
  55.  
  56. \ 3) The Z := Z^2+C function can be broken down into two calculations:
  57. \    one calculates what happens to the real part of Z (call this value az), 
  58. \    the other calculates what happens to the imaginary part (call this bz).
  59. \        NEWaz := az*az - bz*bz + aconst
  60. \        NEWbz := 2*az*bz + bconst
  61. \    Speedup: Notice that az*az and bz*bz are needed to calculate the new az,
  62. \    as well as to check the magnitude of the result.  So I only calculate these
  63. \    squares ONCE in each iteration, storing the results in registers d5 and d6,
  64. \    and use them for both purposes.
  65.  
  66. \ 4) The routine uses 68000 registers ONLY.  No subroutine calls at all.  This
  67. \    makes for very fast arithmetic!  Register's previous contents are stored 
  68. \    at the top of the routine, then restored at the bottom before it exits.
  69.  
  70. anew task_Mandelbrot-ASM
  71.  
  72. \ maximum iterations is hard coded in this version of the program = 350
  73.  
  74. 350 CONSTANT MandelMax
  75.  
  76. \ These variables determine the region generated.
  77. \ Their values represent scaled values, where the scaler = 16384.
  78. \ That means unity is represented by 16384.
  79. \ The smallest value in our scaled arithmetic is obviously "1",
  80. \ which represents the fraction (1 / 16384) = 0.000061
  81. \ So, smallest distance representable between two adjacent pixels is 0.000061
  82. \ ok so it's not double precision ... 
  83.  
  84. VARIABLE acorner  \ region's top left real part ( x value)
  85. VARIABLE bcorner  \ region's top left imaginary part ( y value)
  86. VARIABLE aconst   \ real part of number being tested
  87. VARIABLE bconst   \ imaginary part of number being tested
  88. VARIABLE xgap     \ horizontal distance between adjacent pixels
  89. VARIABLE ygap     \ vertical distance between adjacent pixels
  90.  
  91. \ the following variables are used to store register contents during
  92. \ the execution of the loop.  Registers are restored when this code returns.
  93.  
  94. VARIABLE temp.d3   
  95. VARIABLE temp.d4   
  96. VARIABLE temp.d6   
  97.  
  98. \ Speed notes:
  99. \ 1000 worst-case calls to this routine takes 2.72 sec. with MandelMax = 35
  100. \ (using 0+0i, which is a Mandelbrot member)
  101. \ Old version, using 10000 scaler and DIVS took 4.58 sec!!!
  102.  
  103.  
  104. \ This uses registers only for super speed. 
  105.  
  106. \ d0 = az
  107. \ d1 = bz
  108. \ d2 = loop counter
  109. \ d3 = aconst
  110. \ d4 = bconst
  111. \ d5 = azsqr
  112. \ d6 = bzsqr
  113. \ a1 = mandelmax
  114.  
  115. ASM quick.guts  ( -- iterations)
  116. \ save d3,d4,d6 in variables
  117.     callcfa    temp.d3
  118.     move.l    d3,$0(org,tos.l)
  119.     move.l    (dsp)+,tos
  120.     callcfa    temp.d4
  121.     move.l    d4,$0(org,tos.l)
  122.     move.l    (dsp)+,tos
  123.     callcfa    temp.d6
  124.     move.l    d6,$0(org,tos.l)
  125.     move.l    (dsp)+,tos
  126.  
  127.  
  128. \ ********************** initialize d3, d4 with aconst and bconst *********
  129.  
  130.     callcfa    aconst
  131.     move.l    $0(org,tos.l),d3
  132.     move.l    (dsp)+,tos
  133.     callcfa    bconst
  134.     move.l    $0(org,tos.l),d4
  135.     move.l    (dsp)+,tos
  136.     callcfa    MandelMax
  137.     move.l    tos,a1
  138.     move.l      (dsp)+,tos
  139. \ save d2,d5 in scratch registers a0,a2
  140.     move.l    d2,a0
  141.     move.l    d5,a2
  142.  
  143. \ ************************ more init's ************************************
  144.     
  145.     moveq    #0,d2            \ initialize loop counter
  146.     moveq    #0,d0            \ initialize az
  147.     moveq    #0,d1            \ initialize bz
  148.  
  149. \ ******************  loop starts, calculate azsqr, bzsqr *****************
  150.  
  151. \ Check if operand greater than 32768.  If so, don't bother squaring: exit loop
  152. \ because magnitude will automatically be greater than 4 = 65536.
  153.  
  154. \ AZSQR
  155. 1$: move.l    d0,d5               \ az copied to d5
  156.     tst.l    d5            \ make sure positive
  157.     bpl.s    22$
  158.     neg.l    d5
  159. 22$: cmpi.l    #32768,d5               \ operand >= 32768 ?
  160.     bge        2$            \ if so, exit loop
  161.     muls.l    d5,d5            \ else square it, d5 now contains azsqr
  162.     asr.l    #$8,d5            \ FAST SCALE BY 16384 !!!
  163.     asr.l    #$6,d5
  164.  
  165. \ BZSQR
  166.     move.l    d1,d6            \ bz copied to d6
  167.     tst.l    d6            \ make sure positive
  168.     bpl.s    32$
  169.     neg.l    d6
  170. 32$: cmpi.l    #32768,d6               \ operand >= 32768
  171.     bge        2$            \ if so, exit loop
  172.     muls.l    d6,d6            \ else square it, d6 now contains bzsqr
  173.     asr.l    #8,d6            \ FAST SCALE!!!
  174.     asr.l    #6,d6
  175.  
  176.  
  177. \ ************************** new BZ **************************
  178.  
  179.  
  180.    tst.l d0                \ test sign of d0
  181.    bpl.s 11$                 \ d0 positive, goto 11$
  182.    neg.l d0                \ else negate d0
  183.    tst.l d1                \ and test d1
  184.    bpl.l 31$                \ d1 pos, d0 was neg, goto 31$
  185.    neg.l d1                \ both neg, negate d1 
  186.    bra   21$                \ and goto 21$ 
  187.  
  188. 11$: tst.l d1                \ d0 pos, test d1
  189.      bpl.s 21$                \ both positive, goto 21$
  190.      neg.l d1                \ else negate d1
  191.      bra   31$                \ and goto 31$
  192.  
  193. 21$: muls.l    d1,d0            \ * bz
  194.      asr.l      #8,d0            \ FAST SCALE !!!
  195.      asr.l    #6,d0
  196.      bra     41$            \ leave answer positive
  197.  
  198. 31$: muls.l    d1,d0            \ * bz
  199.      asr.l      #8,d0            \ FAST SCALE !!!
  200.      asr.l    #6,d0
  201.      neg.l    d0            \ make final answer negative
  202.  
  203.  
  204. 41$: add.l    d0,d0            \ 2* 
  205.      add.l    d4,d0            \ + bconst
  206.      move.l    d0,d1            \ new bz now in d1
  207.  
  208. \ ************************** NEW AZ ***********************************
  209.  
  210.     move.l    d6,d0            \ copy bzsqr out to d0
  211.     neg.l    d0            \ negate its copy
  212.     add.l    d5,d0            \ azsqr-bzsqr
  213.     add.l    d3,d0            \ + aconst, new az now in d0 
  214.  
  215.     
  216. \ ******************* TEST FOR PREMATURE EXIT FROM LOOP ******************
  217.  
  218.     add.l    d5,d6            \ azsqr+bzsqr, magnitude is in d6 now
  219.     cmpi.l    #65536,d6        \ check if magnitude > 65536=4 
  220.     bgt        2$            \ if > , exit loop, else... 
  221.     addi.l    #1,d2             \ ... increment loop counter
  222.     cmp.l    a1,d2             \ MandelMax iterations stored in A1
  223.     blt        1$             \ back to loop start
  224. 2$: move.l     tos,-(dsp)
  225.     move.l    d2,tos          \ leave iteration count on tos
  226.  
  227. \ now restore register contents
  228.  
  229.     move.l     a0,d2
  230.     move.l    a2,d5
  231.  
  232.     callcfa    temp.d3
  233.     move.l    $0(org,tos.l),d3     \ restore d3 !!!
  234.     move.l    (dsp)+,tos
  235.     callcfa    temp.d4
  236.     move.l    $0(org,tos.l),d4     \ restore d4 !!!
  237.     move.l    (dsp)+,tos
  238.     callcfa    temp.d6
  239.     move.l    $0(org,tos.l),d6     \ restore d6 !!!
  240.     move.l    (dsp)+,tos
  241. END-CODE
  242.  
  243. \ ************************ end quick guts **********************************
  244.  
  245.