;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Wesley Loewer's attempt at writing his own ; 80x87 assembler implementation of Lyapunov fractals ; based on lyapunov_cycles_in_c() in miscfrac.c ; ; Nicholas Wilt, April 1992, originally wrote an 80x87 assembler ; implementation of Lyapunov fractals, but I couldn't get his ; lyapunov_cycles() to work right with fractint. ; So I'm starting mostly from scratch with credits to Nicholas Wilt's ; code marked with NW. .8086 .8087 .MODEL medium,c ; Increase the overflowcheck value at your own risk. ; Bigger value -> check less often -> faster run time, increased risk of overflow. ; I've had failures with 6 when using "set no87" as emulation can be less forgiving. overflowcheck EQU 5 EXTRN Population:QWORD,Rate:QWORD EXTRN colors:WORD, maxit:WORD, lyaLength:WORD EXTRN lyaRxy:WORD, LogFlag:WORD EXTRN fpu:WORD EXTRN filter_cycles:WORD .DATA ALIGN 2 BigNum DD 100000.0 half DD 0.5 ; for rounding to integers e_10 DQ 22026.4657948 ; e^10 e_neg10 DQ 4.53999297625e-5 ; e^(-10) .CODE ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; BifurcLambda - not used in lyapunov.asm, but used in miscfrac.c ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; PUBLIC BifurcLambda BifurcLambda PROC ; Population = Rate * Population * (1 - Population); ; return (fabs(Population) > BIG); push bp ; Establish stack frame mov bp,sp ; sub sp,2 ; Local for 80x87 flags fld Population ; Get population fld st ; Population Population fld1 ; 1 Population Population fsubr ; 1 - Population Population fmul ; Population * (1 - Population) fmul Rate ; Rate * Population * (1 - Population) fst Population ; Store in Population fabs ; Take absolute value fcomp BigNum ; Compare to 100000.0 fstsw [bp-2] ; Return 1 if greater than 100000.0 mov ax,[bp-2] ; sahf ; ja Return1 ; xor ax,ax ; jmp short Done ; Return1: mov ax,1 ; Done: inc sp ; Deallocate locals inc sp ; pop bp ; Restore stack frame ret ; Return BifurcLambda ENDP ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; OneMinusExpMacro - calculates e^x ; parameter : x in st(0) ; returns : e^x in st(0) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; OneMinusExpMacro MACRO LOCAL not_387, positive_x, was_positive_x LOCAL long_way_386, long_way_286 LOCAL positive_x_long, was_positive_x_long, done ; To take e^x, one can use 2^(x*log_2(e)), however on on 8087/287 ; the valid domain for x is 0 <= x <= 0.5 while the 387/487 allows ; -1.0 <= x <= +1.0. ; I wrote the following 387+ specific code to take advantage of the ; improved domain in the f2xm1 command, but the 287- code seems to work ; about as fast. cmp fpu,387 ; is fpu >= 387 ? jl not_387 ;.386 ; a 387 or better is present so we might as well use .386/7 here ;.387 ; so "fstsw ax" can be used. Note that the same can be accomplished .286 ; with .286/7 which is recognized by more more assemblers. .287 ; (ie: MS QuickAssembler doesn't recognize .386/7) ; |x*log_2(e)| must be <= 1.0 in order to work with 386 or greater. ; It usually is, so it's worth taking these extra steps to check. fldl2e ; log_2(e) x fmul ; x*log_2(e) ;begining of short_way code fld1 ; 1 x*log_2(e) fld st(1) ; x*log_2(e) 1 x*log_2(e) fabs ; |x*log_2(e)| 1 x*log_2(e) fcompp ; x*log_2(e) fstsw ax sahf ja long_way_386 ; do it the long way f2xm1 ; e^x-1=2^(x*log_2(e))-1 fchs ; 1-e^x which is what we wanted anyway jmp done ; done here. long_way_386: ; mostly taken from NW's code fld st ; x x frndint ; i x fsub st(1),st ; i x-i fxch ; x-i i f2xm1 ; e^(x-i)-1 i fld1 ; 1 e^(x-i)-1 i fadd ; e^(x-i) i fscale ; e^x i fstp st(1) ; e^x fld1 ; 1 e^x fsubr ; 1-e^x jmp done not_387: .8086 .8087 ; |x*log_2(e)| must be <= 0.5 in order to work with 286 or less. ; It usually is, so it's worth taking these extra steps to check. fldl2e ; log_2(e) x fmul ; x*log_2(e) ;begining of short_way code fld st ; x*log_2(e) x*log_2(e) fabs ; |x*log_2(e)| x*log_2(e) fcomp half ; x*log_2(e) fstsw status_word fwait mov ax,status_word sahf ja long_way_286 ; do it the long way ; 286 or less requires x>=0, if x is negative, use e^(-x) = 1/e^x ftst ; x fstsw status_word fwait mov ax,status_word sahf jae positive_x fabs ; -x f2xm1 ; e^(-x)-1=2^(-x*log_2(e))-1 fld1 ; 1 e^(-x)-1 fadd st,st(1) ; e^(-x) e^(-x)-1 fdiv ; 1-e^x=(e^(-x)-1)/e^(-x) jmp SHORT done ; done here. positive_x: f2xm1 ; e^x-1=2^(x*log_2(e))-1 fchs ; 1-e^x which is what we wanted anyway jmp SHORT done ; done here. long_way_286: ; mostly taken from NW's code fld st ; x x frndint ; i x fsub st(1),st ; i x-i fxch ; x-i i ; x-i could be pos or neg ; 286 or less requires x>=0, if negative, use e^(-x) = 1/e^x ftst fstsw status_word fwait mov ax,status_word sahf jae positive_x_long_way fabs ; i-x f2xm1 ; e^(i-x)-1 i fld1 ; 1 e^(i-x)-1 i fadd st(1),st ; 1 e^(i-x) i fdivr ; e^(x-i) i fscale ; e^x i fstp st(1) ; e^x fld1 ; 1 e^x fsubr ; 1-e^x jmp SHORT done ; done here. positive_x_long_way: ; x-i f2xm1 ; e^(x-i)-1 i fld1 ; 1 e^(x-i)-1 i fadd ; e^(x-i) i fscale ; e^x i fstp st(1) ; e^x fld1 ; 1 e^x fsubr ; 1-e^x done: ENDM ; OneMinusExpMacro ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; modeled from lyapunov_cycles_in_c() in miscfrac.c ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;int lyapunov_cycles_in_c(int filter_cycles, double a, double b) { ; int color, count, i, lnadjust; ; double lyap, total, temp; PUBLIC lyapunov_cycles lyapunov_cycles PROC USES si di, a:QWORD, b:QWORD LOCAL color:WORD, i:WORD, lnadjust:WORD LOCAL halfmax:WORD, status_word:WORD LOCAL total:QWORD ; for (i=0; iorbitcalc()) { ; overflow = TRUE; ; goto jumpout; ; } ; } ; } ; Popalation is left on the stack during most of this procedure fld Population ; Population mov dx,1 ; using dx for overflowcheck counter xor si, si ; using si for i filter_cycles_top: cmp si, filter_cycles ; si < filter_cycles ? jnl filter_cycles_bottom ; if not, end loop mov cx, lyaLength ; using cx for count mov di,OFFSET lyaRxy ; use di as lyaRxy[] index ; no need to compare cx with lyaLength at top of loop ; since lyaLength is guaranteed to be > 0 filter_cycles_count_top: cmp WORD PTR [di],0 ; lyaRxy[count] ? jz filter_cycles_use_b fld a ; if lyaRxy[count]==true, use a jmp SHORT filter_cycles_used_a filter_cycles_use_b: fld b ; if lyaRxy[count]==false, use b filter_cycles_used_a: ; leave Rate on stack for use in BifurcLambdaMacro ; BifurcLambdaMacro ; returns in flag register fld st(1) ; Population Rate Population fmul st,st ; Population^2 Rate Population fsubp st(2),st ; Rate Population-Population^2 ; =Population*(1-Population) fmul ; Rate*Population*(1-Population) dec dx ; decrement overflowcheck counter jnz filter_cycles_skip_overflow_check ; can we skip check? fld st ; NewPopulation NewPopulation fabs ; Take absolute value fcomp BigNum ; Compare to 100000.0 fstsw status_word ; fwait mov dx,overflowcheck ; reset overflowcheck counter mov ax,status_word ; sahf ; NewPopulation ja overflowed_with_Pop filter_cycles_skip_overflow_check: add di,2 ; di += sizeof(*lyaRxy) loop filter_cycles_count_top ; if --cx > 0 then loop inc si jmp filter_cycles_top ; let's do it again filter_cycles_bottom: ; total = 1.0; ; lnadjust = 0; fld1 fstp total mov lnadjust,0 ; for (i=0; i < maxit/2; i++) { ; for (count = 0; count < lyaLength; count++) { ; Rate = lyaRxy[count] ? a : b; ; if (curfractalspecific->orbitcalc()) { ; overflow = TRUE; ; goto jumpout; ; } ; temp = fabs(Rate-2.0*Rate*Population); ; if ((total *= (temp))==0) { ; overflow = TRUE; ; goto jumpout; ; } ; } ; while (total > 22026.4657948) { ; total *= 0.0000453999297625; ; lnadjust += 10; ; } ; while (total < 0.0000453999297625) { ; total *= 22026.4657948; ; lnadjust -= 10; ; } ; } ; don't forget Population is still on stack mov ax,maxit ; calculate halfmax shr ax,1 mov halfmax,ax xor si, si ; using si for i halfmax_top: cmp si, halfmax ; si < halfmax ? ; too far for a short jump, need a near jump here jl init_halfmax_count ; yes, continue on with loop jmp halfmax_bottom ; if not, end loop init_halfmax_count: mov cx, lyaLength ; using cx for count mov di,OFFSET lyaRxy ; use di as lyaRxy[] index ; no need to compare cx with lyaLength at top of loop ; since lyaLength is guaranteed to be > 0 halfmax_count_top: cmp WORD PTR [di],0 ; lyaRxy[count] ? jz halfmax_use_b fld a ; if lyaRxy[count]==true, use a jmp SHORT halfmax_used_a halfmax_use_b: fld b ; if lyaRxy[count]==false, use b halfmax_used_a: ; save Rate, but leave it on stack for use in BifurcLambdaMacro fst Rate ; save for not_overflowed use ;BifurcLambdaMacro ; returns in flag register fld st(1) ; Population Rate Population fmul st,st ; Population^2 Rate Population fsubp st(2),st ; Rate Population-Population^2 ; =Population*(1-Population) fmul ; Rate*Population*(1-Population) dec dx ; decrement overflowcheck counter jnz not_overflowed ; can we skip overflow check? fld st ; NewPopulation NewPopulation fabs ; Take absolute value fcomp BigNum ; Compare to 100000.0 fstsw status_word ; fwait mov dx,overflowcheck ; reset overflowcheck counter mov ax,status_word ; sahf ; jbe not_overflowed ; putting overflowed_with_Pop: here saves 2 long jumps in inner loops overflowed_with_Pop: fstp Population ; save Population and clear stack jmp color_zero ; not_overflowed: ; Population fld st ; Population Population ; slightly faster _not_ to fld Rate here fmul Rate ; Rate*Population Population fadd st,st ; 2.0*Rate*Population Population fsubr Rate ; Rate-2.0*Rate*Population Population fabs ; fabs(Rate-2.0*Rate*Population) Population fmul total ; total*fabs(Rate-2.0*Rate*Population) Population ftst ; compare to 0 fstp total ; save the new total fstsw status_word ; Population fwait mov ax,status_word sahf jz overflowed_with_Pop ; if total==0, then treat as overflow add di,2 ; di += sizeof(*lyaRxy) loop halfmax_count_top ; if --cx > 0 then loop fld total ; total Population too_big_top: fcom e_10 ; total Population fstsw status_word fwait mov ax,status_word sahf jna too_big_bottom fmul e_neg10 ; total*e_neg10 Population add lnadjust,10 jmp SHORT too_big_top too_big_bottom: too_small_top: ; total Population fcom e_neg10 ; total Population fstsw status_word fwait mov ax,status_word sahf jnb too_small_bottom fmul e_10 ; total*e_10 Population sub lnadjust,10 jmp SHORT too_small_top too_small_bottom: fstp total ; save as total inc si jmp halfmax_top ; let's do it again halfmax_bottom: ; better make another check, just to be sure ; NewPopulation cmp dx,overflowcheck ; was overflow just checked? jl last_overflowcheck ; if not, better check one last time fstp Population ; save Population and clear stack jmp short jumpout ; skip overflowcheck last_overflowcheck: fst Population ; save new Population fabs ; |NewPopulation| fcomp BigNum ; Compare to 100000.0 fstsw status_word ; fwait mov ax,status_word sahf ja color_zero ; overflowed jumpout: ; if (overflow || total <= 0 || (temp = log(total) + lnadjust) > 0) ; color = 0; ; else { ; if (LogFlag) ; lyap = -temp/((double) lyaLength*i); ; else ; lyap = 1 - exp(temp/((double) lyaLength*i)); ; color = 1 + (int)(lyap * (colors-1)); ; } ; return color; ;} ; no use testing for the overflow variable as you ; cannot get here and have overflow be true fld total ; total ftst ; is total <= 0 ? fstsw status_word fwait mov ax,status_word sahf ja total_not_neg ; if not, continue fstp st ; pop total from stack jmp SHORT color_zero ; if so, set color to 0 total_not_neg: ; total is still on stack fldln2 ; ln(2) total fxch ; total ln(2) fyl2x ; ln(total)=ln(2)*log_2(total) fiadd lnadjust ; ln(total)+lnadjust ftst ; compare against 0 fstsw status_word fwait mov ax,status_word sahf jbe not_positive fstp st ; pop temp from stack color_zero: xor ax,ax ; return color 0 jmp thats_all not_positive: ; temp is still on stack fild lyaLength ; lyaLength temp fimul halfmax ; lyaLength*halfmax temp fdiv ; temp/(lyaLength*halfmax) cmp LogFlag,0 ; is LogFlag set? jz LogFlag_not_set ; if !LogFlag goto LogFlag_not_set: fchs ; -temp/(lyaLength*halfmax) jmp calc_color LogFlag_not_set: ; temp/(lyaLength*halfmax) OneMinusExpMacro ; 1-exp(temp/(lyaLength*halfmax)) calc_color: ; what is now left on the stack is what the C code calls lyap ; lyap mov ax,colors ; colors dec ax ; colors-1 mov i,ax ; temporary storage fimul i ; lyap*(colors-1) ; the "half" makes ASM round like C does fsub half ; sub 0.5 for rounding purposes fistp i ; temporary = (int)(lyap*(colors-1)) fwait ; one moment please... mov ax,i ; ax = temporary inc ax ; ax = 1 + (int)(lyap * (colors-1)); thats_all: mov color, ax ret lyapunov_cycles endp END