Artifact b180c91a7c9f4fc60f65583e7ce93705c02a654578fd2facb6b9c76716a26af9:
- File
r34/xmpl/roots.tst
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 28858) [annotate] [blame] [check-ins using] [more...]
% Tests of the root finding package. % Author: Stanley L. Kameny (stan%valley.uucp@rand.org) % Answers are rounded to the value given by rootacc (default = 6) unless % higher accuracy is needed to separate roots. Root order and format % may differ from that given here, but values should agree. % In the following, problems 20), 78) and 82) are time consuming and % been commented out to speed up the test. roots x; % To load roots package. write "This is Roots Package test ", symbolic roots!-mod$ % Simple root finding. showtime; % 1) multiple real and imaginary roots plus two real roots. zz:= (x-3)**2*(100x**2+113)**2*(1000000x-10000111)*(x-1); roots zz; %{X=1.06301*I,X=1.06301*I,X=-1.06301*I,X=-1.06301*I, %X=3.0,X=3.0,X=1,X=10.0001} (rootacc caused rounding to 6 places) % Accuracy is increased whenever necessary to separate distinct roots. % 2) accuracy increase to 7 required for two roots. zz:=(x**2+1)*(x-2)*(1000000x-2000001); roots zz; %{X=2.0,X=I,X= - I,X=2.000001} % 3) accuracy increase to 8 required. zz:= (x-3)*(10000000x-30000001); roots zz; %{X=3.0,X=3.0000001} % 4) accuracy increase required here to separate repeated root from % simple root. zz := (x-3)*(1000000x-3000001)*(x-3)*(1000000x-3241234); roots zz; %{X=3.0,X=3.0,X=3.000001,X=3.24123} % other simple examples % 5) five real roots with widely different spacing. zz:= (x-1)*(10x-11)*(x-1000)*(x-1001)*(x-100000); roots zz; %{X=1,X=1.1,X=1000.0,X=1001.0,X=100000.0} % 6) a cluster of 5 roots in complex plane in vicinity of x=1. zz:= (x-1)*(10000x**2-20000x+10001)*(10000x**2-20000x+9999); roots zz; %{X=1,X=1.01,X=0.99,X=1 + 0.01*I,X=1 - 0.01*I} % 7) four closely spaced real roots. zz := (x-1)*(100x-101)*(100x-102)*(100x-103); roots zz; %{X=1.02,X=1.01,X=1,X=1.03} % 8) five closely spaced roots, 3 real + 1 complex pair. zz := (x-1)*(100x-101)*(100x-102)*(100x**2-200x+101); roots zz; %{X=1.01,X=1,X=1.02,X=1 + 0.1*I,X=1 - 0.1*I} % 9) symmetric cluster of 5 roots, 3 real + 1 complex pair. zz := (x-2)*(10000x**2-40000x+40001)*(10000x**2-40000x+39999); roots zz; %{X=2.0,X=2.01,X=1.99,X=2.0 + 0.01*I,X=2.0 - 0.01*I} % 10) closely spaced real and complex pair. ss:= (x-2)*(100000000x**2-400000000x+400000001); roots ss; %{X=2.0,X=2.0 + 0.0001*I,X=2.0 - 0.0001*I} % 11) Zero roots and multiple roots cause no problem. % Multiple roots are shown when the switch multiroot is on %(normally on.) zz:= x*(x-1)**2*(x-4)**3*(x**2+1); roots zz; %{X=0,X=4.0,X=4.0,X=4.0,X=1,X=1,X=I,X= - I} % 12) nearestroot will find a single root "near" a value, real or % complex. nearestroot(zz,2i); %{X=I} % More difficult examples. % Three examples in which root scaling is needed in the complex % iteration process. % 13) nine roots, 3 real and 3 complex pairs. zz:= x**9-45x-2; roots zz; %{X= - 0.0444444,X=1.61483, % X=0.00555 + 1.60944*I,X=0.00555 - 1.60944*I, % X= - 1.60371, % X=1.14348 + 1.13804*I,X=1.14348 - 1.13804*I, % X=-1.13237 + 1.13805*I,X=-1.13237 - 1.13805*I} % In the next two examples, there are complex roots with extremely % small real parts (new capability in Mod 1.91.) % 14) nine roots, 1 real and 4 complex pairs. zz:= x**9-9999x**2-0.01; roots zz; %{X=5.0E-29 + 0.00100005*I,X=5.0E-29 - 0.00100005*I, %X=3.72754, %X=-0.829456 + 3.63408*I, X=-0.829456 + 3.63408*I, %X=-3.3584 + 1.61732*I, X=-3.3584 - 1.61732*I, %X=2.32408 + 2.91431*I, X=2.32408 - 2.91431*I} % rootacc 7 produces 7 place accuracy; answers will print in bigfloat % format if floating point print >6 digits is not implemented. % 15) nine roots, 1 real and 4 complex pairs. rootacc 7; zz:= x**9-500x**2-0.001; roots zz; %{X=1.6E-26 + 0.001414214*I,X=1.6E-26 - 0.001414214*I, % X=-0.540677 + 2.368861*I, X=-0.540677 - 2.368861*I, % X=-2.189157 + 1.054242*I, X=-2.189157 - 1.054242*I, % X=1.514944 + 1.899679*I, X=1.514944 - 1.899679*I} % the famous Wilkinson "ill-conditioned" polynomial and its family. % 16) W(6) four real roots plus one complex pair. zz:= 10000*(for j:=1:6 product(x+j))+27x**5; roots zz; %{X= - 2.950367,X=-4.452438 + 0.021235*I,X=-4.452438 - 0.021235*I, % X= - 0.9999775,X= - 2.003647,X= - 6.143833} % 17) W(8) 4 real roots plus 2 complex pairs. zz:= 1000*(for j:=1:8 product(x+j))+2x**7; roots zz; %{X=-4.295858 + 0.28151*I,X=-4.295858 - 0.28151*I, % X= - 2.982725, % X=-6.494828 + 1.015417*I,X=-6.494828 - 1.015417*I, % X= - 0.9999996,X= - 2.000356,X= - 8.437546} % 18) W(10) 6 real roots plus 2 complex pairs. zz:=1000*(for j:= 1:10 product (x+j))+x**9; roots zz; %{X= - 4.616444,X=-6.046279 + 1.134321*I,X=-6.046279 - 1.134321*I, % X= - 4.075943,X= - 2.998063, % X=-8.70405 + 1.691061*I,X=-8.70405 - 1.691061*I, % X= -1,X= - 2.000013,X= - 10.80988} % 19) W(12) 6 real roots plus 3 complex pairs. zz:= 10000*(for j:=1:12 product(x+j))+4x**11; roots zz; %{X=-5.985629 + 0.809425*I,X=-5.985629 - 0.809425*I, % X= - 4.880956,X= - 4.007117, % X=-7.953917 + 1.948001*I,X=-7.953917 - 1.948001*I, % X= -1,X=-11.02192 + 2.23956*I,X=-11.02192 - 2.23956*I, % X= - 2.0,X= - 2.999902,X= - 13.1895} % 20) W(20) 10 real roots plus 5 complex pairs. (The original problem) % This example is commented out, since it takes significant time: % zz:= x**19+10**7*for j:=1:20 product (x+j); roots zz; %{X=-10.12155 + 0.6013*I,X=-10.12155 - 0.6013*I, % X= - 8.928803, % X=-11.82101 + 1.59862*I,X=-11.82101 - 1.59862*I, % X= - 8.006075,X= - 6.999746, % X=-14.01105 + 2.44947*I,X=-14.01105 - 2.44947*I, % X=-1,X= - 6.000006, % X=-19.45964 + 1.87436*I,X=-19.45964 - 1.87436*I, % X= - 2.0,X= - 5.0, % X=-16.72504 + 2.73158*I,X=-16.72504 - 2.73158*I, % X= - 3.0,X= - 4.0,X= - 20.78881} rootacc 6; % 21) Finding one of a cluster of 8 roots. zz:= (10**16*(x-1)**8-1); nearestroot(zz,2); %{X=1.01} % 22) Six real roots spaced 0.01 apart. c := 100; zz:= (x-1)*for i:=1:5 product (c*x-(c+i)); roots zz; %{X=1.03,X=1.02,X=1.04,X=1.0,X=1.01,X=1.05} % 23) Six real roots spaced 0.001 apart. c := 1000; zz:= (x-1)*for i:=1:5 product (c*x-(c+i)); roots zz; %{X=1.003,X=1.002,X=1.004,X=1,X=1.001,X=1.005} % 24) Five real roots spaced 0.0001 apart. c := 10000; zz:= (x-1)*for i:=1:4 product (c*x-(c+i)); roots zz; %{X=1.0002,X=1.0003,X=1,X=1.0001,X=1.0004} % 25) A cluster of 9 roots, 5 real, 2 complex pairs; spacing 0.1. zz:= (x-1)*(10**8*(x-1)**8-1); roots zz; %{X=1,X=1.1,X=1 + 0.1*I,X=1 - 0.1*I,X=0.9, % X=0.929289 + 0.070711*I,X=0.929289 - 0.070711*I, % X=1.07071 + 0.07071*I,X=1.07071 - 0.07071*I} % 26) Same, but with spacing 0.01. zz:= (x-1)*(10**16*(x-1)**8-1); roots zz; %{X=1,X=1.01,X=1 + 0.01*I,X=1 - 0.01*I,X=0.99, % X=0.992929 + 0.007071*I,X=0.992929 - 0.007071*I, % X=1.00707 + 0.00707*I,X=1.00707 - 0.00707*I} % 27) Spacing reduced to 0.001. zz:= (x-1)*(10**24*(x-1)**8-1); roots zz; %{X=1,X=1.001,X=1 + 0.001*I,X=1 - 0.001*I,X=0.999, % X=0.999293 + 0.000707*I,X=0.999293 - 0.000707*I, % X=1.00071 + 0.000707*I,X=1.00071 - 0.000707*I} % 28) Eight roots divided into two clusters. zz:= (10**8*(x-1)**4-1)*(10**8*(x+1)**4-1); roots zz; %{X=1.01,X= - 1.01,X=1 + 0.01*I,X=-1 - 0.01*I, % X=-1 + 0.01*I,X=1 - 0.01*I,X=0.99,X= - 0.99} % 29) A cluster of 8 roots in a different configuration. zz:= (10**8*(x-1)**4-1)*(10**8*(100x-102)**4-1); roots zz; %{X=1.01, X=1.0199, X=1.02 + 0.0001*I, X=1.02 - 0.0001*I, % X=1 + 0.01*I, X=1 - 0.01*I, X=0.99, X=1.0201} % 30) A cluster of 8 complex roots. zz:= ((10x-1)**4+1)*((10x+1)**4+1); roots zz; %{X=0.0292893 + 0.0707107*I,X=-0.0292893 - 0.0707107*I, % X=-0.0292893 + 0.0707107*I,X=0.0292893 - 0.0707107*I, % X=0.170711 + 0.070711*I,X=-0.170711 - 0.070711*I, % X=-0.170711 + 0.070711*I,X=0.170711 - 0.070711*I} % In these examples, accuracy increase is required to separate a % repeated root from a simple root. % 31) Using allroots; zz:= (x-4)*(x-3)**2*(1000000x-3000001); roots zz; %{X=3.0,X=3.0,X=3.000001,X=4.0} % 32) Using realroots; realroots zz; %{X=3.0,X=3.0,X=3.000001,X=4.0} comment Tests of new capabilities in mod 1.87 for handling complex polynomials and polynomials with very small imaginary parts or very small real roots. A few real examples are shown, just to demonstrate that these still work; % 33) A trivial complex case (but degrees 1 and 2 are special cases); zz:= x-i; roots zz; %{X=I} % 34) Real case. zz:= y-7; roots zz; %{Y=7.0} % 35) Roots with small imaginary parts (new capability); zz := 10**16*(x**2-2x+1)+1; roots zz; %{X=1 + 0.00000001*I,X=1 - 0.00000001*I} % 36) One real, one complex root. zz:=(x-9)*(x-5i-7); roots zz; %{X=9.0,X=7.0 + 5.0*I} % 37) Three real roots. zz:= (x-1)*(x-2)*(x-3); roots zz; %{X=2.0,X=1,X=3.0} % 38) 2 real + 1 imaginary root. zz:=(x**2-8)*(x-5i); roots zz; %{X=2.82843,X= - 2.82843,X=5.0*I} % 39) 2 complex roots. zz:= (x-1-2i)*(x+2+3i); roots zz; %{X=1 + 2.0*I,X=-2.0 - 3.0*I} % 40) 2 irrational complex roots. zz:= x**2+(3+2i)*x+7i; roots zz; %{X=0.14936 - 2.21259*I,X=-3.14936 + 0.21259*I} % 41) 2 complex roots of very different magnitudes with small imaginary % parts. zz:= x**2+(1000000000+12i)*x-1000000000; roots zz; %{X=1 - 0.000000012*I,X=-1.0E+9 - 12.0*I} % 42) Multiple real and complex roots cause no difficulty, provided % that input is given in integer or rational form and mode is complex. zz :=(x**2-2i*x+5)**3*(x-2i)*(x-11/10)**2; roots zz; %{X=-1.44949*I, X=-1.44949*I, X=-1.44949*I, % X=3.44949*I, X=3.44949*I, X=3.44949*I, X=1.1, X=1.1, X=2.0*I} % 43) 2 real, 2 complex roots. zz:= (x**2-4)*(x**2+3i*x+5i); roots zz; %{X=2.0,X= - 2.0,X=-1.2714 + 0.466333*I,X=1.2714 - 3.46633*I} % 44) 4 complex roots. zz:= x**4+(0.000001i)*x-16; roots zz; %{X=2.0 - 0.0000000625*I,X=-2.0*I,X=-2.0 - 0.0000000625*I,X=2.0*I} % 45) 2 real, 2 complex roots. zz:= (x**2-4)*(x**2+2i*x+8); roots zz; %{X=2.0,X= - 2.0,X=2.0*I,X=-4.0*I} % 46) Using realroots to find only real roots. realroots zz; %{X= - 2.0,X=2.0} % 47) Same example, applying nearestroot to find a single root. zz:= (x**2-4)*(x**2+2i*x+8); nearestroot(zz,1); %{X=2.0} % 48) Same example, but focusing on imaginary point. nearestroot(zz,i); %{X=2.0*I} % 49) The seed parameter can be complex also. nearestroot(zz,1+i); %{X=2.0*I} % 50) One more nearestroot example. Nearest root to real point may be %complex. zz:= (x**2-4)*(x**2-i); roots zz; %{X=2.0,X= - 2.0,X=0.707107 + 0.707107*I,X=-0.707107 - 0.707107*I} nearestroot (zz,1); %{X=0.707107 + 0.707107*I} % 51) 1 real root plus 5 complex roots. zz:=(x**3-3i*x**2-5x+9)*(x**3-8); roots zz; %{X=2.0, X=-1 + 1.73205*I, X=-1 - 1.73205*I, % X=0.981383 - 0.646597*I, X=-2.41613 + 1.19385*I, % X=1.43475 + 2.45274*I} nearestroot(zz,1); %{X=0.981383 - 0.646597*I} % 52) roots can be computed to any accuracy desired, eg. (note that the % imaginary part of the second root is truncated because of its size, % and that the imaginary part of a complex root is never polished away, % even if it is smaller than the accuracy would require.) zz := x**3+10**(-20)*i*x**2+8; rootacc 12; roots zz; rootacc 0; %{X=1 - 1.73205080757*I,X=1 + 1.73205080757*I, % X=-2.0 - 3.33333E-21*I} % 53) Precision increase to 12 required to find small imaginary root, % but standard accuracy can be used. zz := x**2+123456789i*x+1; roots zz; %{X=0.0000000081*I,X=-123457000.0*I} % 54) Small real root is found with root 10*18 times larger(new). zz := (x+1)*(x**2+123456789*x+1); roots zz; %{X= - 0.0000000081,X= - 1.0,X= - 123457000.0} % 55) 2 complex, 3 real irrational roots. ss := (45*x**2+(-10i+12)*x-10i)*(x**3-5x**2+1); roots ss; %{X=0.469832,X= - 0.429174,X=4.95934,X=0.18139 + 0.417083*I, % X=-0.448056 - 0.19486*I} % 56) Complex polynomial with floating coefficients. zz := x**2+1.2i*x+2.3i+6.7; roots zz; %{X=-0.42732 + 2.09121*I,X=0.42732 - 3.29121*I} % multiple roots will be found if coefficients read in exactly. ZZ := X**3 + (1.09 - 2.4*I)*X**2 + (-1.44 - 2.616*I)*X + -1.5696; roots zz; %{X=1.2*I,X=1.2*I,X= - 1.09} % 57) Realroots, isolater and rlrootno accept 1, 2 or 3 arguments: (new) zz:= for j:=-1:3 product (x-j); rlrootno zz; % 5 realroots zz; %{X=0,X= -1,X=1,X=2.0,X=3.0} rlrootno(zz,positive); %positive selects positive, excluding 0. % 3 rlrootno(zz,negative); %negative selects negative, excluding 0. % 1 realroots(zz,positive); %{X=1,X=2.0,X=3.0} rlrootno(zz,-1.5,2); %the format with 3 arguments selects a range. % 4 realroots(zz,-1.5,2); %the range is inclusive, except that: %{X=0,X=-1,X=1,X=2.0} % A specific limit b may be excluded by using exclude b. Also, the % limits infinity and -infinity can be specified. realroots(zz,exclude 0,infinity); % equivalent to realroots(zz,positive). %{X=1,X=2.0,X=3.0} rlrootno(zz,-infinity,exclude 0); % equivalent to rlrootno(zz,negative). % 1 rlrootno(zz,-infinity,0); % 2 rlrootno(zz,infinity,-infinity); %equivalent to rlrootno zz; (order of limits does not matter.) % 5 realroots(zz,1,infinity); % finds all real roots >= 1. %{X=1,X=2.0,X=3.0} realroots(zz,1,positive); % finds all real roots > 1. %{X=2.0,X=3.0} % New capabilities added in mod 1.88. % 58) 3 complex roots, with two separated by very small real difference. zz :=(x+i)*(x+10**8i)*(x+10**8i+1); roots zz; %{X= - I,X=-1.0E+8*I,X=-1 - 1.0E+8*I} % 59) Real polynomial with two complex roots separated by very small % imaginary part. zz:= (10**14x+123456789000000+i)*(10**14x+123456789000000-i); roots zz; %{X=-1.23457 + 1.0E-14*I,X=-1.23457 - 1.0E-14*I} % 60) Real polynomial with two roots extremely close together. zz:= (x+2)*(10**10x+12345678901)*(10**10x+12345678900); roots zz; %{X= - 1.2345678901,X= - 1.23456789,X= - 2.0} % 61) Real polynomial with multiple root extremely close to simple root. zz:= (x-12345678/10000000)*(x-12345679/10000000)**2; roots zz; %{X=1.2345679,X=1.2345679,X=1.2345678} % 62) Similar problem using realroots. zz:=(x-2**30/10**8)**2*(x-(2**30+1)/10**8); realroots zz; %{X=10.73741824,X=10.73741824,X=10.73741825} % 63) Three complex roots with small real separation between two. zz:= (x-i)*(x-1-10**8i)*(x-2-10**8i); roots zz; %{X=I,X=1 + 1.0E+8*I,X=2.0 + 1.0E+8*I} % 64) Use of nearestroot to isolate one of the close roots. nearestroot(zz,10**8i+99/100); % {X=1 + 1.0E+8*I} % 65) Slightly more complicated example with close complex roots. zz:= (x-i)*(10**8x-1234-10**12i)*(10**8x-1233-10**12i); roots zz; %{X=I,X=0.00001233 + 10000.0*I,X=0.00001234 + 10000.0*I} % 66) Four closely spaced real roots with varying spacings. zz:= (x-1+1/10**7)*(x-1+1/10**8)*(x-1)*(x-1-1/10**7); roots zz; %{X=1,X=0.9999999,X=0.99999999,X=1.0000001} % 67) Complex pair plus two close real roots. zz:= (x**2+1)*(x-12345678/10000000)*(x-12345679/10000000); roots zz; %{X=1.2345678,X=1.2345679,X=I,X= - I} % 68) Same problem using realroots to find only real roots. realroots zz; %{X=1.2345678,X=1.2345679} % The switch ratroot causes output to be given in rational form. % 69) Two complex roots with output in rational form. on ratroot,complex; zz:=x**2-(5i+1)*x+1; sss:= roots zz; % 17343 - 93179*I 96531 + 518636*I %SSS := {X=-----------------,X=------------------} % 500000 100000 % With roots in rational form, mkpoly can be used to reconstruct a % polynomial. zz1 := mkpoly sss; % 2 %ZZ1 := 50000000000*X - (49999800000 + 250000100000*I)*X + ( % % 50000120977 + 42099*I) % Finding the roots of the new polynomial zz1. rr:= roots zz1; % 17343 - 93179*I 96531 + 518636*I %RR := {X=-----------------,X=------------------} % 500000 100000 % The roots are stable to the extent that rr=ss, although zz1 and % zz may differ. zz1 - zz; % 2 %49999999999*X - (49999799999 + 250000099995*I)*X + (50000120976 + % % 42099*I) % 70) Same type of problem in which roots are found exactly. zz:=(x-10**8+i)*(x-10**8-i)*(x-10**8+3i/2)*(x-i); rr := roots zz; % 4 3 2 %ZZ := (2*X - (600000000 - I)*X + 60000000000000005*X - ( % % 2000000000000000800000000 + 29999999999999999*I)*X + ( % % 30000000000000003 + 2000000000000000200000000*I))/2 %RR := {X=100000000 + I,X=100000000 - I,X=I,X= % % 200000000 - 3*I % -----------------} % 2 % Reconstructing a polynomial from the roots. ss := mkpoly rr; % 4 3 2 %SS := 2*X - (600000000 - I)*X + 60000000000000005*X - ( % % 2000000000000000800000000 + 29999999999999999*I)*X + ( % % 30000000000000003 + 2000000000000000200000000*I) % In this case, the same polynomial is obtained. ss - num zz; % 0 % 71) Finding one of the complex roots using nearestroot. nearestroot(zz,10**8-2i); % 200000000 - 3*I %{X=-----------------} % 2 % Finding the other complex root using nearestroot. nearestroot(zz,10**8+2i); %{X=100000000 + I} % 72) A realroots problem which requires accuracy increase to avoid % confusion of two roots. zz:=(x+1)*(10000000x-19999999)*(1000000x-2000001)*(x-2); realroots zz; % 19999999 2000001 % {X=-1,X=----------,X=2,X=---------} % 10000000 1000000 % 73) Without the accuracy increase, this example would produce the % obviously incorrect answer 2. realroots(zz,3/2,exclude 2); % 19999999 % {X=----------} % 10000000 % Rlrootno also gives the correct answer in this case. rlrootno(zz,3/2,exclude 2); % 1 % 74) Roots works equally well in this problem. rr := roots zz; % 19999999 2000001 % RR := {X=----------,X=-1,X=2,X=---------} % 10000000 1000000 % 75) The function getroot is convenient for obtaining the value of a % root. rr1 := getroot(1,rr); % 19999999 % RR1 := ---------- % 10000000 % 76) For example, the value can be used as an argument to nearestroot. nearestroot(zz,rr1); % 19999999 % {X=----------} % 10000000 comment New capabilities added to Mod 1.90 for avoiding floating point exceptions and exceeding iteration limits; % 77) This and the next example would previously have aborted because %of exceeding iteration limits: off ratroot; zz := x**16 - 900x**15 -2; roots zz; %{X= - 0.665423,X=0.069527 + 0.661817*I,X=0.069527 - 0.661817*I, % X=0.650944 + 0.138369*I,X=0.650944 - 0.138369*I, % X=-0.205664 + 0.632867*I,X=-0.205664 - 0.632867*I, % X=-0.607902 + 0.270641*I,X=-0.607902 - 0.270641*I, % X=0.332711 + 0.57633*I,X=0.332711 - 0.57633*I, % X=0.538375 + 0.391176*I,X=0.538375 - 0.391176*I, % X=-0.44528 + 0.494497*I,X=-0.44528 - 0.494497*I,X=900.0} % 78) a still harder example. % This example is commented out, since it takes significant time: % z := x**30 - 900x**29 - 2; roots zz; %{X= - 0.810021, % X=-0.04388 + 0.808856*I,X=-0.04388 - 0.808856*I, % X=0.805322 + 0.087587*I,X=0.805322 - 0.087587*I, % X=0.131027 + 0.799383*I,X=0.131027 - 0.799383*I, % X=-0.791085 + 0.174125*I,X=-0.791085 - 0.174125*I, % X=-0.216732 + 0.780507*I,X=-0.216732 - 0.780507*I, % X=0.767663 + 0.258664*I,X=0.767663 - 0.258664*I, % X=0.299811 + 0.752532*I,X=0.299811 - 0.752532*I, % X=-0.735162 + 0.340111*I,X=-0.735162 - 0.340111*I, % X=-0.379447 + 0.715665*I,X=-0.379447 - 0.715665*I, % X=0.694106 + 0.417645*I,X=0.694106 - 0.417645*I, % X=-0.524417 + 0.617362*I,X=-0.524417 - 0.617362*I, % X=0.454578 + 0.67049*I,X=0.454578 - 0.67049*I, % X=-0.644866 + 0.490195*I,X=-0.644866 - 0.490195*I, % X=0.588091 + 0.557094*I,X=0.588091 - 0.557094*I,X=900.0} % 79) this deceptively simple example previously caused floating point % overflows on some systems: aa := x**6 - 4*x**3 + 2; realroots aa; %{X=0.836719,X=1.50579} % 80) a harder problem, which would have failed on almost all systems: rr := x**16 - 90000x**15 - x**2 -2; realroots rr; %{X= - 0.493299,X=90000.0} % 81) this example would have failed because of floating point % exceptions on almost all computer systems. rr := X**30 - 9*10**10*X**29 - 2; realroots rr; %{X= - 0.429188,X=9.0E+10} % 82) a test of allroot on this example. % This example is commented out, since it takes significant time: % roots rr; %{X= - 0.429188, % X=-0.023236 + 0.428559*I,X=-0.023236 - 0.428559*I, % X=0.426672 + 0.046403*I,X=0.426672 - 0.046403*I, % X=0.069435 + 0.423534*I,X=0.069435 - 0.423534*I, % X=-0.419154 + 0.092263*I,X=-0.419154 - 0.092263*I, % X=-0.11482 + 0.413544*I,X=-0.11482 - 0.413544*I, % X=0.406722 + 0.13704*I,X=0.406722 - 0.13704*I, % X=0.158859 + 0.398706*I,X=0.158859 - 0.398706*I, % X=-0.389521 + 0.180211*I,X=-0.389521 - 0.180211*I, % X=-0.201035 + 0.379193*I,X=-0.201035 - 0.379193*I, % X=0.367753 + 0.22127*I,X=0.367753 - 0.22127*I, % X=-0.277851 + 0.327111*I,X=-0.277851 - 0.327111*I, % X=0.240855 + 0.355234*I,X=0.240855 - 0.355234*I, % X=-0.341674 + 0.259734*I,X=-0.341674 - 0.259734*I, % X=0.311589 + 0.295153*I,X=0.311589 - 0.295153*I,X=9.0E+10} % 83) test of starting point for iteration: no convergence if good % real starting point is not found. zz := x**30 -9*10**12x**29 -2; firstroot zz; %{X= - 0.36617} % 84) a case in which there are no real roots and good imaginary % starting point must be used or roots cannot be found. zz:= 9x**16 - x**5 +1; roots zz; %{X=0.182294 + 0.828368*I,X=0.182294 - 0.828368*I, % X=-0.697397 + 0.473355*I,X=-0.697397 - 0.473355*I, % X=0.845617 + 0.142879*I,X=0.845617 - 0.142879*I, % X=-0.161318 + 0.87905*I,X=-0.161318 - 0.87905*I, % X=-0.866594 + 0.193562*I,X=-0.866594 - 0.193562*I, % X=0.459373 + 0.737443*I,X=0.459373 - 0.737443*I, % X=0.748039 + 0.494348*I,X=0.748039 - 0.494348*I, % X=-0.510014 + 0.716449*I,X=-0.510014 - 0.716449*I} % 85) five complex roots. zz := x**5 - x**3 + i; roots zz; %{X=-0.83762*I,X=-0.664702 + 0.636663*I,X=0.664702 + 0.636663*I, % X=1.16695 - 0.21785*I,X=-1.16695 - 0.21785*I} % Additional capabilities in Mod 1.91. % 86) handling of polynomial with huge or infinitesimal coefficients. on rounded; zz:= 1.0e-500x**3+x**2+x; roots zz; off rounded; %{X=0,X=-1,X= - 1.0E+500} % switch roundbf will have been turned on in the last example in most % computer systems. This will inhibit the use of hardware floating % point unless roundbf is turned off. % polynomials which make use of powergcd substitution and cascaded % solutions. % Uncomplicated cases. switch powergcd; % introduced here to verify that same answers are % obtained with and without employing powergcd strategy. Roots are % found faster for applicable cases when !*powergcd=t (default state.) % 87) powergcd done at the top level. zz := x**12-5x**9+1; roots zz; %{X=0.848444,X=-0.424222 + 0.734774*I,X=-0.424222 - 0.734774*I, % X=0.152522 - 0.816316*I, % X=-0.783212 + 0.276071*I,X=0.63069 + 0.540246*I, % X=0.152522 + 0.816316*I, % X=-0.783212 - 0.276071*I,X=0.63069 - 0.540246*I, % X=1.70906,X=-0.85453 + 1.48009*I,X=-0.85453 - 1.48009*I} off powergcd; roots zz; on powergcd; %{X=0.848444,X=-0.424222 + 0.734774*I,X=-0.424222 - 0.734774*I, % X=1.70906,X=-0.783212 + 0.276071*I,X=-0.783212 - 0.276071*I, % X=0.152522 + 0.816316*I,X=0.152522 - 0.816316*I, % X=0.63069 + 0.540246*I,X=0.63069 - 0.540246*I, % X=-0.85453 + 1.48009*I,X=-0.85453 - 1.48009*I} % 88) powergcd done after square free factoring. zz := (x-1)**2*zz; roots zz; %{X=1,X=1, % X=0.848444,X=-0.424222 + 0.734774*I,X=-0.424222 - 0.734774*I, % X=0.152522 - 0.816316*I,X=-0.783212 + 0.276071*I, % X=0.63069 + 0.540246*I, % X=0.152522 + 0.816316*I,X=-0.783212 - 0.276071*I, % X=0.63069 - 0.540246*I, % X=1.70906,X=-0.85453 + 1.48009*I,X=-0.85453 - 1.48009*I} off powergcd; roots zz; on powergcd; %{X=1,X=1, % X=0.848444,X=-0.424222 + 0.734774*I,X=-0.424222 - 0.734774*I, % X=1.70906,X=-0.783212 + 0.276071*I,X=-0.783212 - 0.276071*I, % X=0.152522 + 0.816316*I,X=0.152522 - 0.816316*I, % X=0.63069 + 0.540246*I,X=0.63069 - 0.540246*I, % X=-0.85453 + 1.48009*I,X=-0.85453 - 1.48009*I} % 89) powergcd done after separation into real and complex polynomial. zz := x**5-i*x**4+x**3-i*x**2+x-i; roots zz; %{X=0.5 + 0.866025*I,X=-0.5 - 0.866025*I, % X=-0.5 + 0.866025*I,X=0.5 - 0.866025*I,X=I} off powergcd; roots zz; on powergcd; %{X=-0.5 + 0.866025*I,X=-0.5 - 0.866025*I, % X=0.5 + 0.866025*I,X=0.5 - 0.866025*I,X=I} % Cases where root separation requires accuracy and/or precision % increase. In some examples we get excess accuracy, but it is hard % avoid this and still get all roots separated. % 90) accuracy increase required to separate close roots; let x=y**2; zz:= (x-3)*(100000000x-300000001); roots zz; %{Y=1.732050808,Y= - 1.732050808,Y=1.73205081,Y= - 1.73205081} off powergcd; roots zz; on powergcd; %{Y=1.732050808,Y= - 1.732050808,Y= - 1.73205081,Y=1.73205081} % 91) roots to be separated are on different square free factors. zz:= (x-3)**2*(10000000x-30000001); roots zz; %{Y=1.73205081 ,Y=1.73205081 ,Y= - 1.73205081 ,Y= - 1.73205081 , % Y=1.73205084 ,Y= - 1.73205084} off powergcd; roots zz; on powergcd; %{Y=1.73205081 ,Y=1.73205081 ,Y= - 1.73205081,Y= - 1.73205081, % Y=1.73205084 ,Y= - 1.73205084} % 92) roots must be separated in the complex polynomial factor only. zz :=(y+1)*(x+10**8i)*(x+10**8i+1); roots zz; %{Y=-1, % Y=-7071.067812 + 7071.067812*I,Y=7071.067812 - 7071.067812*I, % Y=-7071.067777 + 7071.067847*I,Y=7071.067777 - 7071.067847*I} % 93) zz := (x-2)**2*(1000000x-2000001)*(y-1); roots zz; %{Y=1.4142136,Y=1.4142136,Y= - 1.4142136,Y= - 1.4142136, % Y=1,Y=1.4142139,Y= - 1.4142139} % 94) zz := (x-2)*(10000000x-20000001); roots zz; %{Y=1.41421356 ,Y= - 1.41421356 ,Y=1.4142136,Y= - 1.4142136} % 95) zz := (x-3)*(10000000x-30000001); roots zz; %{Y=1.73205081 ,Y= - 1.73205081 ,Y=1.73205084 ,Y= - 1.73205084} % 96) zz := (x-9)**2*(1000000x-9000001); roots zz; %{Y=3.0,Y=3.0,Y= - 3.0,Y= - 3.0,Y=3.00000017,Y= - 3.00000017} % 97) zz := (x-3)**2*(1000000x-3000001); roots zz; %{Y=1.7320508,Y=1.7320508,Y= - 1.7320508,Y= - 1.7320508, % Y=1.7320511,Y= - 1.7320511} % 98) the accuracy of the root -sqrt 3 depends upon another close root. on rounded; zz := (x-3)*(y+1.732058); roots zz; %{Y= - 1.732051,Y=1.73205,Y= - 1.732058} zz := (x-3)*(y+1.732051); roots zz; %{Y= - 1.73205081,Y=1.73205,Y= - 1.732051} zz := (x-3)*(y+1.73205081); roots zz; %{Y= - 1.732050808,Y=1.73205,Y= - 1.73205081} % 99) minimum accuracy specified using rootacc. rootacc 10; roots zz; %{Y= - 1.732050808,Y=1.732050808,Y= - 1.73205081} % 100) off rounded; rootacc 6; zz := (y-i)*(x-2)*(1000000x-2000001); roots zz; %{Y=1.4142136,Y= - 1.4142136,Y=1.4142139,Y= - 1.4142139,Y=I} % 101) this example requires accuracy 15. zz:= (y-2)*(100000000000000y-200000000000001); roots zz; %{Y=2.0,Y=2.00000000000001} % 102) still higher precision needed. zz:= (y-2)*(10000000000000000000y-20000000000000000001); roots zz; %{Y=2.000 00000 00000 00000 1,Y=2.0} % 103) increase in precision required for substituted polynomial. zz:= (x-2)*(10000000000x-20000000001); roots zz; %{Y=1.41421356241,Y= - 1.41421356241,Y=1.41421356237, % Y= - 1.41421356237} % 104) still higher precision required for substituted polynomial. zz:= (x-2)*(100000000000000x-200000000000001); roots zz; %{Y=1.414 21356 23730 99,Y= - 1.414 21356 23730 99, % Y=1.414 21356 23730 95,Y= - 1.414 21356 23730 95} % 105) accuracy must be increased to separate root of complex factor % from root of real factor. zz:=(9y-10)*(y-2)*(9y-10-9i/100000000); roots zz; %{Y=1.111111111,Y=2.0,Y=1.111111111 + 0.00000001*I} % 106) realroots does the same accuracy increase for real root based % upon the presence of a close complex root in the same polynomial. % The reason for this might not be obvious unless roots is called. realroots zz; %{Y=1.111111111,Y=2.0} % 107) realroots now uses powergcd logic whenever it is applicable. zz := (x-1)*(x-2)*(x-3); realroots zz; %{Y=-1,Y=1,Y= - 1.41421,Y=1.41421,Y= - 1.73205,Y=1.73205} realroots(zz,exclude 1,2); %{Y=1.41421,Y=1.73205} % 108) root of degree 1 polynomial factor must be evaluated at % precision 18 and accuracy 10 in order to separate it from a root of % another real factor. clear x; zz:=(9x-10)**2*(9x-10-9/100000000)*(x-2); roots zz; %{X=1.111111111,X=1.111111111,X=1.111111121,X=2.0} nearestroot(zz,1); %{X=1.111111111} nearestroot(zz,1.5); %{X=1.111111121} nearestroot(zz,1.65); %{X=2.0} % 109) in this example, precision >=40 is used and two roots need to be % found to accuracy 16 and two to accuracy 14. zz := (9x-10)*(7x-8)*(9x-10-9/10**12)*(7x-8-7/10**14); roots zz; %{X=1.1111111111121,X=1.142 85714 28571 43, % X=1.1111111111111,X=1.142 85714 28571 53} % 110) very small real or imaginary parts of roots require high % precision or exact computations, or they will be lost or incorrectly % found. zz := 1000000*R**18 + 250000000000*R**4 - 1000000*R**2 + 1; roots zz; %{R=0.00141421 + 1.6E-26*I,R=-0.00141421 - 1.6E-26*I, % R=0.00141421 - 1.6E-26*I,R=-0.00141421 + 1.6E-26*I, % R=2.36886 + 0.54068*I,R=-2.36886 - 0.54068*I, % R=-2.36886 + 0.54068*I,R=2.36886 - 0.54068*I, % R=1.05424 + 2.18916*I,R=-1.05424 - 2.18916*I, % R=-1.05424 + 2.18916*I,R=1.05424 - 2.18916*I, % R=2.42978*I,R=-2.42978*I, % R=1.89968 + 1.51494*I,R=-1.89968 - 1.51494*I, % R=-1.89968 + 1.51494*I,R=1.89968 - 1.51494*I} showtime; end;