Sat Jun 29 13:48:36 PDT 1991
REDUCE 3.4, 15-Jul-91 ...
1: 1:
2: 2:
3: 3: % 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;
{X=0}
% To load roots package.
write "This is Roots Package test ", symbolic roots!-mod$
This is Roots Package test Mod 1.92, 10 Oct 1990.
% Simple root finding.
showtime;
Time: 476 ms
% 1) multiple real and imaginary roots plus two real roots.
zz:= (x-3)**2*(100x**2+113)**2*(1000000x-10000111)*(x-1);
8 7 6
ZZ := 10000000000*X - 170001110000*X + 872607770000*X
5 4 3
- 1974219158600*X + 2833796550200*X - 3810512046359*X
2
+ 3119397498913*X - 2030292260385*X + 1149222756231
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}
%{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);
4 3 2
ZZ := 1000000*X - 4000001*X + 5000002*X - 4000001*X + 4000002
roots zz;
{X=2.0,X=I,X= - I,X=2.000001}
%{X=2.0,X=I,X= - I,X=2.000001}
% 3) accuracy increase to 8 required.
zz:= (x-3)*(10000000x-30000001);
2
ZZ := 10000000*X - 60000001*X + 90000003
roots zz;
{X=3.0,X=3.0000001}
%{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);
4 3 2
ZZ := 2*(500000000000*X - 6120617500000*X + 28085557620617*X
- 57256673223702*X + 43756673585553)
roots zz;
{X=3.0,X=3.0,X=3.000001,X=3.24123}
%{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);
5 4 3 2
ZZ := 10*X - 1020031*X + 2013152032*X - 1005224243011*X
+ 2104312111000*X - 1101100000000
roots zz;
{X=1,X=1.1,X=1000.0,X=1001.0,X=100000.0}
%{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);
5 4 3 2
ZZ := 100000000*X - 500000000*X + 1000000000*X - 1000000000*X
+ 499999999*X - 99999999
roots zz;
{X=1,
X=1.01,
X=0.99,
X=1 + 0.01*I,
X=1 - 0.01*I}
%{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);
4 3 2
ZZ := 2*(500000*X - 2030000*X + 3090550*X - 2091103*X + 530553)
roots zz;
{X=1.02,X=1.01,X=1,X=1.03}
%{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);
5 4 3 2
ZZ := 2*(500000*X - 2515000*X + 5065100*X - 5105450*X + 2575601*X
- 520251)
roots zz;
{X=1.01,
X=1,
X=1.02,
X=1 + 0.1*I,
X=1 - 0.1*I}
%{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);
5 4 3 2
ZZ := 100000000*X - 1000000000*X + 4000000000*X - 8000000000*X
+ 7999999999*X - 3199999998
roots zz;
{X=2.0,
X=2.01,
X=1.99,
X=2.0 + 0.01*I,
X=2.0 - 0.01*I}
%{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);
3 2
SS := 100000000*X - 600000000*X + 1200000001*X - 800000002
roots ss;
{X=2.0,X=2.0 + 0.0001*I,X=2.0 - 0.0001*I}
%{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);
7 6 5 4 3 2
ZZ := X*(X - 14*X + 74*X - 186*X + 249*X - 236*X + 176*X - 64)
roots zz;
{X=0,
X=4.0,
X=4.0,
X=4.0,
X=1,
X=1,
X=I,
X= - I}
%{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}
%{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;
9
ZZ := X - 45*X - 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}
%{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;
9 2
100*X - 999900*X - 1
ZZ := ------------------------
100
roots zz;
*** precision increased to 27
{X=5.0E-29 + 0.00100005*I,
X=5.0E-29 - 0.00100005*I,
X=3.72754,
X=-0.82946 + 3.63408*I,
X=-0.82946 - 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}
%{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;
7
zz:= x**9-500x**2-0.001;
9 2
1000*X - 500000*X - 1
ZZ := -------------------------
1000
roots zz;
{X=1.6E-26 + 0.001414214*I,
X=1.6E-26 - 0.001414214*I,
X=2.429781,
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}
%{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;
6 5 4 3 2
ZZ := 10000*X + 210027*X + 1750000*X + 7350000*X + 16240000*X
+ 17640000*X + 7200000
roots zz;
{X= - 2.950367,
X=-4.452438 + 0.021235*I,
X=-4.452438 - 0.021235*I,
X= - 2.003647,
X= - 0.9999775,
X= - 6.143833}
%{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;
8 7 6 5 4
ZZ := 2*(500*X + 18001*X + 273000*X + 2268000*X + 11224500*X
3 2
+ 33642000*X + 59062000*X + 54792000*X + 20160000)
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}
%{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;
10 9 8 7 6
ZZ := 1000*X + 55001*X + 1320000*X + 18150000*X + 157773000*X
5 4 3
+ 902055000*X + 3416930000*X + 8409500000*X
2
+ 12753576000*X + 10628640000*X + 3628800000
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}
%{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;
12 11 10 9
ZZ := 4*(2500*X + 195001*X + 6792500*X + 139425000*X
8 7 6
+ 1873657500*X + 17316585000*X + 112475577500*X
5 4 3
+ 515175375000*X + 1643017090000*X + 3535037220000*X
2
+ 4828898880000*X + 3716107200000*X + 1197504000000)
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}
%{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;
6
% 21) Finding one of a cluster of 8 roots.
zz:= (10**16*(x-1)**8-1);
8 7
ZZ := 10000000000000000*X - 80000000000000000*X
6 5
+ 280000000000000000*X - 560000000000000000*X
4 3
+ 700000000000000000*X - 560000000000000000*X
2
+ 280000000000000000*X - 80000000000000000*X
+ 9999999999999999
nearestroot(zz,2);
{X=1.01}
%{X=1.01}
% 22) Six real roots spaced 0.01 apart.
c := 100;
C := 100
zz:= (x-1)*for i:=1:5 product (c*x-(c+i));
6 5 4
ZZ := 40*(250000000*X - 1537500000*X + 3939625000*X
3 2
- 5383556250*X + 4137919435*X - 1696170123*X + 289681938
)
roots zz;
{X=1.03,X=1.02,X=1.04,X=1.01,X=1,X=1.05}
%{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;
C := 1000
zz:= (x-1)*for i:=1:5 product (c*x-(c+i));
6 5 4
ZZ := 40*(25000000000000*X - 150375000000000*X + 376877125000000*X
3 2
- 503758505625000*X + 378762766881850*X
- 151883516888703*X + 25377130631853)
roots zz;
{X=1.003,X=1.002,X=1.004,X=1.001,X=1,X=1.005}
%{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;
C := 10000
zz:= (x-1)*for i:=1:4 product (c*x-(c+i));
5 4
ZZ := 8*(1250000000000000*X - 6251250000000000*X
3 2
+ 12505000437500000*X - 12507501312562500*X
+ 6255001312625003*X - 1251250437562503)
roots zz;
{X=1.0002,X=1.0003,X=1.0001,X=1,X=1.0004}
%{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);
9 8 7 6
ZZ := 100000000*X - 900000000*X + 3600000000*X - 8400000000*X
5 4 3
+ 12600000000*X - 12600000000*X + 8400000000*X
2
- 3600000000*X + 899999999*X - 99999999
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}
%{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);
9 8
ZZ := 10000000000000000*X - 90000000000000000*X
7 6
+ 360000000000000000*X - 840000000000000000*X
5 4
+ 1260000000000000000*X - 1260000000000000000*X
3 2
+ 840000000000000000*X - 360000000000000000*X
+ 89999999999999999*X - 9999999999999999
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}
%{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);
9 8
ZZ := 1000000000000000000000000*X - 9000000000000000000000000*X
7
+ 36000000000000000000000000*X
6
- 84000000000000000000000000*X
5
+ 126000000000000000000000000*X
4
- 126000000000000000000000000*X
3
+ 84000000000000000000000000*X
2
- 36000000000000000000000000*X + 8999999999999999999999999*X
- 999999999999999999999999
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}
%{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);
8 6
ZZ := 10000000000000000*X - 40000000000000000*X
4 2
+ 59999999800000000*X - 40000001200000000*X
+ 9999999800000001
roots zz;
{X=1.01,
X= - 1.01,
X=0.99,
X= - 0.99,
X=1 + 0.01*I,
X=-1 - 0.01*I,
X=-1 + 0.01*I,
X=1 - 0.01*I}
%{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);
8 7
ZZ := 1000000000000000000000000*X - 8080000000000000000000000*X
6
+ 28562400000000000000000000*X
5
- 57694432000000000000000000*X
4
+ 72836160149999999900000000*X
3
- 58848320599199999600000000*X
2
+ 29716320897575999400000000*X - 8574560597551679600000000*X
+ 1082432149175678300000001
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}
%{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);
8 6 4 2
ZZ := 4*(25000000*X - 1000000*X + 20000*X + 200*X + 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}
%{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);
ZZ :=
4 3 2
1000000*X - 13000001*X + 63000010*X - 135000033*X + 108000036
roots zz;
{X=3.0,X=3.0,X=3.000001,X=4.0}
%{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}
%{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;
ZZ := - I + X
roots zz;
{X=I}
%{X=I}
% 34) Real case.
zz:= y-7;
ZZ := Y - 7
roots zz;
{Y=7.0}
%{Y=7.0}
% 35) Roots with small imaginary parts (new capability);
zz := 10**16*(x**2-2x+1)+1;
2
ZZ := 10000000000000000*X - 20000000000000000*X + 10000000000000001
roots zz;
{X=1 + 0.00000001*I,X=1 - 0.00000001*I}
%{X=1 + 0.00000001*I,X=1 - 0.00000001*I}
% 36) One real, one complex root.
zz:=(x-9)*(x-5i-7);
2
ZZ := - 5*I*X + 45*I + X - 16*X + 63
roots zz;
{X=9.0,X=7.0 + 5.0*I}
%{X=9.0,X=7.0 + 5.0*I}
% 37) Three real roots.
zz:= (x-1)*(x-2)*(x-3);
3 2
ZZ := X - 6*X + 11*X - 6
roots zz;
{X=2.0,X=1,X=3.0}
%{X=2.0,X=1,X=3.0}
% 38) 2 real + 1 imaginary root.
zz:=(x**2-8)*(x-5i);
2 3
ZZ := - 5*I*X + 40*I + X - 8*X
roots zz;
{X=2.82843,X= - 2.82843,X=5.0*I}
%{X=2.82843,X= - 2.82843,X=5.0*I}
% 39) 2 complex roots.
zz:= (x-1-2i)*(x+2+3i);
2
ZZ := I*X - 7*I + X + X + 4
roots zz;
{X=1 + 2.0*I,X=-2.0 - 3.0*I}
%{X=1 + 2.0*I,X=-2.0 - 3.0*I}
% 40) 2 irrational complex roots.
zz:= x**2+(3+2i)*x+7i;
2
ZZ := 2*I*X + 7*I + X + 3*X
roots zz;
{X=0.14936 - 2.21259*I,X=-3.14936 + 0.21259*I}
%{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;
2
ZZ := 12*I*X + X + 1000000000*X - 1000000000
roots zz;
{X=1 - 0.000000012*I,X=-1.0E+9 - 12.0*I}
%{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;
8 7 6 5 4
ZZ := ( - 800*I*X + 1760*I*X - 6768*I*X + 12760*I*X - 25018*I*X
3 2 9
+ 39600*I*X - 46780*I*X + 55000*I*X - 30250*I + 100*X
8 7 6 5 4 3
- 220*X - 779*X + 1980*X - 9989*X + 19580*X - 28269*X
2
+ 38500*X - 21175*X)/100
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}
%{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);
3 2 4 2
ZZ := 3*I*X + 5*I*X - 12*I*X - 20*I + X - 4*X
roots zz;
{X=2.0,
X= - 2.0,
X=-1.2714 + 0.46633*I,
X=1.2714 - 3.46633*I}
%{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;
4
I*X + 1000000*X - 16000000
ZZ := -----------------------------
1000000
roots zz;
{X=2.0 - 0.0000000625*I,
X=-2.0 - 0.0000000625*I,
X=-2.0*I,
X=2.0*I}
%{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);
3 4 2
ZZ := 2*I*X - 8*I*X + X + 4*X - 32
roots zz;
{X=2.0,X= - 2.0,X=2.0*I,X=-4.0*I}
%{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}
%{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);
3 4 2
ZZ := 2*I*X - 8*I*X + X + 4*X - 32
nearestroot(zz,1);
{X=2.0}
%{X=2.0}
% 48) Same example, but focusing on imaginary point.
nearestroot(zz,i);
{X=2.0*I}
%{X=2.0*I}
% 49) The seed parameter can be complex also.
nearestroot(zz,1+i);
{X=2.0*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);
2 4 2
ZZ := - I*X + 4*I + X - 4*X
roots zz;
{X=2.0,
X= - 2.0,
X=0.707107 + 0.707107*I,
X=-0.707107 - 0.707107*I}
%{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}
%{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);
5 2 6 4 3
ZZ := - 3*I*X + 24*I*X + X - 5*X + X + 40*X - 72
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}
%{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}
%{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;
1 2 3
ZZ := -----------------------*I*X + X + 8
100000000000000000000
rootacc 12;
12
roots zz;
{X=1 + 1.73205080757*I,X=1 - 1.73205080757*I,X=-2.0 - 3.33333E-21*I}
rootacc 0;
6
%{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;
2
ZZ := 123456789*I*X + X + 1
roots zz;
{X=0.0000000081*I,X=-123457000.0*I}
%{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);
3 2
ZZ := X + 123456790*X + 123456790*X + 1
roots zz;
{X= - 0.0000000081,X=-1,X= - 123457000.0}
%{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);
4 3 2 5 4
SS := - 10*I*X + 40*I*X + 50*I*X - 10*I*X - 10*I + 45*X - 213*X
3 2
- 60*X + 45*X + 12*X
roots ss;
{X=0.469832,
X= - 0.429174,
X=4.95934,
X=0.18139 + 0.417083*I,
X=-0.448056 - 0.19486*I}
%{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;
2
12*I*X + 23*I + 10*X + 67
ZZ := ----------------------------
10
roots zz;
{X=-0.42732 + 2.09121*I,X=0.42732 - 3.29121*I}
%{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;
2 3 2
- 6000*I*X - 6540*I*X + 2500*X + 2725*X - 3600*X - 3924
ZZ := -------------------------------------------------------------
2500
roots zz;
{X=1.2*I,X=1.2*I,X= - 1.09}
%{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);
4 3 2
ZZ := X*(X - 5*X + 5*X + 5*X - 6)
rlrootno zz;
5
% 5
realroots zz;
{X=0,
X= - 1,
X=1,
X=2.0,
X=3.0}
%{X=0,X= -1,X=1,X=2.0,X=3.0}
rlrootno(zz,positive);
3
%positive selects positive, excluding 0.
% 3
rlrootno(zz,negative);
1
%negative selects negative, excluding 0.
% 1
realroots(zz,positive);
{X=1,X=2.0,X=3.0}
%{X=1,X=2.0,X=3.0}
rlrootno(zz,-1.5,2);
4
%the format with 3 arguments selects a range.
% 4
realroots(zz,-1.5,2);
{X=0,X= - 1,X=1,X=2.0}
%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);
{X=1,X=2.0,X=3.0}
% equivalent to realroots(zz,positive).
%{X=1,X=2.0,X=3.0}
rlrootno(zz,-infinity,exclude 0);
1
% equivalent to rlrootno(zz,negative).
% 1
rlrootno(zz,-infinity,0);
2
% 2
rlrootno(zz,infinity,-infinity);
5
%equivalent to rlrootno zz; (order of limits does not matter.)
% 5
realroots(zz,1,infinity);
{X=1,X=2.0,X=3.0}
% finds all real roots >= 1.
%{X=1,X=2.0,X=3.0}
realroots(zz,1,positive);
{X=2.0,X=3.0}
% 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);
2 3 2
ZZ := 200000001*I*X + 100000001*I*X - 10000000000000000*I + X + X
- 10000000200000000*X - 100000000
roots zz;
{X= - I,X=-1.0E+8*I,X=-1 - 1.0E+8*I}
%{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);
2
ZZ := 10000000000000000000000000000*X
+ 24691357800000000000000000000*X
+ 15241578750190521000000000001
roots zz;
{X=-1.23457 + 1.0E-14*I,X=-1.23457 - 1.0E-14*I}
%{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);
3 2
ZZ := 100*(1000000000000000000*X + 4469135780100000000*X
+ 6462429435342508889*X + 3048315750285017778)
roots zz;
{X= - 1.2345678901,X= - 1.23456789,X= - 2.0}
%{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;
3 2
ZZ := (500000000000000000000*X - 1851851800000000000000*X
+ 2286236726108825000000*X - 940838132549050755399)/
500000000000000000000
roots zz;
{X=1.2345679,X=1.2345679,X=1.2345678}
%{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);
3 2
ZZ := (610351562500000000*X - 19660800006103515625*X
+ 211106232664064000000*X - 755578637962830675968)/
610351562500000000
realroots zz;
{X=10.73741824,X=10.73741824,X=10.73741825}
%{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);
2 3
ZZ := - 200000001*I*X + 300000003*I*X + 9999999999999998*I + X
2
- 3*X - 10000000199999998*X + 300000000
roots zz;
{X=I,X=1 + 1.0E+8*I,X=2.0 + 1.0E+8*I}
%{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}
% {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);
2
ZZ := 2*( - 100005000000000000000*I*X + 1233623350000000*I*X
3
+ 499999999999999999239239*I + 5000000000000000*X
2
- 123350000000*X - 500099999999999999239239*X
+ 1233500000000000)
roots zz;
{X=I,X=0.00001233 + 10000.0*I,X=0.00001234 + 10000.0*I}
%{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);
4 3
ZZ := (10000000000000000000000*X - 39999999900000000000000*X
2
+ 59999999699999900000000*X - 39999999699999800000001*X
+ 9999999899999900000001)/10000000000000000000000
roots zz;
{X=1,X=0.99999999,X=0.9999999,X=1.0000001}
%{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);
4 3 2
ZZ := (50000000000000*X - 123456785000000*X + 126207888812681*X
- 123456785000000*X + 76207888812681)/50000000000000
roots zz;
{X=1.2345678,X=1.2345679,X=I,X= - I}
%{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}
%{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;
2
ZZ := X - (1 + 5*I)*X + 1
sss:= roots zz;
17343 - 93179*I 96531 + 518636*I
SSS := {X=-----------------,X=------------------}
500000 100000
% 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)
% 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
% 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)
% 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);
4 3 2
ZZ := (2*X - (600000000 - I)*X + 60000000000000005*X
- (2000000000000000800000000 + 29999999999999999*I)*X
+ (30000000000000003 + 2000000000000000200000000*I))/2
rr := roots zz;
200000000 - 3*I
RR := {X=100000000 + I,X=100000000 - I,X=I,X=-----------------}
2
% 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)
% 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
% 0
% 71) Finding one of the complex roots using nearestroot.
nearestroot(zz,10**8-2i);
200000000 - 3*I
{X=-----------------}
2
% 200000000 - 3*I
%{X=-----------------}
% 2
% Finding the other complex root using nearestroot.
nearestroot(zz,10**8+2i);
{X=100000000 + I}
%{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);
4 3 2
ZZ := 10000000000000*X - 50000009000000*X + 60000026999999*X
+ 40000000000001*X - 80000035999998
realroots zz;
19999999 2000001
{X=-1,X=----------,X=2,X=---------}
10000000 1000000
% 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
% 19999999
% {X=----------}
% 10000000
% Rlrootno also gives the correct answer in this case.
rlrootno(zz,3/2,exclude 2);
1
% 1
% 74) Roots works equally well in this problem.
rr := roots zz;
19999999 2000001
RR := {X=----------,X=2,X=-1,X=---------}
10000000 1000000
% 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
% 19999999
% RR1 := ----------
% 10000000
% 76) For example, the value can be used as an argument to nearestroot.
nearestroot(zz,rr1);
19999999
{X=----------}
10000000
% 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;
16 15
ZZ := X - 900*X - 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}
%{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;
6 3
AA := X - 4*X + 2
realroots aa;
{X=0.836719,X=1.50579}
%{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;
16 15 2
RR := X - 90000*X - X - 2
realroots rr;
{X= - 0.493299,X=90000.0}
%{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;
30 29
RR := X - 90000000000*X - 2
realroots rr;
{X= - 0.429188,X=9.0E+10}
%{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;
30 29
ZZ := X - 9000000000000*X - 2
firstroot zz;
{X= - 0.36617}
%{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;
16 5
ZZ := 9*X - X + 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}
%{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;
5 3
ZZ := X - X + 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}
%{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;
*** Domain mode COMPLEX changed to COMPLEX-ROUNDED
zz:= 1.0e-500x**3+x**2+x;
2
ZZ := X*(1.0E-500*X + X + 1)
roots zz;
{X=0,X=-1,X= - 1.0E+500}
off rounded;
*** Domain mode COMPLEX-ROUNDED changed to COMPLEX
%{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;
12 9
ZZ := X - 5*X + 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}
%{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;
{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}
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;
14 13 12 11 10 9 2
ZZ := X - 2*X + X - 5*X + 10*X - 5*X + X - 2*X + 1
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}
%{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;
{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}
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;
5 4 3 2
ZZ := X - I*X + X - I*X + 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}
%{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;
{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}
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);
4 2
ZZ := 100000000*Y - 600000001*Y + 900000003
roots zz;
*** Roots precision increase to 22
{Y=1.732050808,Y= - 1.732050808,Y=1.73205081,Y= - 1.73205081}
%{Y=1.732050808,Y= - 1.732050808,Y=1.73205081,Y= - 1.73205081}
off powergcd;
roots zz;
{Y=1.732050808,Y= - 1.732050808,Y= - 1.73205081,Y=1.73205081}
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);
6 4 2
ZZ := 10000000*Y - 90000001*Y + 270000006*Y - 270000009
roots zz;
{Y=1.73205081,
Y=1.73205081,
Y= - 1.73205081,
Y= - 1.73205081,
Y=1.73205084,
Y= - 1.73205084}
%{Y=1.73205081 ,Y=1.73205081 ,Y= - 1.73205081 ,Y= - 1.73205081 ,
% Y=1.73205084 ,Y= - 1.73205084}
off powergcd;
roots zz;
{Y=1.73205081,
Y=1.73205081,
Y= - 1.73205081,
Y= - 1.73205081,
Y=1.73205084,
Y= - 1.73205084}
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);
5 4 3 2
ZZ := Y + Y + (1 + 200000000*I)*Y + (1 + 200000000*I)*Y
- (10000000000000000 - 100000000*I)*Y
- (10000000000000000 - 100000000*I)
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}
%{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);
7 6 5 4 3
ZZ := 1000000*Y - 1000000*Y - 6000001*Y + 6000001*Y + 12000004*Y
2
- 12000004*Y - 8000004*Y + 8000004
roots zz;
{Y=1.4142136,
Y=1.4142136,
Y= - 1.4142136,
Y= - 1.4142136,
Y=1,
Y=1.4142139,
Y= - 1.4142139}
%{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);
4 2
ZZ := 10000000*Y - 40000001*Y + 40000002
roots zz;
{Y=1.41421356,Y= - 1.41421356,Y=1.4142136,Y= - 1.4142136}
%{Y=1.41421356 ,Y= - 1.41421356 ,Y=1.4142136,Y= - 1.4142136}
% 95)
zz := (x-3)*(10000000x-30000001);
4 2
ZZ := 10000000*Y - 60000001*Y + 90000003
roots zz;
{Y=1.73205081,Y= - 1.73205081,Y=1.73205084,Y= - 1.73205084}
%{Y=1.73205081 ,Y= - 1.73205081 ,Y=1.73205084 ,Y= - 1.73205084}
% 96)
zz := (x-9)**2*(1000000x-9000001);
6 4 2
ZZ := 1000000*Y - 27000001*Y + 243000018*Y - 729000081
roots zz;
{Y=3.0,
Y=3.0,
Y= - 3.0,
Y= - 3.0,
Y=3.00000017,
Y= - 3.00000017}
%{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);
6 4 2
ZZ := 1000000*Y - 9000001*Y + 27000006*Y - 27000009
roots zz;
{Y=1.7320508,
Y=1.7320508,
Y= - 1.7320508,
Y= - 1.7320508,
Y=1.7320511,
Y= - 1.7320511}
%{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;
*** Domain mode COMPLEX changed to COMPLEX-ROUNDED
zz := (x-3)*(y+1.732058);
3 2
ZZ := Y + 1.732058*Y - 3*Y - 5.196174
roots zz;
{Y= - 1.732051,Y=1.73205,Y= - 1.732058}
%{Y= - 1.732051,Y=1.73205,Y= - 1.732058}
zz := (x-3)*(y+1.732051);
3 2
ZZ := Y + 1.732051*Y - 3*Y - 5.196153
roots zz;
{Y= - 1.73205081,Y=1.73205,Y= - 1.732051}
%{Y= - 1.73205081,Y=1.73205,Y= - 1.732051}
zz := (x-3)*(y+1.73205081);
3 2
ZZ := Y + 1.73205081*Y - 3*Y - 5.19615243
roots zz;
{Y= - 1.732050808,Y=1.73205,Y= - 1.73205081}
%{Y= - 1.732050808,Y=1.73205,Y= - 1.73205081}
% 99) minimum accuracy specified using rootacc.
rootacc 10;
10
roots zz;
{Y= - 1.732050808,Y=1.732050808,Y= - 1.73205081}
%{Y= - 1.732050808,Y=1.732050808,Y= - 1.73205081}
% 100)
off rounded;
*** Domain mode COMPLEX-ROUNDED changed to COMPLEX
rootacc 6;
6
zz := (y-i)*(x-2)*(1000000x-2000001);
5 4 3 2
ZZ := 1000000*Y - (1000000*I)*Y - 4000001*Y + (4000001*I)*Y
+ 4000002*Y - 4000002*I
roots zz;
{Y=1.4142136,Y= - 1.4142136,Y=1.4142139,Y= - 1.4142139,Y=I}
%{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);
2
ZZ := 100000000000000*Y - 400000000000001*Y + 400000000000002
roots zz;
{Y=2.00000000000001,Y=2.0}
%{Y=2.0,Y=2.00000000000001}
% 102) still higher precision needed.
zz:= (y-2)*(10000000000000000000y-20000000000000000001);
2
ZZ := 10000000000000000000*Y - 40000000000000000001*Y
+ 40000000000000000002
roots zz;
{Y=2.000 00000 00000 00000 1,Y=2.0}
%{Y=2.000 00000 00000 00000 1,Y=2.0}
% 103) increase in precision required for substituted polynomial.
zz:= (x-2)*(10000000000x-20000000001);
4 2
ZZ := 10000000000*Y - 40000000001*Y + 40000000002
roots zz;
{Y=1.41421356241,Y= - 1.41421356241,Y=1.41421356237,Y
= - 1.41421356237}
%{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);
4 2
ZZ := 100000000000000*Y - 400000000000001*Y + 400000000000002
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}
%{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);
3 2
ZZ := (8100000000*Y - (34200000000 + 81*I)*Y
+ (46000000000 + 252*I)*Y - (20000000000 + 180*I))/100000000
roots zz;
{Y=1.111111111,Y=2.0,Y=1.111111111 + 0.00000001*I}
%{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}
%{Y=1.111111111,Y=2.0}
% 107) realroots now uses powergcd logic whenever it is applicable.
zz := (x-1)*(x-2)*(x-3);
6 4 2
ZZ := Y - 6*Y + 11*Y - 6
realroots zz;
{Y= - 1,
Y=1,
Y= - 1.41421,
Y=1.41421,
Y= - 1.73205,
Y=1.73205}
%{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}
%{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);
4 3 2
ZZ := (72900000000*X - 388800000729*X + 756000003078*X
- 640000004140*X + 200000001800)/100000000
roots zz;
{X=1.111111111,X=1.111111111,X=1.111111121,X=2.0}
%{X=1.111111111,X=1.111111111,X=1.111111121,X=2.0}
nearestroot(zz,1);
{X=1.111111111}
%{X=1.111111111}
nearestroot(zz,1.5);
{X=1.111111121}
%{X=1.111111121}
nearestroot(zz,1.65);
{X=2.0}
%{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);
4
ZZ := (396900000000000000000000000000*X
3
- 1789200000000400869000000000000*X
2
+ 3024400000001361556000000003969*X
- 2272000000001541380000000008946*X
+ 640000000000581600000000005040)/100000000000000000000000000
roots zz;
{X=1.1111111111121,X=1.142 85714 28571 43,X=1.1111111111111,X
=1.142 85714 28571 53}
%{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;
18 4 2
ZZ := 1000000*R + 250000000000*R - 1000000*R + 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}
%{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;
Time: 573342 ms
end;
4: 4:
Quitting
Sat Jun 29 13:58:15 PDT 1991