Artifact d76b8e5ace15982d11abb18f7ec7ec8563b6a3b0bc635bd3146923e999544aaa:
- File
r34/xlog/roots.log
— 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: 55621) [annotate] [blame] [check-ins using] [more...]
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