print("D Broadhurst, KP combined test, 8 Jul 2001");
print("Hard coded for Fib(25561) - please hand edit!");

N=fibonacci(25561);
F=prod(k=1,length(lsm),lsm[k]);
G=prod(k=1,length(lsp),lsp[k]);
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");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();

\p1500
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"));