QE on maxima
以下は,現在作成中の QE ツールの出力例です.
・有効桁数が fpprec の多倍長浮動小数点数を使用.
・2つの数値解間の距離が %ez より小ならばそれらは同じ値の根(重根や共通根)に対応し,%ez 以上ならばそれらは異なる値の根に対応する数値解と見做しています.
Maxima 5.39.0 http://maxima.sourceforge.net using Lisp CMU Common Lisp 21b (21B Unicode) Distributed under the GNU Public License. See the file COPYING. Dedicated to the memory of William Schelter. The function bug_report() provides bug reporting information. (%i1) (qvpeds ([all], [c,b,a,x],0,h1,r11,0 ),qe( bfpcad(ext( '( x^4+a*x^2+b*x+c>=0 ) ))) ); Evaluation took 11.3900 seconds (11.4100 elapsed) using 1307.972 MB. (%o1) [[c = root(c,1),b = root(126976*c^3-135*b^4,1), root(8*a*c-9*b^2-2*a^3,1) <= a,true], [root(c,1) < c,b < root(b,1), root(256*c^3-128*a^2*c^2+(144*a*b^2+16*a^4)*c-27*b^4-4*a^3*b^2,2) <= a,true], [root(c,1) < c,b = root(b,1), root(256*c^3-128*a^2*c^2+(144*a*b^2+16*a^4)*c-27*b^4-4*a^3*b^2,1) <= a,true], [root(c,1) < c,root(b,1) < b, root(256*c^3-128*a^2*c^2+(144*a*b^2+16*a^4)*c-27*b^4-4*a^3*b^2,2) <= a,true]] (%i2) (qvpeds ([all], [c,b,a,x],1,h1,r11,0 ),qe( bfpcad(ext( '( x^4+a*x^2+b*x+c>=0 ) ))) ); [fpprec,fpprintprec,%ez,ratepsilon]: [30,30,1.0b-5,1.0b-30] x^4+a*x^2+b*x+c >= 0 x^4+a*x^2+b*x+c >= 0 [[c], [b,256*c^3-27*b^4,1024*c^3-2187*b^4,1024*c^3+3*b^4,4096*c^3+27*b^4, 126976*c^3-135*b^4], [a,8*a*c-9*b^2-2*a^3, 256*c^3-128*a^2*c^2+(144*a*b^2+16*a^4)*c-27*b^4-4*a^3*b^2], [x^4+a*x^2+b*x+c]] 1 multi-roots: 0.0b0 0.0b0 [256*c^3-27*b^4,-27*b^4] dist: 0.0b0 2 multi-roots: 0.0b0 0.0b0 [256*c^3-27*b^4,-27*b^4] dist: 0.0b0 3 multi-roots: 0.0b0 0.0b0 [256*c^3-27*b^4,-27*b^4] dist: 0.0b0 1 multi-roots: 0.0b0 0.0b0 [1024*c^3-2187*b^4,-2187*b^4] dist: 0.0b0 2 multi-roots: 0.0b0 0.0b0 [1024*c^3-2187*b^4,-2187*b^4] dist: 0.0b0 3 multi-roots: 0.0b0 0.0b0 [1024*c^3-2187*b^4,-2187*b^4] dist: 0.0b0 1 multi-roots: 0.0b0 0.0b0 [1024*c^3+3*b^4,3*b^4] dist: 0.0b0 2 multi-roots: 0.0b0 0.0b0 [1024*c^3+3*b^4,3*b^4] dist: 0.0b0 3 multi-roots: 0.0b0 0.0b0 [1024*c^3+3*b^4,3*b^4] dist: 0.0b0 1 multi-roots: 0.0b0 0.0b0 [4096*c^3+27*b^4,27*b^4] dist: 0.0b0 2 multi-roots: 0.0b0 0.0b0 [4096*c^3+27*b^4,27*b^4] dist: 0.0b0 3 multi-roots: 0.0b0 0.0b0 [4096*c^3+27*b^4,27*b^4] dist: 0.0b0 1 multi-roots: 0.0b0 0.0b0 [126976*c^3-135*b^4,-135*b^4] dist: 0.0b0 2 multi-roots: 0.0b0 0.0b0 [126976*c^3-135*b^4,-135*b^4] dist: 0.0b0 3 multi-roots: 0.0b0 0.0b0 [126976*c^3-135*b^4,-135*b^4] dist: 0.0b0 common-roots: 0.0b0 0.0b0 dist: 0.0b0 common-roots: 0.0b0 0.0b0 dist: 0.0b0 common-roots: 0.0b0 0.0b0 dist: 0.0b0 common-roots: 0.0b0 0.0b0 dist: 0.0b0 common-roots: 0.0b0 0.0b0 dist: 0.0b0 1 multi-roots: 4.52295331325015833831979417428b-16*%i-2.44948974278317847258028706083b0 (-4.52295331325015213319900815541b-16*%i)-2.44948974278317772381428108859b0 [256*c^3-128*a^2*c^2+(144*a*b^2+16*a^4)*c-27*b^4-4*a^3*b^2, (-5.0b-1*(16*a^4+6.27069374152493593138504723125b2*a)) -1.74185937264581553649584645312b1*a^3-3.2b1*a^2-5.44b2] dist: 1.17428054512251951790867706837b-15 common-roots: -2.44948974278317809819728407471b0 -2.44948974278317809819728407471b0 dist: 1.57772181044202361082345713057b-30 1 multi-roots: -1.41421356237309504880168872421b0*%i (-1.41421356237309504880168872421b0*%i)-3.86941038573562221171322357299b-43 [256*c^3-128*a^2*c^2+(144*a*b^2+16*a^4)*c-27*b^4-4*a^3*b^2, (-8.0b0*a^4)-3.2b1*a^2-3.2b1] dist: 7.88860905221011805411728660181b-31 2 multi-roots: 1.41421356237309504880261489839b0*%i-2.95418118567534149019656617711b-22 1.41421356237309504880076255003b0*%i+2.9541811856753414902004355875b-22 [256*c^3-128*a^2*c^2+(144*a*b^2+16*a^4)*c-27*b^4-4*a^3*b^2, (-8.0b0*a^4)-3.2b1*a^2-3.2b1] dist: 1.94429469761416284867540568328b-21 common-roots: 0.0b0 0.0b0 dist: 0.0b0 1 multi-roots: 4.52295330147550776509999452971b-16*%i-2.44948974278317847258028887801b0 (-4.52295330147551542569579253843b-16*%i)-2.44948974278317772381427927141b0 [256*c^3-128*a^2*c^2+(144*a*b^2+16*a^4)*c-27*b^4-4*a^3*b^2, (-5.0b-1*(16*a^4+6.27069374152493593138504723125b2*a)) -1.74185937264581553649584645313b1*a^3-3.2b1*a^2 -5.44000000000000000000000000001b2] dist: 1.17428054562583903487602624948b-15 common-roots: -2.44948974278317809819728407471b0 -2.44948974278317809819728407471b0 dist: 1.57772181044202361082345713057b-30 1 multi-roots: 0.0b0 0.0b0 [8*a*c-9*b^2-2*a^3,-2*a^3] dist: 0.0b0 2 multi-roots: 0.0b0 0.0b0 [8*a*c-9*b^2-2*a^3,-2*a^3] dist: 0.0b0 common-roots: 0.0b0 0.0b0 dist: 0.0b0 common-roots: 0.0b0 0.0b0 dist: 0.0b0 1 multi-roots: 8.16496580927726032731894282802b-1-1.6183536594179216924381257046b-22*%i 1.61835365941792169243742042834b-22*%i+8.16496580927726032732961767002b-1 [8*a*c-9*b^2-2*a^3,(-2*a^3)+4.0b0*a-2.17732421580726942061980806641b0] dist: 1.11547535241290012202233114443b-21 1 multi-roots: (-8.47286141428164487440940737737b-20*%i)-1.41421356237309504909493338483b0 8.47286141428164487616706655174b-20*%i-1.41421356237309504850844406359b0 [256*c^3-128*a^2*c^2+(144*a*b^2+16*a^4)*c-27*b^4-4*a^3*b^2, 8.0b0*a^4-3.2b1*a^2+3.2b1] dist: 6.10479709847422546500953724134b-19 2 multi-roots: 1.41421356237309504879618188426b0-1.96859282073234240509586747592b-21*%i 1.96859282073234238751927573217b-21*%i+1.41421356237309504880719556416b0 [256*c^3-128*a^2*c^2+(144*a*b^2+16*a^4)*c-27*b^4-4*a^3*b^2, 8.0b0*a^4-3.2b1*a^2+3.2b1] dist: 1.16962633227785224686827731017b-20 common-roots: 1.41421356237309504880168872421b0 1.41421356237309504880168872421b0 dist: 0.0b0 common-roots: 0.0b0 0.0b0 dist: 0.0b0 common-roots: -1.41421356237309504880168872421b0 -1.41421356237309504880168872421b0 dist: 0.0b0 1 multi-roots: 8.1649658092772603273200008634b-1-1.30109252983440372998384918995b-22*%i 1.30109252983440372998339457472b-22*%i+8.16496580927726032732855963464b-1 [8*a*c-9*b^2-2*a^3,(-2*a^3)+4.0b0*a-2.17732421580726942061980806641b0] dist: 8.94560966410596198552565674457b-22 common-roots: 0.0b0 0.0b0 dist: 0.0b0 1 multi-roots: (-1.16997742006955532323735088816b-17*%i)-3.35546475144558657978639517931b-1 1.16997742006955532102907767889b-17*%i-3.35546475144558652208708812765b-1 [x^4+a*x^2+b*x+c, x^4-4.77860936379455500840332276085b0*x^2 -3.05577241698508990342139722064b0*x-5.0b-1] dist: 2.41004349697564216587243138313b-17 1 multi-roots: (-2.24047760152893241737421755167b-20*%i)-1.0109852688012821415609910299b0 2.24047760152893269859199915596b-20*%i-1.01098526880128214151321562239b0 [x^4+a*x^2+b*x+c, x^4-3.55546677120777806279111679242b0*x^2 -3.05577241698508990342139722064b0*x-5.0b-1] dist: 6.55010344684744407831148213933b-20 1 multi-roots: 2.22935715805258081302014670225b-18*%i-4.14333797761439630281210468432b-1 (-2.22935715805258088784912140497b-18*%i)-4.14333797761439627962993519833b-1 [x^4+a*x^2+b*x+c, x^4-3.42753994633320075795235578053b0*x^2 -2.55577241698508990342139722064b0*x-5.0b-1] dist: 5.02536199426577147169718392348b-18 1 multi-roots: (-7.63751536911083876261312227633b-19*%i)-8.96295167765296091409417104731b-1 7.63751536911083876548501856105b-19*%i-8.96295167765296090286367856344b-1 [x^4+a*x^2+b*x+c, x^4-3.03243266180629550631765220844b0*x^2 -2.55577241698508990342139722064b0*x-5.0b-1] dist: 1.89591805066530344113644713169b-18 1 multi-roots: 4.44329956060120776495651206012b-21*%i-4.7445898461131539492366468499b-1 (-4.44329956060120776491378978904b-21*%i)-4.74458984611315394915833267933b-1 [x^4+a*x^2+b*x+c, x^4-2.89645721365288302103527609247b0*x^2 -2.32127592854140316276702227643b0*x-5.0b-1] dist: 1.18449456333131606287348159343b-20 1 multi-roots: (-5.13518395612155029572934523107b-18*%i)-8.19617912427072101369893178291b-1 5.13518395612155030519077584595b-18*%i-8.19617912427072096157128200993b-1 [x^4+a*x^2+b*x+c, x^4-2.75961902983845013509216731558b0*x^2 -2.32127592854140316276702227643b0*x-5.0b-1] dist: 1.15175247237145207854523711501b-17 1 multi-roots: (-6.2414826969557199979633340867b-13*%i)-6.38943104254668599248510111018b-1 (-6.95918202422128079528203919324b-12*%i)-6.38943104241533885902022292701b-1 [x^4+a*x^2+b*x+c, x^4-2.44948974278317809819728407471b0*x^2 -2.08677944009771642211264733221b0*x-5.0b-1] dist: 1.45826385597869537841530414623b-11 2 multi-roots: (-6.95918202422128079528203919324b-12*%i)-6.38943104241533885902022292701b-1 7.58333029391685279498189770221b-12*%i-6.38943104242614942415515511763b-1 [x^4+a*x^2+b*x+c, x^4-2.44948974278317809819728407471b0*x^2 -2.08677944009771642211264733221b0*x-5.0b-1] dist: 1.4582638571553690760678343437b-11 1 multi-roots: 6.38943104246144034405208996102b-1-4.422502465774425554602186706b-14*%i 1.33349305972666630723236835579b-13*%i+6.38943104246298402619456223129b-1 [x^4+a*x^2+b*x+c, x^4-2.44948974278317809819728407471b0*x^2 +2.08677944009771642211264733221b0*x-5.0b-1] dist: 2.3529170930722626582908906827b-13 2 multi-roots: 1.33349305972666630723236835579b-13*%i+6.38943104246298402619456223129b-1 6.38943104246374990541382696251b-1-8.91242813149223756054934246976b-14*%i [x^4+a*x^2+b*x+c, x^4-2.44948974278317809819728407471b0*x^2 +2.08677944009771642211264733221b0*x-5.0b-1] dist: 2.35287498234870156742355194154b-13 1 multi-roots: 3.72610157565100454685245607331b-20*%i+4.74458984611315394913552119431b-1 4.74458984611315394925945833492b-1-3.72610157565100454690915365651b-20*%i [x^4+a*x^2+b*x+c, x^4-2.89645721365288302103527609247b0*x^2 +2.32127592854140316276702227643b0*x-5.0b-1] dist: 7.55455976814101740503286439967b-20 1 multi-roots: 8.19617912427072098761958869152b-1-3.13943116569295322775812726995b-21*%i 3.13943116569295322775468371035b-21*%i+8.19617912427072098765062510133b-1 [x^4+a*x^2+b*x+c, x^4-2.75961902983845013509216731558b0*x^2 +2.32127592854140316276702227643b0*x-5.0b-1] dist: 7.00404879414403669951788837319b-21 1 multi-roots: 2.46936845160507392345254636383b-17*%i+4.14333797761439612007308924107b-1 4.14333797761439646236895064158b-1-2.46936845160507398188288178659b-17*%i [x^4+a*x^2+b*x+c, x^4-3.42753994633320075795235578053b0*x^2 +2.55577241698508990342139722064b0*x-5.0b-1] dist: 6.00897394505265450935583153395b-17 1 multi-roots: 8.96295167765296090779340009951b-1-8.69963119046618331638821004739b-20*%i 8.69963119046618331598771326227b-20*%i+8.96295167765296090916444951125b-1 [x^4+a*x^2+b*x+c, x^4-3.03243266180629550631765220844b0*x^2 +2.55577241698508990342139722064b0*x-5.0b-1] dist: 2.21520197802208732404660033187b-19 1 multi-roots: 3.35546475144558654573716194067b-1-3.2321873145433364152600510801b-19*%i 3.23218731454333641635946277158b-19*%i+3.35546475144558655613632136628b-1 [x^4+a*x^2+b*x+c, x^4-4.77860936379455500840332276085b0*x^2 +3.05577241698508990342139722064b0*x-5.0b-1] dist: 1.2244617433979626836819761986b-18 1 multi-roots: 1.01098526880128213782781893505b0-2.07015345747187669505738755424b-18*%i 2.07015345747187669030859350914b-18*%i+1.01098526880128214524638771724b0 [x^4+a*x^2+b*x+c, x^4-3.55546677120777806279111679242b0*x^2 +3.05577241698508990342139722064b0*x-5.0b-1] dist: 8.49572269592501254899699170041b-18 1 multi-roots: 1.09731695855008129432850901121b-16*%i-6.29960524947436974960114115801b-1 (-1.09731695855008175020983658478b-16*%i)-6.29960524947436189807096491478b-1 [x^4+a*x^2+b*x+c,x^4-1.19055078897614960606377922946b0*x^2-5.0b-1*x] dist: 8.15248085790726489126442096385b-16 1 multi-roots: 0.0b0 0.0b0 [x^4+a*x^2+b*x+c,x^4-5.0b-1*x^2] dist: 0.0b0 1 multi-roots: 6.29960524947436582381894689913b-1-5.1915207100556332839659272546b-22*%i 5.19152071005563328395652911378b-22*%i+6.29960524947436582385315917366b-1 [x^4+a*x^2+b*x+c,x^4-1.19055078897614960606377922946b0*x^2+5.0b-1*x] dist: 3.57531436008650390541071687971b-21 1 multi-roots: (-8.67770249132643708238919796611b-23*%i)-1.31549862931202979917891670099b0 8.67770249132643708239032323801b-23*%i-1.31549862931202979917841983362b0 [x^4+a*x^2+b*x+c, x^4-4.90268221512206394838567282523b0*x^2 -3.79286913672985525962281904563b0*x+5.0b-1] dist: 5.26306169342786570228478616546b-22 1 multi-roots: 2.61198261984724332021478496099b-1-3.73328701456410231125396194552b-23*%i 3.73328701456410231125394316086b-23*%i+2.61198261984724332021514776701b-1 [x^4+a*x^2+b*x+c, x^4+7.12406849810570318368614385351b0*x^2 -3.79286913672985525962281904563b0*x+5.0b-1] dist: 8.30135823846948556696946435113b-23 1 multi-roots: (-8.35793034501348062623930605099b-22*%i)-1.26839633131681513709955065996b0 8.35793034501348062625014797815b-22*%i-1.26839633131681513709478200116b0 [x^4+a*x^2+b*x+c, x^4-4.51570275957450820467053128978b0*x^2 -3.29286913672985525962281904563b0*x+5.0b-1] dist: 5.05314820570176695705561472995b-21 1 multi-roots: 3.49188104507532910404015377995b-23*%i+2.98842309414307218065285325202b-1 2.98842309414307218065613316054b-1-3.49188104507532910403992032144b-23*%i [x^4+a*x^2+b*x+c, x^4+5.33076228416614959738464811917b0*x^2 -3.29286913672985525962281904563b0*x+5.0b-1] dist: 3.35343542885119089772093274405b-22 1 multi-roots: (-5.10955124220135577792387388171b-22*%i)-1.1495332358829680004419098722b0 5.10955124220135577792837906846b-22*%i-1.14953323588296800043902649925b0 [x^4+a*x^2+b*x+c, x^4-3.58590105561473348013717487389b0*x^2 -2.16812942838935673533957135587b0*x+5.0b-1] dist: 3.05910773620696643491835445611b-21 1 multi-roots: 4.29760429000264513227806165013b-1-4.9825698827127380563723750383b-17*%i 4.9825698827127377393878597465b-17*%i+4.29760429000264633921325792835b-1 [x^4+a*x^2+b*x+c, x^4+2.15309806062091933677807113863b0*x^2 -2.16812942838935673533957135587b0*x+5.0b-1] dist: 1.56516218758945019086233255227b-16 1 multi-roots: (-1.09877730976256428518192469349b-21*%i)-1.00617122801969079931775349078b0 1.09877730976254185519229728272b-21*%i-1.00617122801969079931174334307b0 [x^4+a*x^2+b*x+c, x^4-2.54325618867219212419273854599b0*x^2 -1.04338972004885821105632366611b0*x+5.0b-1] dist: 6.39930635240535273049586094033b-21 1 multi-roots: 6.38943104246272475524586123887b-1-7.84439082707517235056253391037b-20*%i 7.84439082708190432497441062349b-20*%i+6.38943104246272476186112486434b-1 [x^4+a*x^2+b*x+c,x^4-1.04338972004885821105632366611b0*x+5.0b-1] dist: 6.79875661664383746823339492112b-19 1 multi-roots: (-4.8300774831537619911512151596b-22*%i)-9.65994523816328465899519279333b-1 4.8300774707834158989233615011b-22*%i-9.65994523816328465899094849168b-1 [x^4+a*x^2+b*x+c, x^4-2.26361409430139643717899746966b0*x^2 -7.67624175513389459396781893817b-1*x+5.0b-1] dist: 1.05514307182643029247458493579b-21 1 multi-roots: 1.119137819361247989040796719b-15*%i+6.94963770046457076218883279372b-1 6.94963770046457864444214852857b-1-1.11913781936124711403020335568b-15*%i [x^4+a*x^2+b*x+c, x^4-4.13672885217763849890412496005b-1*x^2 -7.67624175513389459396781893817b-1*x+5.0b-1] dist: 2.37301011549313127333803802423b-15 1 multi-roots: 6.33988972265475161381443097681b-22*%i-9.23416931585132456001273979169b-1 (-6.33988972265475161353799507713b-22*%i)-9.23416931585132455992195460933b-1 [x^4+a*x^2+b*x+c, x^4-1.9717229852393250817227538327b0*x^2 -4.91858630977920707737240121529b-1*x+5.0b-1] dist: 9.16663850192252560176095861665b-21 1 multi-roots: 7.49558352273654755133389531399b-1-9.2980021082756567072476690926b-18*%i 9.29800210827546330169426966815b-18*%i+7.49558352273654772974086850138b-1 [x^4+a*x^2+b*x+c, x^4-7.95576487395121227355033572426b-1*x^2 -4.91858630977920707737240121529b-1*x+5.0b-1] dist: 2.57701737215882295851533416416b-17 1 multi-roots: 1.03401156272799309188393045594b-19*%i-8.83251837810700396804617510575b-1 (-1.03401156272799304752829508618b-19*%i)-8.83251837810700395251773718947b-1 [x^4+a*x^2+b*x+c, x^4-1.69948573504594766924564690096b0*x^2 -2.45929315488960353868620060764b-1*x+5.0b-1] dist: 1.56655387320991285239749179312b-18 1 multi-roots: 7.96303936790265302555727618361b-1-3.40419552049028707429902475843b-21*%i 3.40419552049944803968351960399b-21*%i+7.96303936790265302568553128937b-1 [x^4+a*x^2+b*x+c, x^4-1.1137806681726386454165880796b0*x^2 -2.45929315488960353868620060764b-1*x+5.0b-1] dist: 1.45206029523254927419705052618b-20 1 multi-roots: (-5.70609945371498130990553057283b-21*%i)-8.40896415253714543050874201664b-1 5.70609945371498131018153821252b-21*%i-8.40896415253714543011376750803b-1 [x^4+a*x^2+b*x+c,x^4-1.41421356237309504880168872421b0*x^2+5.0b-1] dist: 4.11130989883287052898186383085b-20 2 multi-roots: 8.40896415253714543012849282445b-1-6.53339775168399828589095914175b-21*%i 6.53339775168399828561495150206b-21*%i+8.40896415253714543049401670022b-1 [x^4+a*x^2+b*x+c,x^4-1.41421356237309504880168872421b0*x^2+5.0b-1] dist: 3.88177560186052570543478329165b-20 1 multi-roots: 2.15893555484106906165912452485b-16*%i+8.83251837810700183846822208909b-1 8.83251837810700608209569020613b-1-2.15893555484106797595047500567b-16*%i [x^4+a*x^2+b*x+c, x^4-1.69948573504594766924564690097b0*x^2 +2.45929315488960353868620060765b-1*x+5.0b-1] dist: 6.05412132418777109501344790091b-16 1 multi-roots: 6.37787066397529175282449379547b-20*%i-7.96303936790265302727748521798b-1 (-6.37787066396800919273345251121b-20*%i)-7.963039367902653023965322255b-1 [x^4+a*x^2+b*x+c, x^4-1.1137806681726386454165880796b0*x^2 +2.45929315488960353868620060765b-1*x+5.0b-1] dist: 3.54929751662912295409283425277b-19 1 multi-roots: 9.23416931585132422243653899598b-1-6.27472747421841058511382108183b-18*%i 6.2747274742183885879114162131b-18*%i+9.23416931585132489749815540504b-1 [x^4+a*x^2+b*x+c, x^4-1.9717229852393250817227538327b0*x^2 +4.91858630977920707737240121529b-1*x+5.0b-1] dist: 6.86627313685598614587092899877b-17 1 multi-roots: (-1.99636848170413194768633600558b-15*%i)-7.49558352273655235457557488862b-1 1.99636848170413258293302923552b-15*%i-7.49558352273654292649918892678b-1 [x^4+a*x^2+b*x+c, x^4-7.95576487395121227355033572426b-1*x^2 +4.91858630977920707737240121529b-1*x+5.0b-1] dist: 4.10254003056179335802065997111b-15 1 multi-roots: 9.65994523816328465899290503781b-1-2.66934379792368063978571497834b-22*%i 2.66934379792368063978568047301b-22*%i+9.65994523816328465899323624718b-1 [x^4+a*x^2+b*x+c, x^4-2.26361409430139643717899746966b0*x^2 +7.67624175513389459396781893818b-1*x+5.0b-1] dist: 5.34895175649003922093343053329b-22 1 multi-roots: (-4.29381502396318990315032909949b-19*%i)-6.9496377004645747104791949945b-1 4.29381502396139055973603979073b-19*%i-6.94963770046457469615178632778b-1 [x^4+a*x^2+b*x+c, x^4-4.13672885217763849890412496004b-1*x^2 +7.67624175513389459396781893818b-1*x+5.0b-1] dist: 1.67039524946424430892380307845b-18 1 multi-roots: 1.00617122801969079926553973356b0-1.60075239406241646868279412713b-19*%i 1.60075239406241646862513736191b-19*%i+1.00617122801969079936395710029b0 [x^4+a*x^2+b*x+c, x^4-2.54325618867219212419273854599b0*x^2 +1.04338972004885821105632366611b0*x+5.0b-1] dist: 3.34936273279166823298162816124b-19 1 multi-roots: (-3.77185753862812953984534199095b-24*%i)-6.38943104246272475855369706517b-1 3.77185746499376146262928840184b-24*%i-6.38943104246272475855328903805b-1 [x^4+a*x^2+b*x+c,x^4+1.04338972004885821105632366611b0*x+5.0b-1] dist: 4.14942038334135943047182928317b-23 1 multi-roots: 1.14953323588296799935017434052b0-1.25168566885221419821833160929b-18*%i 1.25168566885221419738370885239b-18*%i+1.14953323588296800153076203092b0 [x^4+a*x^2+b*x+c, x^4-3.58590105561473348013717487389b0*x^2 +2.16812942838935673533957135587b0*x+5.0b-1] dist: 3.31991426545575543023864496728b-18 1 multi-roots: (-2.58295399767902797007565021885b-19*%i)-4.29760429000264573649957309613b-1 2.58295399767902797030577470367b-19*%i-4.29760429000264573499174648236b-1 [x^4+a*x^2+b*x+c, x^4+2.15309806062091933677807113863b0*x^2 +2.16812942838935673533957135587b0*x+5.0b-1] dist: 5.38146323165938018060210022577b-19 1 multi-roots: 1.2683963313168151249994797521b0-1.02128852191271037829509986173b-17*%i 1.02128852191271037157319651058b-17*%i+1.26839633131681514919485290902b0 [x^4+a*x^2+b*x+c, x^4-4.51570275957450820467053128978b0*x^2 +3.29286913672985525962281904563b0*x+5.0b-1] dist: 3.16643045115257576716717302982b-17 1 multi-roots: 8.98733613243950412511068357208b-21*%i-2.98842309414307218067753665713b-1 (-8.98733613243950412511912321165b-21*%i)-2.98842309414307218063144975543b-1 [x^4+a*x^2+b*x+c, x^4+5.33076228416614959738464811916b0*x^2 +3.29286913672985525962281904563b0*x+5.0b-1] dist: 1.85561005632311052538934339445b-20 1 multi-roots: 1.315498629312029794561387534b0-3.5935851244110502148689984252b-18*%i 3.5935851244110502062082304748b-18*%i+1.31549862931202980379594900061b0 [x^4+a*x^2+b*x+c, x^4-4.90268221512206394838567282523b0*x^2 +3.79286913672985525962281904563b0*x+5.0b-1] dist: 1.1701817878693290753327731022b-17 1 multi-roots: 2.54178527994425949722680256301b-23*%i-2.61198261984724332021510649724b-1 (-2.54178527994425949722681244282b-23*%i)-2.61198261984724332021482623075b-1 [x^4+a*x^2+b*x+c, x^4+7.12406849810570318368614385351b0*x^2 +3.79286913672985525962281904563b0*x+5.0b-1] dist: 5.80496512960753920779203206589b-23 [T,F]: [92,175] Evaluation took 13.7400 seconds (13.7800 elapsed) using 1299.854 MB. (%o2) [[c = root(c,1),b = root(126976*c^3-135*b^4,1), root(8*a*c-9*b^2-2*a^3,1) <= a,true], [root(c,1) < c,b < root(b,1), root(256*c^3-128*a^2*c^2+(144*a*b^2+16*a^4)*c-27*b^4-4*a^3*b^2,2) <= a,true], [root(c,1) < c,b = root(b,1), root(256*c^3-128*a^2*c^2+(144*a*b^2+16*a^4)*c-27*b^4-4*a^3*b^2,1) <= a,true], [root(c,1) < c,root(b,1) < b, root(256*c^3-128*a^2*c^2+(144*a*b^2+16*a^4)*c-27*b^4-4*a^3*b^2,2) <= a,true]]
CAD on maxima
SARAG https://github.com/andrejv/maxima/tree/master/share/contrib/sarag を待ちきれない方のために作ってみました.ただし,lifting のネックである根の分離は realroots に任せ,符号判定も近似的です.
使用例(wxmaxima 等で実行すると見易いと思います)
(vpeds ( [x,a,b],0,0,3,0 ),cnnm(cad ( '( a*x^2+b<0 and x>0 )) ) ); (vpeds ( [x,a,b],0,0,11,0 ),cnnm(cad ( '( a*x^2+b<0 and x>0 )) ) ); (vpeds ( [x,a,b],0,0,11,0 ),cnnm(cad ( '( a*x^5+b*x+1<0 and x>0 )) ) ); (vpeds ( [x,a,b],0,0,3,0 ),cnnm(cad ( '( sqrt(x)<a*x+b ) )) );
(%i1) (vpeds ( [x,a],0,0,11,0 ),cnnm(cad ( '( x*(x-2)*(x-3)*(x-5)*(x-7)=a ) )) ); connecting cells. connecting cells. matrix([a < root(3125*a^4+20592*a^3-16391308*a^2+244684800*a+2540160000,1), x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,1)], [a = root(3125*a^4+20592*a^3-16391308*a^2+244684800*a+2540160000,1), x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,1) or x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,2)], [root(3125*a^4+20592*a^3-16391308*a^2+244684800*a+2540160000,1) < a and a < root(3125*a^4+20592*a^3-16391308*a^2+244684800*a+2540160000, 2), x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,1) or x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,2) or x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,3)], [a = root(3125*a^4+20592*a^3-16391308*a^2+244684800*a+2540160000,2), x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,1) or x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,2) or x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,3) or x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,4)], [root(3125*a^4+20592*a^3-16391308*a^2+244684800*a+2540160000,2) < a and a < root(3125*a^4+20592*a^3-16391308*a^2+244684800*a+2540160000, 3), x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,1) or x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,2) or x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,3) or x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,4) or x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,5)], [a = root(3125*a^4+20592*a^3-16391308*a^2+244684800*a+2540160000,3), x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,1) or x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,2) or x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,3) or x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,4)], [root(3125*a^4+20592*a^3-16391308*a^2+244684800*a+2540160000,3) < a and a < root(3125*a^4+20592*a^3-16391308*a^2+244684800*a+2540160000, 4), x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,1) or x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,2) or x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,3)], [a = root(3125*a^4+20592*a^3-16391308*a^2+244684800*a+2540160000,4), x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,1) or x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,2)], [root(3125*a^4+20592*a^3-16391308*a^2+244684800*a+2540160000,4) < a, x = root(x^5-17*x^4+101*x^3-247*x^2+210*x-a,1)]) Evaluation took 2.3600 seconds (2.3600 elapsed) using 540.776 MB.
なお,例によって,コードの大半は結果の簡約に費やされています.
/* 20170324 22:04:46 */ showtime:on$ display2d:false$ prton2():=(prt:print,prt2:print)$ prton():=(prt:print)$ prtoff():=kill(prt,prt2)$ vpeds(v,p,e,d,s):=block([], %vv:v, algexact:true, bsolve:false, dtype:0, if p=1 then (prtoff(),prton()) elseif p=2 then prton2() else (ratprint:off,prtoff()), if e=0 then delv0 (L,vv):=delv0s (L,vv) else delv0 (L,vv):=delv0m (L,vv), if d=0 then getsample (a,b,c):=getsample0 (a,b,c) elseif d=1 then getsample (a,b,c):=getsample1 (a,b,c) elseif d=11 then (dtype:11,getsample (a,b,c):=getsample11 (a,b,c)) elseif d=2 then getsample (a,b,c):=getsample2 (a,b,c) else (getsample (a,b,c):=getsample1 (a,b,c), bsolve:true), if s=1 then dispsample:true else dispsample:false, usersimp:simp, simp:false )$ matchdeclare([aa,bb,cc],true)$ kill("implies")$ infix("implies",60,40)$ (x implies y):=(not(x) or y)$ %P:[]$ getlcP(f):=block([ff:expand(f),d,lc], d:hipow(ff,%v),lc:coeff(ff,%v,d),/*print(%v),*/ if member(%v,listofvars(ff)) then %P:append(%P,[ff]) elseif listofvars(ff)#[] then %C:append(%C,[ff]), if not(constantp(lc)) then (%C:append(%C,[lc]),getlcP(ff-lc*(%v)^d)) )$ gcd2(f,g,v):=block([p1:f,p2:g,p3:1,k], for k:0 unless p3=0 do (p3:remainder (p1,p2,v),p1:p2,p2:p3),factor(content(p1,v)[2]))$ sqfree(f,v):=block([], if member(v,listofvars (f)) then content( divide (f,gcd2 (f,diff(f,v),v),v) [1] ,v) [2] else f )$ sqfreem(L,v):=block([], if is(L=[]) then [] else listify(setify( map(lambda([f],sqfree (f,v)),L) )) )$ resultant1(f,v):=block ([g:sqfree(f,v)], resultant (g,diff (g,v),v ) )$ resultant2(f1,f2,v):=block ( [g1: quotient (f1,gcd2 (f1,f2,v),v), g2: quotient (f2,gcd2 (f1,f2,v),v)], resultant (g1,g2,v ) )$ delv0m(L,v):=block([P:[],PP,i,lp], delv0 (LL,vv):=delv0s (LL,vv), %C:[],%v:v, map(lambda([f],if member(%v,listofvars(f)) then P:append(P,[f]) else %C:append(%C,[f])),L), PP:map("[",P),prt(PP), /* PP:map(lambda([f],(%P:[],getlcP(f),%P)),P),prt(PP), */ map(lambda([f],%C:append(%C,[resultant1(f,%v)])),flatten(PP)), lp:length(PP),/*prt([%C,P]),*/ %C:append(%C,flatten(makelist(makelist(makelist(makelist(resultant2(PP[i][ii],PP[j][jj],%v),ii,1,length(PP[i])),jj,1,length(PP[j])),j,i+1,lp),i,1,lp-1))), listify(setify(delete(false,map(lambda([f],if constantp(f) then false else f),%C)))) )$ delv0s(L,v):=block([P:[],PP,i,lp], %C:[],%v:v, map(lambda([f],if member(%v,listofvars(f)) then P:append(P,[f]) else %C:append(%C,[f])),L), /* PP:map("[",P),print(PP), */ PP:map(lambda([f],(%P:[],getlcP(f),%P)),P),prt(PP), map(lambda([f],%C:append(%C,[resultant1(f,%v)])),flatten(PP)), lp:length(PP),/*prt([%C,P]),*/ %C:append(%C,flatten(makelist(makelist(makelist(makelist(resultant2(PP[i][ii],PP[j][jj],%v),ii,1,length(PP[i])),jj,1,length(PP[j])),j,i+1,lp),i,1,lp-1))), listify(setify(delete(false,map(lambda([f],if constantp(f) then false else f),%C)))) )$ load ("grobner")$ delv(L,vl):=block([p:[sqfreem(L,vl[1])],k], for k:1 thru length(vl)-1 do p:append(p,[sqfreem( delv0(last(p),vl[k]) ,vl[k+1])]), p:append(rest (p,-1), [listify(setify( poly_normalize_list (last (p), [last (vl) ] )))]), reverse(p) )$ factors(f):=block ([g:factor(f)], if atom(g) then g elseif op(g)="+" then g else (if op(g)="/" then g:args(g)[1], if atom(g) then g elseif op(g)="-" then g:-g, if atom(g) then g elseif op(g)="*" then delete(0, map(lambda([e], if constantp(e) then 0 elseif atom(e) then e elseif op(e)="^" then args(e)[1] else e), args(g))) elseif op(g)="^" then args(g)[1] else g ) )$ factorsm(L):=listify(setify(flatten(map(factors,L))))$ delv(L,vl):=block([p:[factorsm(L)],k], for k:1 thru length(vl)-1 do p:append(p,[factorsm( delv0(last(p),vl[k]) )]), p:append(rest (p,-1),[factorsm(last(p))]), reverse(p) )$ realonly:on$ getsample0(sub,p,v):=block([p1,p2,q,k], p1:listify(setify(flatten(map(lambda([f],algsys([f],listofvars(f))),subst(sub[2],p))))), if p1=[] then p2:[0] else (p1:sort(map(rhs,p1),"<"),q:map("[",p1), p2:sort(listify(setify(append(p1,(append([first(p1)-1],p1)+append(p1,[last(p1)+1]))/2))),"<")), makelist([[0.5*k],[v=p2[k]]],k,1,length(p2)) )$ getsample1(sub,p,v):=block([p1,sf,rt,T:[],S:[],p2,q,k], p1:apply(append,map(lambda([f], (sf: subst(sub[2],f), rt:if constantp (sf) then [] else flatten(solve( sf ,v)), rt:if rt=[] then [] else flatten(map(lambda([e],if is(imagpart(rhs(e))=0) then e else [] ), rt )), prt2([p,sf,v,rt]), rt:if rt=[] then [] else (rt:sort (map (rhs,rt),"<" ), makelist([[f,k ],rt[k]],k,1,length (rt) )) )),p)),prt2("p1:"),prt2(p1), if p1=[] then p2:[0] else (map (lambda ( [e], if not(member (e [2],T )) then (T:append(T,[e[2]]),S:append(S,[e])) ),p1 ), prt2 ("S:"), prt2(S), if S=[] then (p1:[],p2:[0]) else ( S:sort (S,lambda ( [u,t], u[2]<t[2] ) ), [q,p1]:[map (first,S),map (last,S)], prt2 ("q:"), prt2(q),prt2 ("p1:"), prt2(p1), p2:sort(listify(setify(append(p1, (append([first(p1)-1],p1)+append(p1,[last(p1)+1]))/2))),"<") ) ), prt2 ("p2:"), prt2(p2), makelist([[ if p1=[] then true elseif is(k=1) then v<first(q) elseif is(k=length (p2) ) then last(q)<v elseif evenp(k) then v=q[k/2] else (q[(k-1)/2]<v and v<q[(k+1)/2]) ],[v=p2[k]]],k,1,length(p2)) )$ load ("grobner") $ getrt(sub,f,v):=block([eqv,subv:listofvars (sub),p1,sf,rt00:[],rt01:[],rt0,rt:[],sl,k],prt2(sub,f,v), p1:makelist ( if atom(sub [1] [k]) then sub[2][k][1] elseif is (op (sub [1] [k] )="=") then args(rhs(sub [1] [k] ))[1] else sub[2][k][1], k,1,length (sub [2] ) ), p1:append (p1, [f] ),prt2 ("p1:",p1), sf:poly_reduced_grobner (p1,append (subv,[v] ) ),prt2 ("sf:",sf), map (lambda ([e],if listofvars (e)=[v] then rt00:append (rt00, [e] ) else rt01:append (rt01, [e] ) ),sf ), if rt00=[] then rt:[] else (rt0:apply(realroots,rt00), if rt0=[] then rt:[] else (sub2:apply("+",map(lambda ( [e], abs (lhs (e[2])-rhs (e[2]) ) ),sub [2] )), /* cf. apply("+",[]) = 0 */ for eqv in rt0 do (sl:algsys(subst (eqv,rt01),subv), ans:apply(min, map (lambda ( [e], subst (e,sub2) ),sl )), if is(ans<10^(-5)) then rt:append (rt, [eqv] ), prt2(eqv,sl,"min:",float(ans)) ) ) ), rt:if rt=[] then [] else (rt:sort (map (rhs,rt),"<" ), makelist([root(f,k),[last(rt00),rt[k]]],k,1,length (rt) )), prt2("rt:",rt),rt )$ rootsepsilon:1.0e-10$ getsample11(sub,p,v):=block([p1,S:[],T:[],p2,p0,q,k], p1:apply(append,map(lambda([f],if member (v,listofvars (f) ) then getrt(sub,f,v) else []),p)), prt2("p1 in 11:",p1), if p1=[] then p2:[0] else (map (lambda ( [e], if not(member (e[2][2],T )) then (S:append(S,[e]),T:append(T,[e[2][2]])) ),p1 ), prt2 ("S:",S), if S=[] then (p1:[],p2:[0]) else ( S:sort (S,lambda ( [u,t], u[2][2]<t[2][2] ) ), [q,p0]:[map(first,S),map(last,S)], p0p:map(first,p0),p1:map(last,p0), prt2 ("q:"), prt2(q),prt2 ("p1:"), prt2(p1), p2:sort(listify(setify(append(p1, (append([first(p1)-1],p1)+append(p1,[last(p1)+1]))/2))),"<") ) ), prt2 ("p2:"), prt2(p2), makelist([[ if p1=[] then true elseif is(k=1) then v<first(q) elseif is(k=length (p2) ) then last(q)<v elseif evenp(k) then v=q[k/2] else (q[(k-1)/2]<v and v<q[(k+1)/2]) ],[ [(if oddp(k) then v-p2[k] else p0p[k/2]), (if oddp(k) then v=p2[k] else v=p1[k/2])] ]],k,1,length(p2)) )$ getsample2(sub,p,v):=block([sub1,rt,q,p1,p2,k], sub1:delete([],map(lambda ( [e], if atom (e) then [] elseif op(e)="=" then e else []), sub[1] ) ), rt:subst(sub1,p), rt:flatten(map (lambda ( [e], if constantp (e) then [] else e ), rt )), if rt=[] then rt:[] else rt:flatten(map (lambda ( [e],solve (e,v) ),rt ) ), if rt=[] then rt:[] else (rt:map (lambda ( [e],[e,ratsimp(subst (sub [2],e )) ] ), map(rhs,rt) ), rt:delete([],map(lambda([e],if imagpart(e [2])=0 then e else [] ), rt )), rt:listify(setify(rt)) ), if rt=[] then (p1:[],p2:[0]) else ( rt:sort (rt,lambda ( [u,t], u[2]<t[2] )),/*print (rt),*/ [q,p1]:[map (first,rt),map (last,rt)], p2:sort(listify(setify(append(p1, (append([first(p1)-1],p1)+append(p1,[last(p1)+1]))/2))),"<") ), /*prt2 ("p2:"), prt2(p2),*/ makelist([[ if p1=[] then true elseif is(k=1) then v<first(q) elseif is(k=length (p2) ) then last(q)<v elseif evenp(k) then v=q[k/2] else (q[(k-1)/2]<v and v<q[(k+1)/2]) ],[v=ratsimp( p2[k] )]],k,1,length(p2)) )$ is2(f):= if atom (f) then f elseif member (op (f), ["or","and"] ) then map (is,f) else is (f)$ kill("Lo","Lc","Go","Gc","Eq")$ infix("Lo",60,40)$ infix("Lc",60,40)$ infix("Go",60,40)$ infix("Gc",60,40)$ infix("Eq",60,40)$ (x Lo y):=if member(extd,listofvars([x,y])) then false else float(10^(-5)+x<y)$ (x Lc y):=if member(extd,listofvars([x,y])) then false else float(x<y+10^(-5))$ (x Go y):=if member(extd,listofvars([x,y])) then false else float(x>y+10^(-5))$ (x Gc y):=if member(extd,listofvars([x,y])) then false else float(10^(-5)+x>y)$ (x Eq y):=if member(extd,listofvars([x,y])) then false else abs(float(x-y))<10^(-5)$ kill("implies")$ infix(implies,60,40)$ "implies"(x,y):=((not(x)) or (y))$ defrule(imp2m, (aa) implies (bb), mor ( mnot((aa)),(bb) )); defrule(noteq2m, (aa) # (bb), mnot (equal((aa),(bb)))); defrule(symneq2m, notequal((aa),(bb)), mnot (equal((aa),(bb)))); cad_monic(F0):=block ( [F,mpp:[],p,sx,sxl,i,k], F:addD(subst(["="=equal,"and"=mand,"or"=mor,"not"=mnot],F0)), F:apply1(F,imp2m,noteq2m,symneq2m), simp:on, for vk in %vv do F:sepaF(F,vk), for vk in %vv do F:sepaS(sepaS(F,vk),vk), for vk in %vv do F:sepaF(F,vk), F:ev(subst([mand="and",mor="or",mnot="not"],F),eval),prt(F), F:scanmap(lambda([f],if atom(f) then f elseif member(op(f),["<","<=",">",">=",equal,notequal]) then (p:dnf(sepalcv(f,%vv[1])),prt(p), if op(p)="or" then (mpp:append (mpp, [args(p)] ), concat(mp,length (mpp) )) else f ) else f),F), print([F,mpp]), if mpp=[] then [F] else ( i:0, sxl:map(lambda([n],(i:i+1,product(concat(sx,i)-k,k,1,n))),map(length,mpp)), sxl:solve(sxl), map(lambda([sxl],subst(ev(subst(mppp=mpp,subst(sxl,makelist(concat(mp,i)=mppp[i][concat(sx,i)],i,1,length(mpp)))),eval),F)),sxl) ) ) $ /* epsilon sgn-det ver. */ cad (F0) :=block ([F,vk,P2:{},pp,vv:reverse(%vv),pn:[[[],[]]],pn0], F:addD(subst(["="=equal,"and"=mand,"or"=mor,"not"=mnot],F0)), F:apply1(F,imp2m,noteq2m,symneq2m), simp:on, for vk in %vv do F:sepaF(F,vk), for vk in %vv do F:sepaS(sepaS(F,vk),vk), for vk in %vv do F:sepaF(F,vk), F1:subst(["<"="Lo",">"="Go","<="="Lc",">="="Gc",equal="Eq"],F), F:ev(subst([mand="and",mor="or",mnot="not"],F),eval), F1:ev(subst([mand="and",mor="or",mnot="not"],F1),eval), prt2(F1), prt(F), /* F:subst(["<"="Lo",">"="Go","<="="Lc",">="="Gc",equal="Eq"],F), */ scanmap (lambda ( [f], if atom (f) then f elseif member (op (f), ["<","<=",">",">=",equal,notequal]) then P2:append (P2, {lhs (f)-rhs (f) } ) else f ),F ), simp:on,P2:listify (P2), /*P2:map(lambda ([f], sqfree (f,(%vv)[1])),P2), */ pp:map(factor,delv(P2,%vv)), prt(pp),%pp:pp, for n:1 thru length (vv) do pn:apply(append,map(lambda([e],map(lambda([t],map(append,e,t)),getsample(e,pp[n],vv[n]))),pn)),prt(length (pn)),pn00:pn, pn: if is(dtype=11) then delete(false,map(lambda([l],if is2( (subst( map(last,l[2]), F1)) ) then (if is(dispsample=true) or is(bsolve=true) then l else l[1])),pn)) else delete(false,map(lambda([l],if is2( (subst(l[2], F1)) ) then (if is(dispsample=true) or is(bsolve=true) then l else l[1])),pn)), /*print(pn),*/ simp:usersimp, if bsolve=true and dispsample=true then cad2s (pn) elseif bsolve=true and dispsample=false then cad2 (pn) else pn )$ cadm(f):=block([usr_display2d], thisout:cad(f), usr_display2d:display2d, set_display('xml), print(apply(matrix,thisout)), display2d:usr_display2d,return(done) )$ addD0(f):=block([DN], if atom(f) then return(f) elseif member(op(f),["<",">","<=",">=",equal]) then (DN:{}, scanmap(lambda([t],if atom(t) then t else (DN:ev(append(DN,{notequal(ratsimp(denom(t)),0)}),simp),t)), f),prt2([f,DN]), DN:ev(substpart("and",DN,0),eval),prt2([f,DN]), if is(DN=false) then return(false) else return( mand((f) , (apply1(DN,symneq2m))) )) else return(f) )$ addD(F):=block([addD0,S], prt2("addD:"), S:scanmap(addD0,F,bottomup),prt2(S), return(S) )$ sepaF0(f,v):=block([p,DN,NU], if atom(f) then return(f) elseif member(op(f),["<",">","<=",">=",equal]) then (p:factor(lhs(f)-rhs(f)), NU:num(p),DN:denom(p),prt2([p,NU,DN]), if member(v,listofvars(DN)) and member(op(f),["<",">"]) then return(apply(op(f),[NU*DN,0])) elseif member(v,listofvars(DN)) and member(op(f),["<=",">="]) then return( mand( (apply(op(f),[NU*DN,0])) , (mnot(equal(DN,0))) ) ) elseif member(v,listofvars(DN)) and op(f)=equal then return( mand( (equal(NU,0)) , (mnot(equal(DN,0))) ) ) else return(f)) else return(f))$ sepaF(F,v):=block([S], prt2("sepaF for",v,":"), S:scanmap(lambda([s],sepaF0(s,v)),F), prt2(S),return(S) )$ sepaSlocal(ff,v):=block([sepaSlocal0,getS,getSqrt,p,Sqrt,A,B,C], sepaSlocal0(f):=block([], if atom(f) or is(getS=true) then return(f) elseif op(f)=sqrt and member(v,listofvars(f)) then (getS:true,getSqrt:f,return(Sqrt)) else return(f)), if atom(ff) then return(ff) elseif member(op(ff),["<",">","<=",">=",equal]) then ( getS:false,p:scanmap(sepaSlocal0,lhs(ff)-rhs(ff)), (if is(getS=true) then (p:subst(getSqrt=Sqrt,p), [A,B]:[coeff(expand(p),Sqrt,1),getSqrt^2],C:expand(A*Sqrt-p), return(mand( [A,B,C,op(ff)] , (getSqrt^2 >= 0) ))) else return(ff)) ) else return(ff))$ defrule(sepaS1,[aa,bb,cc,equal], mand(mor(mand(aa >= 0,cc >= 0),mand(aa < 0,cc <= 0)),/*bb >= 0,*/ equal(cc^2-aa^2*bb,0)) )$ defrule(sepaS2,[aa,bb,cc,"<="], mand(mor(mand(aa < 0,aa^2*bb-cc^2 >= 0), mand(cc >= 0,aa^2*bb-cc^2 <= 0))/*,bb >= 0*/) )$ defrule(sepaS3,[aa,bb,cc,">="], mand(mor(mand(aa > 0,cc^2-aa^2*bb <= 0), mand(cc <= 0,cc^2-aa^2*bb >= 0))/*,bb >= 0*/) )$ defrule(sepaS4,[aa,bb,cc,"<"], mand(mor(mand(aa < 0,cc > 0),mand(aa < 0,aa^2*bb-cc^2 > 0), mand(cc > 0,aa^2*bb-cc^2 < 0))/*,bb >= 0*/) )$ defrule(sepaS5,[aa,bb,cc,">"], mand(mor(mand(aa > 0,cc < 0),mand(aa > 0,cc^2-aa^2*bb < 0), mand(cc < 0,cc^2-aa^2*bb > 0))/*,bb >= 0*/) )$ sepaS(F,v):=block([S], prt2("sepaS:"), S:scanmap(lambda([s],sepaSlocal(s,v)),expand(F)), prt2(S), S:expand(apply1(S,sepaS1,sepaS2,sepaS3,sepaS4,sepaS5)), prt2(S),return(S) )$ sepalcv(f,v):=block([pp,dg,lc], pp:expand(lhs(f)-rhs(f)),dg:hipow(pp,v),lc:ratcoef(pp,v,dg), if atom(f) then return(f) elseif member(op(f),["<",">","<=",">=",equal,notequal]) then (if dg=0 then return(f) else return( (notequal(lc,0) and f) or (equal(lc,0) and sepalcv(apply(op(f),[expand(pp-lc*v^dg),0]),v)) )) else return(f))$ /* cnf dnf */ kill("Or"); kill("And"); infix(Or,60,40)$ infix(And,60,40)$ matchdeclare([aa,bb,cc],true)$ defrule(ru05r, (aa) Or ((bb) And (cc)), ((aa) Or (bb)) And ((aa) Or (cc)))$ defrule(ru08r, ((bb) And (cc)) Or (aa), ((aa) Or (bb)) And ((aa) Or (cc)))$ defrule(ru05a, (aa) And ((bb) Or (cc)), ((aa) And (bb)) Or ((aa) And (cc)))$ defrule(ru08a, ((bb) Or (cc)) And (aa), ((aa) And (bb)) Or ((aa) And (cc)))$ defrule(ru00r, (aa) Or (bb), (aa) or (bb))$ defrule(ru00a, (aa) And (bb), (aa) and (bb))$ orand2orand(f):=block([g],g:f, if atom(g) then return(g), if op(g)="or" then g:xreduce("Or",args(g)), if op(g)="and" then g:xreduce("And",args(g)),g)$ cnf(f):=applyb2(applyb2(scanmap(orand2orand,f),ru05r,ru08r),ru00r,ru00a)$ dnf(f):=applyb2(applyb2(scanmap(orand2orand,f),ru05a,ru08a),ru00r,ru00a)$ cad2(F):=map(lambda([L], (U:[], map(lambda([f], (if atom (f) then f elseif op(f)="=" then (sol:solve(subst(U,[first(rhs(f))]),lhs(f)),prt2(sol), sol:lhs(f)=(sort( delete([],map(lambda([e], if imagpart(subst(L[2],e))=0 then [e,subst(L[2],e)] else []), map(rhs,sol) )) , lambda ( [ss,tt], ss[2]<tt[2] ) ))[last(rhs(f))][1], U:append (U, [sol] ),sol) elseif op(f)="<" and listp(rhs (f)) then (sol:solve(subst(U,[first(rhs(f))]),lhs(f)),prt2(sol), sol:lhs(f)<(sort( delete([],map(lambda([e], if imagpart(subst(L[2],e))=0 then [e,subst(L[2],e)] else []), map(rhs,sol) )) , lambda ( [ss,tt], ss[2]<tt[2] ) ))[last(rhs(f))][1]) elseif op(f)="<" and listp(lhs (f)) then (sol:solve(subst(U,[first(lhs(f))]),rhs(f)),prt2(sol), sol:(sort( delete([],map(lambda([e], if imagpart(subst(L[2],e))=0 then [e,subst(L[2],e)] else []), map(rhs,sol) )) , lambda ( [ss,tt], ss[2]<tt[2] ) ))[last(lhs(f))][1]<rhs(f)) elseif op(f)="and" then (sol1:solve(subst(U,[first(lhs((args(f))[1]))]),rhs((args(f))[1])), sol1:map(rhs,sol1), sol2:solve(subst(U,[first(rhs((args(f))[2]))]),lhs((args(f))[2])), sol2:map(rhs,sol2), prt2([sol1,sol2]), (sort( delete([],map(lambda([e], if imagpart(subst(L[2],e))=0 then [e,subst(L[2],e)] else []), sol1 )) , lambda ( [ss,tt], ss[2]<tt[2] ) ))[last(lhs((args(f))[1]))][1] < rhs((args(f))[1]) and lhs((args(f))[2]) < (sort( delete([],map(lambda([e], if imagpart(subst(L[2],e))=0 then [e,subst(L[2],e)] else []), sol2 )) , lambda ( [ss,tt], ss[2]<tt[2] ) ))[last(rhs((args(f))[2]))][1] ) else f) ),L[1])) ),F)$ cad2s(F):=map(lambda([L], [(U:[], map(lambda([f], (if atom (f) then f elseif op(f)="=" then (sol:solve(subst(U,[first(rhs(f))]),lhs(f)),prt2(sol), sol:lhs(f)=(sort( delete([],map(lambda([e], if imagpart(subst(L[2],e))=0 then [e,subst(L[2],e)] else []), map(rhs,sol) )) , lambda ( [ss,tt], ss[2]<tt[2] ) ))[last(rhs(f))][1], U:append (U, [sol] ),sol) elseif op(f)="<" and listp(rhs (f)) then (sol:solve(subst(U,[first(rhs(f))]),lhs(f)),prt2(sol), sol:lhs(f)<(sort( delete([],map(lambda([e], if imagpart(subst(L[2],e))=0 then [e,subst(L[2],e)] else []), map(rhs,sol) )) , lambda ( [ss,tt], ss[2]<tt[2] ) ))[last(rhs(f))][1]) elseif op(f)="<" and listp(lhs (f)) then (sol:solve(subst(U,[first(lhs(f))]),rhs(f)),prt2(sol), sol:(sort( delete([],map(lambda([e], if imagpart(subst(L[2],e))=0 then [e,subst(L[2],e)] else []), map(rhs,sol) )) , lambda ( [ss,tt], ss[2]<tt[2] ) ))[last(lhs(f))][1]<rhs(f)) elseif op(f)="and" then (sol1:solve(subst(U,[first(lhs((args(f))[1]))]),rhs((args(f))[1])), sol1:map(rhs,sol1), sol2:solve(subst(U,[first(rhs((args(f))[2]))]),lhs((args(f))[2])), sol2:map(rhs,sol2), prt2([sol1,sol2]), (sort( delete([],map(lambda([e], if imagpart(subst(L[2],e))=0 then [e,subst(L[2],e)] else []), sol1 )) , lambda ( [ss,tt], ss[2]<tt[2] ) ))[last(lhs((args(f))[1]))][1] < rhs((args(f))[1]) and lhs((args(f))[2]) < (sort( delete([],map(lambda([e], if imagpart(subst(L[2],e))=0 then [e,subst(L[2],e)] else []), sol2 )) , lambda ( [ss,tt], ss[2]<tt[2] ) ))[last(rhs((args(f))[2]))][1] ) else f) ),L[1])),L [2] ] ),F)$ cnn(LL0):=block([n,k,LL:subst("="=equal,LL0),LA,LB,LB0,s,t,L0,M], print("connecting cells."), if is(dispsample=true) then (n:length(LL[1][1]), for k:n step -1 thru 1 do ( LL:append(LL,[[append(makelist(0,n),[1]),makelist([0,0],n)]]), L0:LL[1][1],t:delete(L0[k],L0),s:false, LL:delete(0,map(lambda([L], (LA:L[1], if is(t=delete(LA[k],LA)) then (s:(s or LA[k]),LB:L[2],0) else (M:subst(L0[k]=s,L0),L0:LA,t:delete(LA[k],LA),s:LA[k],LB0:LB,LB:L[2],append([M],[LB0])) ) ),LL)) ), thisout0:LL, thisout:map(lambda([L],[map(lambda([e],if atom(e) then e elseif op(e)="or" then cnn0(e) else e),L[1]),L[2]] ),LL) ) else (n:length(LL[1]), for k:n step -1 thru 1 do ( LL:append(LL,[append(makelist(0,n),[1])]), L0:LL[1],t:delete(L0[k],L0),s:false, LL:delete(0,map(lambda([L], if is(t=delete(L[k],L)) then (s:(s or L[k]),0) else (M:subst(L0[k]=s,L0),L0:L,t:delete(L[k],L),s:L[k],M) ),LL)) ), thisout0:LL, thisout:map(lambda([L],map(lambda([e],if atom(e) then e elseif op(e)="or" then cnn0(e) else e),L) ),LL) ), thisout:subst([equal="="],thisout) )$ distb(f):=block([], if atom(f) then f elseif member(op(f),["or","and"]) then map(distb,f) elseif member(op(f),[equal,"<",">","<=",">="]) then apply(op(f),[{lhs(f)},{rhs(f)}]) else f )$ cnn0(f):=block( [f1,S:[],v,Sb,T,f01], if atom(f) then f else ( map(lambda([e], (if not(member(lhs(e),S)) then S:append(S,[lhs(e)]), if not(member(rhs(e),S)) then S:append(S,[rhs(e)]))), args(subst("and"="or",f))), v:if is("and"=op((args(f))[1])) then S[2] else S[1], Sb:map("{",delete(v,S)),T:makelist({10^5*k},k,length(Sb)), prt2(S,Sb,T,map("=",Sb,T),f,subst(map("=",Sb,T),distb(f))), f01:ev(subst("{"="(",subst(map("=",Sb,T),distb(f))),eval),prt2(f01), f01:ora(subst(["="=equal],cineqs2(f01))), if atom (f01) then f01 else subst(map("=",append([{v}],T),append([v],delete(v,S))),distb(f01)) ) )$ cnnm(LL):=block([k,L0:LL,L1:cnn(LL),L2,usr_display2d], for k:0 unless is(L0=L1) do (L2:cnn(L1),L0:L1,L1:L2), usr_display2d:display2d, set_display('xml), print(apply(matrix,L1)), display2d:usr_display2d,return(done) )$ infix ("==",50,50)$ nary ("&&",50,50)$ nary ("||",50,50)$ sqrtdispflag:off$ mma(F):=subst (["="="==","and"="&&","or"="||"],ora(subst (["="="==",%i=I],F)))$ sqrtdispflag:on$ ora(f):=block([], if not(listp(f)) then return(f), if atom(f) then return(f) else apply("or",map(lambda([e],apply("and",e)),subst("="=equal,f))) )$ cineqs2(F):=block([lv,v,imagunit,L,M,N,k,K,J], lv:listofvars(F),if length(lv)#1 then return([[is(F)]]) else v:first(lv), L:sort(listify(setify(flatten(scanmap( lambda([e], if atom(e) then e elseif member(op(e),["and","or"]) then substpart("[",e,0) elseif member(op(e),["<",">","<=",">=","=","#",equal]) then map(lambda([t],if member(imagunit,listofvars(t)) then [] else rhs(t)),subst(%i=imagunit,solve(substpart("=",e,0),v))) else e), F)))),"<"), if L=[] then return([[is2(subst(v=0,F))]]), L:sort(append(L,(append([first(L)-1],L)+append(L,[last(L)+1]))/2),"<"), M:map(lambda([e],if member(op(F),["and","or"]) then map(is2,rat(subst(v=e,F))) else is2(rat(subst(v=e,F)))),L), if length(setify(M))=1 then return([[is2(subst(v=0,F))]]), L:append([0,0],L,[0,0]),M:append([true,true],M,[true,true]), K:[],for k:3 thru length(M)-2 do ( if oddp(k) then if part(M,k) then (if not(part(M,k-1)) then K:append(K,[part(L,k-1)<v]) elseif not(part(M,k-2)) then K:append(K,[part(L,k-1)<=v]) else K:K, if not(part(M,k+1)) then K:append(K,[v<part(L,k+1)]) elseif not(part(M,k+2)) then K:append(K,[v<=part(L,k+1)]) else K:K) else K:K else if not(part(M,k-1)) and part(M,k) and not(part(M,k+1)) then K:append(K,[v=part(L,k)]) else K:K), N:[],J:[],for k in K do if rhs(k)=v then J:[k] else (N:append(N,[append(J,[k])]),J:[]), delete([],append(N,[J])) )$
ご冗談でしょう,Maxima さん.
1と2は特別なのでしょうか?
k@k ~ $ rmaxima --init= Maxima 5.39.0 http://maxima.sourceforge.net using Lisp CMU Common Lisp 21b (21B Unicode) Distributed under the GNU Public License. See the file COPYING. Dedicated to the memory of William Schelter. The function bug_report() provides bug reporting information. (%i1) newcontext();facts(); (%o1) context2947 (%o2) [] (%i3) apply(assume,[x > 0,equal(1,sqrt(x))]); (%o3) [x > 0, equal(1, sqrt(x))] (%i4) apply(assume, [equal(-1,-sqrt(x))]); *A2 gc_alloc_new_region failed, nbytes=8. CMUCL has run out of dynamic heap space (512 MB). You can control heap size with the -dynamic-space-size commandline option. Imminent dynamic space overflow has occurred: Only a small amount of dynamic space is available now. Please note that you will be returned to the Top-Level without warning if you run out of space while debugging. Maxima encountered a Lisp error: Heap (dynamic space) overflow Automatically continuing. To enable the Lisp debugger set *debugger-hook* to nil. (%i5) Maxima encountered a Lisp error: Interrupted at #xF7FDB430. Automatically continuing. To enable the Lisp debugger set *debugger-hook* to nil. (%i5) facts(); (%o5) [x > 0, equal(1, sqrt(x))] (%i6) newcontext();facts(); (%o6) context178628 (%o7) [] (%i8) apply(assume,[x > 0,equal(2,sqrt(x))]); (%o8) [x > 0, equal(2, sqrt(x))] (%i9) apply(assume, [equal(-2,-sqrt(x))]); *A2 gc_alloc_new_region failed, nbytes=8. CMUCL has run out of dynamic heap space (512 MB). Returning to top-level. *A2 gc_alloc_new_region failed, nbytes=8. CMUCL has run out of dynamic heap space (512 MB). Returning to top-level. *A2 gc_alloc_new_region failed, nbytes=8. CMUCL has run out of dynamic heap space (512 MB). Returning to top-level. rlwrap: warning: maxima killed by SIGTRAP. rlwrap has not crashed, but for transparency, it will now kill itself (without dumping core)with the same signal Trace/breakpoint trap k@k ~ $ rmaxima --init= Maxima 5.39.0 http://maxima.sourceforge.net using Lisp CMU Common Lisp 21b (21B Unicode) Distributed under the GNU Public License. See the file COPYING. Dedicated to the memory of William Schelter. The function bug_report() provides bug reporting information. (%i1) newcontext();facts(); (%o1) context2947 (%o2) [] (%i3) apply(assume,[x > 0,equal(2,sqrt(x))]); (%o3) [x > 0, equal(2, sqrt(x))] (%i4) apply(assume, [equal(-2,-sqrt(x))]); Maxima encountered a Lisp error: Interrupted at #x1006880B. Automatically continuing. To enable the Lisp debugger set *debugger-hook* to nil. (%i5) facts(); (%o5) [x > 0, equal(2, sqrt(x))] (%i6) newcontext();facts(); (%o6) context101642 (%o7) [] (%i8) apply(assume,[x > 0,equal(3,sqrt(x))]); (%o8) [x > 0, equal(3, sqrt(x))] (%i9) apply(assume, [equal(-3,-sqrt(x))]); (%o9) [redundant] (%i10) facts(); (%o10) [x > 0, equal(3, sqrt(x))] (%i11) newcontext();facts(); (%o11) context101663 (%o12) [] (%i13) apply(assume,[x > 0,equal(4,sqrt(x))]); (%o13) [x > 0, equal(4, sqrt(x))] (%i14) apply(assume, [equal(-4,-sqrt(x))]); (%o14) [redundant]
cineqs4 更新
変更点
・名前が cineqs4g になりました.
・処理の効率を改善.数パーセント高速化しましたが,下記の処理のため結果として遅くなりました(苦).
(%i15) for c:1 thru 100 do (ep(0,0),cineqs4g('( x^3+abs(x+1)^(-1/4)-2>0 )))$ Evaluation took 5.2900 seconds (5.2900 elapsed) using 216.100 MB. (%i16) for c:1 thru 100 do (ep(0,0),cineqs4('( x^3+abs(x+1)^(-1/4)-2>0 )))$ Evaluation took 1.7900 seconds (1.7900 elapsed) using 145.687 MB.
・解集合の境界点の近似値の表示を5桁までにしました.
・より正確な値を知る(というか厳密解を特定する)ための関数 NP() を新設.境界点の近似値,対応する厳密解を根とする整数係数多項式のリストのリストを出力します.
・その多項式を得るために load("grobner") しているので,初回コール時は少しモタツキます.
(%i17) ep(0,0)$cineqs4g('( sqrt(x)+sqrt(3-x)+sqrt(3+x)<22/5 )); Evaluation took 0.0000 seconds (0.0000 elapsed) using 736 bytes. Evaluation took 0.0500 seconds (0.0600 elapsed) using 1.867 MB. (%o18) [[0 <= x,x < 0.96473],[2.9379 < x,x <= 3]] (%i19) NP(); Evaluation took 0.1800 seconds (0.1800 elapsed) using 4.819 MB. (%o19) [[0,x], [0.96473, 9765625*x^4+285875000*x^3+1561490000*x^2-10002150400*x+7930971136], [2.9379, 9765625*x^4+285875000*x^3+1561490000*x^2-10002150400*x+7930971136], [3,3-x]] (%i20) bfloat(realroots( (NP())[2][2] ,1.0e-100)),fpprec:100; Evaluation took 0.1900 seconds (0.1900 elapsed) using 6.588 MB. (%o20) [x = 9.647334915176576305835677154346565758106422109360258068951819981502359285249692172963134559965398747b-1, x = 2.937906436935922003266179147367196899336830443976134830225254296066863763304387097967652650705038142b0]
/* 17.02.07.Tue. 22:19:01 */ showtime:on$ display2d:false$ prtoff():=kill(prt)$ prton():=block([],prt([x]):=apply(print,x),return(done))$ ep(e,p):=( if e=1 then algexact:true else algexact:false, if p=1 then prton() else (ratprint:off,prtoff()), usersimp:simp,simp:off )$ matchdeclare([aa,bb,cc],true)$ kill("implies")$ infix(implies,60,40)$ "implies"(x,y):=((not(x)) or (y))$ Powalgsyssingle2g(f,v):=block([F,eqs,inqs,db,pr0,prk,pr,el], local(Powalgsyssinglelocal2g), Powalgsyssinglelocal2g(t):=block([], if atom(t) then return(t) elseif op(t)=Pow3 then (X:first(t),P:second(t),NU:num(P),DN:denom(P),prt([X,P,NU,DN]), if DN=1 and P>=0 then return(X^P) elseif evenp(DN) then (pslc:pslc+1,prk:concat(pr,pslc), eqs:append(eqs,[prk^DN=X^NU]), inqs:append(inqs,[concat(pr,pslc)>=0]), return(prk)) elseif oddp(DN) then (pslc:pslc+1,prk:concat(pr,pslc), eqs:append(eqs,[prk^DN=X^NU]), return(prk)) else return(X^P) ) else return(t)), prt("Powalgsyssingle2g:"),pslc:0,eqs:[],inqs:[], F:scanmap(Powalgsyssinglelocal2g,f,bottomup),prt(F), if pslc>0 then (F:map(num,factor(map(lambda([p],lhs(p)-rhs(p)), append([F],eqs) ))),prt([F,inqs]), el:[F, append(makelist(concat(pr,k),k,1,pslc),[v])], %PL:append(%PL,[el]), F:apply(algsys,el),prt(F), F:map(lambda([l],subst(l,[append([v],inqs)])),F), F:if is(%rnum_list=[]) then F else map( lambda([l],if listofvars(l)=[] then l else [[]]), F), F:ev(map(lambda([l],[substpart("and",first(l),0)]),F),eval)), prt(F),return(F) )$ kill("Lo","Lc2","Go","Gc","Eq")$ infix(Lo,60,40)$ infix(Lc2,60,40)$ infix(Go,60,40)$ infix(Gc,60,40)$ infix(Eq,60,40)$ (x Lo y):=if member(extd,listofvars([x,y])) then false else float(10^(-5)+x<y)$ (x Lc2 y):=if member(extd,listofvars([x,y])) then false else float(x<y+10^(-5))$ (x Go y):=if member(extd,listofvars([x,y])) then false else float(x>y+10^(-5))$ (x Gc y):=if member(extd,listofvars([x,y])) then false else float(10^(-5)+x>y)$ (x Eq y):=if member(extd,listofvars([x,y])) then false else abs(float(x-y))<10^(-5)$ /* float_approx_equal_tolerance float_approx_equal(,) bfloat_approx_equal(,) */ usersimp:simp$simp:off$ defrule(dev2pow, (aa) / (bb), (aa) * (bb) ^ (-1))$ defrule(sqrt2pow, sqrt((aa)), (aa) ^ (1/2))$ simp:usersimp$ defrule(pow2Pow, (aa) ^ (bb), if rat(aa)=0 then (if rat(bb)>0 then 0 elseif rat(bb)=0 then 1 else Pow3(false0,bb,Pow1)) elseif listofvars(aa)=[] or (denom(rat(bb))=1 and rat(bb)>=0) then (aa) ^ (bb) else Pow3(aa,bb,Pow1))$ defrule(abs2Pow, abs((aa)), Pow3((aa)^2,1/2,Pow1))$ defrule(cineqs4n2N1, (aa) # (bb), Not( equal((aa),(bb)) ))$ defrule(cineqs4n2N2, notequal((aa), (bb)), Not( equal((aa),(bb)) ))$ defrule(defbound, Pow3((aa),(bb),(cc)), (prt("defbound:"), if is(bb<0) then Lnz:append(Lnz,{Not( equal((aa),0) )}) elseif evenp(denom(rat(bb))) then Lnn:append(Lnn,{equal((aa),0)}), Pow3((aa),(bb),(cc))))$ realonly:true$ algexact:true$ /* algepsilon:10^7$ rootsepsilon:1.0e-17$ */ cineqs4g(f):=block([L0,F0,F,lv,v,L,sl,func,M,N,k,K,J,el],local(Pow3), %PL:[], lv:listofvars(f),if length(lv)#1 then (simp:usersimp,return([[(true) and (f)]])) else v:first(lv), F:apply1(f,dev2pow,sqrt2pow,pow2Pow), simp:on, F:subst(["<"="Lo",">"="Go","<="="Lc2",">="="Gc","="="Eq","not"=Not],F), F:apply1(F,abs2Pow,cineqs4n2N1,cineqs4n2N2), prt(F), Lnn:{}, F:if member(Pow0,listofvars(F)) or member(Pow1,listofvars(F)) then scanmap(lambda([e], if atom(e) then e elseif member(op(e),["Lo","Go","Lc2","Gc","Eq",equal]) then (if member(Pow0,listofvars(e)) then false elseif member(Pow1,listofvars(e)) then (Lnz:{},(applyb1(e,defbound)) and (substpart("and",Lnz,0))) else e) else e) ,F,bottomup) else F, if is(listofvars(F)=[]) then (F:ev(subst([Not="not"],F),eval),simp:usersimp, return([[(true) and (F)]])), prt([F,Lnn]), L:sort(delete(true,delete(false,listify(setify(flatten( map( lambda([e], if atom(e) then e elseif member(op(e),["Lo","Go","Lc2","Gc","Eq",equal,notequal]) then (if member(Pow1,listofvars(e)) then (Powalgsyssingle2g(e,v)) else (el:apply(num,[factor(lhs(e)-rhs(e))]),%PL:append(%PL,[[[el],[v]]]), sl:map(rhs,flatten(algsys([el], [v]))), delete(first(append(%rnum_list,[gensym()])),sl))) else e), flatten([subst(["%and"="[","%or"="[","and"="[","or"="[",mand="[",mor="[",Not="[","implies"="["],[F,listify(Lnn)])]) ) ))))),"<"), F:subst([Not="not"],F), prt([L,F]), Pow3(x,y,z):=if (is(x<0) and evenp(denom(y))) or (is(x Eq 0) and is(y<0)) then extd else x^y, func(t):=subst(v=t,F), /*prt(func(gensym("x"))),*/ if L=[] then (simp:usersimp,return([[apply(is,[func(0)])]])), L:sort(append(L,(append([first(L)-1],L)+append(L,[last(L)+1]))/2),"<"), M:map(is,(map(func,rat(L)))),prt([L,M]), if length(setify(M))=1 then (simp:usersimp,return([[apply(is,[func(0)])]])), L:append([0,0],L,[0,0]),M:append([true,true],M,[true,true]), K:[],for k:3 thru length(M)-2 do ( if oddp(k) then if part(M,k) then (if not(part(M,k-1)) then K:append(K,[part(L,k-1)<v]) elseif not(part(M,k-2)) then K:append(K,[part(L,k-1)<=v]) else K:K, if not(part(M,k+1)) then K:append(K,[v<part(L,k+1)]) elseif not(part(M,k+2)) then K:append(K,[v<=part(L,k+1)]) else K:K) else K:K else if not(part(M,k-1)) and part(M,k) and not(part(M,k+1)) then K:append(K,[v=part(L,k)]) else K:K), N:[],J:[],for k in K do if rhs(k)=v then J:[k] else (N:append(N,[append(J,[k])]),J:[]), simp:usersimp, %BL:delete([],append(N,[J])), eval_string(ev(string(%BL),fpprintprec:5)) )$ NP():=block([pl,v,bl,bp], if not ?boundp ('poly_monomial_order) then load("grobner"), v:first(listofvars(%BL)), pl:factor(listify(setify(flatten(map( lambda([l], map(lambda([p],if polynomialp(p,[v]) then p else []),l) ),map(lambda([l],apply(poly_grobner,l)),%PL)))))), pl:map(lambda([p],[p,map(rhs,flatten(algsys([p],[v])))]),pl), prt(pl), bl:delete(v,flatten(subst(["<"="[","<="="[",">"="[",">="="[","="="["],%BL))), [pl,bl]:eval_string(ev(string([pl,bl]),fpprintprec:5)), makelist( flatten(map( lambda([p],if member(bp,second(p)) then [bp,first(p)] else []), pl)), bp,bl) )$
cineqs4 原理
cineqs4 は「真理値に関する中間値定理」により,原始式同士の結合の様子に立ち入ることなく,系を解いています.すなわち,...
一変数半代数系 F(x) と実数 a,b (a,<,>=,<=,=,#}) について,{rel(f(a),0),rel(f(b),0)}={true,false}
例えば,f(a)>0はtrue,f(b)>0はfales(つまり,f(b)<=0はtrue)なので
⇒F(x) のある原始式 rel(f(x),0) (rel∈{>,<,>=,<=,=,#}) について,f(c)=0,a<=c<=b となる実数 c がある
という性質(つまり,連続関数 f の零点に触れなければ F の真理値は変わらないこと)により,系に現れるすべての f のすべての相異なる実数の零点 n (>=0) 個とその間,および,両脇との合計 2n+1 個の代表点についての F の真理値は,各点,各区間における真理値に一致するので,true となった代表点に対する点,区間を表す式の選言が解集合の式となる訳です.
具体例です.
(%i1) ep(1,0)$cineqs4('( ((x-1)*(x^2-3)<=0 and x>0) or x^2-x-1=0 )); Evaluation took 0.0000 seconds (0.0000 elapsed) using 856 bytes. Evaluation took 0.0800 seconds (0.0800 elapsed) using 2.356 MB. (%o2) [[x = -(sqrt(5)-1)/2],[1 <= x,x <= sqrt(3)]]
処理の様子は,コントローラー ep (algExact and Print) の第二引数を 1 にすれば確認できます.
(%i3) ep(1,1)$cineqs4('( ((x-1)*(x^2-3)<=0 and x>0) or x^2-x-1=0 )); Evaluation took 0.0000 seconds (0.0000 elapsed) using 992 bytes. ((x-1)*(x^2-3) Lc2 0) and (x Go 0) or (x^2-x-1 Eq 0) ((x-1)*(x^2-3) Lc2 0) and (x Go 0) or (x^2-x-1 Eq 0) [((x-1)*(x^2-3) Lc2 0) and (x Go 0) or (x^2-x-1 Eq 0),{}] [[-sqrt(3),-(sqrt(5)-1)/2,0,1,(sqrt(5)+1)/2,sqrt(3)], ((x-1)*(x^2-3) Lc2 0) and (x Go 0) or (x^2-x-1 Eq 0)] [[((-2*sqrt(3))-1)/2,-sqrt(3),((-(sqrt(5)-1)/2)-sqrt(3))/2,-(sqrt(5)-1)/2, -(sqrt(5)-1)/4,0,1/2,1,((sqrt(5)+1)/2+1)/2,(sqrt(5)+1)/2, ((sqrt(5)+1)/2+sqrt(3))/2,sqrt(3),(2*sqrt(3)+1)/2], [false,false,false,true,false,false,false,true,true,true,true,true,false]] Evaluation took 0.0300 seconds (0.0400 elapsed) using 2.036 MB. (%o4) [[x = -(sqrt(5)-1)/2],[1 <= x,x <= sqrt(3)]]
この例では多項式なので最初の3行はあまり変化しません(Lc2 は less close,Go は greater open,Eq は equal).
次の [-sqrt(3),-(sqrt(5)-1)/2,0,1,(sqrt(5)+1)/2,sqrt(3)] が (x-1)*(x^2-3),x,x^2-x-1 の零点 6 個であり,最後に表示されたリストが,上で述べた代表点 13 個と対応する ((x-1)*(x^2-3) Lc2 0) and (x Go 0) or (x^2-x-1 Eq 0) の真理値です.
その結果において,左から 4 番目が true なので [x = -(sqrt(5)-1)/2],8から12番目までがtrue なので [1 <= x,x <= sqrt(3)],他は false,ということで,(%o4) が得られます.
cineqs4 使用例
作動確認バージョン.
Maxima 5.39.0 using Lisp CMU Common Lisp 21b (21B Unicode)
高次不等式.
(%i1) ep(1,0)$cineqs4('( (x-1)*(x^2-2)*(x^3-3)*(x^2+4)<=0 )); Evaluation took 0.0000 seconds (0.0000 elapsed) using 856 bytes. Evaluation took 0.0600 seconds (0.0600 elapsed) using 1.337 MB. (%o2) [[-sqrt(2) <= x,x <= 1],[sqrt(2) <= x,x <= 3^(1/3)]]
代数的な表示が出来ない場合は,近似値(内部では小数第6位以下を切り捨てています)で答えます.
(%i3) ep(1,0)$cineqs4('( (x-1)*(x^2-2)*(x^3-3)*(x^2+4)<=1/10 )); Evaluation took 0.0000 seconds (0.0100 elapsed) using 720 bytes. Evaluation took 0.0100 seconds (0.0100 elapsed) using 1.091 MB. (%o4) [[-1.414631932468004 <= x,x <= 1.010331702011963], [1.374072765807135 <= x,x <= 1.473559962228517]]
分数関数で表された不等式.
(%i5) ep(1,0)$cineqs4('( (x-1)/x<2 )); Evaluation took 0.0000 seconds (0.0000 elapsed) using 720 bytes. Evaluation took 0.0000 seconds (0.0100 elapsed) using 520.242 KB. (%o6) [[x < -1],[0 < x]]
連言.
(%i7) ep(1,0)$cineqs4('( (x-1)*(x^2-2)*(x^3-3)<=0 and (x-1)/x<2 )); Evaluation took 0.0000 seconds (0.0000 elapsed) using 720 bytes. Evaluation took 0.0400 seconds (0.0500 elapsed) using 2.571 MB. (%o8) [[-sqrt(2) <= x,x < -1],[0 < x,x <= 1],[sqrt(2) <= x,x <= 3^(1/3)]]
選言.
(%i9) ep(1,0)$cineqs4('( (x-1)*(x^2-2)*(x^3-3)<=0 or (x-1)/x<2 )); Evaluation took 0.0000 seconds (0.0000 elapsed) using 720 bytes. Evaluation took 0.0300 seconds (0.0300 elapsed) using 2.804 MB. (%o10) [[true]]
こんな場合は...
(%i11) ep(1,0)$cineqs4('( x^4-x-1<=0 )); Evaluation took 0.0000 seconds (0.0000 elapsed) using 720 bytes. Evaluation took 0.2500 seconds (0.2500 elapsed) using 24.248 MB. (%o12) [[-(sqrt(sqrt(6)*sqrt((sqrt(849)+9)^(1/3)*(3*sqrt(849)+155)^(1/3) *(2^(8/3)*18^(2/3)*sqrt(849) -9*2^(8/3)*18^(2/3)) +(sqrt(849)+9)^(1/3)*(3*sqrt(849)+155)^(2/3) *(2^(1/3)*18^(2/3)*sqrt(849) -9*2^(1/3)*18^(2/3)) +(sqrt(849)+9)^(1/3) *(32*18^(2/3)*sqrt(849)-16*18^(5/3))) +3*2^(11/3)*(3^(3/2)*sqrt(283)-27)^(1/3) -8*18^(2/3)*(sqrt(3)*sqrt(283)+9)^(1/3)) -sqrt(6)*sqrt(2^(1/3)*(3*sqrt(849)-27)^(1/3)*(3*sqrt(849)+155)^(1/3) -2^(8/3)*(3*sqrt(849)-27)^(1/3))) /24 <= x, x <= (sqrt(sqrt(6)*sqrt((sqrt(849)+9)^(1/3)*(3*sqrt(849)+155)^(1/3) *(2^(8/3)*18^(2/3) *sqrt(849) -9*2^(8/3)*18^(2/3)) +(sqrt(849)+9)^(1/3)*(3*sqrt(849)+155)^(2/3) *(2^(1/3)*18^(2/3) *sqrt(849) -9*2^(1/3)*18^(2/3)) +(sqrt(849)+9)^(1/3) *(32*18^(2/3)*sqrt(849)-16*18^(5/3))) +3*2^(11/3)*(3^(3/2)*sqrt(283)-27)^(1/3) -8*18^(2/3)*(sqrt(3)*sqrt(283)+9)^(1/3)) +sqrt(6)*sqrt(2^(1/3)*(3*sqrt(849)-27)^(1/3) *(3*sqrt(849)+155)^(1/3) -2^(8/3)*(3*sqrt(849)-27)^(1/3))) /24]]
コントローラー ep の第一引数を 0 にすると,近似値で答えます.
(%i13) ep(0,0)$cineqs4('( x^4-x-1<=0 )); Evaluation took 0.0000 seconds (0.0000 elapsed) using 736 bytes. Evaluation took 0.0000 seconds (0.0100 elapsed) using 175.711 KB. (%o14) [[-0.7244919786096257 <= x,x <= 1.220744081172491]]
無理関数で表された不等式.
(%i15) ep(1,0)$cineqs4('( sqrt(x)<x-1 )); Evaluation took 0.0000 seconds (0.0000 elapsed) using 720 bytes. Evaluation took 0.0100 seconds (0.0100 elapsed) using 674.156 KB. (%o16) [[(sqrt(5)+3)/2 < x]]
含意.
(%i17) ep(1,0)$cineqs4('( 1<sqrt(x) implies 1<x )); Evaluation took 0.0000 seconds (0.0000 elapsed) using 720 bytes. Evaluation took 0.0100 seconds (0.0100 elapsed) using 428.797 KB. (%o18) [[true]]
根号が複数あっても.
(%i19) ep(1,0)$cineqs4('( sqrt(x)+sqrt(3-x)<2 )); Evaluation took 0.0000 seconds (0.0000 elapsed) using 720 bytes. Evaluation took 0.0300 seconds (0.0200 elapsed) using 1.524 MB. (%o20) [[0 <= x,x < -(2^(3/2)-3)/2],[(2^(3/2)+3)/2 < x,x <= 3]] (%i21) ep(0,0)$cineqs4('( sqrt(x)+sqrt(3-x)+sqrt(3+x)<22/5 )); Evaluation took 0.0000 seconds (0.0000 elapsed) using 736 bytes. Evaluation took 0.0200 seconds (0.0200 elapsed) using 2.500 MB. (%o22) [[0 <= x,x < 0.9647335423197492],[2.937906564163217 < x,x <= 3]]
ネストされても.
(%i23) ep(1,0)$cineqs4('( 1/sqrt(sqrt(1+x)-x)<2 )); Evaluation took 0.0000 seconds (0.0000 elapsed) using 720 bytes. Evaluation took 0.0700 seconds (0.0800 elapsed) using 4.835 MB. (%o24) [[-1 <= x,x < 5/4]]
有理数冪でも.
(%i25) ep(1,0)$cineqs4('( x^4-x^(1/3)-1>0 )); Evaluation took 0.0000 seconds (0.0000 elapsed) using 720 bytes. Evaluation took 0.0300 seconds (0.0400 elapsed) using 1.039 MB. (%o26) [[x < -0.6196702671972711],[1.198342827550492 < x]]
一変数有理冪半代数系ソルバー cineqs4
cineqs2 を分数関数や有理冪乗(fractional power)関数で表された方程式,不等式まで扱えるように拡張しました.
コードの大半は,Maxima の非論理的挙動を封じるためのものです(笑).
作動確認は Maxima 5.39.0 on Ubuntu 及び windows で行っています.
/* 17.02.05.Sun. 12:17:08 */ showtime:on$ display2d:false$ prtoff():=kill(prt)$ prton():=block([],prt([x]):=apply(print,x),return(done))$ ep(e,p):=( if e=1 then algexact:true else algexact:false, if p=1 then prton() else (ratprint:off,prtoff()), usersimp:simp,simp:off )$ matchdeclare([aa,bb,cc],true)$ kill("implies")$ infix(implies,60,40)$ "implies"(x,y):=((not(x)) or (y))$ powalgsyssingle2(f,v):=block([powalgsyssinglelocal2,F,eqs,inqs,db,paral,pr0,prk,pr], powalgsyssinglelocal2(t):=block([], if atom(t) then return(t) elseif op(t)=Pow3 then (X:first(t),P:second(t),NU:num(P),DN:denom(P),prt([X,P,NU,DN]), if is(DN=1) and is(P>=0) then return(X^P) elseif evenp(DN) then (pslc:pslc+1,prk:concat(pr,pslc), eqs:append(eqs,[prk^DN=X^NU]), inqs:append(inqs,[concat(pr,pslc)>=0]), db:append(db,map(rhs,flatten(algsys([X],[v])))), paral:append(paral,%rnum_list), return(prk)) elseif oddp(DN) then (pslc:pslc+1,prk:concat(pr,pslc), eqs:append(eqs,[prk^DN=X^NU]), db:append(db,map(rhs,flatten(algsys([X],[v])))), return(prk)) else return(X^P) ) else return(t)), prt("powalgsyssingle2:"),pslc:0,eqs:[],inqs:[],db:[],paral:[], F:scanmap(powalgsyssinglelocal2,apply1(f,cineqs4r1),bottomup),prt(F), if is(pslc>0) then (F:map(num,factor(map(lambda([p],lhs(p)-rhs(p)), append([F],eqs) ))),prt([F,inqs,db]), F:algsys(F, append([v],makelist(concat(pr,k),k,1,pslc)) ),prt(F), F:map(lambda([l],subst(l,[append([v],inqs),db])),F), F:if is(%rnum_list=[]) then F else subst( map( lambda([e],e=false), %rnum_list) ,F), F:ev(map(lambda([l],[substpart("and",first(l),0),last(l)]),F),eval)), prt(F),return(F) )$ kill("Lo","Lc2","Go","Gc","Eq")$ infix(Lo,60,40)$ infix(Lc2,60,40)$ infix(Go,60,40)$ infix(Gc,60,40)$ infix(Eq,60,40)$ (x Lo y):=if member(extd,listofvars([x,y])) then false else float(10^(-5)+x<y)$ (x Lc2 y):=if member(extd,listofvars([x,y])) then false else float(x<y+10^(-5))$ (x Go y):=if member(extd,listofvars([x,y])) then false else float(x>y+10^(-5))$ (x Gc y):=if member(extd,listofvars([x,y])) then false else float(10^(-5)+x>y)$ (x Eq y):=if member(extd,listofvars([x,y])) then false else abs(float(x-y))<10^(-5)$ /* float_approx_equal_tolerance float_approx_equal(,) bfloat_approx_equal(,) */ defrule(cineqs4r1, (aa) ^ (bb), if is(listofvars(aa)=[]) or ( is(denom(bb)=1) and is(bb>=0) ) then (aa) ^ (bb) else Pow3(aa,bb,Pow0))$ defrule(cineqs4r2, (aa) / (bb), if is(listofvars(bb)=[]) then (aa) / (bb) else (aa) * Pow3(bb,-1,Pow0))$ defrule(cineqs4r3, (aa) # (bb), Not( equal((aa),(bb)) ))$ defrule(cineqs4r4, notequal((aa), (bb)), Not( equal((aa),(bb)) ))$ defrule(cineqs4r5, Div((aa), 0), Div((aa),0,Div0))$ defrule(cineqs4r6, Div((aa), (bb)), Div((aa),(bb),Div1))$ defrule(cineqs4r7, abs((aa)), Pow3((aa)^2,1/2,Pow0))$ defrule(defbound, Pow3((aa),(bb),(cc)), (L0:append(L0,{equal((aa),0)}),Pow3((aa),(bb),(cc))))$ defrule(defbound2, Div((aa),(bb),(cc)), (prt("defbound2:"),L0:append(L0,{Not( equal((bb),0)) }),(aa)/(bb)))$ defrule(defbound3, Div((aa),(bb)), ((aa)/(bb)))$ realonly:true$ algexact:true$ algepsilon:10^17$ rootsepsilon:1.0e-17$ cineqs4(f):=block([L0,F0,F,lv,v,L,sl,func,M,N,k,K,J],local(Pow3), lv:listofvars(f), if length(lv)#1 then (simp:usersimp,return([[(true) and (f)]])) else v:first(lv), F:subst(["<"="Lo",">"="Go","<="="Lc2",">="="Gc", "="="Eq","not"=Not,"/"=Div,0^0=1],f), F:apply1(F,cineqs4r3,cineqs4r4,cineqs4r5,cineqs4r6,cineqs4r7), prt(F), F0:if member(Div0,listofvars(F)) then scanmap(lambda([e], if atom(e) then e elseif member(op(e),["Lo","Go","Lc2","Gc","Eq",equal]) and member(Div0,listofvars(e)) then false else e) ,F,bottomup) else F, if is(listofvars(F0)=[]) then (F0:ev(subst([Not="not"],F0),eval),simp:usersimp, return([[(true) and (F0)]])), F1:if member(Div1,listofvars(F)) then scanmap(lambda([e], if atom(e) then e elseif member(op(e),["Lo","Go","Lc2","Gc","Eq",equal]) and member(Div1,listofvars(e)) then (L0:{},(applyb1(e,defbound2)) and (substpart("and",L0,0))) else e) ,F0,bottomup) else F0, F:F1, prt(F), simp:on,L0:{}, F:apply1(F,cineqs4r1,defbound), prt([F,L0]), L:sort(delete(false,listify(setify(flatten(scanmap( lambda([e], if atom(e) then e elseif member(op(e),["Lo","Go","Lc2","Gc","Eq",equal,notequal]) then (if member(Pow0,listofvars(e)) then (powalgsyssingle2(e,v)) else (sl:map(rhs,flatten(algsys( [apply(num,[factor(lhs(e)-rhs(e))])], [v]))), delete(first(append(%rnum_list,[gensym()])),sl)) ) elseif member(op(e),["%and","%or","and","or",mand,mor,Not,"implies"]) then substpart("[",e,0) else e), mand(F, substpart(mand,L0,0)) ))))),"<"), F:subst([Not="not"],F), F:apply1(F,defbound3), prt([L,F]), Pow3(x,y,z):=if (is(x<0) and evenp(denom(y))) or (is(x Eq 0) and is(y<0)) then extd else x^y, func(t):=subst(v=t,F), /*prt(func(gensym("x"))),*/ if L=[] then (simp:usersimp,return([[apply(is,[func(0)])]])), L:sort(append(L,(append([first(L)-1],L)+append(L,[last(L)+1]))/2),"<"), M:map(is,ratsimp(map(func,L))),prt([L,M]), /* M:map(is,(map(func,L))),prt([L,M]), fpprec:64, M:map(lambda([e],is(rationalize(bfloat(subst([v=e],F))))),L),prt([L,M]), fpprec:16, M:map(lambda([e],is(subst(v=float(e),F))),L),prt([L,M]), */ if length(setify(M))=1 then (simp:usersimp,return([[apply(is,[func(0)])]])), L:append([0,0],L,[0,0]),M:append([true,true],M,[true,true]), K:[],for k:3 thru length(M)-2 do ( if oddp(k) then if part(M,k) then (if not(part(M,k-1)) then K:append(K,[part(L,k-1)<v]) elseif not(part(M,k-2)) then K:append(K,[part(L,k-1)<=v]) else K:K, if not(part(M,k+1)) then K:append(K,[v<part(L,k+1)]) elseif not(part(M,k+2)) then K:append(K,[v<=part(L,k+1)]) else K:K) else K:K else if not(part(M,k-1)) and part(M,k) and not(part(M,k+1)) then K:append(K,[v=part(L,k)]) else K:K), N:[],J:[],for k in K do if rhs(k)=v then J:[k] else (N:append(N,[append(J,[k])]),J:[]), simp:usersimp, delete([],append(N,[J])) )$