Pedro Fortuny Ayuso on Thu, 02 Mar 2017 15:44:27 +0100


[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]

Re: Mathematica "Reduce" function


Dear Denis,

Thanks a lot. I have tried it and (for the moment) it seems to
be somewhat slower than what Bill suggested.

In any case, I guess I have now more information than what I
imagined I could have and am trying to digest it.

Thanks a lot again,


Pedro.

On Thu, Mar 02, 2017 at 11:37:58AM +0100, Denis Simon wrote:
> Dear Pedro,
> 
> Here is a small gp script that mixes exhaustive search and lifting.
> This should be not too bad if the prime number p is small (typically 2 or 3 as in your examples),
> and in any cases much better than naive exhaustive search.
> 
> Idea of the algo : exhaustive search mod p, then exhaustive lift mod p^2, mod p^3...
> 
> Sincerely,
> Denis SIMON.
> 
> 
> p = 3;  \\ your prime number
> k = 2;  \\ the exponent of p
> eqns = [ x^2 + 3*y^2 - 4, 3*x^3 - 4*y^2 + x*y - 11 ];   \\ list of equations
> 
> {
> vars = variables(eqns);                                 \\ list of variables
> n = #vars;
> sol_modpi = [vector(n,j,0)]; \\ this vector contains the solutions mod p^i, initialiszed for i=0
> pi = 1;                      \\ = p^i
> pi1= p;                      \\ = p^(i+1)
> 
> for( i = 0, k-1,
>   sol_modpi1 = [];             \\ solution mod p^(i+1)
>   for( j = 1, #sol_modpi,
>     sol = sol_modpi[j];      \\ select a solution mod p^i
>     forvec( X = vector(n,j,[0,p-1]),
>       try = sol + pi*X;
>       if(substvec(eqns,vars,try)%pi1==0,sol_modpi1 = concat(sol_modpi1,[try]));
>     );
>   );
>   print("nb of sol mod ",p,"^",i+1," = ",#sol_modpi1);
>   sol_modpi = sol_modpi1;
>   pi = pi1;
>   pi1*=p;
> );
> }
> 
> 
> ----- Mail original -----
> > De: "Bill Allombert" <Bill.Allombert@math.u-bordeaux.fr>
> > À: pari-users@pari.math.u-bordeaux.fr
> > Envoyé: Jeudi 2 Mars 2017 11:26:35
> > Objet: Re: Mathematica "Reduce" function
> > 
> > On Thu, Mar 02, 2017 at 10:00:48AM +0100, Pedro Fortuny Ayuso wrote:
> > > Thanks to all.
> > > 
> > > My specific problem is trying to solve equations like
> > > 
> > > 6x^2 + 12y^2 +20z^2 = 0
> > > 
> > > over Z/(2^k)Z. That is, finding the points of that surface
> > > over that ring.
> > 
> > Solutions of homogenous degree-2 equation in three variables can be
> > parametrized as soon as one solution is known using qfparam:
> > For example [0,1,1] is a (primitive) solution mod 2^5 so set
> > 
> > ? M=qfparam(matdiagonal([6,12,20]),[0,1,1])
> > %9 = [0,-20,0;3,0,10;3,0,-10]
> > ? v = y^2 * M*[1,x/y,(x/y)^2]~
> > %10 = [-20*y*x,10*x^2+3*y^2,-10*x^2+3*y^2]~
> > 
> > so for all x,y, (-20*y*x,10*x^2+3*y^2,-10*x^2+3*y^2) is a solution mod
> > 2^5:
> > 
> > ? content(6*(-20*y*x)^2+12*(10*x^2+3*y^2)^2+20*(-10*x^2+3*y^2)^2)
> > %12 = 32
> > 
> > > Bill's reply of counting
> > > 
> > > length([[x,y,z]|x<-[0..2^k-1];y<-[0..2^k-1];z<-[0..2^k-1],6*x^2+12*y^2+20*z^2==0])
> > 
> > You are missing a reduction mod 2^k at the end.
> > 
> > > is the fastest but it ***looks like*** a lot slower than
> > > Mathematica (but please notice I am working on a system
> > > with pari/gp and my colleague on a different one with Mathematica,
> > > so that it may have nothing to do with pari/Mathematica).
> > 
> > This is quite possible, I do not know what Mathematica is doing.
> > what you can do is to check whether (6*x^2+12*y^2)/-20 is a square
> > instead of iterating over z:
> > 
> > [[x,y,z]|x<-2*[0..2^(k-1)-1];y<-[0..2^k-1],issquare((6*x^2+12*y^2)/-20*Mod(1,2^k),&z)]
> > 
> > Cheers,
> > Bill.
> > 
> > 
> 

-- 
Pedro Fortuny Ayuso
http://pfortuny.net


EPIG, Campus de Viesques, Gijon
Dpto. de Matematicas
Universidad de Oviedo