📄 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<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 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
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -