File: dos\lyapunov.asm
1 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
2 ; 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
3 ; lyapunov_cycles() to work right with fractint.
4 ; So I'm starting mostly from scratch with credits to Nicholas Wilt's
5 ; code marked with NW.
6
7 .8086
8 .8087
9 .MODEL medium,c
10
11 ; Increase the overflowcheck value at your own risk.
12 ; Bigger value -> check less often -> faster run time, increased risk of overflow.
13 ; 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:DWORD
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)
14
15 ; |x*log_2(e)| must be <= 1.0 in order to work with 386 or greater.
16 ; 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
17 fld st ; x x
18 frndint ; i x
19 fsub st(1),st ; i x-i
20 fxch ; x-i i
21 f2xm1 ; e^(x-i)-1 i
22 fld1 ; 1 e^(x-i)-1 i
23 fadd ; e^(x-i) i
24 fscale ; e^x i
25 fstp st(1) ; e^x
26 fld1 ; 1 e^x
27 fsubr ; 1-e^x
28 jmp done
29
30 not_387:
31 .8086
32 .8087
33 ; |x*log_2(e)| must be <= 0.5 in order to work with 286 or less.
34 ; 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
35 fld st ; x x
36 frndint ; i x
37 fsub st(1),st ; i x-i
38 fxch ; x-i i
39 ; x-i could be pos or neg
40 ; 286 or less requires x>=0, if negative, use e^(-x) = 1/e^x
41 ftst
42 fstsw status_word
43 fwait
44 mov ax,status_word
45 sahf
46 jae positive_x_long_way
47 fabs ; i-x
48 f2xm1 ; e^(i-x)-1 i
49 fld1 ; 1 e^(i-x)-1 i
50 fadd st(1),st ; 1 e^(i-x) i
51 fdivr ; e^(x-i) i
52 fscale ; e^x i
53 fstp st(1) ; e^x
54 fld1 ; 1 e^x
55 fsubr ; 1-e^x
56 jmp SHORT done ; done here.
57
58 positive_x_long_way: ; x-i
59 f2xm1 ; e^(x-i)-1 i
60 fld1 ; 1 e^(x-i)-1 i
61 fadd ; e^(x-i) i
62 fscale ; e^x i
63 fstp st(1) ; e^x
64 fld1 ; 1 e^x
65 fsubr ; 1-e^x
66 done:
67 ENDM ; OneMinusExpMacro
68
69 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
70 ; modeled from lyapunov_cycles_in_c() in miscfrac.c
71 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
72 ;int lyapunov_cycles_in_c(int filter_cycles, double a, double b) {
73 ; int color, count, i, lnadjust;
74 ; double lyap, total, temp;
75
76 PUBLIC lyapunov_cycles
77 lyapunov_cycles PROC USES si di, a:QWORD, b:QWORD
78 LOCAL color:WORD, i:WORD, lnadjust:WORD
79 LOCAL halfmax:DWORD, status_word:WORD
80 LOCAL total:QWORD
81
82
83 ; for (i=0; i 84 ; for (count=0; count 85 ; Rate = lyaRxy[count] ? a : b;
86 ; if (curfractalspecific->orbitcalc()) {
87 ; overflow = TRUE;
88 ; goto jumpout;
89 ; }
90 ; }
91 ; }
92
93 ; Popalation is left on the stack during most of this procedure
94 fld Population ; Population
95 mov bx,1 ; using bx for overflowcheck counter
96 mov ax,word ptr filter_cycles
97 mov dx,word ptr filter_cycles+2
98 add ax,1
99 adc dx,0 ; add 1 because dec at beginning
100 mov si, ax ; using si for low part of i
101 mov i, dx
102 filter_cycles_top:
103 dec si
104 jne hop_skip_and_jump ; si > 0 :
105 cmp i,0
106 je filter_cycles_bottom
107 dec i
108 hop_skip_and_jump:
109 mov cx, lyaLength ; using cx for count
110 mov di,OFFSET lyaRxy ; use di as lyaRxy[] index
111 ; no need to compare cx with lyaLength at top of loop
112 ; since lyaLength is guaranteed to be > 0
113
114 filter_cycles_count_top:
115 cmp WORD PTR [di],0 ; lyaRxy[count] ?
116 jz filter_cycles_use_b
117 fld a ; if lyaRxy[count]==true, use a
118 jmp SHORT filter_cycles_used_a
119 filter_cycles_use_b:
120 fld b ; if lyaRxy[count]==false, use b
121 filter_cycles_used_a:
122 ; leave Rate on stack for use in BifurcLambdaMacro
123 ; BifurcLambdaMacro ; returns in flag register
124
125 fld st(1) ; Population Rate Population
126 fmul st,st ; Population^2 Rate Population
127 fsubp st(2),st ; Rate Population-Population^2
128 ; =Population*(1-Population)
129 fmul ; Rate*Population*(1-Population)
130 dec bx ; decrement overflowcheck counter
131 jnz filter_cycles_skip_overflow_check ; can we skip check?
132 fld st ; NewPopulation NewPopulation
133 fabs ; Take absolute value
134 fcomp BigNum ; Compare to 100000.0
135 fstsw status_word ;
136 fwait
137 mov bx,overflowcheck ; reset overflowcheck counter
138 mov ax,status_word ;
139 sahf ; NewPopulation
140
141 ; ja overflowed_with_Pop
142 jna filter_cycles_skip_overflow_check
143 jmp overflowed_with_Pop
144 filter_cycles_skip_overflow_check:
145 add di,2 ; di += sizeof(*lyaRxy)
146 loop filter_cycles_count_top ; if --cx > 0 then loop
147
148 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
149 mov ax,word ptr maxit ; calculate halfmax
150 mov dx,word ptr maxit+2 ; calculate halfmax
151 shr dx,1
152 rcr ax,1
153 mov word ptr halfmax,ax
154 mov word ptr halfmax+2,dx
155 add ax,1
156 adc dx,0 ; add 1 because dec at beginning
157 mov i,dx ; count down from halfmax
158 mov si,ax ; using si for low part of i
159 halfmax_top:
160 dec si
161 jne init_halfmax_count ; si > 0 ?
162 cmp i,0
163 je step_to_halfmax_bottom
164 dec i
165 jmp short init_halfmax_count ; yes, continue on with loop
166 step_to_halfmax_bottom:
167 jmp halfmax_bottom ; if not, end loop
168 init_halfmax_count:
169 mov cx, lyaLength ; using cx for count
170 mov di,OFFSET lyaRxy ; use di as lyaRxy[] index
171 ; no need to compare cx with lyaLength at top of loop
172 ; since lyaLength is guaranteed to be > 0
173
174 halfmax_count_top:
175 cmp WORD PTR [di],0 ; lyaRxy[count] ?
176 jz halfmax_use_b
177 fld a ; if lyaRxy[count]==true, use a
178 jmp SHORT halfmax_used_a
179 halfmax_use_b:
180 fld b ; if lyaRxy[count]==false, use b
181 halfmax_used_a:
182 ; save Rate, but leave it on stack for use in BifurcLambdaMacro
183 fst Rate ; save for not_overflowed use
184 ;BifurcLambdaMacro ; returns in flag register
185
186 fld st(1) ; Population Rate Population
187 fmul st,st ; Population^2 Rate Population
188 fsubp st(2),st ; Rate Population-Population^2
189 ; =Population*(1-Population)
190 fmul ; Rate*Population*(1-Population)
191 dec bx ; decrement overflowcheck counter
192 jnz not_overflowed ; can we skip overflow check?
193 fld st ; NewPopulation NewPopulation
194 fabs ; Take absolute value
195 fcomp BigNum ; Compare to 100000.0
196 fstsw status_word ;
197 fwait
198 mov bx,overflowcheck ; reset overflowcheck counter
199 mov ax,status_word ;
200 sahf ;
201
202 jbe not_overflowed
203
204 ; putting overflowed_with_Pop: here saves 2 long jumps in inner loops
205 overflowed_with_Pop:
206 fstp Population ; save Population and clear stack
207 jmp color_zero ;
208
209 not_overflowed: ; Population
210 fld st ; Population Population
211 ; slightly faster _not_ to fld Rate here
212 fmul Rate ; Rate*Population Population
213 fadd st,st ; 2.0*Rate*Population Population
214 fsubr Rate ; Rate-2.0*Rate*Population Population
215 fabs ; fabs(Rate-2.0*Rate*Population) Population
216 fmul total ; total*fabs(Rate-2.0*Rate*Population) Population
217 ftst ; compare to 0
218 fstp total ; save the new total
219 fstsw status_word ; Population
220 fwait
221 mov ax,status_word
222 sahf
223 jz overflowed_with_Pop ; if total==0, then treat as overflow
224
225 add di,2 ; di += sizeof(*lyaRxy)
226 loop halfmax_count_top ; if --cx > 0 then loop
227
228 fld total ; total Population
229 too_big_top:
230 fcom e_10 ; total Population
231 fstsw status_word
232 fwait
233 mov ax,status_word
234 sahf
235 jna too_big_bottom
236 fmul e_neg10 ; total*e_neg10 Population
237 add lnadjust,10
238 jmp SHORT too_big_top
239 too_big_bottom:
240
241 too_small_top: ; total Population
242 fcom e_neg10 ; total Population
243 fstsw status_word
244 fwait
245 mov ax,status_word
246 sahf
247 jnb too_small_bottom
248 fmul e_10 ; total*e_10 Population
249 sub lnadjust,10
250 jmp SHORT too_small_top
251 too_small_bottom:
252
253 fstp total ; save as total
254 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 word ptr LogFlag,0 ; is LogFlag set?
jz LogFlag_not_set ; if !LogFlag goto LogFlag_not_set:
cmp word ptr LogFlag+2,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