本日のC.A.D.
今回は，第一種 Chebyshev 多項式 $T_4(x)=64x^7-112x^5+56x^3-7x$ に関する
$\exists a\ \exists b\ \exists c\ \forall x\ (\, -1\le x\le1 \to 64 x^7+a x^5+b x^3+c x\le m\,)\quad \cdots\quad (*)$
が $1\le m $ と等価であるという性質です．
まず
$\forall x\ (\, -1\le x\le1 \to 64 x^7+a x^5+b x^3+c x\le 1\,)$
が $a=-112\ \land\ b=56\ \land\ c=-7$ と等価であることは
? tst12([all],[c,a,b,x],imp,"x^2<=1,64*x^7+a*x^5+b*x^3+c*x<=1");Ans(1); *** using Lazard's method (MPP17). [x,3] [b,3] [a,6] [c,25] time = 14,301 ms. 151 4461(14,161) 55129(267,4271) 11(36373,18580) *** combined adjacent 10 cells. 1[c = [c+7,1],a = [a+112,1],b = [b-56,1],true] time = 2min, 31,786 ms.
のように判るので，$1\le m $ ならば $(*)$ を得ます．そして
$\forall x\ (\, -1\le x\le1 \to 64 x^7+a x^5+b x^3+c x< 1\,)$
が false と等価であることも
? tst12([all],[c,a,b,x],imp,"x^2<=1,64*x^7+a*x^5+b*x^3+c*x<1");Ans(); *** using Lazard's method (MPP17). [x,3] [b,3] [a,6] [c,25] time = 14,440 ms. 151 4461(14,161) 55129(267,4271) 0(36373,18580) 1[false] time = 2min, 30,091 ms.
のように判るので，$(*)$ ならば $1\le m $ を得ます．また，前半の結果を用いるなら，後半は
? tst12([ex],[v,x],and,"x^2<=1,64*x^7-112*x^5+56*x^3-7*x==v");Ans(); *** using Lazard's method (MPP17). [x,3] [v,2] time = 42 ms. 5 3(12,4) *** combined adjacent 2 cells. 1[[v+1,1] <= v <= [v-1,1],true] time = 21 ms.
のように $1$ が値域に属することから容易に得られます．
Wolfram Language 12.3.0 Engine for Linux x86 (64-bit)
Copyright 1988-2021 Wolfram Research, Inc.
In[1]:= Resolve[ForAll[x,x^2<=1,64*x^7+a*x^5+b*x^3+c*x<=1],Reals,WorkingPrecision->50]//Timing
Resolve::cadpr: The cylindrical algebraic decomposition algorithm used by Resolve
failed due to a too low WorkingPrecision. Increasing the value of WorkingPrecision may allow the algorithm
to succeed.
Out[1]= {22.957144, ForAll[{x}, x^2 <= 1, c*x + b*x^3 + a*x^5 + 64*x^7 <= 1]}
In[2]:= Resolve[ForAll[x,x^2<=1,64*x^7+a*x^5+b*x^3+c*x<=1],Reals,WorkingPrecision->100]//Timing
Out[2]= {77.871964, a == -112 && b == 56 && c == -7}
本日のC.A.D.
A.W.Strzebonski "Cylindrical Algebraic Decomposition using validated numerics" の 7.2. Randomly generated systems で使用された入力です．元の code では CylindricalDecomposition ですが，変数順序を heuristic に委ねるため，無指定の Reduce で実行しました．
Wolfram Language 12.3.0 Engine for Linux x86 (64-bit)
Copyright 1988-2021 Wolfram Research, Inc.
In[1]:= rana={-363 + 177*x1 - 125*x1*x2^2 + 444*x2^3 - 766*x2*x3 + 477*x1*x2*x3 + 447*x3^2 == 0 && 864 - 55*x1^3*
x2 - 491*x2^2 + 959*x1*x2^3 - 465*x1*x3^2 > 0,566 + 180*x1^2*x2^2 + 264*x1^2*x3 + 919*x2^2*x3 + 941*x2^3*x3 + 481
*x1*x2*x3^2 - 519*x2*x3^3 >= 0,634 - 508*x3 + 896*x3^3 - 342*x1*x2*x4 - 462*x2*x3*x4 + 144*x2*x4^2 >= 0,-120 + 38
*x2^3*x4 + 207*x3*x4 + 872*x3^2*x4 + 165*x2*x3^2*x4 + 672*x2*x3*x4^2 + 87*x3^2*x4^2 == 0,-394 + 902*x2 - 282*x1*x
2 + 347*x3 + 389*x2*x3 - 209*x1*x2*x3 + 987*x2*x3^2 == 0 && -202 + 292*x2^3 + 326*x2^2*x3 + 728*x2*x3^2 - 80*x3^3
> 0 && -866 - 1024*x1^4 - 886*x2*x3 - 736*x1*x2^2*x3 + 763*x2*x3^3 >= 0};
In[2]:= ranb={446 + 511*x3 - 190*x1^3*x3 + 132*x2^2*x3 + 968*x1^2*x2*x4 + 226*x1*x2*x3*x4 + 645*x2*x3*x4^2 > 0,-4
47 - 312*x3 + 623*x3^2 - 912*x2*x4 + 555*x4^2 == 0 && -730 - 978*x2*x3 + 557*x2*x3^2 + 487*x1^2*x4 - 701*x2*x3*x4
+ 470*x2*x4^2 + 297*x3*x4^2 >= 0,-894 - 884*x2*x3 - 290*x3^2 + 144*x1*x3^2 - 887*x2*x3^2 + 20*x3^3 == 0 && -100
+ 571*x1*x3 - 314*x2^3*x3 + 398*x1*x3^3 - 522*x2*x3^3 >= 0,505 + 834*x2^3*x3 + 887*x2*x3^2 + 151*x3^3 - 243*x2*x3
^3 == 0 && -323 - 408*x1^2*x2 + 484*x1*x2^2 + 556*x3 + 705*x1*x2*x3 > 0 && 422 - 225*x1^4 + 903*x1^2*x3^2 - 944*x
2^2*x3^2 + 996*x3^4 >= 0,-263 - 172*x2 + 78*x2^2 + 41*x2*x3 + 864*x4 + 623*x2*x4 + 293*x3*x4 > 0 && 781 + 948*x1^
2 - x3 - 946*x1^3*x4 - 606*x1*x3*x4 >= 0};
In[3]:= ranc={-417 + 159*x1^3 + 210*x1^3*x2 + 358*x1*x2^2*x3 + 605*x2*x3^3 == 0 && -1016 + 1010*x2 + 516*x1*x2 -
960*x3 + 571*x1*x3 - x3^2 > 0 && -4 + 627*x1^2 - 100*x1*x2 - 670*x3 + 670*x1*x3 >= 0,-502 - 1013*x1*x2*x3 + 328*x
2^2*x3 + 330*x1*x3^2 + 884*x3*x4 + 502*x1*x3*x4 == 0 && 313 - 636*x2 - 997*x1*x3 - 573*x3^2 + 591*x3*x4 >= 0,790
+ 446*x2^3 + 189*x1^2*x3 + 399*x1*x2*x3 - 110*x3^2 - 19*x2*x3^2 > 0 && 60 - 531*x1^2 - 691*x1*x2 + 908*x1^2*x3 -
935*x1*x2*x3 + 431*x2^2*x3 - 747*x3^3 >= 0,275 + 130*x3^2 - 279*x2*x4 - 113*x3*x4 + 620*x4^2 > 0 && 542 + 680*x1*
x3 + 75*x2^2*x3 + 550*x1*x3^2 + 621*x1*x3*x4 >= 0,-218 + 955*x1^2*x2*x3 + 520*x2*x3^3 + 1017*x3^4 + 250*x1*x2*x4
- 391*x2^3*x4 - 525*x1*x3*x4^2 > 0};
In[4]:= First@Timing@Reduce[#,Reals]&/@rana//Timing
Out[4]= {1.426756, {0.467609, 0.346341, 0.038994, 0.43843, 0.134579}}
In[5]:= First@Timing@Reduce[#,Reals]&/@ranb//Timing
Out[5]= {7.619839, {0.886294, 0.062116, 0.02397, 4.808107, 1.835422}}
In[6]:= First@Timing@Reduce[#,Reals]&/@ranc//Timing
Out[6]= {17.547247, {3.909086, 3.40725, 4.59773, 0.081719, 5.538123}}同じ入力に対する tst12 の実行結果です．
GP/PARI CALCULATOR Version 2.14.0 (development git-b0a4238356)
amd64 running linux (x86-64/GMP-6.2.1 kernel) 64-bit version
compiled: Feb 21 2021, gcc version 9.3.0 (Ubuntu 9.3.0-17ubuntu1~20.04)
threading engine: single
(readline v8.0 enabled, extended help enabled)
Copyright (C) 2000-2021 The PARI Group
PARI/GP is free software, covered by the GNU General Public License, and comes WITHOUT ANY WARRANTY WHATSOEVER.
Type ? for help, \q to quit.
Type ?17 for how to get moral (and possibly technical) support.
parisizemax = 20000002048, primelimit = 500000
? wt=getwalltime();rana7(7*23);print("total = ",strtime(getwalltime()-wt),".");print();
[[[x1,3,6],[x2,3,7],[x3,2,4]]]
*** reordering by degree.
*** using the sum of squares projection.
[x3,2]
[x1,4]
[x2,11]
time = 261 ms.
45 743(16,30) 704(324,304)
*** combined adjacent 173 cells.
[[[x1,2,3],[x2,3,5],[x3,3,5]]]
*** reordering by degree.
*** using the sum of squares projection.
[x1,1]
[x3,2]
[x2,3]
time = 66 ms.
13 109(8,2) 277(70,0)
*** combined adjacent 184 cells.
[[[x1,1,1],[x2,1,3],[x3,3,3],[x4,2,3]]]
*** reordering by degree.
*** using the sum of squares projection.
[x1,1]
[x2,2]
[x4,2]
[x3,2]
time = 115 ms.
5 23(0,2) 69(10,10) 98(0,0)
*** combined adjacent 44 cells.
[[[x1,0,0],[x2,3,3],[x3,2,5],[x4,2,6]]]
*** reordering by degree.
*** using the sum of squares projection.
[x1]
[x2,1]
[x3,1]
[x4,3]
time = 74 ms.
11 47(8,0) 78(38,0) 78(0,0)
*** combined adjacent 13 cells.
[[[x1,4,4],[x2,3,11],[x3,3,10]]]
*** reordering by degree.
*** using the sum of squares projection.
[x1,2]
[x3,5]
[x2,10]
time = 1,151 ms.
33 143(16,23) 4(58,6)
*** combined adjacent 2 cells.
total = 3,205 ms.
? wt=getwalltime();rana7(7*29);print("total = ",strtime(getwalltime()-wt),".");print();
[[[x1,3,4,6],[x2,3,4,7],[x3,2,3,4]]]
*** reordering by d+sotd+o.
*** using the sum of squares projection.
[x3,2]
[x1,4]
[x2,11]
time = 230 ms.
45 743(16,30) 704(324,304)
*** combined adjacent 173 cells.
[[[x1,2,4,3],[x2,3,4,5],[x3,3,4,5]]]
*** reordering by d+sotd+o.
*** using the sum of squares projection.
[x1,1]
[x3,2]
[x2,3]
time = 66 ms.
13 109(8,2) 277(70,0)
*** combined adjacent 184 cells.
[[[x1,1,3,1],[x2,1,3,3],[x3,3,3,3],[x4,2,3,3]]]
*** reordering by d+sotd+o.
*** using the sum of squares projection.
[x1,1]
[x2,2]
[x4,2]
[x3,2]
time = 107 ms.
5 23(0,2) 69(10,10) 98(0,0)
*** combined adjacent 44 cells.
[[[x1,0,0,0],[x2,3,4,3],[x3,2,4,5],[x4,2,4,6]]]
*** reordering by d+sotd+o.
*** using the sum of squares projection.
[x1]
[x2,1]
[x3,1]
[x4,3]
time = 78 ms.
11 47(8,0) 78(38,0) 78(0,0)
*** combined adjacent 13 cells.
[[[x1,4,4,4],[x2,3,4,11],[x3,3,4,10]]]
*** reordering by d+sotd+o.
*** using the sum of squares projection.
[x1,2]
[x3,5]
[x2,10]
time = 1,169 ms.
33 143(16,23) 4(58,6)
*** combined adjacent 2 cells.
total = 3,205 ms.
? wt=getwalltime();rana7(7*31);print("total = ",strtime(getwalltime()-wt),".");print();
[[[x1,8,13,18,6],[x2,13,17,21,7],[x3,6,10,10,4]]]
*** reordering by sum(d+td in terms)+o.
*** using the sum of squares projection.
[x3,2]
[x1,4]
[x2,11]
time = 222 ms.
45 743(16,30) 704(324,304)
*** combined adjacent 173 cells.
[[[x1,5,8,11,3],[x2,9,11,19,5],[x3,8,11,18,5]]]
*** reordering by sum(d+td in terms)+o.
*** using the sum of squares projection.
[x1,1]
[x3,2]
[x2,3]
time = 66 ms.
13 109(8,2) 277(70,0)
*** combined adjacent 184 cells.
[[[x1,1,2,3,1],[x2,3,2,9,3],[x3,5,6,7,3],[x4,4,5,9,3]]]
*** reordering by sum(d+td in terms)+o.
*** using the sum of squares projection.
[x1,1]
[x2,2]
[x4,2]
[x3,2]
time = 108 ms.
5 23(0,2) 69(10,10) 98(0,0)
*** combined adjacent 44 cells.
[[[x1,0,0,0,0],[x2,5,11,12,3],[x3,8,8,17,5],[x4,8,8,21,6]]]
*** reordering by sum(d+td in terms)+o.
*** using the sum of squares projection.
[x1]
[x2,1]
[x3,1]
[x4,3]
time = 72 ms.
11 47(8,0) 78(38,0) 78(0,0)
*** combined adjacent 13 cells.
[[[x1,7,14,13,4],[x2,15,16,30,11],[x3,16,22,28,10]]]
*** reordering by sum(d+td in terms)+o.
*** using the sum of squares projection.
[x1,2]
[x2,5]
[x3,6]
time = 767 ms.
15 55(12,11) 6(26,6)
*** combined adjacent 1 cells.
total = 2,406 ms.
? wt=getwalltime();ranb7(7*23);print("total = ",strtime(getwalltime()-wt),".");print();
[[[x1,3,3],[x2,2,4],[x3,1,5],[x4,2,3]]]
*** using the sum of squares projection.
[x4,1]
[x3,3]
[x2,5]
[x1,8]
time = 174 ms.
21 175(12,15) 915(102,86) 1192(390,0)
*** combined adjacent 44 cells.
[[[x1,2,1],[x2,1,5],[x3,2,6],[x4,2,6]]]
*** reordering by degree.
*** using the sum of squares projection.
[x1,1]
[x2,2]
[x4,5]
[x3,8]
time = 145 ms.
25 195(16,18) 180(0,60) 270(64,0)
*** combined adjacent 178 cells.
[[[x1,1,3],[x2,3,4],[x3,3,9]]]
*** reordering by degree.
*** using the sum of squares projection.
[x1,2]
[x2,2]
[x3,3]
time = 134 ms.
7 31(4,0) 21(0,24)
*** combined adjacent 6 cells.
[[[x1,4,5],[x2,3,7],[x3,4,9]]]
*** reordering by degree.
*** using the sum of squares projection.
[x1,2]
[x2,6]
[x3,16]
time = 1,358 ms.
95 152(16,87) 251(14,8)
*** combined adjacent 151 cells.
[[[x1,3,3],[x2,2,4],[x3,1,4],[x4,1,5]]]
*** reordering by degree.
*** using the sum of squares projection.
[x3,2]
[x4,5]
[x2,8]
[x1,21]
time = 311 ms.
75 1681(14,80) 13707(546,668) 13910(0,5540)
*** combined adjacent 7087 cells.
total = 10,369 ms.
? wt=getwalltime();ranb7(7*29);print("total = ",strtime(getwalltime()-wt),".");print();
[[[x1,3,4,3],[x2,2,4,4],[x3,1,4,5],[x4,2,4,3]]]
*** using the sum of squares projection.
[x4,1]
[x3,3]
[x2,5]
[x1,8]
time = 174 ms.
21 175(12,15) 915(102,86) 1192(390,0)
*** combined adjacent 44 cells.
[[[x1,2,3,1],[x2,1,3,5],[x3,2,3,6],[x4,2,3,6]]]
*** reordering by d+sotd+o.
*** using the sum of squares projection.
[x1,1]
[x2,2]
[x4,5]
[x3,8]
time = 145 ms.
25 195(16,18) 180(0,60) 270(64,0)
*** combined adjacent 178 cells.
[[[x1,1,4,3],[x2,3,4,4],[x3,3,4,9]]]
*** reordering by d+sotd+o.
*** using the sum of squares projection.
[x1,2]
[x2,2]
[x3,3]
time = 108 ms.
7 31(4,0) 21(0,24)
*** combined adjacent 6 cells.
[[[x1,4,4,5],[x2,3,4,7],[x3,4,4,9]]]
*** reordering by d+sotd+o.
*** using the sum of squares projection.
[x1,2]
[x2,6]
[x3,16]
time = 1,371 ms.
95 152(16,87) 251(14,8)
*** combined adjacent 151 cells.
[[[x1,3,4,3],[x2,2,2,4],[x3,1,3,4],[x4,1,4,5]]]
*** reordering by d+sotd+o.
*** using the sum of squares projection.
[x3,2]
[x2,3]
[x4,6]
[x1,12]
time = 216 ms.
49 709(4,26) 4191(300,386) 5002(0,2064)
*** combined adjacent 2701 cells.
total = 6,579 ms.
? wt=getwalltime();ranb7(7*31);print("total = ",strtime(getwalltime()-wt),".");print();
[[[x1,6,11,12,3],[x2,5,7,15,4],[x3,5,3,16,5],[x4,4,8,12,3]]]
*** using the sum of squares projection.
[x4,1]
[x3,3]
[x2,5]
[x1,8]
time = 174 ms.
21 175(12,15) 915(102,86) 1192(390,0)
*** combined adjacent 44 cells.
[[[x1,2,5,3,1],[x2,5,3,13,5],[x3,8,7,14,6],[x4,9,7,16,6]]]
*** reordering by sum(d+td in terms)+o.
*** using the sum of squares projection.
[x1,1]
[x2,2]
[x3,4]
[x4,6]
time = 158 ms.
19 139(20,14) 136(0,52) 218(56,0)
*** combined adjacent 124 cells.
[[[x1,3,5,9,3],[x2,6,13,13,4],[x3,18,15,27,9]]]
*** reordering by sum(d+td in terms)+o.
*** using the sum of squares projection.
[x1,2]
[x2,2]
[x3,3]
time = 113 ms.
7 31(4,0) 21(0,24)
*** combined adjacent 6 cells.
[[[x1,10,17,17,5],[x2,11,24,24,7],[x3,19,25,30,9]]]
*** reordering by sum(d+td in terms)+o.
*** using the sum of squares projection.
[x1,2]
[x2,6]
[x3,16]
time = 1,371 ms.
95 152(16,87) 251(14,8)
*** combined adjacent 151 cells.
[[[x1,6,8,9,3],[x2,5,2,7,4],[x3,4,3,8,4],[x4,5,4,12,5]]]
*** reordering by sum(d+td in terms)+o.
*** using the sum of squares projection.
[x2,1]
[x3,2]
[x4,4]
[x1,6]
time = 123 ms.
23 199(4,8) 556(92,48) 1058(228,0)
*** combined adjacent 11 cells.
total = 4,510 ms.
? wt=getwalltime();ranc7(7*23);print("total = ",strtime(getwalltime()-wt),".");print();
[[[x1,3,8],[x2,2,6],[x3,3,7]]]
*** reordering by degree.
*** using the sum of squares projection.
[x2,3]
[x3,7]
[x1,22]
time = 385 ms.
87 1795(42,71) 561(314,1257)
*** combined adjacent 244 cells.
[[[x1,1,4],[x2,2,3],[x3,2,8],[x4,1,3]]]
*** reordering by degree.
*** using the sum of squares projection.
[x4,2]
[x2,3]
[x1,4]
[x3,10]
time = 244 ms.
41 245(22,18) 901(276,140) 667(0,506)
*** combined adjacent 268 cells.
[[[x1,2,6],[x2,3,6],[x3,3,8]]]
*** reordering by degree.
*** using the sum of squares projection.
[x1,2]
[x2,4]
[x3,16]
time = 304 ms.
59 675(30,29) 1162(326,306)
*** combined adjacent 521 cells.
[[[x1,1,3],[x2,2,2],[x3,2,6],[x4,2,4]]]
*** reordering by degree.
*** using the sum of squares projection.
[x2,2]
[x1,2]
[x4,4]
[x3,4]
time = 161 ms.
11 53(4,7) 189(0,0) 317(76,64)
*** combined adjacent 189 cells.
[[[x1,2,3],[x2,3,4],[x3,4,4],[x4,2,3]]]
*** reordering by degree.
*** using the sum of squares projection.
[x4,1]
[x1,3]
[x2,4]
[x3,6]
time = 896 ms.
37 277(21,24) 1677(182,84) 1725(850,0)
*** combined adjacent 70 cells.
total = 9,355 ms.
? wt=getwalltime();ranc7(7*29);print("total = ",strtime(getwalltime()-wt),".");print();
[[[x1,3,4,8],[x2,2,4,6],[x3,3,4,7]]]
*** reordering by d+sotd+o.
*** using the sum of squares projection.
[x2,3]
[x3,7]
[x1,22]
time = 379 ms.
87 1795(42,71) 561(314,1257)
*** combined adjacent 244 cells.
[[[x1,1,3,4],[x2,2,3,3],[x3,2,3,8],[x4,1,3,3]]]
*** reordering by d+sotd+o.
*** using the sum of squares projection.
[x4,2]
[x2,3]
[x1,4]
[x3,10]
time = 217 ms.
41 245(22,18) 901(276,140) 667(0,506)
*** combined adjacent 268 cells.
[[[x1,2,3,6],[x2,3,3,6],[x3,3,3,8]]]
*** reordering by d+sotd+o.
*** using the sum of squares projection.
[x1,2]
[x2,4]
[x3,16]
time = 273 ms.
59 675(30,29) 1162(326,306)
*** combined adjacent 521 cells.
[[[x1,1,3,3],[x2,2,3,2],[x3,2,3,6],[x4,2,3,4]]]
*** reordering by d+sotd+o.
*** using the sum of squares projection.
[x2,2]
[x1,2]
[x4,4]
[x3,4]
time = 163 ms.
11 53(4,7) 189(0,0) 317(76,64)
*** combined adjacent 189 cells.
[[[x1,2,4,3],[x2,3,4,4],[x3,4,4,4],[x4,2,4,3]]]
*** reordering by d+sotd+o.
*** using the sum of squares projection.
[x4,1]
[x1,3]
[x2,4]
[x3,6]
time = 886 ms.
37 277(21,24) 1677(182,84) 1725(850,0)
*** combined adjacent 70 cells.
total = 9,202 ms.
? wt=getwalltime();ranc7(7*31);print("total = ",strtime(getwalltime()-wt),".");
[[[x1,13,14,21,8],[x2,7,10,17,6],[x3,10,14,16,7]]]
*** reordering by sum(d+td in terms)+o.
*** using the sum of squares projection.
[x2,3]
[x3,7]
[x1,22]
time = 390 ms.
87 1795(42,71) 561(314,1257)
*** combined adjacent 244 cells.
[[[x1,4,3,11,4],[x2,4,5,7,3],[x3,10,7,20,8],[x4,3,3,7,3]]]
*** reordering by sum(d+td in terms)+o.
*** using the sum of squares projection.
[x4,2]
[x2,3]
[x1,4]
[x3,10]
time = 222 ms.
41 245(22,18) 901(276,140) 667(0,506)
*** combined adjacent 268 cells.
[[[x1,9,10,16,6],[x2,9,11,17,6],[x3,12,11,23,8]]]
*** reordering by sum(d+td in terms)+o.
*** using the sum of squares projection.
[x1,2]
[x2,4]
[x3,16]
time = 321 ms.
59 675(30,29) 1162(326,306)
*** combined adjacent 521 cells.
[[[x1,3,2,8,3],[x2,3,6,5,2],[x3,8,7,15,6],[x4,5,4,9,4]]]
*** reordering by sum(d+td in terms)+o.
*** using the sum of squares projection.
[x2,2]
[x1,2]
[x4,4]
[x3,4]
time = 163 ms.
11 53(4,7) 189(0,0) 317(76,64)
*** combined adjacent 189 cells.
[[[x1,4,8,11,3],[x2,6,11,15,4],[x3,9,12,16,4],[x4,4,8,11,3]]]
*** reordering by sum(d+td in terms)+o.
*** using the sum of squares projection.
[x4,1]
[x1,3]
[x2,4]
[x3,6]
time = 874 ms.
37 277(21,24) 1677(182,84) 1725(850,0)
*** combined adjacent 70 cells.
total = 9,311 ms.
本日のC.A.D.
CAD において，変数順序の選定は非常に重要です． 例えば
? projs([x1,x2,x3,x4],[-218 + 955*x1^2*x2*x3 + 520*x2*x3^3 + 1017*x3^4 + 250*x1*x2*x4 - 391*x2^3*x4 - 525*x1*x3*x4^2]); *** using Lazard's method (MPP17). [x4,1] [x3,3] [x2,4] [x1,4] time = 3min, 10,999 ms.
Random-c-5, A.W.Strzebonski "Cylindrical Algebraic Decomposition using validated numerics"
ですが，変数順序を変えると
? projs([x3,x2,x1,x4],[-218 + 955*x1^2*x2*x3 + 520*x2*x3^3 + 1017*x3^4 + 250*x1*x2*x4 - 391*x2^3*x4 - 525*x1*x3*x4^2]); *** using Lazard's method (MPP17). [x4,1] [x1,3] [x2,4] [x3,6] time = 734 ms. ? projs([x3,x2,x4,x1],[-218 + 955*x1^2*x2*x3 + 520*x2*x3^3 + 1017*x3^4 + 250*x1*x2*x4 - 391*x2^3*x4 - 525*x1*x3*x4^2]); *** using Lazard's method (MPP17). [x1,1] [x4,2] [x2,4] [x3,7] time = 686 ms.
といった結果を得ます．
変数（消去）順序については
１．次数が低い変数から
２．それを含む項の全次数の最大値が小さい変数から
３．出現が少ない変数から
とする heuristic (C.W.Brown "Companion to the Tutorial Cylindrical Algebraic Decomposition") が有名で，提唱者の Brown 博士の tarski ( https://github.com/chriswestbrown/tarski ) にも suggest-var-order として実装されています．
> (suggest-var-order [-218 + 955*x1^2*x2*x3 + 520*x2*x3^3 + 1017*x3^4 + 250*x1*x2*x4 - 391*x2^3*x4 - 525*x1*x3*x4^2 > 0]) ( x3:sym x2:sym x4:sym x1:sym )
上記と同じ ISSAC 2004 ( http://www3.risc.jku.at/conferences/issac2004/program.html ) で発表された Redlog の heuristic ( https://www.researchgate.net/publication/2899124_Efficient_Projection_Orders_for_CAD ) は，各変数に対して，それを主変数とする projection を 1 step だけ求め，得られた全ての項の全次数の和が最小の変数を各 level 毎に選んでゆくものです．
Redlog Revision 5715 of 2021-03-11, 12:26:36Z
(c) 1992-2021 T. Sturm and A. Dolzmann (www.redlog.eu)
6: rlcadporder(-218 + 955*x1^2*x2*x3 + 520*x2*x3^3 + 1017*x3^4 + 250*x1*x2*x4 - 391*x2^3*x4 - 525*x1*x3*x4^2 > 0);
+++ Optimizing projection order.
+ input order by blocks: -> (x4 x3 x2 x1)
+ current input block: (x4 x3 x2 x1)
[x4:50] [x3:312] [x2:122] [x1:48] choose x1
[x4:309] [x3:880] [x2:428] choose x4
[x3:2146] [x2:1398] choose x2
+ reordered block: (x1 x4 x2 x3)
+ optimized order: -> (x1 x4 x2 x3)
{x1,x4,x2,x3}$
Time: 201330 ms plus GC time: 18470 ms比較：
https://www.cl.cam.ac.uk/~lp15/papers/Arith/Huang-3heuristics.pdf
ML の導入例：
https://arxiv.org/pdf/1904.11061.pdf
Taxicab number $\operatorname{Ta}(4)$
今回は「 $x^3+y^3=a,\ 1\le x\le y$ を満たす格子点 $(x,\ y)$ の個数が $4$ である」ような $a$ の最小値，および，それに対応する $4$ 個の格子点を求める問題です．
処理時間とメモリーとのバランスから，適当な上界を与え，出現をカウントしつつ，$x^3+y^3$ を生成する方針を採りました．利用した Map は，いわゆる array や dictionary にあたるオブジェクトですが，個数ではなく，格子点の $x$ 座標を係数とする多項式（リストを用いると 30sec 以上処理時間が延びます）を対応させることでワンパスを実現しています．
GP/PARI CALCULATOR Version 2.13.1 (released)
amd64 running mingw (x86-64/GMP-6.1.2 kernel) 64-bit version
compiled: Jan 16 2021, gcc version 8.3-posix 20190406 (GCC)
threading engine: single
(readline v8.0 enabled, extended help enabled)
Copyright (C) 2000-2020 The PARI Group
PARI/GP is free software, covered by the GNU General Public License, and comes WITHOUT ANY WARRANTY
WHATSOEVER.
Type ? for help, \q to quit.
Type ?17 for how to get moral (and possibly technical) support.
parisize = 30000000000, primelimit = 500000
(18:56) gp > #
timer = 1 (on)
(18:56) gp > n=7*10^12;d=4;
(18:56) gp > A=Map();e=d-1;for(i=1,floor((n/2)^(1/3)),i3=i^3;for(j=i,floor((n-i^3)^(1/3)),a=i3+j^3;if(mapisdefined(A,a,&p),mapput(A,a,p*x+i);P=mapget(A,a);if(poldegree(P)==e,print([a,Vec(P)])),mapput(A,a,i))));
[6963472309248, [2421, 5436, 10200, 13322]]
time = 11min, 57,031 ms.
(19:09) gp > #A
%3 = 161516889
(19:10) gp > 27792341+16331743+47676463+37255436+32460906
%4 = 161516889また，領域 $\ n_1< x^3+y^3\le n_2\ $ により分割して
ta(d,n1,n2)={my(A=Map(),e=d-1);
for(i=1,(n1/2)^(1/3),for(j=floor((n1-i^3)^(1/3))+1,(n2-i^3)^(1/3),a=i^3+j^3;if(mapisdefined(A,a,&p),mapput(A,a,p*x+i);P=mapget(A,a);if(poldegree(P)==e,print([a,Vec(P)])),mapput(A,a,i))));
for(i=floor((n1/2)^(1/3))+1,(n2/2)^(1/3),for(j=i,(n2-i^3)^(1/3),a=i^3+j^3;if(mapisdefined(A,a,&p),mapput(A,a,p*x+i);P=mapget(A,a);if(poldegree(P)==e,print([a,Vec(P)])),mapput(A,a,i))));
#A;};のように生成すれば
parisize = 10000000000 (18:34) gp > ta(4,1,5*10^11) time = 2min, 9,063 ms. %2 = 27792341 (18:35) gp > ta(4,5*10^11,1*10^12) time = 1min, 20,875 ms. %2 = 16331743 (18:35) gp > ta(4,1*10^12,3*10^12) time = 3min, 57,625 ms. %2 = 47676463 (18:35) gp > ta(4,3*10^12,5*10^12) time = 3min, 12,968 ms. %2 = 37255436 (18:35) gp > ta(4,5*10^12,7*10^12) [6963472309248, [2421, 5436, 10200, 13322]] time = 2min, 48,703 ms. %2 = 32460906
といった並列化も可能です．
なお，今回の記事は jurupapa さんの記事 ( https://maxima.hatenablog.jp/entry/2021/06/10/003011 ) にインスパイアされたもので，以下は，最初に書いた maxima 版．
C:\maxima-5.45.0\bin>maxima.bat Maxima 5.45.0 https://maxima.sourceforge.io using Lisp SBCL 2.0.0 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) (display2d:off,showtime:on)$ Evaluation took 0.0000 seconds (0.0000 elapsed) using 0 bytes. (%i2) (n:9*10^7,d:3)$ /* 87,539,319 */ Evaluation took 0.0000 seconds (0.0000 elapsed) using 0 bytes. (%i3) /* method 1: using a generating function */ (q:subst([a=1,x=1],p:args(rat(expand(sum(a^i*x^(i^3)*sum(x^(j^3),j,i,(n-i^3)^(1/3)),i,1,(n/2)^(1/3)))))), for i in q do if d=pop(q) then disp(pop(p)) else pop(p))$ (a^255+a^228+a^167)*x^87539319 Evaluation took 8.4690 seconds (8.4800 elapsed) using 919.738 MB. (%i4) /* method 2: using a counting array */ (kill(C),C[i]:=[],for i:1 thru (n/2)^(1/3) do for j:i thru (n-i^3)^(1/3) do (k:i^3+j^3,if d=length(C[k]:append(C[k],[[i,j]])) then disp(k,C[k])))$ 87539319 [[167,436],[228,423],[255,414]] Evaluation took 2.6870 seconds (2.6810 elapsed) using 598.288 MB. (%i5) /* method 3: find adjacent occurrences in the sorted list */ (A:sort(A0:create_list(i^3+j^3,i,1,floor((n/2)^(1/3)),j,i,floor((n-i^3)^(1/3)))), a:pop(A),len:1,for i in A do if a=b:pop(A) then len:len+1 else (if d=len then (disp(a),return()),a:b,len:1), pos:sublist_indices(A0,lambda([x],x=a)),ij:create_list([i,j],i,1,floor((n/2)^(1/3)),j,i,floor((n-i^3)^(1/3))),disp(map(lambda([k],ij[k]),pos)))$ 87539319 [[167,436],[228,423],[255,414]] Evaluation took 0.9840 seconds (0.9860 elapsed) using 250.403 MB.
このうち，個数が $4$ の場合に適用できたのは method 3 のみでした．
C:\maxima-5.45.0\bin>maxima.bat -X "--dynamic-space-size 30000" Lisp options: (--dynamic-space-size 30000) Maxima 5.45.0 https://maxima.sourceforge.io using Lisp SBCL 2.0.0 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) (display2d:off,showtime:on)$ Evaluation took 0.0000 seconds (0.0000 elapsed) using 0 bytes. (%i2) (n:7*10^12,d:4)$ Evaluation took 0.0000 seconds (0.0000 elapsed) using 0 bytes. (%i3) (A:sort(A0:create_list(i^3+j^3,i,1,floor((n/2)^(1/3)),j,i,floor((n-i^3)^(1/3)))),a:pop(A),len:1,for i in A do if a=b:pop(A) then len:len+1 else (if d=len then (disp(a),return()),a:b,len:1),pos:sublist_indices(A0,lambda([x],x=a)),ii:create_list(i,i,1,floor((n/2)^(1/3)),j,i,floor((n-i^3)^(1/3))),disp(map(lambda([k],[iik:ii[k],(a-iik^3)^(1/3)]),pos)))$ 6963472309248 [[2421,19083],[5436,18948],[10200,18072],[13322,16630]] Evaluation took 1295.1720 seconds (1297.1910 elapsed) using 320617.394 MB.
様々な言語によるコード https://rosettacode.org/wiki/Taxicab_numbers
本日のC.A.D.
Quantified Formula: QE(forall a b c d) ((a^2*c + -1*a*c^2 + -1*a*d^2 + b^2*c + -1*a^3*c^2 + -1*a^3*d^2 + a^2*c^3 + a^2*c*d^2 + -1*a*b^2*c^2 + -1*a*b^2*d^2 + b^2*c^3 + b^2*c*d^2 = 0 and -1*a^2*d + -1*b^2*d + b*c^2 + b*d^2 + -1*a^2*b*c^2 + -1*a^2*b*d^2 + a^2*c^2*d + a^2*d^3 + -1*b^3*c^2 + -1*b^3*d^2 + b^2*c^2*d + b^2*d^3 = 0 and -1*b < 0 and -1*d < 0) => (-1*a + c = 0 and -1*b + d = 0))
https://github.com/ths-rwth/smtrat/blob/JLAMP-CDCAC/benchmarks/31.smt2
そのまま実行すると
? tst12([all,all,all,all],[b,a,c,d],(f1,f2,f3,f4,f5,f6)->imp(f1*f2*f3*f4,f5*f6),"a^2*c + -1*a*c^2 + -1*a*d^2 + b^2*c + -1*a^3*c^2 + -1*a^3*d^2 + a^2*c^3 + a^2*c*d^2 + -1*a*b^2*c^2 + -1*a*b^2*d^2 + b^2*c^3 + b^2*c*d^2 = 0, -1*a^2*d + -1*b^2*d + b*c^2 + b*d^2 + -1*a^2*b*c^2 + -1*a^2*b*d^2 + a^2*c^2*d + a^2*d^3 + -1*b^3*c^2 + -1*b^3*d^2 + b^2*c^2*d + b^2*d^3 = 0, -1*b < 0, -1*d < 0, -1*a + c = 0, -1*b + d = 0",7);Ans() *** using the sum of squares projection. [d,4] [c,9] [a,24] [b,28] time = 3,877 ms. 111 4121(150,127) 63963(2022,4503) 643841(28592,52664) *** combined adjacent 643840 cells. 1[true,true,true,true] time = 2min, 25,341 ms.
となり，否定が false を見る方針なら
? tst12([ex,ex,ex,ex],[b,a,c,d],(f1,f2,f3,f4,f5,f6)->f1*f2*f3*f4*not(f5*f6),"a^2*c + -1*a*c^2 + -1*a*d^2 + b^2*c + -1*a^3*c^2 + -1*a^3*d^2 + a^2*c^3 + a^2*c*d^2 + -1*a*b^2*c^2 + -1*a*b^2*d^2 + b^2*c^3 + b^2*c*d^2 = 0, -1*a^2*d + -1*b^2*d + b*c^2 + b*d^2 + -1*a^2*b*c^2 + -1*a^2*b*d^2 + a^2*c^2*d + a^2*d^3 + -1*b^3*c^2 + -1*b^3*d^2 + b^2*c^2*d + b^2*d^3 = 0, -1*b < 0, -1*d < 0, -1*a + c = 0, -1*b + d = 0");Ans() *** using Lazard's method (MPP17). [d,4] [c,9] [a,24] [b,28] time = 1,924 ms. 55 2057(32,44) 31935(982,2225) 0(14262,26188) 1[false] time = 42,614 ms.
となりますが，$p\to[q\land r]$ は $[p\to q]\land[p\to r]$，更に $[p\to q]\land[[p\land q]\to r]$ と変形できるので，分割して実行すると
? tst12([all,all,all,all],[b,a,c,d],(f1,f2,f3,f4,f5)->imp(f1*f2*f3*f4,f5),"a^2*c + -1*a*c^2 + -1*a*d^2 + b^2*c + -1*a^3*c^2 + -1*a^3*d^2 + a^2*c^3 + a^2*c*d^2 + -1*a*b^2*c^2 + -1*a*b^2*d^2 + b^2*c^3 + b^2*c*d^2 = 0, -1*a^2*d + -1*b^2*d + b*c^2 + b*d^2 + -1*a^2*b*c^2 + -1*a^2*b*d^2 + a^2*c^2*d + a^2*d^3 + -1*b^3*c^2 + -1*b^3*d^2 + b^2*c^2*d + b^2*d^3 = 0, -1*b < 0, -1*d < 0, -1*a + c = 0");Ans()
*** using Lazard's method (MPP17).
[d,3]
[c,7]
[a,18]
[b,10]
time = 1,016 ms.
35 1189(118,60) 15615(574,1077) 129501(8480,10604)
*** combined adjacent 129500 cells.
1[true,true,true,true]
time = 22,563 ms.
? tst12([all,all,all],[b,a,d],(f1,f2,f3,f4,f5)->imp(f1*f2*f3*f4,f5),strrep("a^2*c + -1*a*c^2 + -1*a*d^2 + b^2*c + -1*a^3*c^2 + -1*a^3*d^2 + a^2*c^3 + a^2*c*d^2 + -1*a*b^2*c^2 + -1*a*b^2*d^2 + b^2*c^3 + b^2*c*d^2 = 0, -1*a^2*d + -1*b^2*d + b*c^2 + b*d^2 + -1*a^2*b*c^2 + -1*a^2*b*d^2 + a^2*c^2*d + a^2*d^3 + -1*b^3*c^2 + -1*b^3*d^2 + b^2*c^2*d + b^2*d^3 = 0, -1*b < 0, -1*d < 0, -1*b + d = 0",["c"],["a"]));Ans()
*** using Lazard's method (MPP17).
[d,4]
[a,5]
[b,4]
time = 83 ms.
11 125(20,14) 1111(44,119)
*** combined adjacent 1110 cells.
1[true,true,true]
time = 133 ms.となり，否定が false を見る方針でも
? tst12([ex,ex,ex,ex],[b,a,c,d],andx,"a^2*c + -1*a*c^2 + -1*a*d^2 + b^2*c + -1*a^3*c^2 + -1*a^3*d^2 + a^2*c^3 + a^2*c*d^2 + -1*a*b^2*c^2 + -1*a*b^2*d^2 + b^2*c^3 + b^2*c*d^2 = 0, -1*a^2*d + -1*b^2*d + b*c^2 + b*d^2 + -1*a^2*b*c^2 + -1*a^2*b*d^2 + a^2*c^2*d + a^2*d^3 + -1*b^3*c^2 + -1*b^3*d^2 + b^2*c^2*d + b^2*d^3 = 0, -1*b < 0, -1*d < 0, -1*a + c != 0");Ans()
*** using Lazard's method (MPP17).
[d,3]
[c,7]
[a,18]
[b,10]
time = 999 ms.
17 591(26,14) 7174(258,524) 0(4094,3930)
1[false]
time = 8,024 ms.
? tst12([ex,ex,ex],[b,a,d],andx,strrep("a^2*c + -1*a*c^2 + -1*a*d^2 + b^2*c + -1*a^3*c^2 + -1*a^3*d^2 + a^2*c^3 + a^2*c*d^2 + -1*a*b^2*c^2 + -1*a*b^2*d^2 + b^2*c^3 + b^2*c*d^2 = 0, -1*a^2*d + -1*b^2*d + b*c^2 + b*d^2 + -1*a^2*b*c^2 + -1*a^2*b*d^2 + a^2*c^2*d + a^2*d^3 + -1*b^3*c^2 + -1*b^3*d^2 + b^2*c^2*d + b^2*d^3 = 0, -1*b < 0, -1*d < 0, -1*b + d != 0",["c"],["a"]),7);Ans()
*** using the sum of squares projection.
[d,4]
[a,5]
[b,4]
time = 125 ms.
5 59(6,2) 0(20,48)
1[false]
time = 49 ms.のように高速化できます．
ところでこの問題，入力の 2 つの等式
$\begin{eqnarray}
&& a^3 c^2+a^3 d^2-a^2 c^3-a^2 c d^2-a^2 c+a b^2 c^2+a b^2 d^2+a c^2+a d^2-b^2 c^3-b^2 c d^2-b^2 c=0,\\
&& a^2 b c^2+a^2 b d^2-a^2 c^2 d-a^2 d^3+a^2 d+b^3 c^2+b^3 d^2-b^2 c^2 d-b^2 d^3+b^2 d-b c^2-b d^2=0
\end{eqnarray}$
は，それぞれ
$\begin{eqnarray}
&& a(a^2+b^2+1)(c^2+d^2)=c(c^2+d^2+1)(a^2+b^2),\\
&& b(a^2+b^2-1)(c^2+d^2)=d(c^2+d^2-1)(a^2+b^2)
\end{eqnarray}$
と整理できるので，辺々を $(a^2+b^2)(c^2+d^2)$ で割りたくなりますが，$0< b\ \land\ 0< d$ により，それが可能であり
$\begin{eqnarray}
&& a+\frac{a}{a^2+b^2}=c+\frac{c}{c^2+d^2},\\
&& b-\frac{b}{a^2+b^2}=d-\frac{d}{c^2+d^2}
\end{eqnarray}$
そして，$s=a+bi,\ t=c+di$ とおくと
${\displaystyle s+\frac{1}{s}=t+\frac{1}{t}}$
つまり
${\displaystyle \frac{(s-t)(st-1)}{st}=\frac{1}{s}(s-t) \left( s-\frac{\bar{{t}}}{|t|^2}\right)=0}$
ここで，再び，$\mathrm{Im}{\, s}=b>0\ \land\ \mathrm{Im}{\, \bar{t}}=-d<0$ により，$s=t$ なので
$\begin{eqnarray}
\forall a\ \forall b\ \forall c\ \forall d &&[\, a(a^2+b^2+1)(c^2+d^2)=c(c^2+d^2+1)(a^2+b^2)\\
&& \land\ b(a^2+b^2-1)(c^2+d^2)=d(c^2+d^2-1)(a^2+b^2)\\
&& \land\ 0< b\ \land\ 0< d\ \to\ [\,a=c\ \land\ b=d\,]\,]
\end{eqnarray}$
を得ます．是即ち「上半平面を定義域とする Joukowsky 変換は単射である」です．
また，上記から入力の 2 つの等式の連言は $s=t\ \lor\ st=1\ \lor s=0\ \lor\ t=0$ と等価，つまり
$\begin{eqnarray}
\forall a\ \forall b\ \forall c\ \forall d &&[\,[\, a(a^2+b^2+1)(c^2+d^2)=c(c^2+d^2+1)(a^2+b^2)\\
&& \land\ b(a^2+b^2-1)(c^2+d^2)=d(c^2+d^2-1)(a^2+b^2)\,]\\
&& \leftrightarrow\ [\,[\,a=c\ \land\ b=d\,]\ \lor\ [\,ac-bd=1\ \land\ ad+bc=0\,]\\
&& \lor\ [\,a=0\ \land\ b=0\,]\ \lor\ [\,c=0\ \land\ d=0\,]\,]\,]
\end{eqnarray}$
であることも解かるので，実行してみると
? mulp=5; /* -> */ ? tst12([ex,ex,ex,ex],[b,a,c,d],(f1,f2,f3,f4,f5,f6,f7,f8,f9,f10)->f1*f2*not(or(or(or(f3*f4,f5*f6),f7*f8),f9*f10)),"a^2*c + -1*a*c^2 + -1*a*d^2 + b^2*c + -1*a^3*c^2 + -1*a^3*d^2 + a^2*c^3 + a^2*c*d^2 + -1*a*b^2*c^2 + -1*a*b^2*d^2 + b^2*c^3 + b^2*c*d^2 = 0, -1*a^2*d + -1*b^2*d + b*c^2 + b*d^2 + -1*a^2*b*c^2 + -1*a^2*b*d^2 + a^2*c^2*d + a^2*d^3 + -1*b^3*c^2 + -1*b^3*d^2 + b^2*c^2*d + b^2*d^3 = 0, -1*a + c = 0, -1*b + d = 0, a*c - b*d = 1, a*d + b*c = 0, a = 0, b = 0, c = 0, d = 0");Ans() *** using Lazard's method (MPP17). [d,6] [c,13] [a,45] [b,209] time = 17,627 ms. 987 82344(314,704) 2147780(24596,93892) 0(638480,2355568) 1[false] time = 1h, 44min, 23,022 ms. /* <- */ ? tst12([ex,ex,ex,ex],[b,a,c,d],(f1,f2,f3,f4,f5,f6,f7,f8,f9,f10)->not(f1*f2)*or(or(or(f3*f4,f5*f6),f7*f8),f9*f10),"a^2*c + -1*a*c^2 + -1*a*d^2 + b^2*c + -1*a^3*c^2 + -1*a^3*d^2 + a^2*c^3 + a^2*c*d^2 + -1*a*b^2*c^2 + -1*a*b^2*d^2 + b^2*c^3 + b^2*c*d^2 = 0, -1*a^2*d + -1*b^2*d + b*c^2 + b*d^2 + -1*a^2*b*c^2 + -1*a^2*b*d^2 + a^2*c^2*d + a^2*d^3 + -1*b^3*c^2 + -1*b^3*d^2 + b^2*c^2*d + b^2*d^3 = 0, -1*a + c = 0, -1*b + d = 0, a*c - b*d = 1, a*d + b*c = 0, a = 0, b = 0, c = 0, d = 0");Ans() *** using Lazard's method (MPP17). [d,6] [c,13] [a,45] [b,209] time = 17,751 ms. 987 82345(314,704) 2141738(24602,93897) 0(634478,2355040) 1[false] time = 1h, 43min, 47,907 ms.
といった結果を得ます．
本日のC.A.D.
$\exists a\ \exists b\ \exists c\ \forall x\ \forall y\ [\ \neg[a=0 \land b=0] \land\, [\,ax+by+c=0\ \to\ dx^2+exy+fy^2=1\,]\,]$
C.W.Brown "Simple CAD Construction and its Applications"
? tst12([ex,ex,ex,all,all],[d,e,f,a,b,c,x,y],(f1,f2,f3,f4)->not(f1*f2)*imp(f3,f4),"a==0,b==0,a*x+b*y+c==0,d*x^2+e*x*y+f*y^2==1",7);Ans() *** using the sum of squares projection. [y,2] [x,4] [c,7] [b,15] [a,1] [f,27] [e,1] [d,1] time = 2,808 ms. 3 9(0,0) 115(102,69) 345(0,0) 1640(2712,1450) 7840(1840,692) 47180(7340,5152) 80(1712,24) *** combined adjacent 64 cells. 1[d = [d,1],e = [e,1],[f,1] < f,true,true,true,true,true] 2[[d,1] < d,e < [e,1],f = [4*f*d-e^2,1],true,true,true,true,true] 3[[d,1] < d,e = [e,1],f = [f,1],true,true,true,true,true] 4[[d,1] < d,[e,1] < e,f = [4*f*d-e^2,1],true,true,true,true,true] time = 19,025 ms. ? tst12([ex,ex,ex,all,all],[f,d,e,a,b,c,x,y],(f1,f2,f3,f4)->not(f1*f2)*imp(f3,f4),"a==0,b==0,a*x+b*y+c==0,d*x^2+e*x*y+f*y^2==1",7);Ans() *** using the sum of squares projection. [y,2] [x,4] [c,7] [b,15] [a,1] [e,27] [d,1] [f,1] time = 2,257 ms. 3 9(0,0) 115(482,107) 345(0,0) 1640(2834,1498) 7840(1832,660) 47280(7452,5092) 88(1712,24) *** combined adjacent 72 cells. 1[f = [f,1],[d,1] < d,e = [e,1],true,true,true,true,true] 2[[f,1] < f,d = [d,1],e = [e,1],true,true,true,true,true] 3[[f,1] < f,[d,1] < d,e = [-4*f*d+e^2,1] or e = [-4*f*d+e^2,2],true,true,true,true,true] time = 15,817 ms.
Wolfram Language 12.3.0 Engine for Linux x86 (64-bit)
Copyright 1988-2021 Wolfram Research, Inc.
In[1]:= SetSystemOptions["InequalitySolvingOptions"->"CADSortVariables"->False];
In[2]:= CylindricalDecomposition[Exists[{a,b,c},Not[a==0 && b==0],ForAll[{x,y},a*x+b*y+c==0,d*x^2+e*x*y+f*y^2==1]
],{d,e,f}]// Timing
Out[2]= {6.796608, (d == 0 && e == 0 && f > 0) || (d > 0 && f == e^2/(4*d))}
In[3]:= CylindricalDecomposition[Exists[{a,b,c},Not[a==0 && b==0],ForAll[{x,y},a*x+b*y+c==0,d*x^2+e*x*y+f*y^2==1]
],{f,d,e}]// Timing
Out[3]=
{6.390617, (f == 0 && d > 0 && e == 0) ||
(f > 0 && ((d == 0 && e == 0) || (d > 0 && (e == -2*Sqrt[d*f] || e == 2*Sqrt[d*f]))))}
本日のC.A.D.
a, b, c がある 3 角形の 3 辺の長さ，ma, mb, mc がその各辺の中点と対頂点との距離（つまり，中線の長さ）ならば
${\displaystyle\frac{3}{4}}(a+b+c) < m_a+m_b+m_c < a+b+c$
8.1 Bottema, et. al. "GEOMETRIC INEQUALITIES"
等式制約は中線定理，非退化 3 角形の存在条件は WLOG で 0＜a, a≦b, b≦c, c＜a+b，そして，a+b+c=1，従って，比の値は ma+mb+mc=x とおけるので，c, mc を消去したのち，parisizemax = 20G で実行すると
? tst12([ex,ex,ex,ex],[x,a,b,ma,mb],andx,"0<a,a<=b,b<=(1-a-b),(1-a-b)<a+b,0<ma,0<mb,0<x-ma-mb,a^2+b^2==2*((x-ma-mb)^2+((1-a-b)/2)^2),b^2+(1-a-b)^2==2*(ma^2+(a/2)^2),(1-a-b)^2+a^2==2*(mb^2+(b/2)^2),0<x");Ans() *** using Lazard's method (MPP17). [mb,4] [ma,6] [b,15] [a,68] [x,1332] time = 3min, 24,355 ms. 5771 984695(168,6317) 869098(224436,1848995) 869111(307784,1653038) 503(0,229566) *** combined adjacent 502 cells. 1[[4*x-3,1] < x < [x-1,1],true,true,true,true] time = 3h, 13min, 15,870 ms.
のようにちょっと大変です．
こうした cell の爆発の回避には，やはり下位の lifting で真理値の定まる制約の採用が効果的なので，求むべきは，x の必要条件です．
ところで，今回の問題の難しさは，2 次式で定義された ma, mb, mc の和を扱うことにありますが，2 乗和であれば
? tst12([ex,ex,ex,ex,ex],[y,a,b,ma,mb,mc],andx,"0<a,a<=b,b<=(1-a-b),(1-a-b)<a+b,0<ma,0<mb,0<mc,a^2+b^2==2*(mc^2+((1-a-b)/2)^2),b^2+(1-a-b)^2==2*(ma^2+(a/2)^2),(1-a-b)^2+a^2==2*(mb^2+(b/2)^2),ma^2+mb^2+mc^2==y,0<y");Ans() *** using Lazard's method (MPP17). [mc,3] [mb,4] [ma,6] [b,14] [a,38] [y,242] time = 9,897 ms. 727 47811(48,910) 30824(13820,74979) 30824(4108,56040) 30824(3392,7928) 114(1292,1748) *** combined adjacent 113 cells. 1[[4*y-1,1] <= y < [8*y-3,1],true,true,true,true,true] time = 4min, 13,126 ms.
のように難しくなく，これにより x の値域の拡大集合が得られます．つまり，この 2 乗和の条件と一般的な 2 乗和と和との関係から
? tst12([ex,ex,ex,ex],[x,y,ma,mb,mc],andx,"0<ma,0<mb,0<mc,ma^2+mb^2+mc^2==y,1/4<=y,y<3/8,ma+mb+mc==x,0<x");Ans() *** using Lazard's method (MPP17). [mc,3] [mb,4] [ma,5] [y,6] [x,7] time = 67 ms. 11 28(0,12) 228(8,16) 924(108,222) 7(62,146) *** combined adjacent 6 cells. 1[[2*x-1,1] < x < [8*x^2-9,2],true,true,true,true] time = 342 ms.
という必要条件が得られ，これを 0＜x の代わりの制約として与えれば
? tst12([ex,ex,ex,ex],[x,a,b,ma,mb],andx,"0<a,a<=b,b<=(1-a-b),(1-a-b)<a+b,0<ma,0<mb,0<x-ma-mb,a^2+b^2==2*((x-ma-mb)^2+((1-a-b)/2)^2),b^2+(1-a-b)^2==2*(ma^2+(a/2)^2),(1-a-b)^2+a^2==2*(mb^2+(b/2)^2),1/2<x,8*x^2-9<0");Ans() *** using Lazard's method (MPP17). [mb,4] [ma,6] [b,15] [a,68] [x,1333] time = 3min, 24,448 ms. 1279 223043(41,1481) 130080(50064,411377) 130080(58236,277036) 503(0,36682) *** combined adjacent 502 cells. 1[[4*x-3,1] < x < [x-1,1],true,true,true,true] time = 37min, 17,195 ms.
といった結果を得ます．
この方法は問題の形に依存していますが，一般に適用可能なものとして，束縛変数のみの制約に対する解集合の各束縛変数の空間への射影を用いる方法が考えられます．例えば，最下位の束縛変数 a の束縛を外し，自由変数とした式を QE すれば
? tst12([ex,ex,ex,ex],[a,b,ma,mb,mc],andx,"0<a,a<=b,b<=(1-a-b),(1-a-b)<a+b,0<ma,0<mb,0<mc,a^2+b^2==2*(mc^2+((1-a-b)/2)^2),b^2+(1-a-b)^2==2*(ma^2+(a/2)^2),(1-a-b)^2+a^2==2*(mb^2+(b/2)^2)");Ans() *** using Lazard's method (MPP17). [mc,2] [mb,2] [ma,2] [b,6] [a,6] time = 97 ms. 9 8(6,9) 8(0,0) 8(0,0) 4(0,0) *** combined adjacent 3 cells. 1[[a,1] < a <= [3*a-1,1],true,true,true,true] time = 23 ms.
という必要条件が得られ，これを 0＜a の代わりの制約として与えれば
? tst12([ex,ex,ex,ex],[x,a,b,ma,mb],andx,"0<a,a<=b,b<=(1-a-b),(1-a-b)<a+b,0<ma,0<mb,0<x-ma-mb,a^2+b^2==2*((x-ma-mb)^2+((1-a-b)/2)^2),b^2+(1-a-b)^2==2*(ma^2+(a/2)^2),(1-a-b)^2+a^2==2*(mb^2+(b/2)^2),1/2<x,8*x^2-9<0,3*a-1<0");Ans() *** using Lazard's method (MPP17). [mb,4] [ma,6] [b,15] [a,68] [x,1333] time = 3min, 23,477 ms. 1279 27839(41,1481) 128801(1220,37978) 128801(55678,269357) 503(0,36678) *** combined adjacent 502 cells. 1[[4*x-3,1] < x < [x-1,1],true,true,true,true] time = 15min, 25,279 ms.
といった結果を得ます．
さて，こうした入出力の組合わせ，書き換え，段階的な QE には，総合環境を目指して開発されている tarski ( https://www.usna.edu/Users/cs/wcbrown/tarski/index.html ) が適しており，以下のような前処理により，処理時間の短縮が可能です．
$ tarski > (version) "tarski version 1.29, Fri Mar 19 16:06:57 EDT 2021":str
まず，qfr ( quantified formula rewriting, https://www.usna.edu/Users/cs/cstech/tr/reports/2007-01.orig.pdf ) を適用して，1 次の等式制約により，mc, c を消去，次に，QEPCADB を段階的に適用して，mb, b を消去します．
> (def F [0<a /\ a<=b /\ b<=c /\ c<a+b /\ 0<ma /\ 0<mb /\ 0<mc /\ a^2+b^2=2(mc^2+(c/2)^2) /\ b^2+c^2=2(ma^2+(a/2)^2) /\ c^2+a^2=2(mb^2+(b/2)^2) /\ ma+mb+mc=x /\ 0<x /\ a+b+c=1]) :void > (def F (qfr[ex c,mc[$F]])) :void > (def F (qepcad-qe[ex mb[$F]])) :void > (def F (qepcad-qe[ex b[$F]])) F :void [a > 0 /\ ma > 0 /\ 3 a - 1 <= 0 /\ x > 0 /\ 4 ma^2 + 2 a - 1 >= 0 /\ 4 ma^2 - 9 a^2 + 8 a - 2 <= 0 /\ 2 ma + a - 1 < 0 /\ 8 x^2 - 16 ma x + 4 ma^2 - 9 a^2 > 0 /\ 16 x^4 - 64 ma x^3 + 80 ma^2 x^2 - 36 a^2 x^2 - 32 ma^3 x + 72 a^2 ma x - 72 a ma^2 + 36 ma^2 + 18 a^3 - 45 a^2 + 36 a - 9 = 0]:tar
ここから先は，+N, +L オプションを付けても "Qepcad failure!":err となります．一方，必要条件は
> (def G [0<a /\ a<=b /\ b<c+a /\ c<a+b /\ 0<ma /\ 0<mb /\ 0<mc /\ a^2+b^2=2(mc^2+(c/2)^2) /\ b^2+c^2=2(ma^2+(a/2)^2) /\ c^2+a^2=2(mb^2+(b/2)^2) /\ ma^2+mb^2+mc^2=y /\ 0<y /\ a+b+c=1]) :void > (def G (qfr[ex c[$G]])) :void > (def G (qepcad-qe[ex a,b,ma,mb,mc[$G]])) G :void [4 y - 1 >= 0 /\ 8 y - 3 < 0]:tar > (def G (qepcad-qe(qfr[ex y,ma,mb,mc[$G /\ 0<ma /\ 0<mb /\ 0<mc /\ ma^2+mb^2+mc^2=y /\ ma+mb+mc=x /\ 0<x]]))) G :void [x > 0 /\ 8 x^2 - 9 < 0 /\ 2 x - 1 > 0]:tar
そして
> (def H [0<a /\ a<=b /\ b<=c /\ c<a+b /\ 0<ma /\ 0<mb /\ 0<mc /\ a^2+b^2=2(mc^2+(c/2)^2) /\ b^2+c^2=2(ma^2+(a/2)^2) /\ c^2+a^2=2(mb^2+(b/2)^2) /\ a+b+c=1]) :void > (def H (qepcad-qe(qfr(exclose H '(a))))) H :void [a > 0 /\ 3 a - 1 <= 0]:tar
となるので，これらの連言を整理，書式の近い redlog に翻訳して
> (syntax 'redlog (normalize [$F /\ $G /\ $H])) "and(a > 0, ma > 0, 3*a - 1 <= 0, 4*ma^2 + 2*a - 1 >= 0, 4*ma^2 - 9*a^2 + 8*a - 2 <= 0, 2*ma + a - 1 < 0, 8*x^2 - 16*ma*x + 4*ma^2 - 9*a^2 > 0, 16*x^4 - 64*ma*x^3 + 80*ma^2*x^2 - 36*a^2*x^2 - 32*ma^3*x + 72*a^2*ma*x - 72*a*ma^2 + 36*ma^2 + 18*a^3 - 45*a^2 + 36*a - 9 = 0, 8*x^2 - 9 < 0, 2*x - 1 > 0)":str
原子式の列をコピペ，tst12 で実行すると
? tst12([ex,ex],[x,a,ma],andx,"a > 0, ma > 0, 3*a - 1 <= 0, 4*ma^2 + 2*a - 1 >= 0, 4*ma^2 - 9*a^2 + 8*a - 2 <= 0, 2*ma + a - 1 < 0, 8*x^2 - 16*ma*x + 4*ma^2 - 9*a^2 > 0, 16*x^4 - 64*ma*x^3 + 80*ma^2*x^2 - 36*a^2*x^2 - 32*ma^3*x + 72*a^2*ma*x - 72*a*ma^2 + 36*ma^2 + 18*a^3 - 45*a^2 + 36*a - 9 = 0, 8*x^2 - 9 < 0, 2*x - 1 > 0");Ans() *** using Lazard's method (MPP17). [ma,6] [a,20] [x,102] time = 854 ms. 103 1508(18,124) 53(8,1123) *** combined adjacent 52 cells. 1[[4*x-3,1] < x < [x-1,1],true,true] time = 3,578 ms.
G, H を付さない場合でも
? tst12([ex,ex],[x,a,ma],andx,"a > 0, ma > 0, 3*a - 1 <= 0, x > 0, 4*ma^2 + 2*a - 1 >= 0, 4*ma^2 - 9*a^2 + 8*a - 2 <= 0, 2*ma + a - 1 < 0, 8*x^2 - 16*ma*x + 4*ma^2 - 9*a^2 > 0, 16*x^4 - 64*ma*x^3 + 80*ma^2*x^2 - 36*a^2*x^2 - 32*ma^3*x + 72*a^2*ma*x - 72*a*ma^2 + 36*ma^2 + 18*a^3 - 45*a^2 + 36*a - 9 = 0");Ans() *** using Lazard's method (MPP17). [ma,6] [a,20] [x,102] time = 849 ms. 317 4184(30,357) 53(60,4430) *** combined adjacent 52 cells. 1[[4*x-3,1] < x < [x-1,1],true,true] time = 11,089 ms.
のように大幅に高速化，「1 人では出来ない（こともないけれど可也しんどい）ことも 2 人居れば何でもない」わけです．
なお，不等式の証明のみであれば，1 人でも何でもありません．
? tst12([ex,ex,ex,ex,ex],[a,b,c,ma,mb,mc],andx,"0<a,a<=b,b<=c,c<a+b,0<ma,0<mb,0<mc,a^2+b^2==2*(mc^2+(c/2)^2),b^2+c^2==2*(ma^2+(a/2)^2),c^2+a^2==2*(mb^2+(b/2)^2),ma+mb+mc<=3/4*(a+b+c)");Ans() *** using Lazard's method (MPP17). [mc,3] [mb,4] [ma,6] [c,18] [b,98] [a,1] time = 1,656 ms. 1 52(0,0) 192(12,80) 192(104,488) 192(0,0) 0(0,0) 1[false] time = 1,358 ms. ? tst12([ex,ex,ex,ex,ex],[a,b,c,ma,mb,mc],andx,"0<a,a<=b,b<=c,c<a+b,0<ma,0<mb,0<mc,a^2+b^2==2*(mc^2+(c/2)^2),b^2+c^2==2*(ma^2+(a/2)^2),c^2+a^2==2*(mb^2+(b/2)^2),ma+mb+mc>=1*(a+b+c)");Ans() *** using Lazard's method (MPP17). [mc,3] [mb,4] [ma,6] [c,16] [b,97] [a,1] time = 1,660 ms. 1 76(0,0) 276(14,100) 276(152,704) 276(0,0) 0(0,0) 1[false] time = 2,007 ms.
