Последний раз редактировалось dmd 28.09.2020, 09:58, всего редактировалось 1 раз.
Если кому будет интересно, предлагаю поучаствовать в шлифовке gp-кода решателя биквадратных уравнений по алгоритму из статьи maxal. (pari/gp код)
Код: partsol(A,B,C)= { my(N, no, xo, yo, zo, k0, k1, l, S= []);
N= iferr(bnfisnorm(bnfinit('x^2 + A*B, 1), -A*C), Err, 0); if(N, if(abs(N[2])==1, no= lift(N[1]); k0= polcoef(no, 0); k1= polcoef(no, 1); if(k0 & k1, xo= k0/A; yo= k1; zo= 1; l= lcm(denominator(xo), denominator(yo)); xo= abs(l*xo); yo= abs(l*yo); zo= abs(l*zo); if(A*xo^2+B*yo^2+C*zo^2==0, S= concat(S, [[xo,yo,zo]])) ) ));
N= iferr(bnfisnorm(bnfinit('x^2 + A*C, 1), -A*B), Err, 0); if(N, if(abs(N[2])==1, no= lift(N[1]); k0= polcoef(no, 0); k1= polcoef(no, 1); if(k0 & k1, xo= k0/A; yo=1; zo= k1; l= lcm(denominator(xo), denominator(zo)); xo= abs(l*xo); yo= abs(l*yo); zo= abs(l*zo); if(A*xo^2+B*yo^2+C*zo^2==0, S= concat(S, [[xo,yo,zo]])) ) ));
N= iferr(bnfisnorm(bnfinit('x^2 + B*C, 1), -B*A), Err, 0); if(N, if(abs(N[2])==1, no= lift(N[1]); k0= polcoef(no, 0); k1= polcoef(no, 1); if(k0 & k1, yo= yo/B; xo= 1; zo= k1; l= lcm(denominator(zo), denominator(yo)); xo= abs(l*xo); yo= abs(l*yo); zo= abs(l*zo); if(A*xo^2+B*yo^2+C*zo^2==0, S= concat(S, [[xo,yo,zo]])) ) ));
N= iferr(bnfisnorm(bnfinit('x^2 + B*A, 1), -B*C), Err, 0); if(N, if(abs(N[2])==1, no= lift(N[1]); k0= polcoef(no, 0); k1= polcoef(no, 1); if(k0 & k1, yo= yo/B; zo= 1; xo= k1; l= lcm(denominator(xo), denominator(yo)); xo= abs(l*xo); yo= abs(l*yo); zo= abs(l*zo); if(A*xo^2+B*yo^2+C*zo^2==0, S= concat(S, [[xo,yo,zo]])) ) ));
N= iferr(bnfisnorm(bnfinit('x^2 + C*B, 1), -C*A), Err, 0); if(N, if(abs(N[2])==1, no= lift(N[1]); k0= polcoef(no, 0); k1= polcoef(no, 1); if(k0 & k1, zo= zo/C; yo= k1; xo= 1; l= lcm(denominator(zo), denominator(yo)); xo= abs(l*xo); yo= abs(l*yo); zo= abs(l*zo); if(A*xo^2+B*yo^2+C*zo^2==0, S= concat(S, [[xo,yo,zo]])) ) ));
N= iferr(bnfisnorm(bnfinit('x^2 + C*A, 1), -C*B), Err, 0); if(N, if(abs(N[2])==1, no= lift(N[1]); k0= polcoef(no, 0); k1= polcoef(no, 1); if(k0 & k1, zo= zo/C; yo= 1; xo= k1; l= lcm(denominator(xo), denominator(zo)); xo= abs(l*xo); yo= abs(l*yo); zo= abs(l*zo); if(A*xo^2+B*yo^2+C*zo^2==0, S= concat(S, [[xo,yo,zo]])) ) ));
return(S) };
biquadratic(P)= { SX= Set(); a= polcoef(P, 4); b= polcoef(P, 2); c= polcoef(P, 0); k3= polcoef(P, 3); k1= polcoef(P, 1);
if(k3||k1, print("Polynomial is not biquadratic, exit."); return());
if(b^2-4*a*c==0, print("Discriminant("P") = 0, exit."); return());
A= b^2-4*a*c; B= 4*c; C= -1;
\\ print("(A,B,C) = ("A", "B", "C")");
S1= partsol(A,B,C);
\\ print(S1);
for(i1=1, #S1, xo= S1[i1][1]; yo= S1[i1][2]; zo= S1[i1][3];
Q= abs(2*lcm(A,B)*C*zo);
\\ print("Q = "Q);
Q= divisors(Q);
Px= xo*A*m^2 + 2*yo*B*m*n - xo*B*n^2; Py= -yo*A*m^2 + 2*xo*A*m*n + yo*B*n^2; Pz= zo*A*m^2 + zo*B*n^2;
\\ print("Px = "Px"\nPy = "Py"\nPz = "Pz);
E= 2*c*Px/(Pz - b*Px);
T= subst(E, n, 1); P1= numerator(T); P2= denominator(T);
\\ print("P1 = "P1"\nP2 = "P2); a1= polcoef(P1, 2); b1= polcoef(P1, 1); c1= polcoef(P1, 0);
if(b1^2-4*a1*c1==0, print("Discriminant("P1") = 0, exit."); break());
a2= polcoef(P2, 2); b2= polcoef(P2, 1); c2= polcoef(P2, 0);
T= [a1, 0, a2, 0; b1, a1, b2, a2; c1, b1, c2, b2; 0, c1, 0, c2];
T1= matsolve(T, [1, 0, 0, 0]~); T2= matsolve(T, [0, 0, 0, 1]~);
G= 1; for(i=1, #T1, G= lcm(G, denominator(T1))); for(i=1, #T2, G= lcm(G, denominator(T2)));
\\ print("G = "G);
if(b1==0, r= a1; s= c1; t= G , if(c1!=0, s= 1/4/c1; r= (4*a1*c1-b1^2)*s; Q1= x; Q2= b1*x+2*c1*y; g2= gcd(polcoef(Q2, 1, x), polcoef(Q2, 1, y)); g1= 1; Q2= Q2/g2; s= s*g2^2; ); if(a1!=0, r= 1/4/a1; s= (4*a1*c1-b1^2)*r; Q1= 2*a1*x+b1*y; Q2= y; g1= gcd(polcoef(Q1, 1, x), polcoef(Q1, 1, y)); g2= 1; Q1= Q1/g1; r= r*g1^2; ); if(a1==0 & c1==0, r= s= b1/4; Q1= x+y; Q2= x-y); );
g12= lcm(denominator(r), denominator(s));
r= r*g12; s= s*g12; Q3= z;
D= divisors(G);
for(j=1, #D,
d= D[j];
if(issquarefree(d),
forstep(signg=-1, 1, 2,
g= signg*d;
\\ print("g = "g);
t= g*g12;
if(r<0 & s<0, r= -r; s= -s; t= -t);
\\ print("(r,s,t) = ("r", "s", "t")"); \\ print("(Q1,Q2) = ("Q1", "Q2")"); \\ print("Equation (6): "r" * ("Q1")^2 + "s" * ("Q2")^2 + "t" * ("Q3")^2 = 0");
S2= partsol(r,s,t);
\\ print(S1);
for(i2=1, #S2, ro= S2[i2][1]; so= S2[i2][2]; to= S2[i2][3];
Qx= ro*r*m^2 + 2*so*s*m*n - ro*s*n^2; Qy= -so*r*m^2 + 2*ro*r*m*n + so*s*n^2; Qz= to*r*m^2 + to*s*n^2;
\\ print("Qx = "Qx"\nQy = "Qy"\nQz = "Qz);
if(Q1=='x, Y= (Qy - polcoef(Q2, 1, 'y)*Qx)/polcoef(Q2, 1, 'x); X= Qx); if(Q2=='y, X= (Qx - polcoef(Q1, 1, 'y)*Qy)/polcoef(Q1, 1, 'x); Y= Qy);
\\ print("(X,Y) = ("X", "Y")");
PT= 1*(polcoef(P2, 2)*X^2 + polcoef(P2, 1)*X*Y + polcoef(P2, 0)*Y^2);
S= subst(PT, n, 1);
lg= lcm(denominator(polcoef(S, 4)), denominator(polcoef(S, 3))); lg= lcm(denominator(polcoef(S, 2)), lg); lg= lcm(denominator(polcoef(S, 1)), lg); lg= lcm(denominator(polcoef(S, 0)), lg);
\\ print("lg = "lg);
PT= lg*PT;
\\ print("Thue polynomial: "PT);
Glg= divisors(g*lg);
for(k=1, #Q,
q= Q[k];
for(v=1, #Glg,
p2= Glg[v];
if(gcd(q,p2)==1 & issquare(p2),
\\ print("Thue polynomial: "PT); \\ print("RHS of Thue: ("g" * "lg" * "q"^2 / "p2") = "g*lg*q^2/p2);
Th= thue(subst(PT, n, 1);, g*lg*q^2/p2);
\\ print("Solutions Thue equation:\n"Th);
MN= apply(V -> [subst(subst(X,m,V[1]),n,V[2]), subst(subst(Y,m,V[1]),n,V[2])], Th);
\\MN= select(V -> gcd(V[1],V[2])==1, MN);
\\ print("(m,n) = \n"MN);
X2= apply(V -> [subst(subst(E,m,V[1]),n,V[2])], MN);
for(i=1, #X2, SX= setunion(SX, [X2[i][1]]));
\\ print(SX);
) ) ) ) ) ) ) );
XY= []; for(i=1, #SX, if(issquare(SX[i]), xi= sqrtint(SX[i]); y2= a*xi^4+b*xi^2+c; if(issquare(y2), yi= sqrtint(y2); XY= concat(XY, [[xi,yi]]) ) ) );
\\ print("Solutions of "P" = y^2:\n"XY);
return(XY)
};
Пока у меня такие возникли проблемы. Частное нетривиальное решение уравнения получаю так: приравниваю например , преобразую к уравнению и решаю его при помощи bnfisnorm. Перебираю все возможные варианты такого преобразования, и всё равно при некоторых сочетаниях коэффициентов иногда решение не находится. Как более правильно находить при помощи функции bnfisnorm? Еще почему-то требуемое условие отсекает правильные решения, а без него они находятся. Сейчас пример biquadratic(2*x^4-1) правильно находит решения [[1, 1], [13, 239]], но на biquadratic(2*x^4-7) вычисления сильно подвисают. Если раскоментарить принты Код: \\ print("Thue polynomial: "PT); \\ print("RHS of Thue: ("g" * "lg" * "q"^2 / "p2") = "g*lg*q^2/p2);
то видно, что иногда Туе-полином получается невразумительно огромным. Видимо какие-то варианты надо исключить из рассмотрения, но не пойму какие. Для всех вычислений в коде есть соответствующий комментированный принт, убирая коменты можно отследить в деталях работу алгоритма.
|