home *** CD-ROM | disk | FTP | other *** search
- ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
- ; 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:DWORD, lyaLength:WORD
- EXTRN lyaRxy:WORD, LogFlag:WORD
- EXTRN fpu:WORD
- EXTRN filter_cycles:DWORD
-
- .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:DWORD, status_word:WORD
- LOCAL total:QWORD
-
-
- ; for (i=0; i<filter_cycles; i++) {
- ; for (count=0; count<lyaLength; count++) {
- ; Rate = lyaRxy[count] ? a : b;
- ; if (curfractalspecific->orbitcalc()) {
- ; overflow = TRUE;
- ; goto jumpout;
- ; }
- ; }
- ; }
-
- ; Popalation is left on the stack during most of this procedure
- fld Population ; Population
- mov bx,1 ; using bx for overflowcheck counter
- mov ax,word ptr filter_cycles
- mov dx,word ptr filter_cycles+2
- add ax,1
- adc dx,0 ; add 1 because dec at beginning
- mov si, ax ; using si for low part of i
- mov i, dx
- filter_cycles_top:
- dec si
- jne hop_skip_and_jump ; si > 0 :
- cmp i,0
- je filter_cycles_bottom
- dec i
- hop_skip_and_jump:
- 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 bx ; 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 bx,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
-
- 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,word ptr maxit ; calculate halfmax
- mov dx,word ptr maxit+2 ; calculate halfmax
- shr dx,1
- rcr ax,1
- mov word ptr halfmax,ax
- mov word ptr halfmax+2,dx
- add ax,1
- adc dx,0 ; add 1 because dec at beginning
- mov i,dx ; count down from halfmax
- mov si,ax ; using si for low part of i
- halfmax_top:
- dec si
- jne init_halfmax_count ; si > 0 ?
- cmp i,0
- je step_to_halfmax_bottom
- dec i
- jmp short init_halfmax_count ; yes, continue on with loop
- step_to_halfmax_bottom:
- 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 bx ; 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 bx,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
- jmp halfmax_top ; let's do it again
-
- halfmax_bottom:
- ; better make another check, just to be sure
- ; NewPopulation
- cmp bx,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
-
