果树问题最优解大全

mathe

我们的配置是要求点四点共线，所以约束方程实际上就是由这些点的三点共线的条件构成的，每个三点共线条件构成一个二次方程，所有这些约束构成一个方程组。
当然，由于射影变换可以任意选择，并不影响解的存在性，我会先选择三行，其中一行和另外两行的交点必须有树，把这一行映射为无穷远直线，另外两行分别映射为横坐标和纵坐标，并且继续上这两行上选择单位点（部分算法会采用其它指定方案）。
通常经过这样指定后，如果存在解，那么会只有有限组，但是这些有限解之间并不一定会射影等价。
由于我们总是指定四个点为有理点，所以如果这时求出来的解中包含一个点的坐标不是有理数（实数），那么这个解必然无法通过射影变换变换成全部是有理点（实数点）。
而我们会发现有很多配置的方程求解结果最后变成三次方程，其中一个根是实数根，另外两个是虚数根。这说明了实数根对应的解可以让所有点都是实数点，但是虚数根对应的结果无论怎么变换，必然有部分点坐标不是实数，所以它们必然无法通过射影变换相互转化，这说明了同一种配置对应的约束方程组解出来的解之间是很可能无法通过射影变换相互转化的。
但是另外一方面，用同一个方程得到的多个解之间有可能是可以通过射影变换相互转化的，比如16棵树15行两个实数解之间应该是等价的：
https://bbs.emath.ac.cn/forum.ph ... d=703&pid=11188

mathe

而给定一个配置的两个解，判断它们之间是否可以通过射影变换相互转化，可以选择其中一个解中四个点，要求任意三点不共线，然后穷举在另外一个解中四个点（任意三点不同线），然后求出射影变换这四个点映射到第一个解中给定的四个点，然后判断余下点之间是否可以在射影变换后坐标完全匹配。
另外一种稍微优化的方法可以计算任意一个直线上四个点构成的交比（相同的交比可以根据点的排列顺序不同，有六种不同取值），利用交比不变性，可以更好的选择点之间的配对关系

wayne

我目前计算射影变换矩阵出了问题，就是选定的变换后的目标四个点的坐标不对，正在排查问题，不知道我理解的对不对：
1）给定三个点的射影坐标，其三点 共线应该就是判定矩阵$M_{3\times3}$的秩为2吧，如果是$1$就是存在点重合。
2）然后四个点确定一个射影变换是不是就是按照https://blog.emath.ac.cn/orchard-planting-problem/#i-3 所列的$12\times 13$矩阵处理。

mathe

https://bbs.emath.ac.cn/forum.ph ... d=703&pid=82731
的代码可以计算出所有的自同构群。
然后对原图中任意选择四个没有三点共线的点A,B,C,D
我们可以找出所有的自同构群以后，对于每种自同构（包括所有点不变),A->A',B->B',C->C',D->D'
假设存在一个射影变换，变换前所有点使用第一组坐标，变换后点A',B',C',D'会变换为第二组坐标中的A,B,C,D,然后查看是否所有其它点的坐标也都匹配。
这样计算的算法会效率很高

