/*   Approximation of a real by a rational of the form a/b^2. 
*/






 /** Calcula todas las raices de $x^2 \equiv a$ \pmod m menores o iguales a m/2.
    m puede ser compuesto. 10 Mayo 2000. Revisada 16 Febrero 2006       **/
global(f); /* en f esta la factorizacion de Q */

{
sqrtmod(a,m)=
local(l,p,e,ii,s,aa,ll,pp,ppp,pe,pes,ra,ai,r=[],rr,r1,r2=[0,0],j2,nada=[]);
if(m==1,return([0]));

a%=m;
\\f=factor(m);
l=matsize(f)[1];

for(i=1,l,p=f[i,1];e=f[i,2];      /* a=aa*p^s, calculo de raices mod p^e */
  pe=p^e;aa=a%pe;s=0;if(aa!=0,while(aa%p==0,s++;aa/=p));

  /* caso 1, aa=0 */
  if(aa==0,
    pp=p^floor((e+1)/2);ll=p^(e\2);rr=vector(ll);
    ppp=0;for(j=1,ll,rr[j]=Mod(ppp,pe);ppp+=pp)
  ,

    /* caso 2, s impar */
    if(s%2==1,return(nada));  

    /* caso 3, s par */

    /* p impar */
    if(p!=2,
      if(issquare(Mod(aa,p)),ra=lift(sqrt(Mod(aa,p))),return(nada));
      /* Si hay raices, 'Hensel lift' */
      ll=p^(s/2);rr=vector(2*ll);ai=bezout(2*ra,p)[1];pes=p^(e-s);
      for(j=2,e-s,ra=(ra-ai*(ra*ra-aa))%pes );ra*=p^(s/2);
      /* las raices de a son +/- ra*p^(s/2) mod p^(e-s) */
      ppp=p^(e-s/2);for(j=1,ll,rr[j]=Mod(ra,pe);rr[j+ll]=Mod(-ra,pe);ra+=ppp)
    ,

      /* p = 2 */

      if(aa%8!=1,return(nada)); /* solo hay raices si aa=1 mod 8 */

      if((e-s)==1,r1=[1],
        if((e-s)==2,r1=[1,3],

          if((e-s)>=3,r1=[1,3,5,7];j2=2^3;

            for(j=4,e-s,j2*=2;           /* Lift a mod 2^e */
              ii=1;for(k=1,4,if(r1[k]*r1[k]%j2==aa%j2,r2[ii]=r1[k];ii++));
              r1=[r2[1],r2[1]+2^(j-1),r2[2],r2[2]+2^(j-1)];
   
            )
          )
        )
      );
    
      /* a partir de r1 preparamos rr */
      if(s==0,rr=Mod(r1,pe)
      ,
      
        ll=p^(s/2);ppp=p^(e-s/2);
        if(ll>1,
          if((e-s)==1,ra=p^(s/2);rr=vector(ll);
            for(j=1,ll,rr[j]=Mod(ra,pe);ra+=ppp)
          ,
            if((e-s)==2,
              rr=vector(2*ll);for(j=1,2,r1[j]*=p^(s/2));
              for(j=1,ll,rr[j]=Mod(r1[1],pe);rr[j+ll]=Mod(r1[2],pe);for(k=1,2,r1[k]+=ppp) );
            ,
              if((e-s)>=3,
                rr=vector(4*ll);for(j=1,4,r1[j]*=p^(s/2));
                for(j=1,ll,rr[j]=Mod(r1[1],pe);rr[j+ll]=Mod(r1[2],pe);
                  rr[j+2*ll]=Mod(r1[3],pe);rr[j+3*ll]=Mod(r1[4],pe);for(k=1,4,r1[k]+=ppp))
              )
            )
          )
        ) 
      )
    )  
  );
  r=concat(r,[rr]);
);


/* en r estan las colecciones de raices mod p_i */

/** resto chino de todas las posibles combinaciones **/
l=#r;e=1;for(i=1,l,e*=#r[i]);rr=vector(e); /* En rr ponemos las raices mod m */
e=1;
forvec(X=vector(l,X,[1,#r[X]]),
  rr[e]=lift(chinese(vector(l,Y,r[Y][X[Y]])));e++;
);

/* Se toman solo las raices de 0 a m/2 */

\\return(vecextract(vecsort(rr),concat("1..",Str( (length(rr)+1)\2 + (a==0 & m%4==0)  ) ) ) )

return(rr);

}

sqrtmod2(a,m)=R=[];for(i=0,m\2,if(i^2%m==a%m,R=concat(R,i)));R;









{
approxab2(x,appini=1,appfin=50)=
local(b,e,teta,iter,p,p1,p2,q,q1,q2,pi,R,ep,aa,bb,dat,T=[],brange);




write("approxab2.dat","{[");
gettime();

/* Primeros terminos por fuerza bruta si appini==1 */

print("First terms by brute force ..."); 
print("   b     error=b^- ");
if(appini==1,
  for(b=2,10^3,aa=round(x*b*b);e=abs(x-aa/b^2);teta=-log(e)/log(b);
    if(teta>=3,
      dat=[b,teta];
      print(dat);write1("approxab2.dat",dat);write("approxab2.dat",",");
    );
  );
);




 /* Desarrollo de la FC de x */

b=floor(x);e=1./x; /* condiciones iniciales */
p1=b;p2=1;q1=1;q2=0;



iter=0;
print;print(iter" Convergent "p1"/"q1);
until(iter == appfin,

  /* siguiente convergente */
  e=1./e-b;
  b=floor(1./e);
  p=b*p1+p2;p2=p1;p1=p;
  q=b*q1+q2;q2=q1;q1=q;
  iter++;
  print;print(iter" Convergent "p"/"q);


  brange=ceil(q^.75); 

/* Busqueda del aproximante a/b^2 a partir de p/q */
  
  if(iter>=appini & q > 1,

    print("   b     error=b^-");
    pi=bezout(p,q)[1]; /* inverso de p mod q */
    f=factor(q);
    d=0;
    alfaini=1;alfafin=round(q^.35);
    for(alfa=alfaini,alfafin,
      d=pi*alfa%q;
      R=sqrtmod(d,q);
      for(j=1,#R,
        bb=R[j];
        if(bb<=brange,
          aa=round(bb*bb*x);
          if(bb!=1,
            ep=abs(x-aa/bb^2);teta=-log(ep)/log(bb);
            if(teta >= 3., 
              dat=[bb,teta];
              print(dat);write1("approxab2.dat",dat);write("approxab2.dat",",");
            );
          );
        );
      );

      R=sqrtmod(-d,q);
      for(j=1,#R,
        bb=R[j];
        if(bb<=brange,        
          aa=round(bb*bb*x);
          if(bb!=1,
            ep=abs(x-aa/bb^2);teta=-log(ep)/log(bb);
            if(teta >= 3., 
              dat=[bb,teta];
              print(dat);write1("approxab2.dat",dat);write("approxab2.dat",",");
            );
          );
        );
      );

    );    

  );

); 
}

{
print("\nUsage:  \"approxab2(x,ini=1,end=50)\".\n
       By example, \"approxab2(Pi,1,20)\".\n
       Data generated are written to file \"approxab2.dat\". \n");
 
}