Category : C Source Code
Archive   : FRASRC18.ZIP
Filename : LYAPUNOV.ASM

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; 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; i ; for (count=0; 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 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