以https://bbs.emath.ac.cn/forum.ph ... =3953&pid=82288 中给出的
ADGJBEIJCDHKAFIKCEGLBFHLCJMODINODLMPAHNPGKOPBGMQFJNQAEOQEHMRBKNRCFPRILQRABCSDEFSGHITJKLTMNST
为例子，有8个自同构
ABCDEFGHIJKLMNOPQRST 00 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 20 21 22
GHIJKLABCDEFNMOQPRTS 00 02 01 04 03 05 07 06 12 11 13 09 08 10 15 14 17 16 20 21 18 19 22
EDFBACIHGJLKNMQROPST 01 00 05 04 03 02 12 11 15 14 17 07 06 13 09 08 16 10 19 18 20 21 22
CBADFEKJLHGINMPORQST 02 05 00 04 03 01 09 08 07 06 10 15 14 16 12 11 13 17 18 19 21 20 22
IHGJLKEDFBACMNQORPTS 01 05 00 03 04 02 11 12 06 07 13 14 15 17 08 09 10 16 20 21 19 18 22
KJLHGICBADFEMNPROQTS 02 00 05 03 04 01 08 09 14 15 16 06 07 10 11 12 17 13 21 20 18 19 22
FDEBCALJKHIGMNRQPOST 05 02 01 03 04 00 14 15 11 12 17 08 09 16 06 07 13 10 19 18 21 20 22
LJKHIGFDEBCANMRPQOTS 05 01 02 04 03 00 15 14 09 08 16 12 11 17 07 06 10 13 21 20 19 18 22
解满足参数-1+4*t+1*t^2=0， 所以用-1/t替换t可以得到另外一组解。
        A[+1 ,0 , 0]
        B[0 ,+1 , +1]
        C[-1 ,+1 , +1]
        D[+1 ,-1/2-1/2*t , 0]
        E[0 ,+1/2+1/2*t , +1]
        F[+1 ,0 , +1]
        G[+1 ,-1/2+1/2*t , 0]
        H[+3/2+1/2*t ,-1/2-1/2*t , +1]
        I[0 ,0 , +1]
        J[0 ,+1 , 0]
        K[+1/2+1/2*t ,0 , +1]
        L[+1/2+1/2*t ,+1/2-1/2*t , +1]
        M[-1 ,+3/2-1/2*t , +1]
        N[+1 ,-1/2-1/2*t , +1]
        O[-1 ,+1/2+1/2*t , +1]
        P[+2+1*t ,-1/2-1/2*t , +1]
        Q[+1 ,+1/2+1/2*t , +1]
        R[+2/5+1/5*t ,+3/10-1/10*t , +1]
        S[-1/2-1/2*t ,+1 , +1]
        T[+1/2+1/2*t ,-1*t , +1]

比如选择四个点为ABCD, 然后我们需要穷举所有的8中自同构，比如其中第二个自同构变换A->G,B->H,C->I,D->J
于是选择A,B,C,D的新坐标为[1,0,0], [0,1,1], [-1,1,1], [1,-1/2+1/(2t),0], 分别变换为G,H,I,J的旧坐标[1,-1/2+1/2t,0], [3/2+1/2 t,-1/2-1/2 t,1], [0,0,1], [0,1,0]
求出射影变换然后检查其它点的坐标是否符合变换要求即可（也就是可以全部从新坐标变换为旧坐标了）

mathe

上面有一点选错了，我们应该选择A,B,D,E四个点，因为A,B,C三点共线。
现在利用 https://bbs.emath.ac.cn/forum.ph ... =3953&pid=82399 的代码计算，使用
恒等变换 ABCDEFGHIJKLMNOPQRST
于是使用代码

1.         AA=[+1 ,0 , 0];
   
2.         BB=[0 ,+1 , +1];
   
3.         CC=[-1 ,+1 , +1];
   
4.         DD=[+1 ,-1/2-1/2*t , 0];
   
5.         EE=[0 ,+1/2+1/2*t , +1];
   
6.         FF=[+1 ,0 , +1];
   
7.         GG=[+1 ,-1/2+1/2*t , 0];
   
8.         HH=[+3/2+1/2*t ,-1/2-1/2*t , +1];
   
9.         II=[0 ,0 , +1];
   
10.         JJ=[0 ,+1 , 0];
    
11.         KK=[+1/2+1/2*t ,0 , +1];
    
12.         LL=[+1/2+1/2*t ,+1/2-1/2*t , +1];
    
13.         MM=[-1 ,+3/2-1/2*t , +1];
    
14.         NN=[+1 ,-1/2-1/2*t , +1];
    
15.         OO=[-1 ,+1/2+1/2*t , +1];
    
16.         PP=[+2+1*t ,-1/2-1/2*t , +1];
    
17.         QQ=[+1 ,+1/2+1/2*t , +1];
    
18.         RR=[+2/5+1/5*t ,+3/10-1/10*t , +1];
    
19.         SS=[-1/2-1/2*t ,+1 , +1];
    
