print("D Broadhurst, KP combined test, 8 Jul 2001");
print("Hard coded by TJW for Phi(13681,26371) - please hand edit!");

\p7000
Phi(a)=subst(polcyclo(a),x,26371);
N=(26371^13681-1)/26370;
F=Phi(13680)*Phi(570)*Phi(240)*Phi(228)*Phi(190)*Phi(144)*Phi(120)*Phi(114)*Phi(90)*Phi(80)*Phi(76)*Phi(72)*Phi(60)*Phi(57)*Phi(45);
G=2;
R=(N-1)/F;
S=(N+1)/G;

{if(denominator(R)==1
  &&denominator(S)==1
  &&gcd(F,R)==1
  &&gcd(F,G)==1
  &&F^(1/3)/6>G-1
  &&G-1>3*N/F^(10/3),
  ok=1,
  print("Fail sanity check");ok=0);}

c1=R%F;
c4=(R-c1)/F;
tG=c4%G;
print(tG);

print();

if(issquare((c1+tG*F)^2+4*(tG-c4)),print("Fail "tG);ok=0);

b=contfrac(c1/F);
v=0;vn=1;u=1;un=b[1];n=1;
{while(F^(1/3)>vn,n=n+1;bn=b[n];uo=u;u=un;vo=v;v=vn;
 un=bn*un+uo;vn=bn*vn+vo);}

d=round(c4*v/F);
zs=[v,u*F-c1*v,c4*v-d*F+u,-d];
q=0;
for(k=1,4,print(zs[k]);q=q+zs[k]*x^(4-k));

print();

ps=polroots(q);

{for(k=1,3,r=floor(real(ps[k]));print(r);
 x=r;q1=eval(q);x=r+1;q2=eval(q);
 if(q1*q2>=0,print("Fail root ",k);ok=0));}

if(ok==1,print("OK"),print("not yet"));
quit();