/*
*    DIOF
*
*    Library "DIOF"
*
*    Author: 	Rosario Turco 
*    Copyright:	(C) 2007-2020
*    Version: 	2.1                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                         .1
*    Date:	28/06/10
*
*    Product: 	PARI/GP  version 2.4.2
*
*    Blog:	http://MATHBuildingBlock.blogspot.com
*
*    Site:	www.gruppoeratostene.com
*    Papers: 	Database CNR Solar
*
*
*   
*    SOFTWARE WITHOUT ANY WARRANTY WHATSOEVER.
*
*
*/




debug=0;


/*
** Eq. diofantea di primo grado
**
*/

DiofPG(a,b,c) = local(); { 
 if( Mod(c,gcd(a,b))!=0, error("Equazione diofantea NON risolvibile!"););
 v = bezout(a,b);
 print("Equazione Diofantea di primo grado ",a,"*x+",b,"*y=",c, " risolvibile");
 print("Soluzione :");
 print("x = ",c*v[1],"+",b,"*t");
 print("y = ",c*v[2],"-",a,"*t");
 return;
}

/*
*  Congruenza Lineare
*/
EqMod(a,b,n) = local(k=0,x=0,ciclo=1); { 
 if( b%(gcd(a,n))!= 0, error("Congruenza lineare non risolvibile..."););
 print("Congruenza lineare ",a,"*x=",b,"mod",n, " risolvibile");
 while(ciclo==1,
       k=k+1;
	   if( debug, print(k,"\t",b+k*n,"\t",(b+k*n)%a););
	   if((b+k*n)%a==0, ciclo=0;);
 );
 x = (b+k*n)/a;
 print("Soluzione x : ", x);
 return(x);
}

/*
*  Sistema di due congruenze lineari
*/
Sist2EqMod(a,b,n1,c,d,n2)= local(x=0,y=0,z=0); { 
 if( gcd(n1,n2)!= 1, error("Teor. cinese del resto: sistema di 2 congruenze lineare non risolvibile..."););
 print("Sistema di 2 congruenze lineari :");
 print(a,"*x=",b," mod",n1);
 print(c,"*x=",d," mod",n2);
 print("da cui : ");
 x = EqMod(a,b,n1);
 print("x=",x," mod",n1);
 ciclo=1;
 y = EqMod(c,d,n2);
 print("x=",y," mod",n2);
 N=n1*n2;
 bz=bezout(n1,n2);
 print("bezout : ", bz);
 z=x*bz[2]*n2+y*bz[1]*n1;
 print("Soluzione del sistema x : ", z," mod",N);
}

/*
*  Pell's Equation x^2-dy^2=+1
*  simple example with contfrac funtion  
*
*/

/*
*  We must increase the precision to achieve a major expansion of contfrac
*  with \p <number> and allocatemem(x*1024*1024) if we obtain an error
*  on the index of the array or on the dimension of the stack 
*/
   
DiofPell(d,c=1) = local(di=0,k=0); { 

 if( c==-1, if( error("Equazione di Pell NON trattata...");););
 print("Equazione di Pell x^2-dy^2=1 risolvibile");

 vf = contfrac(sqrt(d));
 di=getdimfrac(vf,c);
 
 if(debug,print("dim : ", di););
 
 p=vector(3);
 q=vector(3);
 
 p[1]=1; p[2]=vf[1];q[1]=0; q[2]=1;
 if(debug, print("p[",1,"]=",p[1],"\tq[",1,"]=",q[1]););
 if(debug, print("p[",2,"]=",p[2],"\tq[",2,"]=",q[2]);); 

 precind=2;  
 
 for(i=3,di+1,
 
   k=getIndex(precind);
   if(debug,print("k : ", k););

   precind=k;
   if(k==1,a1=3;a2=2;);
   if(k==2,a1=1;a2=3;);
   if(k==3,a1=2;a2=1;);
 
   p[k]=vf[i-1]*p[a1]+p[a2];    
   q[k]=vf[i-1]*q[a1]+q[a2];

   if( debug, print("vf[",i-1,"]=",vf[i-1]););
   if(debug, print("p[",k,"]=",p[k],"\tq[",k,"]=",q[k]););
 );
 print("x = ", p[k]);
 print("y = ", q[k]);
 print(p[k],"^2-",d,"*",q[k],"^2=",p[k]^2-d*(q[k])^2);
 return;
}

getdimfrac(v,c) = local(dim=0,da1=0,dda1=0); { 
 if(debug, print("v : ", v););
 len = length(v);
 if(debug, print("len : ", len););
 da1 = 2*v[1];
 dda1 = 2*da1;
 if(debug, print("da1 da cercare: ", da1););
 
 for(i=2,len,
     if(debug, print("i : ", i););
     if( v[i] == da1,
	     if(debug,print("da1 trovato"););
	     if( (i-1)%2 == 0 & i<len & c==1,
            dim=i-1;		 
 	        if(debug,print("dim (c=1): ", i-1););		
            i=len;			
		 );
		 if( (i-1)%2 != 0 & i<len & c==1,
 	        if(debug,print("proseguo al prossimo da1 (c=1)"););			
			if( dda1 > len, error("aumentare la precisione a 300 o maggiore e rieseguire"););
		 );
	 );   
 );
 return(dim);
}

getIndex(prec)= local(); { 
 if(debug, print(" prec : ", prec););
 prec++;
 if( prec>3, prec=1;);
 return(prec);
}