20.         TT=[+1/2+1/2*t ,-1*t , +1];
    
21. modp = t^2+4*t-1;
    
22. 
    
23. Q=Mod(subst(AA,t, -4-t),modp);
    
24. V=Mod(subst(BB,t, -4-t), modp);
    
25. U=Mod(subst(DD,t, -4-t), modp);
    
26. N=Mod(subst(EE,t, -4-t), modp);
    
27. 
    
28. T=matrix(12,13);
    
29. for(u=1,3, T[1,u]=Q[u]); T[1,10]=Mod(AA[1],modp);
    
30. for(u=1,3,T[2,u+3]=Q[u]); T[2,10]=Mod(AA[2],modp);
    
31. for(u=1,3,T[3,u+6]=Q[u]); T[3,10]=Mod(AA[3],modp);
    
32. for(u=1,3, T[4,u]=V[u]); T[4,11]=Mod(BB[1],modp);
    
33. for(u=1,3,T[5,u+3]=V[u]); T[5,11]=Mod(BB[2],modp);
    
34. for(u=1,3,T[6,u+6]=V[u]); T[6,11]=Mod(BB[3],modp);
    
35. for(u=1,3, T[7,u]=U[u]); T[7,12]=Mod(DD[1],modp);
    
36. for(u=1,3,T[8,u+3]=U[u]); T[8,12]=Mod(DD[2],modp);
    
37. for(u=1,3,T[9,u+6]=U[u]); T[9,12]=Mod(DD[3],modp);
    
38. for(u=1,3, T[10,u]=N[u]); T[10,13]=Mod(EE[1],modp);
    
39. for(u=1,3,T[11,u+3]=N[u]); T[11,13]=Mod(EE[2],modp);
    
40. for(u=1,3,T[12,u+6]=N[u]); T[12,13]=Mod(EE[3],modp);
    
41. R=matker(T)~[1,1..9];
    
42. S=matrix(3,3);
    
43. for(u=1,3,S[1,u]=R[u];S[2,u]=R[u+3];S[3,u]=R[u+6]);
    
44. S

复制代码

求得变换阵S,并且计算S*II~结果和II坐标不同，淘汰。

于是我们继续采用第二种同构变换
GHIJKLABCDEFNMOQPRTS
即A->G,B->H,D->J, E->K
变换后S*II~同变换前的C不同，还是淘汰

然后采用同构
EDFBACIHGJLKNMQROPST
对应代码

1.         AA=[+1 ,0 , 0];
   
2.         BB=[0 ,+1 , +1];
   
3.         CC=[-1 ,+1 , +1];
   
4.         DD=[+1 ,-1/2-1/2*t , 0];
   
5.         EE=[0 ,+1/2+1/2*t , +1];
   
6.         FF=[+1 ,0 , +1];
   
7.         GG=[+1 ,-1/2+1/2*t , 0];
   
8.         HH=[+3/2+1/2*t ,-1/2-1/2*t , +1];
   
9.         II=[0 ,0 , +1];
   
10.         JJ=[0 ,+1 , 0];
    
11.         KK=[+1/2+1/2*t ,0 , +1];
    
12.         LL=[+1/2+1/2*t ,+1/2-1/2*t , +1];
    
13.         MM=[-1 ,+3/2-1/2*t , +1];
    
14.         NN=[+1 ,-1/2-1/2*t , +1];
    
15.         OO=[-1 ,+1/2+1/2*t , +1];
    
16.         PP=[+2+1*t ,-1/2-1/2*t , +1];
    
17.         QQ=[+1 ,+1/2+1/2*t , +1];
    
18.         RR=[+2/5+1/5*t ,+3/10-1/10*t , +1];
    
19.         SS=[-1/2-1/2*t ,+1 , +1];
    
20.         TT=[+1/2+1/2*t ,-1*t , +1];
    
21. modp = t^2+4*t-1;
    
22. 
    
23. Q=Mod(subst(EE,t, -4-t),modp);
    
24. V=Mod(subst(DD,t, -4-t), modp);
    
