(( Program for iterating simultaneously towards all roots of a polynomial with complex coefficients Andy Korsak 29jun2000 - started out from: GCLOCK.SEQ A simple Graphic Clock program by Tom Zimmer & Robert Smith Fixed bug causing divergence for the polynomials x^n - 1 when n is even -- roles of I & J were reversed in the double loop ajk 05jul2000 )) only forth also definitions 1280 value screen-mwidth 1024 value screen-mheight 400 to screen-width 300 to screen-height \ --------------------------------------------------------------- \ Define the BIT-WINDOW global drawing functions \ --------------------------------------------------------------- Windc demo-dc 2 value bit-originx \ we have a two pixel border around the bitmap 2 value bit-originy 0 value VGA-X \ VGA x coordinate in pixels 0 value VGA-Y \ VGA y coordinate in pixels -1 value prev-x -1 value prev-y : moveto ( x y -- ) 0max screen-height 4 - min swap 0max screen-width 4 - min swap bit-originy + swap bit-originx + swap over prev-x = over prev-y = and IF 2drop ELSE 2dup to prev-y to prev-x MoveTo: demo-dc THEN ; : lineto ( x y -- ) 0max screen-height 4 - min swap 0max screen-width 4 - min swap bit-originy + swap bit-originx + swap over prev-x = over prev-y = and IF 2drop ELSE 2dup to prev-y to prev-x LineTo: demo-dc then ; : line ( x1 y1 x2 y2 -- ) 2swap moveto lineto ; : line-color ( color_object -- ) LineColor: demo-dc ; \ --------------------------------------------------------------- \ Define the BIT-WINDOW window class \ --------------------------------------------------------------- :Class bit-window screenx ( n1 -- n2 ) screen-width 480 */ ; : >screeny ( n1 -- n2 ) screen-width 480 */ ; : xxy-scale ( 1deg scale -- x1 y1 ) >r dup 0 d>f fsin f# 240.0 f* f>d drop >screenx dup r@ center-x */ swap 1 and + center-x + swap 15 + 0 d>f fsin f# 240.0 f* f>d drop >screeny dup r> scale-y */ swap 1 and + negate center-y + ; : xy-scale ( 6deg scale -- x1 y1 ) >r dup sintbl +CELLS @ >screenx dup r@ center-x */ swap 1 and + center-x + swap 15 + sintbl +CELLS @ >screeny dup r> scale-y */ swap 1 and + negate center-y + ; 1 value delay-ms 16 value cdiam 0 value ccolor create colors (( DKGRAY , RED , LTRED , GREEN , LTGREEN , BLUE , LTBLUE , YELLOW , LTYELLOW , MAGENTA , LTMAGENTA , CYAN , LTCYAN , GRAY , WHITE , LTGRAY , : >color ( n1 -- color_object ) 15 and colors +cells @ ; )) WHITE , LTRED , LTYELLOW , LTGREEN , LTCYAN , LTMAGENTA , LTBLUE , BLUE , YELLOW , LTGRAY , LTMAGENTA , CYAN , DKGRAY , GRAY , WHITE , GREEN , : >color ( n1 -- color_object ) 15 and colors +cells @ ; (( : show-circle ( -- ) 1 +TO ccolor ccolor >color line-color 60 0 DO I cdiam xy-scale 2dup 1+ swap 1+ swap line LOOP 5 +TO cdiam cdiam center-x 30 - > IF 16 TO cdiam THEN ; )) : show-border ( -- ) 60 0 \ for minutes on a clock do white line-color i center-x 1 - xy-scale i 1+ center-x 1 - xy-scale line i center-x 12 - xy-scale i 1+ center-x 12 - xy-scale line i 5 mod if ltcyan line-color \ 1 second markers i center-x 12 - xy-scale i center-x 1 - xy-scale line else yellow line-color \ 5 second markers i center-x 20 - xy-scale i center-x 1 - xy-scale line then loop ; 15 value num-roots 1.3e0 fvalue init-magnitude 3 value init-plot-scale init-plot-scale value plot-scale 3.0e0 fvalue root-error 1.0e-6 fvalue root-tolerance 1.0e-14 fvalue float-zero-tolerance B/FLOAT 2* constant COMPLEX-FLOAT create trial-roots COMPLEX-FLOAT num-roots * allot create next-trials COMPLEX-FLOAT num-roots * allot create poly-coefs COMPLEX-FLOAT num-roots 1+ * allot : trial-root ( n -- adr) COMPLEX-FLOAT * trial-roots + ; : next-trial ( n -- adr) COMPLEX-FLOAT * next-trials + ; : poly-coef ( n -- adr) COMPLEX-FLOAT * poly-coefs + ; 3.1415928e0 2.0e0 f* fvalue 2pi 1.1e0 fvalue scalefac 200 value coef-fudge-factor \ a percentage 1000 value angle-tweak \ this is just a resolution level 1000 value ratio-tweak \ " " " " : real-part> ( adr -- | fs: -- r ) f@ ; : imag-part> ( adr -- | fs: -- r ) B/FLOAT + f@ ; : >real-part ( fs: r -- | adr -- ) f! ; : >imag-part ( fs: r -- | adr -- ) B/FLOAT + f! ; : init-poly-coefs poly-coefs num-roots COMPLEX-FLOAT 1+ * erase 1.0e0 num-roots poly-coef >real-part 1.0e0 \ inital scale factor relative to the A0 coeff. num-roots 0 DO init-magnitude f* \ build up scale factor coef-fudge-factor dup random swap 2/ - \ +- percentage 100 + ?dup if s>f f* 1.0e2 f/ endif \ fudged scale factor 100 random 50 > if -1.0e0 f* endif \ flip sign fdup num-roots 1- I - poly-coef >real-part coef-fudge-factor dup random swap 2/ - \ +- percentage 100 + 100 random 50 > if -1.0e0 f* endif \ flip sign ?dup if s>f f* 1.0e2 f/ endif \ fudged scale factor fdup num-roots 1- I - poly-coef >imag-part LOOP fdrop ; : init-trial-roots ( -- ) num-roots 0 DO I s>f 2pi f* num-roots s>f f/ ( angle ) fdup fcos init-magnitude f* I trial-root >real-part fsin init-magnitude f* I trial-root >imag-part LOOP ; : random-angle-tweak ( fs: -- r ) angle-tweak random 1+ s>f angle-tweak 2/ s>f f/ ; : random-ratio-tweak ( fs: -- r ) ratio-tweak random 1+ s>f ratio-tweak 2/ s>f f/ ; : init-next-roots ( -- ) num-roots 0 DO I s>f 2pi f* random-angle-tweak f* num-roots s>f f/ 1.3e0 f* 0.4e0 f+ ( modified angle ) fdup fcos init-magnitude random-ratio-tweak f* f* I trial-root >real-part fsin init-magnitude random-ratio-tweak f* f* I trial-root >imag-part LOOP ; : complex* { adr1 adr2 \ adr1 adr2 -- | fs: -- f_real f_imag } adr1 real-part> adr2 real-part> f* adr1 imag-part> adr2 imag-part> f* f- adr1 real-part> adr2 imag-part> f* adr1 imag-part> adr2 real-part> f* f+ ; 0.0e0 fvalue ftemp0 : complex/ { adr1 adr2 \ adr1 adr2 -- | fs: f_real f_imag } adr2 real-part> fdup f* adr2 imag-part> fdup f* f+ fto ftemp0 adr1 real-part> adr2 real-part> f* adr1 imag-part> adr2 imag-part> f* f+ ftemp0 f/ adr2 real-part> adr1 imag-part> f* adr2 imag-part> adr1 real-part> f* f- ftemp0 f/ ; : complex+ { adr1 adr2 \ adr1 adr2 -- | fs: -- f_real f_imag } adr1 real-part> adr2 real-part> f+ adr1 imag-part> adr2 imag-part> f+ ; : complex- { adr1 adr2 \ adr1 adr2 -- | fs: -- f_real f_imag } adr1 real-part> adr2 real-part> f- adr1 imag-part> adr2 imag-part> f- ; create Zval COMPLEX-FLOAT allot create Pval COMPLEX-FLOAT allot : compute-polynomial ( [complex value at Zval] -- [computed value at Pval] ) Zval Pval COMPLEX-FLOAT cmove \ storages for f-values num-roots 1 DO \ result so far is at P : compute P = P + coeff[N - I] num-roots I - poly-coef ( adr1 ) Pval ( adr1 adr2 ) complex+ ( f_real f_imag ) Pval >imag-part Pval >real-part \ now multiply the result so far by Z Pval Zval complex* ( f_real f_imag ) Pval >imag-part Pval >real-part LOOP 0 poly-coef ( adr1 ) Pval complex+ Pval >imag-part Pval >real-part ; create temp1 COMPLEX-FLOAT allot create temp2 COMPLEX-FLOAT allot create complex_1.0 COMPLEX-FLOAT allot 1.0e0 complex_1.0 >real-part 0.0e0 complex_1.0 >imag-part : compute-next-trial-roots num-roots 0 DO I trial-root ( adr ) Zval COMPLEX-FLOAT cmove compute-polynomial \ Compute product of (Ztrial_i - Ztrial_j), j <> i complex_1.0 temp1 COMPLEX-FLOAT cmove num-roots 0 DO J I = 0= if J trial-root I trial-root complex- temp2 >imag-part temp2 >real-part temp1 temp2 complex* temp1 >imag-part temp1 >real-part endif LOOP \ Compute Ztrial_next_i = Ztrial_i - P(Ztrial_i) / (above product) temp1 real-part> fabs float-zero-tolerance f< temp1 imag-part> fabs float-zero-tolerance f< or if float-zero-tolerance fdup temp1 f! temp1 B/FLOAT + f! endif Pval temp1 complex/ temp2 >imag-part temp2 >real-part I trial-root temp2 complex- I next-trial >imag-part I next-trial >real-part LOOP ; 0 value max-x 0 value max-y : plot-roots 0 TO ccolor ccolor >color line-color num-roots 0 DO screen-width 2/ s>f I trial-root real-part> f* plot-scale s>f f/ f>s screen-width 2/ + \ counting pixels up from the left 4 max screen-width 4 - min \ x coordinate screen-height 2/ s>f I trial-root imag-part> f* plot-scale s>f f/ f>s screen-height 2/ + screen-height swap - \ counting pixels up from the bottom 4 max screen-height 4 - min \ y coordinate moveto screen-width 2/ s>f I next-trial real-part> f* plot-scale s>f f/ f>s screen-width 2/ + \ counting pixels up from the left dup max-x > if dup TO max-x endif dup negate max-x > if dup negate TO max-x endif 4 max screen-width 4 - min \ x coordinate screen-height 2/ s>f I next-trial imag-part> f* plot-scale s>f f/ f>s screen-height 2/ + screen-height swap - \ counting pixels up from the bottom dup max-y > if dup TO max-y endif dup negate max-y > if dup negate TO max-y endif 4 max screen-height 4 - min \ y coordinate lineto 1 +TO ccolor ccolor >color line-color LOOP ; : dump-coefs CR num-roots 0 DO CR I 2 .r space I poly-coef real-part> f. I poly-coef imag-part> f. LOOP ; : dump-roots CR num-roots 0 DO I trial-root real-part> f. I trial-root imag-part> f. LOOP ." err=" root-error f. ; : accum-error ( fs: r -- ) fabs root-error f+ FTO root-error ; : update-trial-roots next-trials trial-roots num-roots COMPLEX-FLOAT * cmove ; : plot-trial-roots screen-width 2 / 1- TO center-x screen-height 2 / 1- TO center-y \ calibrate screen center center-x center-x center-y */ TO scale-y \ calibrate aspect ratio plot-roots Refresh: PolyRoots WINPAUSE ; : ptr plot-trial-roots ; : erase-plot 0 0 screen-width screen-height WHITE PrinterFillArea: ThePrinter 0 0 screen-width screen-height BLACK FillArea: demo-dc ; : init-calc init-poly-coefs init-trial-roots init-next-roots ; : converged? ( -- f ) root-error root-tolerance f< ; : compute-error 0.0e0 FTO root-error num-roots 0 DO I trial-root I next-trial complex- ( real imag ) accum-error accum-error LOOP ; : ?adjust-scale ( -- ) max-x plot-scale > dup if max-x to plot-scale endif max-y plot-scale > dup if max-y to plot-scale endif or if erase-plot ." downscaled " endif ; \ --------------------------------------------------------------- \ Top Level program starts here \ --------------------------------------------------------------- 0 value c-height 0 value c-width : setup-plot init-plot-scale to plot-scale Start: POLYROOTS erase-plot RANDOM-INIT \ initialize random numbers screen-width 2 / 1 - TO center-x screen-height 2 / 1 - TO center-y \ calibrate screen center white line-color \ default color=white screen-width to c-width screen-height to c-height ( show-circle ) show-border ; 0 value iter-num : PolynomialRoots \ { \ c-width c-height -- } setup-plot init-calc compute-error dump-coefs 0 TO iter-num begin key? 0= converged? 0= and while 1 +TO iter-num compute-next-trial-roots compute-error plot-trial-roots cr ." Iter#" iter-num 3 .r ." error = " root-error f. update-trial-roots c-width c-height screen-width screen-height d= 0= if 1 to delay-ms 16 to cdiam screen-width to c-width screen-height to c-height ( show-circle ) show-border then \ ?adjust-scale repeat plot-trial-roots dump-roots ; : pr PolynomialRoots ; : prs begin pr 200 ms key 27 = until ; ' Prs turnkey PolyRoots \ build an application on disk