123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530153115321533153415351536153715381539154015411542154315441545154615471548154915501551155215531554155515561557155815591560156115621563156415651566156715681569157015711572157315741575157615771578157915801581158215831584158515861587158815891590159115921593159415951596159715981599160016011602160316041605160616071608160916101611161216131614161516161617161816191620162116221623162416251626162716281629163016311632163316341635163616371638163916401641164216431644164516461647164816491650165116521653165416551656165716581659166016611662166316641665166616671668166916701671167216731674167516761677167816791680168116821683168416851686168716881689169016911692169316941695169616971698169917001701170217031704170517061707170817091710171117121713171417151716171717181719172017211722172317241725172617271728172917301731173217331734173517361737173817391740174117421743174417451746174717481749175017511752175317541755175617571758175917601761176217631764176517661767176817691770177117721773177417751776177717781779178017811782178317841785178617871788178917901791179217931794179517961797179817991800180118021803180418051806180718081809181018111812181318141815181618171818181918201821182218231824182518261827182818291830183118321833183418351836183718381839184018411842184318441845184618471848184918501851185218531854185518561857185818591860186118621863186418651866186718681869187018711872187318741875187618771878187918801881188218831884188518861887188818891890189118921893189418951896189718981899190019011902190319041905190619071908190919101911191219131914191519161917191819191920192119221923192419251926192719281929193019311932193319341935193619371938193919401941194219431944194519461947194819491950195119521953195419551956195719581959196019611962196319641965196619671968196919701971197219731974197519761977197819791980198119821983198419851986198719881989199019911992199319941995199619971998199920002001200220032004200520062007200820092010201120122013201420152016201720182019202020212022202320242025202620272028202920302031203220332034203520362037203820392040204120422043204420452046204720482049205020512052205320542055205620572058205920602061206220632064206520662067206820692070207120722073207420752076207720782079208020812082208320842085208620872088208920902091209220932094209520962097209820992100210121022103210421052106210721082109211021112112211321142115211621172118211921202121212221232124212521262127212821292130213121322133213421352136213721382139214021412142214321442145214621472148214921502151215221532154215521562157215821592160216121622163216421652166216721682169217021712172217321742175217621772178217921802181218221832184218521862187218821892190219121922193219421952196219721982199220022012202220322042205220622072208220922102211221222132214221522162217221822192220222122222223222422252226222722282229223022312232223322342235223622372238223922402241224222432244224522462247224822492250225122522253225422552256225722582259226022612262226322642265226622672268226922702271227222732274227522762277227822792280228122822283228422852286228722882289229022912292229322942295229622972298229923002301230223032304230523062307230823092310231123122313231423152316231723182319232023212322232323242325232623272328232923302331233223332334233523362337233823392340234123422343234423452346234723482349235023512352235323542355235623572358235923602361236223632364236523662367236823692370237123722373237423752376237723782379238023812382238323842385238623872388238923902391239223932394239523962397239823992400240124022403240424052406240724082409241024112412241324142415241624172418241924202421242224232424242524262427242824292430243124322433243424352436243724382439244024412442244324442445244624472448244924502451245224532454245524562457245824592460246124622463246424652466246724682469247024712472247324742475247624772478247924802481248224832484248524862487248824892490249124922493249424952496249724982499250025012502250325042505250625072508250925102511251225132514251525162517251825192520252125222523252425252526252725282529253025312532253325342535253625372538253925402541254225432544254525462547254825492550255125522553255425552556255725582559256025612562256325642565256625672568256925702571257225732574257525762577257825792580258125822583258425852586258725882589259025912592259325942595259625972598259926002601260226032604260526062607260826092610261126122613261426152616261726182619262026212622262326242625262626272628262926302631263226332634263526362637263826392640264126422643264426452646264726482649265026512652265326542655265626572658265926602661266226632664266526662667266826692670267126722673267426752676267726782679268026812682268326842685268626872688268926902691269226932694269526962697269826992700270127022703270427052706270727082709271027112712271327142715271627172718271927202721272227232724272527262727272827292730273127322733273427352736273727382739274027412742274327442745274627472748274927502751275227532754275527562757275827592760276127622763276427652766276727682769277027712772277327742775277627772778277927802781278227832784278527862787278827892790279127922793279427952796279727982799280028012802280328042805280628072808280928102811281228132814281528162817281828192820282128222823282428252826282728282829283028312832283328342835283628372838283928402841284228432844284528462847284828492850285128522853285428552856285728582859286028612862286328642865286628672868286928702871287228732874287528762877287828792880288128822883288428852886288728882889289028912892289328942895289628972898289929002901290229032904290529062907290829092910291129122913291429152916291729182919292029212922292329242925292629272928292929302931293229332934293529362937293829392940294129422943294429452946294729482949295029512952295329542955295629572958295929602961296229632964296529662967296829692970297129722973297429752976297729782979298029812982298329842985298629872988298929902991299229932994299529962997299829993000300130023003300430053006300730083009301030113012301330143015301630173018301930203021302230233024302530263027302830293030303130323033303430353036303730383039304030413042304330443045304630473048304930503051305230533054305530563057305830593060306130623063306430653066306730683069307030713072307330743075307630773078307930803081308230833084308530863087308830893090309130923093309430953096309730983099310031013102310331043105310631073108310931103111311231133114311531163117311831193120312131223123312431253126312731283129313031313132313331343135313631373138313931403141314231433144314531463147314831493150315131523153315431553156315731583159316031613162316331643165316631673168316931703171317231733174317531763177317831793180318131823183318431853186318731883189319031913192319331943195319631973198319932003201320232033204320532063207320832093210321132123213321432153216321732183219322032213222322332243225322632273228322932303231323232333234323532363237323832393240324132423243324432453246324732483249325032513252325332543255325632573258325932603261326232633264326532663267326832693270327132723273327432753276327732783279328032813282328332843285328632873288328932903291329232933294329532963297329832993300330133023303330433053306330733083309331033113312331333143315331633173318331933203321332233233324332533263327332833293330333133323333333433353336333733383339334033413342334333443345334633473348334933503351335233533354335533563357335833593360336133623363336433653366336733683369337033713372337333743375337633773378337933803381338233833384338533863387338833893390339133923393339433953396339733983399340034013402340334043405340634073408340934103411341234133414341534163417341834193420342134223423342434253426342734283429343034313432343334343435343634373438343934403441344234433444344534463447344834493450345134523453 |
- Codemist Standard Lisp 3.54 for DEC Alpha: May 23 1994
- Dump file created: Mon May 23 10:39:11 1994
- REDUCE 3.5, 15-Oct-93 ...
- Memory allocation: 6023424 bytes
- +++ About to read file ndotest.red
- % Tests of the root finding package.
- % Author: Stanley L. Kameny (stan%valley.uucp@rand.org)
- comment This test file works only with Reduce version 3.5 and later
- and contains examples all of which are solved by roots mod 1.94.
- Answers are rounded to the value given by rootacc (default = 6)
- unless higher accuracy is needed to separate roots. Format may differ
- from that given here, but root order and values should agree exactly.
- (Although the function ROOTS may obtain the roots in a diffferent
- order, they are sorted into a standard order in mod 1.94.)
- In the following, problems 20) and 82) are time consuming and
- have been commented out to speed up the test.
- The hard examples 111) through 115) almost double the test time
- but are necessary to test some logical paths.
- A new "hardest" example has been added as example 116). It is
- commented out, since it is time consuming, but it is solved by roots
- mod 1.94. The time needed to run the three commented-out examples is
- almost exactly equal to the time for the rest of the test. Users of
- fast computers can uncomment the lines marked with %**%. The three
- examples by themselves are contained in the test file rootsxtr.tst.
- When answers are produced which require precision increase for
- printing out or input of roots, roots functions cause precision
- increase to occur. If the precision is already higher than the
- default value, a message is printed out warning of the the precision
- normally needed for input of those values.$
- roots x;
- {x=0}
- % To load roots package.
- write "This is Roots Package test ", symbolic roots!-mod$
- This is Roots Package test Mod 1.94, 28 May 1993.
- % Simple root finding.
- showtime;
- Time: 33 ms plus GC time: 134 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=i,
- x= - i,
- x=2.0,
- x=2.000001}
- %{x=i,x= -i,x=2.0,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=1.0e+5}
- %{x=1,x=1.1,x=1000.0,x=1001.0,x=1.0E+5}
- % 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=0.99,
- x=1,
- x=1 + 0.01*i,
- x=1 - 0.01*i,
- x=1.01}
- %{x=0.99,x=1,x=1 + 0.01*i,x=1 - 0.01*i,x=1.01}
- % 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,
- x=1.01,
- x=1.02,
- x=1.03}
- %{x=1,x=1.01,x=1.02,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,
- x=1 + 0.1*i,
- x=1 - 0.1*i,
- x=1.01,
- x=1.02}
- %{x=1,x=1 + 0.1*i,x=1 - 0.1*i,x=1.01,x=1.02}
- % 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=1.99,
- x=2.0,
- x=2.0 + 0.01*i,
- x=2.0 - 0.01*i,
- x=2.01}
- %{x=1.99,x=2.0,x=2.0 + 0.01*i,x=2.0 - 0.01*i,x=2.01}
- % 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= - 1.60371,
- x=-1.13237 + 1.13805*i,
- x=-1.13237 - 1.13805*i,
- x= - 0.0444444,
- x=0.00555357 + 1.60944*i,
- x=0.00555357 - 1.60944*i,
- x=1.14348 + 1.13804*i,
- x=1.14348 - 1.13804*i,
- x=1.61483}
- %{x= - 1.60371,x=-1.13237 + 1.13805*i,x=-1.13237 - 1.13805*i,
- % x= - 0.0444444,x=0.00555357 + 1.60944*i,x=0.00555357 - 1.60944*i,
- % x=1.14348 + 1.13804*i,x=1.14348 - 1.13804*i,x=1.61483}
- comment 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;
- {x=-3.3584 + 1.61732*i,
- x=-3.3584 - 1.61732*i,
- x=-0.829456 + 3.63408*i,
- x=-0.829456 - 3.63408*i,
- x=5.0025e-29 + 0.00100005*i,
- x=5.0025e-29 - 0.00100005*i,
- x=2.32408 + 2.91431*i,
- x=2.32408 - 2.91431*i,
- x=3.72754}
- %{x=-3.3584 + 1.61732*i,x=-3.3584 - 1.61732*i,
- % x=-0.829456 + 3.63408*i,x=-0.829456 - 3.63408*i,
- % x=5.0025E-29 + 0.00100005*i,x=5.0025E-29 - 0.00100005*i,
- % x=2.32408 + 2.91431*i,x=2.32408 - 2.91431*i,x=3.72754}
- comment 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=-2.189157 + 1.054242*i,
- x=-2.189157 - 1.054242*i,
- x=-0.5406772 + 2.368861*i,
- x=-0.5406772 - 2.368861*i,
- x=1.6e-26 + 0.001414214*i,
- x=1.6e-26 - 0.001414214*i,
- x=1.514944 + 1.899679*i,
- x=1.514944 - 1.899679*i,
- x=2.429781}
- %{x=-2.189157 + 1.054242*i,x=-2.189157 - 1.054242*i,
- % x=-0.5406772 + 2.368861*i,x=-0.5406772 - 2.368861*i,
- % x=1.6E-26 + 0.001414214*i,x=1.6E-26 - 0.001414214*i,
- % x=1.514944 + 1.899679*i,x=1.514944 - 1.899679*i,x=2.429781}
- % 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= - 6.143833,
- x=-4.452438 + 0.02123455*i,
- x=-4.452438 - 0.02123455*i,
- x= - 2.950367,
- x= - 2.003647,
- x= - 0.9999775}
- %{x= - 6.143833,x=-4.452438 + 0.02123455*i,x=-4.452438 - 0.02123455*i,
- % x= - 2.950367,x= - 2.003647,x= - 0.9999775}
- % 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= - 8.437546,
- x=-6.494828 + 1.015417*i,
- x=-6.494828 - 1.015417*i,
- x=-4.295858 + 0.2815097*i,
- x=-4.295858 - 0.2815097*i,
- x= - 2.982725,
- x= - 2.000356,
- x= - 0.9999996}
- %{x= - 8.437546,x=-6.494828 + 1.015417*i,x=-6.494828 - 1.015417*i,
- % x=-4.295858 + 0.2815097*i,x=-4.295858 - 0.2815097*i,
- % x= - 2.982725,x= - 2.000356,x= - 0.9999996}
- % 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= - 10.80988,
- x=-8.70405 + 1.691061*i,
- x=-8.70405 - 1.691061*i,
- x=-6.046279 + 1.134321*i,
- x=-6.046279 - 1.134321*i,
- x= - 4.616444,
- x= - 4.075943,
- x= - 2.998063,
- x= - 2.000013,
- x= - 1}
- %{x= - 10.80988,x=-8.70405 + 1.691061*i,x=-8.70405 - 1.691061*i,
- % x=-6.046279 + 1.134321*i,x=-6.046279 - 1.134321*i,x= - 4.616444,
- % x= - 4.075943,x= - 2.998063,x= - 2.000013,x= - 1}
- % 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= - 13.1895,
- x=-11.02192 + 2.23956*i,
- x=-11.02192 - 2.23956*i,
- x=-7.953917 + 1.948001*i,
- x=-7.953917 - 1.948001*i,
- x=-5.985629 + 0.8094247*i,
- x=-5.985629 - 0.8094247*i,
- x= - 4.880956,
- x= - 4.007117,
- x= - 2.999902,
- x= - 2.0,
- x= - 1}
- %{x= - 13.1895,x=-11.02192 + 2.23956*i,x=-11.02192 - 2.23956*i,
- % x=-7.953917 + 1.948001*i,x=-7.953917 - 1.948001*i,
- % x=-5.985629 + 0.8094247*i,x=-5.985629 - 0.8094247*i,
- % x= - 4.880956,x= - 4.007117,x= - 2.999902,x= - 2.0,x= - 1}
- % 20) W(20) 10 real roots plus 5 complex pairs. (The original problem)
- % This example is commented out, since it takes significant time without
- % being particularly difficult or checking out new paths:
- %**% zz:= x**19+10**7*for j:=1:20 product (x+j); roots zz;
- %{x= - 20.78881,x=-19.45964 + 1.874357*i,x=-19.45964 - 1.874357*i,
- % x=-16.72504 + 2.731577*i,x=-16.72504 - 2.731577*i,
- % x=-14.01105 + 2.449466*i,x=-14.01105 - 2.449466*i,
- % x=-11.82101 + 1.598621*i,x=-11.82101 - 1.598621*i,
- % x=-10.12155 + 0.6012977*i,x=-10.12155 - 0.6012977*i,
- % x= - 8.928803,x= - 8.006075,x= - 6.999746,x= - 6.000006,
- % x= - 5.0,x= - 4.0,x= - 3.0,x= - 2.0,x= - 1}
- 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,
- x=1.01,
- x=1.02,
- x=1.03,
- x=1.04,
- x=1.05}
- %{x=1,x=1.01,x=1.02,x=1.03,x=1.04,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,
- x=1.001,
- x=1.002,
- x=1.003,
- x=1.004,
- x=1.005}
- %{x=1,x=1.001,x=1.002,x=1.003,x=1.004,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,
- x=1.0001,
- x=1.0002,
- x=1.0003,
- x=1.0004}
- %{x=1,x=1.0001,x=1.0002,x=1.0003,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=0.9,
- x=0.929289 + 0.0707107*i,
- x=0.929289 - 0.0707107*i,
- x=1,
- x=1 + 0.1*i,
- x=1 - 0.1*i,
- x=1.07071 + 0.0707107*i,
- x=1.07071 - 0.0707107*i,
- x=1.1}
- %{x=0.9,x=0.929289 + 0.0707107*i,x=0.929289 - 0.0707107*i,
- % x=1,x=1 + 0.1*i,x=1 - 0.1*i,
- % x=1.07071 + 0.0707107*i,x=1.07071 - 0.0707107*i,x=1.1}
- % 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=0.99,
- x=0.992929 + 0.00707107*i,
- x=0.992929 - 0.00707107*i,
- x=1,
- x=1 + 0.01*i,
- x=1 - 0.01*i,
- x=1.00707 + 0.00707107*i,
- x=1.00707 - 0.00707107*i,
- x=1.01}
- %{x=0.99,x=0.992929 + 0.00707107*i,x=0.992929 - 0.00707107*i,
- % x=1,x=1 + 0.01*i,x=1 - 0.01*i,
- % x=1.00707 + 0.00707107*i,x=1.00707 - 0.00707107*i,x=1.01}
- % 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=0.999,
- x=0.999293 + 0.000707107*i,
- x=0.999293 - 0.000707107*i,
- x=1,
- x=1 + 0.001*i,
- x=1 - 0.001*i,
- x=1.00071 + 0.000707107*i,
- x=1.00071 - 0.000707107*i,
- x=1.001}
- %{x=0.999,x=0.999293 + 0.000707107*i,x=0.999293 - 0.000707107*i,
- % x=1,x=1 + 0.001*i,x=1 - 0.001*i,
- % x=1.00071 + 0.000707107*i,x=1.00071 - 0.000707107*i,x=1.001}
- % 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= - 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= - 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}
- % 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=0.99,
- x=1 + 0.01*i,
- x=1 - 0.01*i,
- x=1.01,
- x=1.0199,
- x=1.02 + 0.0001*i,
- x=1.02 - 0.0001*i,
- x=1.0201}
- %{x=0.99,x=1 + 0.01*i,x=1 - 0.01*i,x=1.01,
- % x=1.0199,x=1.02 + 0.0001*i,x=1.02 - 0.0001*i,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.0707107*i,
- x=0.170711 + 0.0707107*i,
- x=-0.170711 + 0.0707107*i,
- x=0.170711 - 0.0707107*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.0707107*i,x=0.170711 + 0.0707107*i,
- % x=-0.170711 + 0.0707107*i,x=0.170711 - 0.0707107*i}
- comment 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=1,x=2.0,x=3.0}
- %{x=1,x=2.0,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=-2.0 - 3.0*i,x=1 + 2.0*i}
- %{x=-2.0 - 3.0*i,x=1 + 2.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=-3.14936 + 0.212593*i,x=0.149358 - 2.21259*i}
- %{x=-3.14936 + 0.21259*i,x=0.149358 - 2.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.0e+9 - 12.0*i,x=1 - 0.000000012*i}
- %{x=-1.0E+9 - 12.0*i,x=1 - 0.000000012*i}
- % 42) Multiple real and complex roots cause no difficulty, provided
- % that input is given in integer or rational form, (or if in decimal
- % fraction format, with switch rounded off or adjprec on and
- % coefficients input explicitly,) so that polynomial is stored exactly.
- 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}
- % 42a) would have failed in roots Mod 1.93 and previously (bug)
- realroots zz;
- {x=1.1,x=1.1}
- %{x=1.1,x=1.1}
- % 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.466333*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*i,
- x=2.0*i,
- 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}
- % 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=-4.0*i,
- x=2.0*i}
- %{x= - 2.0,x=2.0,x=-4.0*i,x=2.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=-1 + 1.73205*i,
- x=-1 - 1.73205*i,
- x=2.0,
- x=-2.41613 + 1.19385*i,
- x=0.981383 - 0.646597*i,
- x=1.43475 + 2.45274*i}
- %{x=-1 + 1.73205*i,x=-1 - 1.73205*i,x=2.0,
- % x=-2.41613 + 1.19385*i,x=0.981383 - 0.646597*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;
- 2 3
- i*x + 100000000000000000000*x + 800000000000000000000
- zz := ---------------------------------------------------------
- 100000000000000000000
- rootacc 12;
- 12
- roots zz;
- {x=-2.0 - 3.33333333333e-21*i,x=1 - 1.73205080757*i,x
- =1 + 1.73205080757*i}
- rootacc 6;
- 6
- %{x=-2.0 - 3.33333333333E-21*i,x=1 - 1.73205080757*i,
- % x=1 + 1.73205080757*i}
- % 53) Precision of 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=-1.23457e+8*i,x=0.0000000081*i}
- %{x=-1.23457E+8*i,x=0.0000000081*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= - 1.23457e+8,x= - 1,x= - 0.0000000081}
- %{x= - 1.23457E+8,x= - 1,x= - 0.0000000081}
- % 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.429174,
- x=0.469832,
- x=4.95934,
- x=-0.448056 - 0.19486*i,
- x=0.18139 + 0.417083*i}
- %{x= - 0.429174,x=0.469832,x=4.95934,
- % x=-0.448056 - 0.19486*i,x=0.18139 + 0.417083*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.427317 + 2.09121*i,x=0.427317 - 3.29121*i}
- %{x=-0.427317 + 2.09121*i,x=0.427317 - 3.29121*i}
- % 56a) multiple roots will be found if coefficients read in exactly.
- % Exact read-in will occur unless dmode is rounded or complex-rounded.
- 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}
- % 57a) Bug corrected in mod 1.94. (handling of rational limits)
- zz := (x-1/3)*(x-1/5)*(x-1/7)*(x-1/11);
- 4 3 2
- 1155*x - 886*x + 236*x - 26*x + 1
- zz := --------------------------------------
- 1155
- realroots(zz,1/11,exclude(1/3));
- {x=0.0909091,x=0.142857,x=0.2}
- %{x=0.0909091,x=0.142857,x=0.2}
- realroots(zz,exclude(1/11),1/3);
- {x=0.142857,x=0.2,x=0.333333}
- %{x=0.142857,x=0.2,x=0.333333}
- % 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=-1 - 1.0e+8*i,x=-1.0e+8*i,x= - i}
- %{x=-1 - 1.0E+8*i,x=-1.0E+8*i,x= - 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= - 2.0,x= - 1.2345678901,x= - 1.23456789}
- %{x= - 2.0,x= - 1.2345678901,x= - 1.23456789}
- % 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=0.9999999,
- x=0.99999999,
- x=1,
- x=1.0000001}
- %{x=0.9999999,x=0.99999999,x=1,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=i,
- x= - i,
- x=1.2345678,
- x=1.2345679}
- %{x=i,x= - i,x=1.2345678,x=1.2345679}
- % 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;
- 346859 - 1863580*i 482657 + 2593180*i
- sss := {x=--------------------,x=--------------------}
- 10000000 500000
-
- % 346859 - 1863580*i 482657 + 2593180*i
- %sss := {x=--------------------,x=--------------------}
- % 10000000 500000
- % With roots in rational form, mkpoly can be used to reconstruct a
- % polynomial.
- zz1 := mkpoly sss;
- 2
- zz1 := 5000000000000*x - (4999999500000 + 25000010000000*i)*x
- + 5000012308763 - 2110440*i
- % 2
- %zz1 := 5000000000000*x - (4999999500000 + 25000010000000*i)*x
- %
- % + 5000012308763 - 2110440*i
- % Finding the roots of the new polynomial zz1.
- rr:= roots zz1;
- 346859 - 1863580*i 482657 + 2593180*i
- rr := {x=--------------------,x=--------------------}
- 10000000 500000
-
- % 346859 - 1863580*i 482657 + 2593180*i
- %rr := {x=--------------------,x=--------------------}
- % 10000000 500000
- % The roots are stable to the extent that rr=ss, although zz1 and
- % zz may differ.
- zz1 - zz;
- 2
- 4999999999999*x - (4999999499999 + 25000009999995*i)*x
- + 5000012308762 - 2110440*i
- % 2
- %4999999999999*x - (4999999499999 + 25000009999995*i)*x
- %
- % + 5000012308762 - 2110440*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;
- rr := {x=100000000 + i,
- x=100000000 - i,
- x=i,
- 200000000 - 3*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,
- %
- % 200000000 - 3*i
- % x=-----------------}
- % 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;
- {x= - 1,
- 19999999
- x=----------,
- 10000000
- x=2,
- 2000001
- x=---------}
- 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;
- rr := {x= - 1,
- 19999999
- x=----------,
- 10000000
- x=2,
- 2000001
- x=---------}
- 1000000
- % 19999999 2000001
- %rr := {x= - 1,x=----------,x=2,x=---------}
- % 10000000 1000000
- % 75) The function getroot is convenient for obtaining the value of a
- % root.
- rr1 := getroot(1,rr);
- rr1 := -1
- % 19999999
- % rr1 := ----------
- % 10000000
- % 76) For example, the value can be used as an argument to nearestroot.
- nearestroot(zz,rr1);
- {x= - 1}
-
- % 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.607902 + 0.270641*i,
- x=-0.607902 - 0.270641*i,
- x=-0.44528 + 0.494497*i,
- x=-0.44528 - 0.494497*i,
- x=-0.205664 + 0.632867*i,
- x=-0.205664 - 0.632867*i,
- x=0.069527 + 0.661817*i,
- x=0.069527 - 0.661817*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.650944 + 0.138369*i,
- x=0.650944 - 0.138369*i,
- x=900.0}
- %{x= - 0.665423,x=-0.607902 + 0.270641*i,x=-0.607902 - 0.270641*i,
- % x=-0.44528 + 0.494497*i, x=-0.44528 - 0.494497*i,
- % x=-0.205664 + 0.632867*i,x=-0.205664 - 0.632867*i,
- % x=0.069527 + 0.661817*i,x=0.069527 - 0.661817*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.650944 + 0.138369*i,x=0.650944 - 0.138369*i,x=900.0}
- % 78) a still harder example.
- zz := x**30 - 900x**29 - 2;
- 30 29
- zz := x - 900*x - 2
- roots zz;
- {x= - 0.810021,
- x=-0.791085 + 0.174125*i,
- x=-0.791085 - 0.174125*i,
- x=-0.735162 + 0.340111*i,
- x=-0.735162 - 0.340111*i,
- x=-0.644866 + 0.490195*i,
- x=-0.644866 - 0.490195*i,
- x=-0.524417 + 0.617362*i,
- x=-0.524417 - 0.617362*i,
- x=-0.379447 + 0.715665*i,
- x=-0.379447 - 0.715665*i,
- x=-0.216732 + 0.780507*i,
- x=-0.216732 - 0.780507*i,
- x=-0.04388 + 0.808856*i,
- x=-0.04388 - 0.808856*i,
- x=0.131027 + 0.799383*i,
- x=0.131027 - 0.799383*i,
- x=0.299811 + 0.752532*i,
- x=0.299811 - 0.752532*i,
- x=0.454578 + 0.67049*i,
- x=0.454578 - 0.67049*i,
- x=0.588091 + 0.557094*i,
- x=0.588091 - 0.557094*i,
- x=0.694106 + 0.417645*i,
- x=0.694106 - 0.417645*i,
- x=0.767663 + 0.258664*i,
- x=0.767663 - 0.258664*i,
- x=0.805322 + 0.0875868*i,
- x=0.805322 - 0.0875868*i,
- x=900.0}
- %{x= - 0.810021,x=-0.791085 + 0.174125*i,x=-0.791085 - 0.174125*i,
- % x=-0.735162 + 0.340111*i,x=-0.735162 - 0.340111*i,
- % x=-0.644866 + 0.490195*i,x=-0.644866 - 0.490195*i,
- % x=-0.524417 + 0.617362*i,x=-0.524417 - 0.617362*i,
- % x=-0.379447 + 0.715665*i,x=-0.379447 - 0.715665*i,
- % x=-0.216732 + 0.780507*i,x=-0.216732 - 0.780507*i,
- % x=-0.04388 + 0.808856*i,x=-0.04388 - 0.808856*i,
- % x=0.131027 + 0.799383*i,x=0.131027 - 0.799383*i,
- % x=0.299811 + 0.752532*i,x=0.299811 - 0.752532*i,
- % x=0.454578 + 0.67049*i,x=0.454578 - 0.67049*i,
- % x=0.588091 + 0.557094*i,x=0.588091 - 0.557094*i,
- % x=0.694106 + 0.417645*i,x=0.694106 - 0.417645*i,
- % x=0.767663 + 0.258664*i,x=0.767663 - 0.258664*i,
- % x=0.805322 + 0.0875868*i,x=0.805322 - 0.0875868*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 because it takes significant time
- % without breaking new ground.
- %**% roots rr;
- %{x= - 0.429188,
- % x=-0.419154 + 0.092263*i,x=-0.419154 - 0.092263*i,
- % x=-0.389521 + 0.180211*i,x=-0.389521 - 0.180211*i,
- % x=-0.341674 + 0.259734*i,x=-0.341674 - 0.259734*i,
- % x=-0.277851 + 0.327111*i,x=-0.277851 - 0.327111*i,
- % x=-0.201035 + 0.379193*i,x=-0.201035 - 0.379193*i,
- % x=-0.11482 + 0.413544*i,x=-0.11482 - 0.413544*i,
- % x=-0.0232358 + 0.428559*i,x=-0.0232358 - 0.428559*i,
- % x=0.0694349 + 0.423534*i,x=0.0694349 - 0.423534*i,
- % x=0.158859 + 0.398706*i,x=0.158859 - 0.398706*i,
- % x=0.240855 + 0.355234*i,x=0.240855 - 0.355234*i,
- % x=0.311589 + 0.295153*i,x=0.311589 - 0.295153*i,
- % x=0.367753 + 0.22127*i,x=0.367753 - 0.22127*i,
- % x=0.406722 + 0.13704*i,x=0.406722 - 0.13704*i,
- % x=0.426672 + 0.0464034*i,x=0.426672 - 0.0464034*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.866594 + 0.193562*i,
- x=-0.866594 - 0.193562*i,
- x=-0.697397 + 0.473355*i,
- x=-0.697397 - 0.473355*i,
- x=-0.510014 + 0.716449*i,
- x=-0.510014 - 0.716449*i,
- x=-0.161318 + 0.87905*i,
- x=-0.161318 - 0.87905*i,
- x=0.182294 + 0.828368*i,
- x=0.182294 - 0.828368*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.845617 + 0.142879*i,
- x=0.845617 - 0.142879*i}
- %{x=-0.866594 + 0.193562*i,x=-0.866594 - 0.193562*i,
- % x=-0.697397 + 0.473355*i,x=-0.697397 - 0.473355*i,
- % x=-0.510014 + 0.716449*i,x=-0.510014 - 0.716449*i,
- % x=-0.161318 + 0.87905*i,x=-0.161318 - 0.87905*i,
- % x=0.182294 + 0.828368*i,x=0.182294 - 0.828368*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.845617 + 0.142879*i,x=0.845617 - 0.142879*i}
- % 85) five complex roots.
- zz := x**5 - x**3 + i;
- 5 3
- zz := x - x + i
- roots zz;
- {x=-1.16695 - 0.217853*i,
- x=-0.664702 + 0.636663*i,
- x=-0.83762*i,
- x=0.664702 + 0.636663*i,
- x=1.16695 - 0.217853*i}
- %{x=-1.16695 - 0.217853*i,x=-0.664702 + 0.636663*i,x=-0.83762*i,
- % x=0.664702 + 0.636663*i,x=1.16695 - 0.217853*i}
- % Additional capabilities in Mod 1.91.
- % 86) handling of polynomial with huge or infinitesimal coefficients.
- precision reset;
- 12
- on rounded;
- *** Domain mode complex changed to complex-rounded
- precision reset;
- 12
- % so that the system will start this example in floating point. Rounded
- % is on so that the polynomial won't fill the page!
- zz:= 1.0e-500x**3+x**2+x;
- *** ROUNDBF turned on to increase accuracy
- 2
- zz := x*(1.0e-500*x + x + 1)
- roots zz;
- {x=0,x= - 1.0e+500,x= - 1}
- off rounded;
- *** Domain mode complex-rounded changed to complex
- % rounded not normally needed for roots.
- %{x=0,x= - 1.0E+500,x= - 1}
- off roundbf;
- comment 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.783212 + 0.276071*i,
- x=0.152522 - 0.816316*i,
- x=0.63069 + 0.540246*i,
- x=-0.783212 - 0.276071*i,
- x=0.152522 + 0.816316*i,
- x=0.63069 - 0.540246*i,
- x=-0.424222 + 0.734774*i,
- x=-0.424222 - 0.734774*i,
- x=0.848444,
- x=-0.85453 + 1.48009*i,
- x=-0.85453 - 1.48009*i,
- x=1.70906}
- %{x=-0.783212 + 0.276071*i,x=0.152522 - 0.816316*i,
- % x=0.63069 + 0.540246*i,x=-0.783212 - 0.276071*i,
- % x=0.152522 + 0.816316*i,x=0.63069 - 0.540246*i,
- % x=-0.424222 + 0.734774*i,x=-0.424222 - 0.734774*i,x=0.848444,
- % x=-0.85453 + 1.48009*i,x=-0.85453 - 1.48009*i,x=1.70906}
- off powergcd;
- roots zz;
- {x=-0.85453 + 1.48009*i,
- x=-0.85453 - 1.48009*i,
- x=-0.783212 + 0.276071*i,
- x=-0.783212 - 0.276071*i,
- x=-0.424222 + 0.734774*i,
- x=-0.424222 - 0.734774*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.848444,
- x=1.70906}
- on powergcd;
- %{x=-0.85453 + 1.48009*i,x=-0.85453 - 1.48009*i,
- % x=-0.783212 + 0.276071*i,x=-0.783212 - 0.276071*i,
- % x=-0.424222 + 0.734774*i,x=-0.424222 - 0.734774*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.848444,x=1.70906}
- % 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.783212 + 0.276071*i,
- x=0.152522 - 0.816316*i,
- x=0.63069 + 0.540246*i,
- x=-0.783212 - 0.276071*i,
- x=0.152522 + 0.816316*i,
- x=0.63069 - 0.540246*i,
- x=-0.424222 + 0.734774*i,
- x=-0.424222 - 0.734774*i,
- x=0.848444,
- x=-0.85453 + 1.48009*i,
- x=-0.85453 - 1.48009*i,
- x=1.70906}
- %{x=1,x=1,
- % x=-0.783212 + 0.276071*i,x=0.152522 - 0.816316*i,
- % x=0.63069 + 0.540246*i,x=-0.783212 - 0.276071*i,
- % x=0.152522 + 0.816316*i,x=0.63069 - 0.540246*i,
- % x=-0.424222 + 0.734774*i,x=-0.424222 - 0.734774*i,x=0.848444,
- % x=-0.85453 + 1.48009*i,x=-0.85453 - 1.48009*i,x=1.70906}
- off powergcd;
- roots zz;
- {x=1,
- x=1,
- x=-0.85453 + 1.48009*i,
- x=-0.85453 - 1.48009*i,
- x=-0.783212 + 0.276071*i,
- x=-0.783212 - 0.276071*i,
- x=-0.424222 + 0.734774*i,
- x=-0.424222 - 0.734774*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.848444,
- x=1.70906}
- on powergcd;
- %{x=1,x=1,
- % x=-0.85453 + 1.48009*i,x=-0.85453 - 1.48009*i,
- % x=-0.783212 + 0.276071*i,x=-0.783212 - 0.276071*i,
- % x=-0.424222 + 0.734774*i,x=-0.424222 - 0.734774*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.848444,x=1.70906}
- % 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;
- {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.73205081,
- y= - 1.732050808,
- y=1.732050808,
- y=1.73205081}
- on powergcd;
- %{y= - 1.73205081,y= - 1.732050808,y=1.732050808,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}
- % 91a) A new capability for nearestroot:
- nearestroot(zz,1.800000000001);
- *** precision increased to 13
- {y=1.732050836436}
- % should find the root to 13 places.
- %{y=1.732050836436}
- % 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.067777 + 7071.067847*i,
- y=7071.067777 - 7071.067847*i,
- y=-7071.067812 + 7071.067812*i,
- y=7071.067812 - 7071.067812*i}
- %{y= - 1,
- % y=-7071.067777 + 7071.067847*i,y=7071.067777 - 7071.067847*i,
- % y=-7071.067812 + 7071.067812*i,y=7071.067812 - 7071.067812*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.4142139,
- y=1,
- y=1.4142139}
- %{y= - 1.4142136,y= - 1.4142136,y=1.4142136,y=1.4142136,
- % y= - 1.4142139,y=1,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 5 depends upon another close root.
- % Although one of the factors is given in decimal notation, it is not
- % necessary to turn rounded on.
- rootacc 10;
- 10
- % using rootacc to specify the minumum desired accuracy.
- zz := (y^2-5)*(y-2.2360679775);
- 3 2
- 400000000*y - 894427191*y - 2000000000*y + 4472135955
- zz := ---------------------------------------------------------
- 400000000
- % in this case, adding one place to the root near sqrt 5 causes a
- % required increase of 4 places in accuracy of the root at sqrt 5.
- roots zz;
- *** precision increased to 14
- {y= - 2.236067977,y=2.2360679774998,y=2.2360679775}
- %{y= - 2.236067977,y=2.2360679774998,y=2.2360679775}
- realroots zz;
- {y= - 2.236067977,y=2.2360679774998,y=2.2360679775}
- % should get the same answer from realroots.
- %{y= - 2.2360679775,y=2.2360679774998,y=2.2360679775}
- % 99) The same thing also happens when the root near sqrt 5 is on a
- % different square-free factor.
- zz := (y^2-5)^2*(y-2.2360679775);
- 5 4 3 2
- zz := (400000000*y - 894427191*y - 4000000000*y + 8944271910*y
- + 10000000000*y - 22360679775)/400000000
- roots zz;
- {y= - 2.236067977,
- y= - 2.236067977,
- y=2.2360679774998,
- y=2.2360679774998,
- y=2.2360679775}
- %{y= - 2.236067977,y= - 2.236067977,y=2.2360679774998,
- % y=2.2360679774998,y=2.2360679775}
- realroots zz;
- {y= - 2.236067977,
- y= - 2.236067977,
- y=2.2360679774998,
- y=2.2360679774998,
- y=2.2360679775}
- % realroots handles this case also.
- %{y= - 2.236067977,y= - 2.236067977,y=2.2360679774998,y=2.2360679774998,
- % y=2.2360679775}
- % 100)
- 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;
- *** precision increased to 15
- {y=2.0,y=2.00000000000001}
- %{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;
- *** precision increased to 20
- {y=2.0,y=2.0000000000000000001}
- %{y=2.0,y=2.0000000000000000001}
- % 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.41421356237,
- y=1.41421356237,
- y= - 1.41421356241,
- y=1.41421356241}
- %{y= - 1.41421356237,y=1.41421356237,y= - 1.41421356241,y=1.41421356241}
- % 104) still higher precision required for substituted polynomial.
- zz:= (x-2)*(100000000000000x-200000000000001);
- 4 2
- zz := 100000000000000*y - 400000000000001*y + 400000000000002
- roots zz;
- *** input of these values may require precision >= 16
- {y= - 1.414213562373095,
- y=1.414213562373095,
- y= - 1.414213562373099,
- y=1.414213562373099}
- %{y= - 1.414213562373095,y=1.414213562373095,
- % y= - 1.414213562373099,y=1.414213562373099}
- % 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}
- % 108a) new cability in mod 1.94.
- realroots 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}
- % 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;
- *** input of these values may require precision >= 16
- {x=1.1111111111111,
- x=1.1111111111121,
- x=1.142857142857143,
- x=1.142857142857153}
- %{x=1.1111111111111,x=1.1111111111121,
- % x=1.142857142857143,x=1.142857142857153}
- % 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=2.42978*i,
- r=-2.42978*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=-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=-1.89968 - 1.51494*i,
- r=1.89968 + 1.51494*i,
- r=-1.89968 + 1.51494*i,
- r=1.89968 - 1.51494*i,
- r=-2.36886 - 0.540677*i,
- r=2.36886 + 0.540677*i,
- r=-2.36886 + 0.540677*i,
- r=2.36886 - 0.540677*i}
- %{r=2.42978*i,r=-2.42978*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=-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=-1.89968 - 1.51494*i,r=1.89968 + 1.51494*i,
- % r=-1.89968 + 1.51494*i,r=1.89968 - 1.51494*i,
- % r=-2.36886 - 0.540677*i,r=2.36886 + 0.540677*i,
- % r=-2.36886 + 0.540677*i,r=2.36886 - 0.540677*i}
- comment These five examples are very difficult root finding problems
- for automatic root finding (not employing problem-specific
- procedures.) They require extremely high precision and high accuracy
- to separate almost multiple roots (multiplicity broken by a small high
- order perturbation.) The examples are roughly in ascending order of
- difficulty.;
- % 111) Two simple complex roots with extremely small real separation.
- c := 10^-6;
- 1
- c := ---------
- 1000000
- zz:=(x-3c^2)^2+i*c*x^7;
- 7 2
- zz := (1000000000000000000*i*x + 1000000000000000000000000*x
- - 6000000000000*x + 9)/1000000000000000000000000
- roots zz;
- *** precision increased to 33
- {x=-15.0732 + 4.89759*i,
- x=-9.31577 - 12.8221*i,
- x=-1.2e-12 + 15.8489*i,
- x=2.99999999999999999999999999999997e-12
- + 3.3068111527572904325663335008527e-44*i,
- x=3.00000000000000000000000000000003e-12
- - 3.30681115275729043256633350085321e-44*i,
- x=9.31577 - 12.8221*i,
- x=15.0732 + 4.89759*i}
- %{x=-15.0732 + 4.89759*i,x=-9.31577 - 12.8221*i,x=-1.2E-12 + 15.8489*i,
- % x=2.99999999999999999999999999999997E-12
- % + 3.3068111527572904325663335008527E-44*i,
- % x=3.00000000000000000000000000000003E-12
- % - 3.30681115275729043256633350085321E-44*i,
- % x=9.31577 - 12.8221*i,x=15.0732 + 4.89759*i}
- % 112) Four simple complex roots in two close sets.
- c := 10^-4;
- 1
- c := -------
- 10000
- zz:=(x^2-3c^2)^2+i*c^2*x^9;
- 9 4 2
- 100000000*i*x + 10000000000000000*x - 600000000*x + 9
- zz := ----------------------------------------------------------
- 10000000000000000
- roots zz;
- *** input of these values may require precision >= 15
- {x=-37.8622 + 12.3022*i,
- x=-23.4002 - 32.2075*i,
- x=-0.00017320508075689 - 2.41778234660324e-18*i,
- x=-0.000173205080756885 + 2.4177823466027e-18*i,
- x=39.8107*i,
- x=0.000173205080756885 + 2.4177823466027e-18*i,
- x=0.00017320508075689 - 2.41778234660324e-18*i,
- x=23.4002 - 32.2075*i,
- x=37.8622 + 12.3022*i}
- %{x=-37.8622 + 12.3022*i,x=-23.4002 - 32.2075*i,
- % x=-0.00017320508075689 - 2.41778234660324E-18*i,
- % x=-0.000173205080756885 + 2.4177823466027E-18*i,
- % x=39.8107*i,
- % x=0.000173205080756885 + 2.4177823466027E-18*i,
- % x=0.00017320508075689 - 2.41778234660324E-18*i,
- % x=23.4002 - 32.2075*i,x=37.8622 + 12.3022*i}
- % 113) Same example, but with higher minimum root accuracy specified.
- rootacc 20;
- 20
- roots zz;
- {x=-37.862241873586290526 + 12.302188128448775345*i,
- x=-23.400152368145827118 - 32.207546656274351069*i,
- x=-0.00017320508075689014714 - 2.417782346603239319e-18*i,
- x=-0.00017320508075688531157 + 2.417782346602699319e-18*i,
- x=39.810717055651151449*i,
- x=0.00017320508075688531157 + 2.417782346602699319e-18*i,
- x=0.00017320508075689014714 - 2.417782346603239319e-18*i,
- x=23.400152368145827118 - 32.207546656274351069*i,
- x=37.862241873586290526 + 12.302188128448775345*i}
- %{x=-37.862241873586290526 + 12.302188128448775345*i,
- % x=-23.400152368145827118 - 32.207546656274351069*i,
- % x=-0.00017320508075689014714 - 2.417782346603239319E-18*i,
- % x=-0.00017320508075688531157 + 2.417782346602699319E-18*i,
- % x=39.810717055651151449*i,
- % x=0.00017320508075688531157 + 2.417782346602699319E-18*i,
- % x=0.00017320508075689014714 - 2.417782346603239319E-18*i,
- % x=23.400152368145827118 - 32.207546656274351069*i,
- % x=37.862241873586290526 + 12.302188128448775345*i}
- precision reset;
- 33
- % This resets precision and rootacc to nominal.
- % 114) Two extremely close real roots plus a complex pair with extremely
- % small imaginary part.
- c := 10^6;
- c := 1000000
- zz:=(c^2*x^2-3)^2+c^2*x^9;
- 9 4
- zz := 1000000000000*x + 1000000000000000000000000*x
- 2
- - 6000000000000*x + 9
- roots zz;
- *** precision increased to 22
- {x= - 251.189,
- x=-77.6216 + 238.895*i,
- x=-77.6216 - 238.895*i,
- x= - 0.000001732050807568877293531,
- x= - 0.000001732050807568877293524,
- x=0.00000173205 + 3.41926e-27*i,
- x=0.00000173205 - 3.41926e-27*i,
- x=203.216 + 147.645*i,
- x=203.216 - 147.645*i}
- %{x= - 251.189,x=-77.6216 + 238.895*i,x=-77.6216 - 238.895*i,
- % x= - 0.000001732050807568877293531,
- % x= - 0.000001732050807568877293524,
- % x=0.00000173205 + 3.41926E-27*i,x=0.00000173205 - 3.41926E-27*i,
- % x=203.216 + 147.645*i,x=203.216 - 147.645*i}
- % 114a) this example is a critical test for realroots as well.
- realroots zz;
- {x= - 251.189,x= - 0.000001732050807568877293531,x
- = - 0.000001732050807568877293524}
- %{x= - 251.189,x= - 0.000001732050807568877293531,
- % x= - 0.000001732050807568877293524}
- % 115) Four simple complex roots in two extremely close sets.
- c := 10^6;
- c := 1000000
- zz:=(c^2*x^2-3)^2+i*c^2*x^9;
- 9 4
- zz := 1000000000000*i*x + 1000000000000000000000000*x
- 2
- - 6000000000000*x + 9
- roots zz;
- {x=-238.895 + 77.6216*i,
- x=-147.645 - 203.216*i,
- x=-0.00000173205080756887729353 - 2.417782346602969319022e-27*i,
- x=-0.000001732050807568877293525 + 2.417782346602969318968e-27*i,
- x=251.189*i,
- x=0.000001732050807568877293525 + 2.417782346602969318968e-27*i,
- x=0.00000173205080756887729353 - 2.417782346602969319022e-27*i,
- x=147.645 - 203.216*i,
- x=238.895 + 77.6216*i}
- %{x=-238.895 + 77.6216*i,x=-147.645 - 203.216*i,
- % x=-0.00000173205080756887729353 - 2.417782346602969319022E-27*i,
- % x=-0.000001732050807568877293525 + 2.417782346602969318968E-27*i,
- % x=251.189*i,
- % x=0.000001732050807568877293525 + 2.417782346602969318968E-27*i,
- % x=0.00000173205080756887729353 - 2.417782346602969319022E-27*i,
- % x=147.645 - 203.216*i,x=238.895 + 77.6216*i}
- % 116) A new "hardest example" type. This polynomial has two sets of
- % extremely close real roots and two sets of extremely close conjugate
- % complex roots, both large and small, with the maximum accuracy and
- % precision required for the largest roots. Three restarts are
- % required, at progressively higher precision, to find all roots.
- % (to run this example, uncomment the following two lines.)
- %**% zz1:= (10^12x^2-sqrt 2)^2+x^7$ zz2:= (10^12x^2+sqrt 2)^2+x^7$
- %**% zzzz := zz1*zz2$ roots zzzz;
- %{x= - 1.00000000000000000000000000009E+8,
- % x= - 9.99999999999999999999999999906E+7,
- % x= - 0.0000011892071150027210667183,
- % x= - 0.0000011892071150027210667167,
- % x=-5.4525386633262882960501E-28 + 0.000001189207115002721066718*i,
- % x=-5.4525386633262882960501E-28 - 0.000001189207115002721066718*i,
- % x=5.4525386633262882960201E-28 + 0.000001189207115002721066717*i,
- % x=5.4525386633262882960201E-28 - 0.000001189207115002721066717*i,
- % x=0.00000118921 + 7.71105E-28*i,
- % x=0.00000118921 - 7.71105E-28*i,
- % x=4.99999999999999999999999999953E+7
- % + 8.66025403784438646763723170835E+7*i,
- % x=4.99999999999999999999999999953E+7
- % - 8.66025403784438646763723170835E+7*i,
- % x=5.00000000000000000000000000047E+7
- % + 8.66025403784438646763723170671E+7*i,
- % x=5.00000000000000000000000000047E+7
- % - 8.66025403784438646763723170671E+7*i}
- % Realroots strategy on this example is different, but determining the
- % necessary precision and accuracy is tricky.
- %**% realroots zzzz;
- %{x= - 1.00000000000000000000000000009E+8,
- % x= - 9.9999999999999999999999999991E+7,
- % x= - 0.0000011892071150027210667183,
- % x= - 0.0000011892071150027210667167}
- showtime;
- Time: 124512 ms plus GC time: 11399 ms
- end;
- (TIME: roots 124545 136078)
- End of Lisp run after 124.56+12.18 seconds
|