25. U=Mod(subst(BB,t, -4-t), modp);
    
26. N=Mod(subst(AA,t, -4-t), modp);
    
27. 
    
28. T=matrix(12,13);
    
29. for(u=1,3, T[1,u]=Q[u]); T[1,10]=Mod(AA[1],modp);
    
30. for(u=1,3,T[2,u+3]=Q[u]); T[2,10]=Mod(AA[2],modp);
    
31. for(u=1,3,T[3,u+6]=Q[u]); T[3,10]=Mod(AA[3],modp);
    
32. for(u=1,3, T[4,u]=V[u]); T[4,11]=Mod(BB[1],modp);
    
33. for(u=1,3,T[5,u+3]=V[u]); T[5,11]=Mod(BB[2],modp);
    
34. for(u=1,3,T[6,u+6]=V[u]); T[6,11]=Mod(BB[3],modp);
    
35. for(u=1,3, T[7,u]=U[u]); T[7,12]=Mod(DD[1],modp);
    
36. for(u=1,3,T[8,u+3]=U[u]); T[8,12]=Mod(DD[2],modp);
    
37. for(u=1,3,T[9,u+6]=U[u]); T[9,12]=Mod(DD[3],modp);
    
38. for(u=1,3, T[10,u]=N[u]); T[10,13]=Mod(EE[1],modp);
    
39. for(u=1,3,T[11,u+3]=N[u]); T[11,13]=Mod(EE[2],modp);
    
40. for(u=1,3,T[12,u+6]=N[u]); T[12,13]=Mod(EE[3],modp);
    
41. R=matker(T)~[1,1..9];
    
42. S=matrix(3,3);
    
43. for(u=1,3,S[1,u]=R[u];S[2,u]=R[u+3];S[3,u]=R[u+6]);
    
44. S;
    
45. S*II~

复制代码

结果S*II~的结果为[Mod(1, t^2 + 4*t - 1), Mod(1/2*t - 1/2, t^2 + 4*t - 1), Mod(0, t^2 + 4*t - 1)]～也就是G的坐标，正好这个变换G->I,完全匹配。
由此我们找到两者之间的变换阵为这时的S, 而且证明了方程两个根对应的解之间可以通过射影变换相互转化，所以它们等价

mathe

发现还有一个问题，上面是以二次方程为例子，比如方程t^2+4t-1=0,那么给定一个根t以后，那么另外一个根根据韦达定理必然是-1/t或-4-t，可以轻松用一个根的多项式来表达另外一个根，但是如果对于三次或以上的方程，那么这种用一个根的多项式精确的表达另外一个根就没有那么容易了。
比如124#引用链接中另外一个配置对应的方程为三次方程1-t-2t^2+t^3=0，于是给定方程的一个根t, 那么在扩张域Q[t]中，我们可以先计算多项式
(X^3-2X^2-X+1)-(t^3-2t^2-t+1)=X^2+(t-2)X+(t-1)^2=0, 于是另外两个根就是这个方程在Q[t]中的根，pari/gp好像没有做这种因子分解的能力，不知道Mathematica能否在
有限域Q[t]中分解上面多项式X^2+(t-2)X+(t-1)^2. 另外我们还可以根据求根公式来求上面的X,于是关键要对判别式(t-2)^2-4(t-1)^2=(t-2+2t-2)(t-2-2t+2)=-t(3t-4)求平方根，pari/gp同样无法处理这种情况，会给出一个神奇的不知道意义的结果：
%44 = Mod(-3*t^2 + 4*t, t^3 - 2*t^2 - t + 1)
? sqrt(%44)
%45 = [2.2665090645181026229368720711415588721*I, 1.1383746936926226201220102497162464711, 2.4816997408423621980266295382639513129*I]~
不知道这种有理数域的扩张域里面该如何求开平方。

