home *** CD-ROM | disk | FTP | other *** search
/ Computerworld 1996 March / Computerworld_1996-03_cd.bin / idg_cd3 / grafika / fraktaly / frasr192 / lyapunov.asm < prev    next >
Assembly Source File  |  1994-04-09  |  21KB  |  511 lines
  1. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  2. ; Wesley Loewer's attempt at writing his own
  3. ; 80x87 assembler implementation of Lyapunov fractals
  4. ; based on lyapunov_cycles_in_c() in miscfrac.c
  5. ;
  6. ; Nicholas Wilt, April 1992, originally wrote an 80x87 assembler
  7. ; implementation of Lyapunov fractals, but I couldn't get his
  8. ; lyapunov_cycles() to work right with fractint.
  9. ; So I'm starting mostly from scratch with credits to Nicholas Wilt's
  10. ; code marked with NW.
  11.  
  12. .8086
  13. .8087
  14. .MODEL medium,c
  15.  
  16. ; Increase the overflowcheck value at your own risk.
  17. ; Bigger value -> check less often -> faster run time, increased risk of overflow.
  18. ; I've had failures with 6 when using "set no87" as emulation can be less forgiving.
  19. overflowcheck EQU 5
  20.  
  21. EXTRN   Population:QWORD,Rate:QWORD
  22. EXTRN   colors:WORD, maxit:DWORD, lyaLength:WORD
  23. EXTRN   lyaRxy:WORD, LogFlag:WORD
  24. EXTRN   fpu:WORD
  25. EXTRN   filter_cycles:DWORD
  26.  
  27. .DATA
  28. ALIGN 2
  29. BigNum      DD      100000.0
  30. half        DD      0.5                 ; for rounding to integers
  31. e_10        DQ      22026.4657948       ; e^10
  32. e_neg10     DQ      4.53999297625e-5    ; e^(-10)
  33.  
  34. .CODE
  35.  
  36. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  37. ; BifurcLambda - not used in lyapunov.asm, but used in miscfrac.c
  38. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  39. PUBLIC  BifurcLambda
  40. BifurcLambda    PROC
  41. ;   Population = Rate * Population * (1 - Population);
  42. ;   return (fabs(Population) > BIG);
  43.         push    bp              ; Establish stack frame
  44.         mov     bp,sp           ;
  45.         sub     sp,2            ; Local for 80x87 flags
  46.         fld     Population      ; Get population
  47.         fld     st              ; Population Population
  48.         fld1                    ; 1 Population Population
  49.         fsubr                   ; 1 - Population Population
  50.         fmul                    ; Population * (1 - Population)
  51.         fmul    Rate            ; Rate * Population * (1 - Population)
  52.         fst     Population      ; Store in Population
  53.         fabs                    ; Take absolute value
  54.         fcomp   BigNum          ; Compare to 100000.0
  55.         fstsw   [bp-2]          ; Return 1 if greater than 100000.0
  56.         mov     ax,[bp-2]       ;
  57.         sahf                    ;
  58.         ja      Return1         ;
  59.         xor     ax,ax           ;
  60.         jmp     short Done      ;
  61. Return1:
  62.         mov     ax,1            ;
  63. Done:   inc     sp              ; Deallocate locals
  64.         inc     sp              ;
  65.         pop     bp              ; Restore stack frame
  66.         ret                     ; Return
  67. BifurcLambda    ENDP
  68.  
  69. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  70. ; OneMinusExpMacro - calculates e^x
  71. ; parameter : x in st(0)
  72. ; returns   : e^x in st(0)
  73. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  74. OneMinusExpMacro   MACRO
  75. LOCAL not_387, positive_x, was_positive_x
  76. LOCAL long_way_386, long_way_286
  77. LOCAL positive_x_long, was_positive_x_long, done
  78.  
  79. ; To take e^x, one can use 2^(x*log_2(e)), however on on 8087/287
  80. ; the valid domain for x is 0 <= x <= 0.5  while the 387/487 allows
  81. ; -1.0 <= x <= +1.0.
  82.  
  83. ; I wrote the following 387+ specific code to take advantage of the
  84. ; improved domain in the f2xm1 command, but the 287- code seems to work
  85. ; about as fast.
  86.  
  87.         cmp     fpu,387                 ; is fpu >= 387 ?
  88.         jl      not_387
  89.  
  90. ;.386   ; a 387 or better is present so we might as well use .386/7 here
  91. ;.387   ; so "fstsw ax" can be used.  Note that the same can be accomplished 
  92. .286    ; with .286/7 which is recognized by more more assemblers.
  93. .287    ; (ie: MS QuickAssembler doesn't recognize .386/7)
  94.  
  95. ; |x*log_2(e)| must be <= 1.0 in order to work with 386 or greater.
  96. ; It usually is, so it's worth taking these extra steps to check.
  97.         fldl2e                          ; log_2(e) x
  98.         fmul                            ; x*log_2(e)
  99. ;begining of short_way code
  100.         fld1                            ; 1 x*log_2(e)
  101.         fld     st(1)                   ; x*log_2(e) 1 x*log_2(e)
  102.         fabs                            ; |x*log_2(e)| 1 x*log_2(e)
  103.         fcompp                          ; x*log_2(e)
  104.         fstsw   ax
  105.         sahf
  106.         ja      long_way_386            ; do it the long way
  107.         f2xm1                           ; e^x-1=2^(x*log_2(e))-1
  108.         fchs                            ; 1-e^x which is what we wanted anyway
  109.     jmp     done                    ; done here.
  110.  
  111. long_way_386:
  112. ; mostly taken from NW's code
  113.         fld     st                      ; x x
  114.         frndint                         ; i x
  115.         fsub    st(1),st                ; i x-i
  116.         fxch                            ; x-i i
  117.         f2xm1                           ; e^(x-i)-1 i
  118.         fld1                            ; 1 e^(x-i)-1 i
  119.         fadd                            ; e^(x-i) i
  120.         fscale                          ; e^x i
  121.         fstp    st(1)                   ; e^x
  122.         fld1                            ; 1 e^x
  123.         fsubr                           ; 1-e^x
  124.         jmp     done
  125.  
  126. not_387:
  127. .8086
  128. .8087
  129. ; |x*log_2(e)| must be <= 0.5 in order to work with 286 or less.
  130. ; It usually is, so it's worth taking these extra steps to check.
  131.         fldl2e                          ; log_2(e) x
  132.         fmul                            ; x*log_2(e)
  133.  
  134. ;begining of short_way code
  135.         fld     st                      ; x*log_2(e) x*log_2(e)
  136.         fabs                            ; |x*log_2(e)| x*log_2(e)
  137.         fcomp   half                    ; x*log_2(e)
  138.         fstsw   status_word
  139.         fwait
  140.         mov     ax,status_word
  141.         sahf
  142.         ja      long_way_286            ; do it the long way
  143.  
  144. ; 286 or less requires x>=0, if x is negative, use e^(-x) = 1/e^x
  145.         ftst                            ; x
  146.         fstsw   status_word
  147.         fwait
  148.         mov     ax,status_word
  149.         sahf
  150.         jae     positive_x
  151.         fabs                            ; -x
  152.         f2xm1                           ; e^(-x)-1=2^(-x*log_2(e))-1
  153.         fld1                            ; 1 e^(-x)-1
  154.         fadd    st,st(1)                ; e^(-x) e^(-x)-1
  155.         fdiv                            ; 1-e^x=(e^(-x)-1)/e^(-x)
  156.         jmp     SHORT done              ; done here.
  157.  
  158. positive_x:
  159.         f2xm1                           ; e^x-1=2^(x*log_2(e))-1
  160.         fchs                            ; 1-e^x which is what we wanted anyway
  161.         jmp     SHORT done              ; done here.
  162.  
  163. long_way_286:
  164. ; mostly taken from NW's code
  165.         fld     st                      ; x x
  166.         frndint                         ; i x
  167.         fsub    st(1),st                ; i x-i
  168.         fxch                            ; x-i i
  169. ; x-i could be pos or neg
  170. ; 286 or less requires x>=0, if negative, use e^(-x) = 1/e^x
  171.         ftst
  172.         fstsw   status_word
  173.         fwait
  174.         mov     ax,status_word
  175.         sahf
  176.         jae     positive_x_long_way
  177.         fabs                            ; i-x
  178.         f2xm1                           ; e^(i-x)-1 i
  179.         fld1                            ; 1 e^(i-x)-1 i
  180.         fadd    st(1),st                ; 1 e^(i-x) i
  181.         fdivr                           ; e^(x-i) i
  182.         fscale                          ; e^x i
  183.         fstp    st(1)                   ; e^x
  184.         fld1                            ; 1 e^x
  185.         fsubr                           ; 1-e^x
  186.         jmp     SHORT done              ; done here.
  187.  
  188. positive_x_long_way:                    ; x-i
  189.         f2xm1                           ; e^(x-i)-1 i
  190.         fld1                            ; 1 e^(x-i)-1 i
  191.         fadd                            ; e^(x-i) i
  192.         fscale                          ; e^x i
  193.         fstp    st(1)                   ; e^x
  194.         fld1                            ; 1 e^x
  195.         fsubr                           ; 1-e^x
  196. done:
  197. ENDM    ; OneMinusExpMacro
  198.  
  199. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  200. ; modeled from lyapunov_cycles_in_c() in miscfrac.c
  201. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  202. ;int lyapunov_cycles_in_c(int filter_cycles, double a, double b) {
  203. ;    int color, count, i, lnadjust;
  204. ;    double lyap, total, temp;
  205.  
  206. PUBLIC lyapunov_cycles
  207. lyapunov_cycles  PROC  USES si di, a:QWORD, b:QWORD
  208.     LOCAL color:WORD, i:WORD, lnadjust:WORD
  209.     LOCAL halfmax:DWORD, status_word:WORD
  210.     LOCAL total:QWORD
  211.  
  212.  
  213. ;    for (i=0; i<filter_cycles; i++) {
  214. ;        for (count=0; count<lyaLength; count++) {
  215. ;            Rate = lyaRxy[count] ? a : b;
  216. ;            if (curfractalspecific->orbitcalc()) {
  217. ;                overflow = TRUE;
  218. ;                goto jumpout;
  219. ;                }
  220. ;            }
  221. ;        }
  222.  
  223. ; Popalation is left on the stack during most of this procedure
  224.         fld     Population              ; Population
  225.         mov     bx,1                    ; using bx for overflowcheck counter
  226.         mov     ax,word ptr filter_cycles
  227.         mov     dx,word ptr filter_cycles+2
  228.         add     ax,1
  229.         adc     dx,0                    ; add 1 because dec at beginning
  230.         mov     si, ax                  ; using si for low part of i
  231.         mov     i, dx
  232. filter_cycles_top:
  233.         dec     si
  234.         jne     hop_skip_and_jump       ; si > 0 :
  235.         cmp     i,0
  236.         je      filter_cycles_bottom
  237.         dec     i
  238. hop_skip_and_jump:
  239.         mov     cx, lyaLength           ; using cx for count
  240.         mov     di,OFFSET lyaRxy        ; use di as lyaRxy[] index
  241.         ; no need to compare cx with lyaLength at top of loop
  242.         ; since lyaLength is guaranteed to be > 0
  243.  
  244. filter_cycles_count_top:
  245.         cmp     WORD PTR [di],0         ; lyaRxy[count] ?
  246.         jz      filter_cycles_use_b
  247.         fld     a                       ; if lyaRxy[count]==true,  use a
  248.         jmp     SHORT filter_cycles_used_a
  249. filter_cycles_use_b:
  250.         fld     b                       ; if lyaRxy[count]==false, use b
  251. filter_cycles_used_a:
  252.         ; leave Rate on stack for use in BifurcLambdaMacro
  253.         ; BifurcLambdaMacro               ; returns in flag register
  254.     
  255.         fld     st(1)                   ; Population Rate Population
  256.         fmul    st,st                   ; Population^2 Rate Population
  257.         fsubp   st(2),st                ; Rate Population-Population^2
  258.                                         ;      =Population*(1-Population)
  259.         fmul                            ; Rate*Population*(1-Population)
  260.         dec     bx                      ; decrement overflowcheck counter
  261.         jnz     filter_cycles_skip_overflow_check ; can we skip check?
  262.         fld     st                      ; NewPopulation NewPopulation
  263.         fabs                            ; Take absolute value
  264.         fcomp   BigNum                  ; Compare to 100000.0
  265.         fstsw   status_word             ;
  266.         fwait
  267.         mov     bx,overflowcheck        ; reset overflowcheck counter
  268.         mov     ax,status_word          ;
  269.         sahf                            ;  NewPopulation
  270.  
  271.         ja      overflowed_with_Pop
  272. filter_cycles_skip_overflow_check:
  273.         add     di,2                    ; di += sizeof(*lyaRxy)
  274.         loop    filter_cycles_count_top ; if --cx > 0 then loop
  275.  
  276.         jmp     filter_cycles_top       ; let's do it again
  277.  
  278. filter_cycles_bottom:
  279.  
  280.  
  281. ;    total = 1.0;
  282. ;    lnadjust = 0;
  283.  
  284.         fld1
  285.         fstp    total
  286.         mov     lnadjust,0
  287.  
  288. ;    for (i=0; i < maxit/2; i++) {
  289. ;        for (count = 0; count < lyaLength; count++) {
  290. ;            Rate = lyaRxy[count] ? a : b;
  291. ;            if (curfractalspecific->orbitcalc()) {
  292. ;                overflow = TRUE;
  293. ;                goto jumpout;
  294. ;                }
  295. ;            temp = fabs(Rate-2.0*Rate*Population);
  296. ;                if ((total *= (temp))==0) {
  297. ;                overflow = TRUE;
  298. ;                goto jumpout;
  299. ;                }
  300. ;            }
  301. ;        while (total > 22026.4657948) {
  302. ;            total *= 0.0000453999297625;
  303. ;            lnadjust += 10;
  304. ;            }
  305. ;        while (total < 0.0000453999297625) {
  306. ;            total *= 22026.4657948;
  307. ;            lnadjust -= 10;
  308. ;            }
  309. ;        }
  310.  
  311.         ; don't forget Population is still on stack
  312.         mov     ax,word ptr maxit       ; calculate halfmax
  313.         mov     dx,word ptr maxit+2     ; calculate halfmax
  314.         shr     dx,1
  315.         rcr     ax,1
  316.         mov     word ptr halfmax,ax
  317.         mov     word ptr halfmax+2,dx
  318.         add     ax,1
  319.         adc     dx,0                   ; add 1 because dec at beginning
  320.         mov     i,dx                   ; count down from halfmax
  321.         mov     si,ax                  ; using si for low part of i
  322. halfmax_top:
  323.         dec     si
  324.         jne     init_halfmax_count     ; si > 0 ?
  325.         cmp     i,0
  326.         je      step_to_halfmax_bottom
  327.         dec     i
  328.         jmp     short init_halfmax_count      ; yes, continue on with loop
  329. step_to_halfmax_bottom:
  330.         jmp     halfmax_bottom          ; if not, end loop
  331. init_halfmax_count:
  332.         mov     cx, lyaLength           ; using cx for count
  333.         mov     di,OFFSET lyaRxy        ; use di as lyaRxy[] index
  334.         ; no need to compare cx with lyaLength at top of loop
  335.         ; since lyaLength is guaranteed to be > 0
  336.  
  337. halfmax_count_top:
  338.         cmp     WORD PTR [di],0         ; lyaRxy[count] ?
  339.         jz      halfmax_use_b
  340.         fld     a                       ; if lyaRxy[count]==true,  use a
  341.         jmp     SHORT halfmax_used_a
  342. halfmax_use_b:
  343.         fld     b                       ; if lyaRxy[count]==false, use b
  344. halfmax_used_a:
  345.         ; save Rate, but leave it on stack for use in BifurcLambdaMacro
  346.         fst     Rate                    ; save for not_overflowed use
  347.         ;BifurcLambdaMacro               ; returns in flag register
  348.  
  349.         fld     st(1)                   ; Population Rate Population
  350.         fmul    st,st                   ; Population^2 Rate Population
  351.         fsubp   st(2),st                ; Rate Population-Population^2
  352.                                         ;      =Population*(1-Population)
  353.         fmul                            ; Rate*Population*(1-Population)
  354.         dec     bx                      ; decrement overflowcheck counter
  355.         jnz     not_overflowed          ; can we skip overflow check?
  356.         fld     st                      ; NewPopulation NewPopulation
  357.         fabs                            ; Take absolute value
  358.         fcomp   BigNum                  ; Compare to 100000.0
  359.         fstsw   status_word             ;
  360.         fwait
  361.         mov     bx,overflowcheck        ; reset overflowcheck counter
  362.         mov     ax,status_word          ;
  363.         sahf                            ;
  364.  
  365.         jbe     not_overflowed
  366.  
  367.         ; putting overflowed_with_Pop: here saves 2 long jumps in inner loops
  368. overflowed_with_Pop:
  369.         fstp    Population              ; save Population and clear stack
  370.         jmp     color_zero              ;
  371.  
  372. not_overflowed:                         ; Population
  373.         fld     st                      ; Population Population
  374.         ; slightly faster _not_ to fld Rate here
  375.         fmul    Rate                    ; Rate*Population Population
  376.         fadd    st,st                   ; 2.0*Rate*Population Population
  377.         fsubr   Rate                    ; Rate-2.0*Rate*Population Population
  378.         fabs                            ; fabs(Rate-2.0*Rate*Population) Population
  379.         fmul    total                   ; total*fabs(Rate-2.0*Rate*Population) Population
  380.         ftst                            ; compare to 0
  381.         fstp    total                   ; save the new total
  382.         fstsw   status_word             ; Population
  383.         fwait
  384.         mov     ax,status_word
  385.         sahf
  386.         jz      overflowed_with_Pop     ; if total==0, then treat as overflow
  387.  
  388.         add     di,2                    ; di += sizeof(*lyaRxy)
  389.         loop    halfmax_count_top       ; if --cx > 0 then loop
  390.  
  391.         fld     total                   ; total Population
  392. too_big_top:
  393.         fcom    e_10                    ; total Population
  394.         fstsw   status_word
  395.         fwait
  396.         mov     ax,status_word
  397.         sahf
  398.         jna     too_big_bottom
  399.         fmul    e_neg10                 ; total*e_neg10 Population
  400.         add     lnadjust,10
  401.         jmp     SHORT too_big_top
  402. too_big_bottom:
  403.  
  404. too_small_top:                          ; total Population
  405.         fcom    e_neg10                 ; total Population
  406.         fstsw   status_word
  407.         fwait
  408.         mov     ax,status_word
  409.         sahf
  410.         jnb     too_small_bottom
  411.         fmul    e_10                    ; total*e_10 Population
  412.         sub     lnadjust,10
  413.         jmp     SHORT too_small_top
  414. too_small_bottom:
  415.  
  416.         fstp    total                   ; save as total
  417.         jmp     halfmax_top             ; let's do it again
  418.  
  419. halfmax_bottom:
  420. ; better make another check, just to be sure
  421.                                         ; NewPopulation
  422.         cmp     bx,overflowcheck        ; was overflow just checked?
  423.         jl      last_overflowcheck      ; if not, better check one last time
  424.         fstp    Population              ; save Population and clear stack
  425.         jmp     short jumpout           ; skip overflowcheck
  426.  
  427. last_overflowcheck:
  428.         fst     Population              ; save new Population
  429.         fabs                            ; |NewPopulation|
  430.         fcomp   BigNum                  ; Compare to 100000.0
  431.         fstsw   status_word             ;
  432.         fwait
  433.         mov     ax,status_word
  434.         sahf
  435.         ja      color_zero              ; overflowed
  436.  
  437. jumpout:
  438.  
  439. ;    if (overflow || total <= 0 || (temp = log(total) + lnadjust) > 0)
  440. ;        color = 0;
  441. ;    else {
  442. ;        if (LogFlag)
  443. ;            lyap = -temp/((double) lyaLength*i);
  444. ;        else
  445. ;            lyap = 1 - exp(temp/((double) lyaLength*i));
  446. ;        color = 1 + (int)(lyap * (colors-1));
  447. ;        }
  448. ;    return color;
  449. ;}
  450.  
  451. ; no use testing for the overflow variable as you
  452. ; cannot get here and have overflow be true
  453.         fld     total                   ; total
  454.         ftst                            ; is total <= 0 ?
  455.         fstsw   status_word
  456.         fwait
  457.         mov     ax,status_word
  458.         sahf
  459.         ja      total_not_neg           ; if not, continue
  460.         fstp    st                      ; pop total from stack
  461.         jmp     SHORT color_zero        ; if so, set color to 0
  462.  
  463. total_not_neg:                          ; total is still on stack
  464.         fldln2                          ; ln(2) total
  465.         fxch                            ; total ln(2)
  466.         fyl2x                           ; ln(total)=ln(2)*log_2(total)
  467.         fiadd   lnadjust                ; ln(total)+lnadjust
  468.         ftst                            ; compare against 0
  469.         fstsw   status_word
  470.         fwait
  471.         mov     ax,status_word
  472.         sahf
  473.         jbe     not_positive
  474.         fstp    st                      ; pop temp from stack
  475.  
  476. color_zero:
  477.         xor     ax,ax                   ; return color 0
  478.         jmp     thats_all
  479.  
  480. not_positive:                           ; temp is still on stack
  481.         fild    lyaLength               ; lyaLength temp
  482.         fimul   halfmax                 ; lyaLength*halfmax temp
  483.         fdiv                            ; temp/(lyaLength*halfmax)
  484.         cmp     LogFlag,0               ; is LogFlag set?
  485.         jz      LogFlag_not_set         ; if !LogFlag goto LogFlag_not_set:
  486.         fchs                            ; -temp/(lyaLength*halfmax)
  487.         jmp     calc_color
  488.  
  489. LogFlag_not_set:                        ; temp/(lyaLength*halfmax)
  490.         OneMinusExpMacro                ; 1-exp(temp/(lyaLength*halfmax))
  491. calc_color:
  492.         ; what is now left on the stack is what the C code calls lyap
  493.                                         ; lyap
  494.         mov     ax,colors               ; colors
  495.         dec     ax                      ; colors-1
  496.         mov     i,ax                    ; temporary storage
  497.         fimul   i                       ; lyap*(colors-1)
  498.         ; the "half" makes ASM round like C does
  499.         fsub    half                    ; sub 0.5 for rounding purposes
  500.         fistp   i                       ; temporary = (int)(lyap*(colors-1))
  501.         fwait                           ; one moment please...
  502.         mov     ax,i                    ; ax = temporary
  503.         inc     ax                      ; ax = 1 + (int)(lyap * (colors-1));
  504.  
  505. thats_all:
  506.         mov     color, ax
  507.         ret
  508. lyapunov_cycles  endp
  509.  
  510. END
  511.