; CALCMAND.ASM - Mandelbrot/Julia Set calculation Routines ; This module runs as part of an overlay with calcfrac.c. ; It must not be called from anywhere other than calcfrac. ; The routines in this code perform Mandelbrot and Julia set ; calculations using 32-bit integer math as opposed to the ; "traditional" floating-point approach. ; This code relies on several tricks to run as quickly as it does. ; One can fake floating point arithmetic by using integer ; arithmetic and keeping track of the implied decimal point ; if things are reasonable -- and in this case, they are. ; I replaced code that looked like: z = x*y with code that ; looks like: ; ix = x * ifudge (outside the loops) ; iy = y * ifudge ; .... ; iz = (ix * iy) / ifudge (inside the loops) ; (and keep remembering that all the integers are "ifudged" bigger) ; The 386 has native 32-bit integer arithmetic, and (briefly) keeps ; 64-bit values around after 32-bit multiplies. If the result is ; divided down right away, you've got 64-bit arithmetic. You just ; have to ensure that the result after the divide is <= 32 bits long. ; CPUs predating the 386 have to emulate 32-bit arithmetic using ; 16-bit arithmetic, which is significantly slower. ; Dividing is slow -- but shifting is fast, and we can select our ; "fudge factor" to be a power of two, permitting us to use that ; method instead. In addition, the 386 can perform 32-bit wide ; shifting -- and even 64-bit shifts with the following logic: ; shdr eax,edx,cl ; shr edx,cl ; so we make sure that our "fudge factor" is a power of 2 and shift ; it down that way. ; Calcmand is hardcoded for a fudge factor of 2**29. ; Bert Tyler ; History since Fractint 16.0 ; (See comments with CJLT in them) ; CJLT=Chris Lusby Taylor who has... ; ; 1. Speeded up 16 bit on 16 bit CPU ; Minor changes, notably prescaling to fg14 before multiplying ; instead of scaling the answer. ; Also, I added overflow detection after adding linit, since it ; seems this could overflow. ; Overall effect is about 10% faster on 386 with debugflag=8088 ; 2. Speeded up 32 bit on 16 bit CPU ; The macro `square' is totally rewritten, as is the logic for 2xy, ; by prescaling x and y to fg31, not fg29. This allows us to do a ; 32 bit multiply in 3, not 4, 16 bit chunks while retaining full ; fg29 accuracy. ; Also, I removed lots of well-meaning but ineffective code handling ; special cases of zeros and tidied up the handling of negative numbers, ; so the routine is quite a bit shorter now and overall throughput of ; Mandel is now over 40% faster on a 386 with debugflag=8088. ; By the way, I was tempted to go the whole hog and replace x*x-y*y ; by (x+y)*(x-y) to reduce 4 16-bit multiplys to 3, but it makes ; escape detection a bit trickier. Another time, maybe. ; ; required for compatibility if Turbo ASM IFDEF ??version MASM51 QUIRKS ENDIF .MODEL medium,c DGROUP group _DATA,_DATA2 .8086 ; these must NOT be in any segment!! ; this get's rid of TURBO-C fixup errors extrn keypressed:far ; this routine is in 'general.asm' extrn getakey:far ; this routine is in 'general.asm' extrn iplot_orbit:far ; this routine is in 'calcfrac.c' extrn scrub_orbit:far ; this routine is in 'calcfrac.c' _DATA2 segment DWORD PUBLIC 'DATA' FUDGEFACTOR equ 29 ; default (non-potential) fudgefactor ; ************************ External variables ***************************** extrn fractype:word ; == 0 if Mandelbrot set, else Julia extrn inside:word ; "inside" color, normally 1 (blue) extrn outside:word ; "outside" color, normally -1 (iter) extrn creal:dword, cimag:dword ; Julia Set Constant extrn delmin:dword ; min increment - precision required extrn maxit:word ; maximum iterations extrn lm:dword ; magnitude bailout limit extrn row:word, col:word ; current pixel to calc extrn color:word ; color calculated for the pixel extrn realcolor:word ; color before inside,etc adjustments extrn reset_periodicity:word ; nonzero if to be reset extrn kbdcount:word ; keyboard counter extrn cpu:word ; cpu type: 86, 186, 286, or 386 extrn dotmode:word extrn show_orbit:word ; "show-orbit" flag extrn orbit_ptr:word ; "orbit pointer" flag extrn periodicitycheck:word ; no periodicity if zero public linitx,linity ; caller sets these public savedmask ; caller sets this ; ************************ Internal variables ***************************** align 4 x dd 0 ; temp value: x y dd 0 ; temp value: y absx dd 0 ; temp value: abs(x) linitx dd 0 ; initial value, set by calcfrac linity dd 0 ; initial value, set by calcfrac savedmask dd 0 ; saved values mask savedx dd 0 ; saved values of X and Y iterations savedy dd 0 ; (for periodicity checks) k dw 0 ; iteration countdown counter oldcolor dw 0 ; prior pixel's escape time k value savedand dw 0 ; AND value for periodicity checks savedincr dw 0 ; flag for incrementing AND value period db 0 ; periodicity, if in the lake _DATA2 ends .CODE ; ***************** Function calcmandasm() ********************************** public calcmandasm FRAME MACRO regs push bp mov bp, sp IRP reg, push reg ENDM ENDM UNFRAME MACRO regs IRP reg, pop reg ENDM pop bp ENDM calcmandasm proc FRAME ; std frame, for TC++ overlays sub ax,ax ; clear ax cmp periodicitycheck,ax ; periodicity checking disabled? je initoldcolor ; yup, set oldcolor 0 to disable it cmp reset_periodicity,ax ; periodicity reset? je short initparms ; inherit oldcolor from prior invocation mov ax,maxit ; yup. reset oldcolor to maxit-250 sub ax,250 ; (avoids slowness at high maxits) initoldcolor: mov oldcolor,ax ; reset oldcolor initparms: mov ax,word ptr creal ; initialize x == creal mov dx,word ptr creal+2 ; ... mov word ptr x,ax ; ... mov word ptr x+2,dx ; ... mov ax,word ptr cimag ; initialize y == cimag mov dx,word ptr cimag+2 ; ... mov word ptr y,ax ; ... mov word ptr y+2,dx ; ... mov ax,maxit ; setup k = maxit inc ax ; (+ 1) mov k,ax ; (decrementing to 0 is faster) cmp fractype,1 ; julia or mandelbrot set? je short dojulia ; julia set - go there ; (Tim wants this code changed so that, for the Mandelbrot, ; Z(1) = (x + iy) + (a + ib). Affects only "fudged" Mandelbrots. ; (for the "normal" case, a = b = 0, and this works, too) ; cmp word ptr x,0 ; Mandelbrot shortcut: ; jne short doeither ; if creal = cimag = 0, ; cmp word ptr x+2,0 ; the first iteration can be emulated. ; jne short doeither ; ... ; cmp word ptr y,0 ; ... ; jne short doeither ; ... ; cmp word ptr y+2,0 ; ... ; jne short doeither ; ... ; dec k ; we know the first iteration passed ; mov dx,word ptr linitx+2 ; copy x = linitx ; mov ax,word ptr linitx ; ... ; mov word ptr x+2,dx ; ... ; mov word ptr x,ax ; ... ; mov dx,word ptr linity+2 ; copy y = linity ; mov ax,word ptr linity ; ... ; mov word ptr y+2,dx ; ... ; mov word ptr y,ax ; ... dec k ; we know the first iteration passed mov dx,word ptr linitx+2 ; add x += linitx mov ax,word ptr linitx ; ... add word ptr x,ax ; ... adc word ptr x+2,dx ; ... mov dx,word ptr linity+2 ; add y += linity mov ax,word ptr linity ; ... add word ptr y,ax ; ... adc word ptr y+2,dx ; ... jmp short doeither ; branch around the julia switch dojulia: ; Julia Set initialization ; "fudge" Mandelbrot start-up values mov ax,word ptr x ; switch x with linitx mov dx,word ptr x+2 ; ... mov bx,word ptr linitx ; ... mov cx,word ptr linitx+2 ; ... mov word ptr x,bx ; ... mov word ptr x+2,cx ; ... mov word ptr linitx,ax ; ... mov word ptr linitx+2,dx ; ... mov ax,word ptr y ; switch y with linity mov dx,word ptr y+2 ; ... mov bx,word ptr linity ; ... mov cx,word ptr linity+2 ; ... mov word ptr y,bx ; ... mov word ptr y+2,cx ; ... mov word ptr linity,ax ; ... mov word ptr linity+2,dx ; ... doeither: ; common Mandelbrot, Julia set code mov period,0 ; claim periodicity of 1 mov savedand,1 ; initial periodicity check mov savedincr,1 ; flag for incrementing periodicity mov word ptr savedx+2,0ffffh; impossible value of "old" x mov word ptr savedy+2,0ffffh; impossible value of "old" y mov orbit_ptr,0 ; clear orbits dec kbdcount ; decrement the keyboard counter jns short nokey ; skip keyboard test if still positive mov kbdcount,10 ; stuff in a low kbd count cmp show_orbit,0 ; are we showing orbits? jne quickkbd ; yup. leave it that way. mov kbdcount,5000 ; else, stuff an appropriate count val cmp cpu,386 ; ("appropriate" to the CPU) je short kbddiskadj ; ... ;; cmp word ptr delmin+2,1 ; is 16-bit math good enough? cmp word ptr delmin+2,8 ; is 16-bit math good enough? ja kbddiskadj ; yes. test less often mov kbdcount,500 ; no. test more often kbddiskadj: cmp dotmode,11 ; disk video? jne quickkbd ; no, leave as is shr kbdcount,1 ; yes, reduce count shr kbdcount,1 ; ... quickkbd: call far ptr keypressed ; has a key been pressed? cmp ax,0 ; ... je nokey ; nope. proceed mov kbdcount,0 ; make sure it goes negative again cmp ax,'o' ; orbit toggle hit? je orbitkey ; yup. show orbits cmp ax,'O' ; orbit toggle hit? jne keyhit ; nope. normal key. orbitkey: call far ptr getakey ; read the key for real mov ax,1 ; reset orbittoggle = 1 - orbittoggle sub ax,show_orbit ; ... mov show_orbit,ax ; ... jmp short nokey ; pretend no key was hit keyhit: mov ax,-1 ; return with -1 mov color,ax ; set color to -1 UNFRAME ; pop stack frame ret ; bail out! nokey: cmp show_orbit,0 ; is orbiting on? jne no16bitcode ; yup. slow down. cmp cpu,386 ; are we on a 386? je short code386bit ; YAY!! 386-class speed! ;; cmp word ptr delmin+2,1 ; OK, we're desperate. 16 bits OK? cmp word ptr delmin+2,8 ; OK, we're desperate. 16 bits OK? ja yes16bitcode ; YAY! 16-bit speed! no16bitcode: call near ptr code32bit ; BOO!! nap time. Full 32 bit math jmp kloopend ; bypass the 386-specific code. yes16bitcode: call near ptr code16bit ; invoke the 16-bit version jmp kloopend ; bypass the 386-specific code. .386 ; 386-specific code starts here code386bit: ;; cmp word ptr delmin+2,3 ; is 16-bit math good enough? cmp word ptr delmin+2,8 ; is 16-bit math good enough? jbe code386_32 ; nope, go do 32 bit stuff IFDEF ??version jmp code386_32 ; TASM screws up IMUL EBX,EBX!! ENDIF ; 16 bit on 386, now we are really gonna move movsx esi,word ptr x+2 ; use SI for X movsx edi,word ptr y+2 ; use DI for Y push ebp mov ebp,-1 shl ebp,FUDGEFACTOR-1 mov cx,FUDGEFACTOR-16 kloop386_16: ; cx=bitshift-16, ebp=overflow.mask mov ebx,esi ; compute (x * x) imul ebx,ebx ; ... test ebx,ebp ; jnz short end386_16 ; (oops. We done.) shr ebx,cl ; get result down to 16 bits mov edx,edi ; compute (y * y) imul edx,edx ; ... test edx,ebp ; say, did we overflow? jnz short end386_16 ; (oops. We done.) shr edx,cl ; get result down to 16 bits mov ax,bx ; compute (x*x - y*y) / fudge sub bx,dx ; for the next iteration add ax,dx ; compute (x*x + y*y) / fudge cmp ax,word ptr lm+2 ; while (xx+yy < lm) jae short end386_16 ; ... imul edi,esi ; compute (y * x) shl edi,1 ; ( * 2 / fudge) sar edi,cl add di,word ptr linity+2 ; (2*y*x) / fudge + linity movsx edi,di ; save as y add bx,word ptr linitx+2 ; (from above) (x*x - y*y)/fudge + linitx movsx esi,bx ; save as x mov ax,oldcolor ; recall the old color cmp ax,k ; check it against this iter jge short chkpd386_16 ; yup. do periodicity check. nonmax386_16: dec k ; while (k < maxit) jnz short kloop386_16 ; try, try again end386_16: pop ebp jmp kloopend ; we done chkpd386_16: mov ax,k ; set up to test for save-time test ax,savedand ; save on 0, check on anything else jz short chksv386_16 ; time to save a new "old" value mov ax,si ; load up x xor ax,word ptr savedx+2 ; does X match? test ax,word ptr savedmask+2 ; truncate to appropriate precision jne short nonmax386_16 ; nope. forget it. mov ax,di ; now test y xor ax,word ptr savedy+2 ; does Y match? test ax,word ptr savedmask+2 ; truncate to appropriate precision jne short nonmax386_16 ; nope. forget it. mov period,1 ; note that we have found periodicity mov k,0 ; pretend maxit reached jmp short end386_16 chksv386_16: mov word ptr savedx+2,si ; save x mov word ptr savedy+2,di ; save y dec savedincr ; time to change the periodicity? jnz short nonmax386_16 ; nope. shl savedand,1 ; well then, let's try this one! inc savedand ; (2**n -1) mov savedincr,4 ; and reset the increment flag jmp short nonmax386_16 ; 32bit on 386: code386_32: mov esi,x ; use ESI for X mov edi,y ; use EDI for Y ; This is the main processing loop. Here, every T-state counts... kloop: ; for (k = 0; k <= maxit; k++) mov eax,esi ; compute (x * x) imul esi ; ... shrd eax,edx,FUDGEFACTOR ; ( / fudge) shr edx,FUDGEFACTOR-1 ; (complete 64-bit shift and check jne short kloopend1 ; bail out if too high mov ebx,eax ; save this for below mov eax,edi ; compute (y * y) imul edi ; ... shrd eax,edx,FUDGEFACTOR ; ( / fudge) shr edx,FUDGEFACTOR-1 ; (complete 64-bit shift and check jne short kloopend1 ; bail out if too high mov ecx,ebx ; compute (x*x - y*y) / fudge sub ebx,eax ; for the next iteration add ecx,eax ; compute (x*x + y*y) / fudge cmp ecx,lm ; while (lr < lm) jae short kloopend1 ; ... mov eax,edi ; compute (y * x) imul esi ; ... shrd eax,edx,FUDGEFACTOR-1 ; ( * 2 / fudge) add eax,linity ; (above) + linity mov edi,eax ; save this as y ; (from the earlier code) ; compute (x*x - y*y) / fudge add ebx,linitx ; + linitx mov esi,ebx ; save this as x mov ax,oldcolor ; recall the old color cmp ax,k ; check it against this iter jge short chkperiod1 nonmax1: dec k ; while (k < maxit) (dec to 0 is faster) jnz short kloop ; while (k < maxit) ... kloopend1: jmp short kloopend32 ; we done. chkperiod1: mov eax,esi xor eax,savedx test eax,savedmask jnz short chksave1 mov eax,edi xor eax,savedy test eax,savedmask jnz short chksave1 mov period,1 ; note that we have found periodicity mov k,0 ; pretend maxit reached jmp short kloopend1 chksave1: mov ax,k test ax,savedand jne short nonmax1 mov savedx,esi mov savedy,edi dec savedincr ; time to change the periodicity? jnz short nonmax1 ; nope. shl savedand,1 ; well then, let's try this one! inc savedand ; (2**n -1) mov savedincr,4 ; and reset the increment flag jmp short nonmax1 kloopend32: .8086 ; 386-specific code ends here kloopend: cmp orbit_ptr,0 ; any orbits to clear? je noorbit2 ; nope. call far ptr scrub_orbit ; clear out any old orbits noorbit2: mov ax,k ; set old color sub ax,10 ; minus 10, for safety mov oldcolor,ax ; and save it as the "old" color mov ax,maxit ; compute color sub ax,k ; (first, re-compute "k") sub kbdcount,ax ; adjust the keyboard count cmp ax,1 ; convert any "outlier" region jge short coloradjust1 ; (where abs(x) > 2 or abs(y) > 2) mov ax,1 ; to look like we ran through coloradjust1: ; at least one loop. mov realcolor,ax ; result before adjustments cmp ax,maxit ; did we max out on iterations? jne short notmax ; nope. mov oldcolor,ax ; set "oldcolor" to maximum cmp inside,0 ; is "inside" >= 0? jl wedone ; nope. leave it at "maxit" mov ax,inside ; reset max-out color to default cmp periodicitycheck,0 ; show periodicity matches? jge wedone ; nope. mov al,period ; reset color to periodicity flag jmp short wedone notmax: cmp outside,0 ; is "outside" >= 0? jl wedone ; nope. leave as realcolor mov ax, outside ; reset to "outside" color wedone: ; mov color,ax ; save the color result UNFRAME ; pop stack frame ret ; and return with color calcmandasm endp ; ******************** Function code16bit() ***************************** ; ; Performs "short-cut" 16-bit math where we can get away with it. ; CJLT has modified it, mostly by preshifting x and y to fg30 from fg29 ; or, since we ignore the lower 16 bits, fg14 from fg13. ; If this shift overflows we are outside x*x+y*y=2, so have escaped. ; Also, he commented out several conditional jumps which he calculated could ; never be taken (e.g. mov ax,si / imul si ;cannot overflow). code16bit proc near mov si,word ptr x+2 ; use SI for X fg13 mov di,word ptr y+2 ; use DI for Y fg13 start16bit: add si,si ;CJLT-Convert to fg14 jo end16bit ;overflows if <-2 or >2 mov ax,si ; compute (x * x) imul si ; Answer is fg14+14-16=fg12 ; cmp dx,0 ;CJLT commented out- ; jl end16bit ;CJLT- imul CANNOT overflow ; mov cx,32-FUDGEFACTOR ;CJLT. FUDGEFACTOR=29 is hard coded loop16bit1: shl ax,1 ; ... rcl dx,1 ; ... jo end16bit ; (oops. overflow) ; loop loop16bit1 ;CJLT...do it once only. dx now fg13. mov bx,dx ; save this for a tad ;ditto for y*y... add di,di ;CJLT-Convert to fg14 jo end16bit ;overflows if <-2 or >2 mov ax,di ; compute (y * y) imul di ; ... ; cmp dx,0 ; say, did we overflow? ; jl end16bit ; (oops. We done.) ; mov cx,32-FUDGEFACTOR ; ( / fudge) ;loop16bit2: shl ax,1 ; ... rcl dx,1 ; ... jo end16bit ; (oops. overflow) ; loop loop16bit2 ; ... mov cx,bx ; compute (x*x - y*y) / fudge sub bx,dx ; for the next iteration add cx,dx ; compute (x*x + y*y) / fudge jo end16bit ; bail out if too high ; js end16bit ; ... cmp cx,word ptr lm+2 ; while (xx+yy < lm) jae end16bit ; ... dec k ; while (k < maxit) jz end16bit ; we done. mov ax,di ; compute (y * x) fg14+14=fg28 imul si ; ... ; mov cx,33-FUDGEFACTOR-2 ; ( * 2 / fudge) ;loop16bit3: shl ax,1 ; ... rcl dx,1 ; ... shl ax,1 ; shift two bits rcl dx,1 ; cannot overflow as |x|<=2, |y|<=2 ; loop loop16bit3 ; ... add dx,word ptr linity+2 ; (2*y*x) / fudge + linity jo end16bit ; bail out if too high mov di,dx ; save as y add bx,word ptr linitx+2 ; (from above) (x*x - y*y)/fudge + linitx jo end16bit ; bail out if too high mov si,bx ; save as x mov ax,oldcolor ; recall the old color cmp ax,k ; check it against this iter jle short nonmax3 ; nope. bypass periodicity check. mov word ptr x+2,si ; save x for periodicity check mov word ptr y+2,di ; save y for periodicity check call checkperiod ; check for periodicity nonmax3: jmp start16bit ; try, try again. end16bit: ; we done. ret code16bit endp ; The following routine checks for periodic loops (a known ; method of decay inside "Mandelbrot Lake", and an easy way to ; bail out of "lake" points quickly). For speed, only the ; high-order sixteen bits of X and Y are checked for periodicity. ; For accuracy, this routine is only fired up if the previous pixel ; was in the lake (which means that the FIRST "wet" pixel was ; detected by the dull-normal maximum iteration process). checkperiod proc near ; periodicity check mov ax,k ; set up to test for save-time test ax,savedand ; save on 0, check on anything else jz checksave ; time to save a new "old" value mov dx,word ptr x+2 ; load up x and dx,word ptr savedmask+2 ; truncate to appropriate precision cmp dx,word ptr savedx+2 ; does X match? jne checkdone ; nope. forget it. mov ax,word ptr x ; load up x and ax,word ptr savedmask ; truncate to appropriate precision cmp ax,word ptr savedx ; does X match? jne checkdone ; nope. forget it. mov dx,word ptr y+2 ; now test y and dx,word ptr savedmask+2 ; truncate to appropriate precision cmp dx,word ptr savedy+2 ; does Y match? jne checkdone ; nope. forget it. mov ax,word ptr y ; load up y and ax,word ptr savedmask ; truncate to appropriate precision cmp ax,word ptr savedy ; does Y match? jne checkdone ; nope. forget it. mov period,1 ; note that we have found periodicity mov k,1 ; pretend maxit reached checksave: mov dx,word ptr x+2 ; load up x and dx,word ptr savedmask+2 ; truncate to appropriate precision mov word ptr savedx+2,dx ; and save it mov ax,word ptr x ; load up x and ax,word ptr savedmask ; truncate to appropriate precision mov word ptr savedx,ax ; and save it mov dx,word ptr y+2 ; load up y and dx,word ptr savedmask+2 ; truncate to appropriate precision mov word ptr savedy+2,dx ; and save it mov ax,word ptr y ; load up y and ax,word ptr savedmask ; truncate to appropriate precision mov word ptr savedy,ax ; and save it dec savedincr ; time to change the periodicity? jnz checkdone ; nope. shl savedand,1 ; well then, let's try this one! inc savedand ; (2**n -1) mov savedincr,4 ; and reset the increment flag checkdone: ret ; we done. checkperiod endp ; ******************** Function code32bit() ***************************** ; ; Perform the 32-bit logic required using 16-bit logic ; ; New twice as fast logic, ; Courtesy of Bill Townsend and Mike Gelvin (CIS:73337,520) ; Even newer, faster still by Chris Lusby Taylor ; who noted that we needn't square the low word if we first multiply ; by 4, since we only need 29 places after the point and this will ; give 30. (We divide answer by two to give 29 bit shift in answer) ; Also, he removed all testing for special cases where a word of operand ; happens to be 0, since testing 65536 times costs more than the saving ; 1 time out of 65536! (I benchmarked it. Just removing the tests speeds ; us up by 3%.) ; ;Note that square returns DI,AX squared in DX,AX now. ; DI,AX is first converted to unsigned fg31 form. ; (For its square to be representable in fg29 (range -4..+3.999) ; DI:AX must be in the range 0..+1.999 which fits neatly into unsigned fg31.) ; This allows us to ignore the part of the answer corresponding to AX*AX as it ; is less than half a least significant bit of the final answer. ; I thought you'd like that. ; ; As we prescaled DI:AX, we need to shift the answer 1 bit to the right to ; end up in fg29 form since 29=(29+2)+(29+2)-32-1 ; However, the mid term AX*DI is needed twice, so the shifts cancel. ; ; Since abs(x) and abs(y) in fg31 form will be needed in calculating 2*X*Y ; we push them onto the stack for later use. ; Note that square does nor affect bl,si,bp ; and leaves highword of argument in di ; but destroys bh,cx square MACRO donepops LOCAL notneg shl ax,1 ;Multiply by 2 to convert to fg30 rcl di,1 ;If this overflows DI:AX was negative jnc notneg not ax ; so negate it not di ; ... add ax,1 ; ... adc di,0 ; ... not bl ; change negswt notneg: shl ax,1 ;Multiply by 2 again to give fg31 rcl di,1 ;If this gives a carry then DI:AX was >=2.0 ;If its high bit is set then DI:AX was >=1.0 ;This is OK, but note that this means that ;DI:AX must now be treated as unsigned. jc donepops push di ; save y or x (in fg31 form) on stack push ax ; ... mul di ;GET MIDDLE PART - 2*A*B mov bh,ah ;Miraculously, it needs no shifting! mov cx,dx mov ax,di mul ax ;SQUARE HIGH HWORD - A*A shl bh,1 ;See if we round up adc ax,1 ;Anyway, add 1 to round up/down accurately adc dx,0 shr dx,1 ;This needs shifting one bit rcr ax,1 add ax,cx ;Add in the 2*A*B term adc dx,0 ENDM ;#EM code32bit proc near ; ; BL IS USED FOR THE "NEGSWT" FLAG ; NEGSWT IS ONE IF EITHER "X" OR "Y" ARE NEGATIVE, BUT NOT BOTH NEGATIVE ; push bp xor bl,bl ; NEGSWT STARTS OFF ZERO ; iteration loop nextit: mov ax,word ptr y ; ax=low(y) mov di,word ptr y+2 ; di=high(y) square done0 ;square y and quit via done0 if it overflows mov si,ax ; square returns results in dx,ax mov bp,dx ; save y*y in bp,si mov ax,word ptr x mov di,word ptr x+2 square done2 ; square x and quit via done2 if it overflows mov cx,ax ; Save low answer in cx. ADD ax,si ; calc y*y + x*x mov ax,bp ADC ax,dx ; ... jno nextxy ; overflow? ;NOTE: The original code tests against lm ;here, but as lm=4<<29 this is the same ;as testing for signed overflow done4: add sp,4 ; discard saved value of |x| fg 31 done2: add sp,4 ; discard saved value of |y| fg 31 done0: pop bp ; restore saved bp ret ;--------------------------------------------------------------------------- nextxy: dec k ; while (k < maxit) jz done4 ; we done. sub cx,si ; subtract y*y from x*x sbb dx,bp ; ... add cx,word ptr linitx ; add "A" adc dx,word ptr linitx+2 ; ... jo done4 ;CJLT-Must detect overflow here ; but increment loop count first mov word ptr x,cx ; store new x = x*x-y*y+a mov word ptr x+2,dx ; ... ; now calculate x*y ; ;More CJLT tricks here. We use the pushed abs(x) and abs(y) in fg31 form ;which, when multiplied, give x*y in fg30, which, luckily, is the same as... ;2*x*y fg29. ;As with squaring, we can ignore the product of the low order words, and still ;be more accurate than the original algorithm. ; pop bp ;Xlow pop di ;Xhigh (already there, actually) pop ax ;Ylow mul di ;Xhigh * Ylow mov bh,ah ;Discard lowest 8 bits of answer mov cx,dx pop ax ;Yhigh mov si,ax ; copy it mul bp ;Xlow * Yhigh xor bp,bp ;Clear answer add bh,ah adc cx,dx adc bp,0 mov ax,si ;Yhigh mul di ;Xhigh * Yhigh shl bh,1 ;round up/down adc ax,cx ;Answer-low adc dx,bp ;Answer-high ;NOTE: The answer is 0..3.9999 in fg29 js done0 ;Overflow if high bit set or bl,bl ; ZERO IF NONE OR BOTH X , Y NEG jz signok ; ONE IF ONLY ONE OF X OR Y IS NEG not ax ; negate result not dx ; ... add ax,1 ; ... adc dx,0 ; ... xor bl,bl ;Clear negswt signok: add ax,word ptr linity adc dx,word ptr linity+2 ; dx,ax = 2(X*Y)+B jo done0 mov word ptr y,ax ; save the new value of y mov word ptr y+2,dx ; ... mov ax,oldcolor ; recall the old color cmp ax,k ; check it against this iter jle short chkmaxit ; nope. bypass periodicity check. call checkperiod ; check for periodicity chkmaxit: cmp show_orbit,0 ; orbiting on? jne horbit ; yep. jmp nextit ;go around again jmpdone0: jmp done0 ; DOES [(X*X)-(Y*Y)+P] BEFORE THE DEC. horbit: push bx ; save my flags mov ax,-1 ; color for plot orbit push ax ; ... push word ptr y+2 ; co-ordinates for plot orbit push word ptr y ; ... push word ptr x+2 ; ... push word ptr x ; ... call far ptr iplot_orbit ; display the orbit add sp,5*2 ; clear out the parameters pop bx ; restore flags jmp nextit ; go around again code32bit endp end