另外一种可行的取巧方案是可以在取一个大素数模p, 比如$p=12373$ （注意选择的素数要充分大，还要保持 t^3 - 2*t^2 - t + 1无法分解）,然后让上面多项式在模p然后在有限域内开平方，这种方法我们知道是存在比较有效的算法。
然后对于开平方结果（这时得到一个t的二次方程，系数都在p阶有限域内），每个系数去尝试和一个分子分母都不是太大的分数去匹配，通常会有很大的概率得到正确的解。
但是我在pari/gp中试验了一下，发现还是无法计算（可能我没有找到合适的函数）
? Mod(Mod(1, 12373)*t^2 + (Mod(1, 12373)*X + Mod(12371, 12373))*t + (Mod(1, 12373)*X^2 + Mod(12371, 12373)*X + Mod(12372, 12373)), Mod(1,12373)*t^3-Mod(2,12373)*t^2-Mod(1,12373)*t+Mod(1,12373))
%63 = Mod(Mod(1, 12373)*t^2 + (Mod(1, 12373)*X + Mod(12371, 12373))*t + (Mod(1, 12373)*X^2 + Mod(12371, 12373)*X + Mod(12372, 12373)), Mod(1, 12373)*t^3 + Mod(12371, 12373)*t^2 + Mod(12372, 12373)*t + Mod(1, 12373))
? sqrt(%63)
  ***   at top-level: sqrt(%63)
  ***                 ^---------
  *** sqrt: incorrect type in roots (t_INTMOD).
  ***   Break loop: type 'break' to go back to GP prompt

mathe

啊哈，果然是方法不对，至少最后一种方法是可用的
? t=ffgen(x^3-2*x^2-x+Mod(1,12373))
%1 = x
? X^2+(t-2)*X+(t-1)^2
%2 = X^2 + (x + 12371)*X + (x^2 + 12371*x + 1)
? factor(%2)
%3 =
[X + (6039*x^2 + 6776*x + 8922) 1]

[X + (6334*x^2 + 5598*x + 3449) 1]

mathe

上面需要分解的多项式是$X^2+(t-2)*X+t^2-2t-1$，前面手工计算错了，所以结果改为
[X + (x^2 + 12372*x + 12371) 1]
[      X + (12372*x^2 + 2*x) 1]
也就是另外两根分别为 (-t^2+t+2)和(t^2-2t)
数值检验方程的三个根-0.80193773580483825247220463901489010233，0.55495813208737119142219487100641048107，2.2469796037174670610500097680084796213
也得可以通过上面两个表达式相互转换，所以计算没错。

然后经计算得出对于变换
BCAEFDIJKLGHMPNORSQT
可以达到目标（将根t转为为-t^2+t+2),变换阵为S=
[   Mod(-2*t^2 + 1, t^3 - 2*t^2 - t + 1)              Mod(0, t^3 - 2*t^2 - t + 1)   Mod(0, t^3 - 2*t^2 - t + 1)]
[Mod(4*t^2 + t - 2, t^3 - 2*t^2 - t + 1) Mod(-2*t^2 - t + 1, t^3 - 2*t^2 - t + 1)   Mod(0, t^3 - 2*t^2 - t + 1)]
[         Mod(-t^2, t^3 - 2*t^2 - t + 1)           Mod(-t^2, t^3 - 2*t^2 - t + 1) Mod(t^2, t^3 - 2*t^2 - t + 1)]

数学星空

我们可以利用待定系数求解：

设\(x^3+ax^2+bx+c=0\)的三个根分别为\(t,m_2t^2+m_1t+m_0,n_2t^2+n_1t+n_0\)

显然有\(x^3+ax^2+bx+c=(x-t)(-x+m_2t^2+m_1t+m_0)(-x+n_2t^2+n_1t+n_0)\)  （1）

利用\(t^3+at^2+bt+c=0\)将（1）关于t的次幂大于等于3全部化为低于3的次幂，令关于x的次幂项系数为0

\((-a^3m_2n_2 + a^2m_1n_2 + a^2m_2n_1 + 2abm_2n_2 - am_0n_2 - am_1n_1 - am_2n_0 - bm_1n_2 - bm_2n_1 - cm_2n_2 + m_0n_1 + m_1n_0)t^2 + (-a^2bm_2n_2 + abm_1n_2 + abm_2n_1 + acm_2n_2 + b^2m_2n_2 - bm_0n_2 - bm_1n_1 - bm_2n_0 - cm_1n_2 - cm_2n_1 + m_0n_0)t - a^2cm_2n_2 + acm_1n_2 + acm_2n_1 + bcm_2n_2 - cm_0n_2 - cm_1n_1 - cm_2n_0 + c=0\)

\((-a^2m_2n_2 + am_1n_2 + am_2n_1 + bm_2n_2 + am_2 + an_2 - m_0n_2 - m_1n_1 - m_2n_0 - m_1 - n_1)t^2 + (-abm_2n_2 + bm_1n_2 + bm_2n_1 + cm_2n_2 + bm_2 + bn_2 - m_0n_1 - m_1n_0 - m_0 - n_0)t - acm_2n_2 + cm_1n_2 + cm_2n_1 + cm_2 + cn_2 - m_0n_0 + b=0\)

\((m_2 + n_2)t^2 + (m_1 + n_1 + 1)t + a + m_0 + n_0=0\)

再分别令上面关于t的次幂项系数为0得到：

[-a^2*c*m2*n2 + a*c*m1*n2 + a*c*m2*n1 + b*c*m2*n2 - c*m0*n2 - c*m1*n1 - c*m2*n0 + c, -a^2*b*m2*n2 + a*b*m1*n2 + a*b*m2*n1 + a*c*m2*n2 + b^2*m2*n2 - b*m0*n2 - b*m1*n1 - b*m2*n0 - c*m1*n2 - c*m2*n1 + m0*n0, -a^3*m2*n2 + a^2*m1*n2 + a^2*m2*n1 + 2*a*b*m2*n2 - a*m0*n2 - a*m1*n1 - a*m2*n0 - b*m1*n2 - b*m2*n1 - c*m2*n2 + m0*n1 + m1*n0, -a*c*m2*n2 + c*m1*n2 + c*m2*n1 + c*m2 + c*n2 - m0*n0 + b, -a*b*m2*n2 + b*m1*n2 + b*m2*n1 + c*m2*n2 + b*m2 + b*n2 - m0*n1 - m1*n0 - m0 - n0, -a^2*m2*n2 + a*m1*n2 + a*m2*n1 + b*m2*n2 + a*m2 + a*n2 - m0*n2 - m1*n1 - m2*n0 - m1 - n1, a + m0 + n0, m1 + n1 + 1, m2 + n2]

求解上面方程得到：

下面记  \(t_0^2=-\frac{1}{4a^3c - a^2b^2 - 18abc + 4b^3 + 27c^2}\)

\(m_0=\frac{(4a^4ct_0 - a^3b^2t_0- 18a^2bct_0 + 4ab^3t_0+ 27ac^2t_0 - a^2b - 3ac + 4b^2)t_0}{2}\)

\(m_1=\frac{-1+(-2a^3 + 7ab - 9c)t_0}{2}\)

\(m_2=-t_0(a^2 - 3b)\)

\(n_0=\frac{(4a^4ct_0 - a^3b^2t_0 - 18a^2bct_0+ 4ab^3t_0 + 27ac^2t_0 + a^2b + 3ac - 4b^2)t_0}{2}\)

\(n_1=\frac{-1+(2a^3 - 7ab + 9c)t_0}{2}\)

\(n_2=t_0(a^2 - 3b)\)

例：对于\(x^3-2x^2-x+1=0\)即\(a=-2,b=-1,c=1\)

可以得到：\({m_0 = 2, m_1 = 1, m_2 = -1, n_0 = 0, n_1 = -2, n_2 = 1,t_0=\frac{1}{7}}\)

mathe

对于每行3棵树的果树问题，大部分最优解都可以让所有点落在一条三次曲线上。
而对于每行四棵树的果树问题，我们到现在为止看到所有最优解对应的代数方程都不超过四次，不知道是否所有这些最优解都可以让所有点落在某条四次曲线上？
参考数据: https://blog.emath.ac.cn/shared/