# source:git/Singular/LIB/ffsolve.lib@d43600

spielwiese
Last change on this file since d43600 was d43600, checked in by Hans Schoenemann <hannes@…>, 9 years ago
format
• Property mode set to `100644`
File size: 26.9 KB
Line
1/////////////////////////////////////////////////////////////////////
2version="version ffsolve.lib 4.0.0.0 Jun_2013 "; // \$Id\$
3category="Symbolic-numerical solving";
4info="
5LIBRARY: ffsolve.lib        multivariate equation solving over finite fields
6AUTHOR: Gergo Gyula Borus, borisz@borisz.net
7KEYWORDS: multivariate equations; finite field
8
9PROCEDURES:
10ffsolve();         finite field solving using heuristically chosen method
11PEsolve();         solve system of multivariate equations over finite field
12simplesolver();    solver using modified exhausting search
13GBsolve();         multivariate solver using Groebner-basis
14XLsolve();         multivariate polynomial solver using linearization
15ZZsolve();         solve system of multivariate equations over finite field
16";
17
18LIB "presolve.lib";
19LIB "general.lib";
20LIB "ring.lib";
21LIB "standard.lib";
22LIB "matrix.lib";
23
24////////////////////////////////////////////////////////////////////
25proc ffsolve(ideal equations, list #)
26"USAGE:         ffsolve(I[, L]); I ideal, L list of strings
27RETURN:         list L, the common roots of I as ideal
28ASSUME:         basering is a finite field of type (p^n,a)
29"
30{
31  list solutions, lSolvers, tempsols;
32  int i,j, k,n, R, found;
33  ideal factors, linfacs;
34  poly lp;
35  // check assumptions
36  if(npars(basering)>1)
37  {
38    ERROR("Basering must have at most one parameter");
39  }
40  if(char(basering)==0)
41  {
42    ERROR("Basering must have finite characteristic");
43  }
44  if(hasGFCoefficient(basering))
45  {
46    ERROR("not implemented for Galois fields");
47  }
48
49  if(size(#))
50  {
51    if(size(#)==1 and typeof(#[1])=="list")
52    { lSolvers = #[1]; }
53    else
54    { lSolvers = #; }
55  }
56  else
57  {
58    if(deg(equations) == 2)
59    {
60      lSolvers = "XLsolve", "PEsolve", "simplesolver", "GBsolve", "ZZsolve";
61    }
62    else
63    {
64      lSolvers = "PEsolve", "simplesolver", "GBsolve", "ZZsolve", "XLsolve";
65    }
66    if(deg(equations) == 1)
67    {
68      lSolvers = "GBsolve";
69    }
70  }
71  n = size(lSolvers);
72  R = random(1, n*(3*n+1) div 2);
73  string solver;
74  for(i=1;i<n+1;i++)
75  {
76    if(R<=(2*n+1-i))
77    {
78      solver = lSolvers[i];
79    }
80    else
81    {
82      R=R-(2*n+1-i);
83    }
84  }
85
86  if(nvars(basering)==1)
87  {
88    return(facstd(equations));
89  }
90
91  // search for the first univariate polynomial
92  found = 0;
93  for(i=1; i<=ncols(equations); i++)
94  {
95    if(univariate(equations[i])>0)
96    {
97      factors=factorize(equations[i],1);
98      for(j=1; j<=ncols(factors); j++)
99      {
100        if(deg(factors[j])==1)
101        {
102          linfacs[size(linfacs)+1] = factors[j];
103        }
104      }
105      if(deg(linfacs[1])>0)
106      {
107        found=1;
108        break;
109      }
110    }
111  }
112  //   if there is, collect its the linear factors
113  if(found)
114  {
115    // substitute the root and call recursively
116    ideal neweqs, invmapideal, ti;
117    map invmap;
118    for(k=1; k<=ncols(linfacs); k++)
119    {
120      lp = linfacs[k];
121      neweqs = reduce(equations, lp);
122
124      def original_ring = basering;
125      def newRing = clonering(nvars(original_ring)-1);
126      setring newRing;
127      ideal mappingIdeal;
128      j=1;
129      for(i=1; i<=size(varexp); i++)
130      {
131        if(varexp[i])
132        {
133          mappingIdeal[i] = 0;
134        }
135        else
136        {
137          mappingIdeal[i] = var(j);
138          j++;
139        }
140      }
141      map recmap = original_ring, mappingIdeal;
142      list tsols = ffsolve(recmap(neweqs), lSolvers);
143      if(size(tsols)==0)
144      {
145        tsols = list(ideal(1));
146      }
147      setring original_ring;
148      j=1;
149      for(i=1;i<=size(varexp);i++)
150      {
151        if(varexp[i]==0)
152        {
153          invmapideal[j] = var(i);
154          j++;
155        }
156      }
157      invmap = newRing, invmapideal;
158      tempsols = invmap(tsols);
159
160      // combine the solutions
161      for(j=1; j<=size(tempsols); j++)
162      {
163        ti =  std(tempsols[j]+lp);
164        if(deg(ti)>0)
165        {
166          solutions = insert(solutions,ti);
167        }
168      }
169    }
170  }
171  else
172  {
173    execute("solutions="+solver+"(equations);") ;
174  }
175  return(solutions);
176}
177example
178{
179  "EXAMPLE:";echo=2;
180  ring R = (2,a),x(1..3),lp;
181  minpoly=a2+a+1;
182  ideal I;
183  I[1]=x(1)^2*x(2)+(a)*x(1)*x(2)^2+(a+1);
184  I[2]=x(1)^2*x(2)*x(3)^2+(a)*x(1);
185  I[3]=(a+1)*x(1)*x(3)+(a+1)*x(1);
186  ffsolve(I);
187}
188////////////////////////////////////////////////////////////////////
189proc PEsolve(ideal L, list #)
190"USAGE:         PEsolve(I[, i]); I ideal, i optional integer
191                solve I (system of multivariate equations) over a
192                finite field using an equvalence property when i is
193                not given or set to 2, otherwise if i is set to 0
194                then check whether common roots exists
195RETURN:         list if optional parameter is not given or set to 2,
196                integer if optional is set to 0
197ASSUME:         basering is a finite field of type (p^n,a)
198NOTE:           When the optional parameter is set to 0, speoff only
199                checks if I has common roots, then return 1, otherwise
200                return 0.
201"
202{
203  int mode, i,j;
204  list results, rs, start;
205  poly g;
206  // check assumptions
207  if(npars(basering)>1)
208  {
209    ERROR("Basering must have at most one parameter");
210  }
211  if(char(basering)==0)
212  {
213    ERROR("Basering must have finite characteristic");
214  }
215  if(hasGFCoefficient(basering))
216  {
217    ERROR("not implemented for Galois fields");
218  }
219
220  if( size(#) > 0 )
221  {
222    mode = #[1];
223  }
224  else
225  {
226    mode = 2;
227  }
228  L = simplify(L,15);
229  g = productOfEqs( L );
230
231  if(g == 0)
232  {
233    if(mode==0)
234    {
235      return(0);
236    }
237    return( list() );
238  }
239  if(g == 1)
240  {
241    list vectors = every_vector();
242    for(j=1; j<=size(vectors); j++)
243    {
244      ideal res;
245      for(i=1; i<=nvars(basering); i++)
246      {
247        res[i] = var(i)-vectors[j][i];
248      }
249      results[size(results)+1] = std(res);
250    }
251    return( results );
252  }
253
254  if( mode == 0 )
255  {
256    return( 1 );
257  }
258  else
259  {
260    for(i=1; i<=nvars(basering); i++)
261    {
262      start[i] = 0:order_of_extension();
263    }
264
265    if( mode == 1)
266    {
267      results[size(results)+1] = melyseg(g, start);
268    }
269    else
270    {
271      while(1)
272      {
273        start = melyseg(g, start);
274        if( size(start) > 0 )
275        {
276          ideal res;
277          for(i=1; i<=nvars(basering); i++)
278          {
279            res[i] = var(i)-vec2elm(start[i]);
280          }
281          results[size(results)+1] = std(res);
282          start = increment(start);
283        }else{
284          break;
285        }
286      }
287    }
288  }
289  return(results);
290}
291example
292{
293  "EXAMPLE:";echo=2;
294  ring R = (2,a),x(1..3),lp;
295  minpoly=a2+a+1;
296  ideal I;
297  I[1]=x(1)^2*x(2)+(a)*x(1)*x(2)^2+(a+1);
298  I[2]=x(1)^2*x(2)*x(3)^2+(a)*x(1);
299  I[3]=(a+1)*x(1)*x(3)+(a+1)*x(1);
300  PEsolve(I);
301}
302////////////////////////////////////////////////////////////////////
303proc simplesolver(ideal E)
304"USAGE:         simplesolver(I); I ideal
305                solve I (system of multivariate equations) over a
306                finite field by exhausting search
307RETURN:         list L, the common roots of I as ideal
308ASSUME:         basering is a finite field of type (p^n,a)
309"
310{
311  int i,j,k,t, correct;
312  list solutions = list(std(ideal()));
313  list partial_solutions;
314  ideal partial_system, curr_sol, curr_sys, factors;
315  poly univar_poly;
316  E = E+defaultIdeal();
317  // check assumptions
318  if(npars(basering)>1)
319  {
320    ERROR("Basering must have at most one parameter");
321  }
322  if(char(basering)==0)
323  {
324    ERROR("Basering must have finite characteristic");
325  }
326  if(hasGFCoefficient(basering))
327  {
328    ERROR("not implemented for Galois fields");
329  }
330  for(k=1; k<=nvars(basering); k++)
331  {
332    partial_solutions = list();
333    for(i=1; i<=size(solutions); i++)
334    {
335      partial_system = reduce(E, solutions[i]);
336      for(j=1; j<=ncols(partial_system); j++)
337      {
338        if(univariate(partial_system[j])>0)
339        {
340          univar_poly = partial_system[j];
341          break;
342        }
343      }
344      factors = factorize(univar_poly,1);
345      for(j=1; j<=ncols(factors); j++)
346      {
347        if(deg(factors[j])==1)
348        {
349          curr_sol = std(solutions[i]+ideal(factors[j]));
350          curr_sys = reduce(E, curr_sol);
351          correct = 1;
352          for(t=1; t<=ncols(curr_sys); t++)
353          {
354            if(deg(curr_sys[t])==0)
355            {
356              correct = 0;
357              break;
358            }
359          }
360          if(correct)
361          {
362            partial_solutions = insert(partial_solutions, curr_sol);
363          }
364        }
365      }
366    }
367    solutions = partial_solutions;
368  }
369  return(solutions);
370}
371example
372{
373  "EXAMPLE:";echo=2;
374  ring R = (2,a),x(1..3),lp;
375  minpoly=a2+a+1;
376  ideal I;
377  I[1]=x(1)^2*x(2)+(a)*x(1)*x(2)^2+(a+1);
378  I[2]=x(1)^2*x(2)*x(3)^2+(a)*x(1);
379  I[3]=(a+1)*x(1)*x(3)+(a+1)*x(1);
380  simplesolver(I);
381}
382////////////////////////////////////////////////////////////////////
383proc GBsolve(ideal equation_system)
384"USAGE:         GBsolve(I); I ideal
385                solve I (system of multivariate equations) over an
386                extension of Z/p by Groebner basis methods
387RETURN:         list L, the common roots of I as ideal
388ASSUME:         basering is a finite field of type (p^n,a)
389"
390{
391  int i,j, prop, newelement, number_new_vars;
392  ideal ls;
393  list results, slvbl, linsol, ctrl, new_sols, varinfo;
394  ideal I, linear_solution, unsolved_part, univar_part, multivar_part, unsolved_vars;
395  intvec unsolved_var_nums;
396  string new_vars;
397  // check assumptions
398  if(npars(basering)>1)
399  {
400    ERROR("Basering must have at most one parameter");
401  }
402  if(char(basering)==0)
403  {
404    ERROR("Basering must have finite characteristic");
405  }
406  if(hasGFCoefficient(basering))
407  {
408    ERROR("not implemented for Galois fields");
409  }
410
411  def original_ring = basering;
412  if(npars(basering)==1)
413  {
414    int prime_coeff_field=0;
415    string minpolystr = "minpoly="+
416    get_minpoly_str(size(original_ring),parstr(original_ring,1))+";" ;
417  }
418  else
419  {
420    int prime_coeff_field=1;
421  }
422
423  option(redSB);
424
425  equation_system = simplify(equation_system,15);
426
427  ideal standard_basis = std(equation_system);
428  list basis_factors = facstd(standard_basis);
429  if( basis_factors[1][1] == 1)
430  {
431    return(results)
432  };
433
434  for(i=1; i<= size(basis_factors); i++)
435  {
436    prop = 0;
437    for(j=1; j<=size(basis_factors[i]); j++)
438    {
439      if( univariate(basis_factors[i][j])>0 and deg(basis_factors[i][j])>1)
440      {
441        prop =1;
442        break;
443      }
444    }
445    if(prop == 0)
446    {
447      ls = solvelinearpart( basis_factors[i] );
448      if(ncols(ls) == nvars(basering) )
449      {
450        ctrl, newelement = add_if_new(ctrl, ls);
451        if(newelement)
452        {
453          results = insert(results, ls);
454        }
455      }
456      else
457      {
458        slvbl = insert(slvbl, list(basis_factors[i],ls) );
459      }
460    }
461  }
462  if(size(slvbl)<>0)
463  {
464    for(int E = 1; E<= size(slvbl); E++)
465    {
466      I = slvbl[E][1];
467      linear_solution = slvbl[E][2];
468      attrib(I,"isSB",1);
469      unsolved_part = reduce(I,linear_solution);
470      univar_part = ideal();
471      multivar_part = ideal();
472      for(i=1; i<=ncols(I); i++)
473      {
474        if(univariate(I[i])>0)
475        {
476          univar_part = univar_part+I[i];
477        }
478        else
479        {
480          multivar_part = multivar_part+I[i];
481        }
482      }
483      varinfo = varaibles(univar_part);
484      unsolved_vars = varinfo[3];
485      unsolved_var_nums = varinfo[4];
486      number_new_vars = ncols(unsolved_vars);
487
488      new_vars = "@y(1.."+string(number_new_vars)+")";
489      def R_new = changevar(new_vars, original_ring);
490      setring R_new;
491      if( !prime_coeff_field )
492      {
493        execute(minpolystr);
494      }
495
496      ideal mapping_ideal;
497      for(i=1; i<=size(unsolved_var_nums); i++)
498      {
499        mapping_ideal[unsolved_var_nums[i]] = var(i);
500      }
501
502      map F = original_ring, mapping_ideal;
503      ideal I_new = F( multivar_part );
504
505      list sol_new;
506      int unsolvable = 0;
507      sol_new = simplesolver(I_new);
508      if( size(sol_new) == 0)
509      {
510        unsolvable = 1;
511      }
512      setring original_ring;
513      if(unsolvable)
514      {
515        list sol_old = list();
516      }
517      else
518      {
519        map G = R_new, unsolved_vars;
520        new_sols = G(sol_new);
521        for(i=1; i<=size(new_sols); i++)
522        {
523          ideal sol = new_sols[i]+linear_solution;
524          sol = std(sol);
525          ctrl, newelement = add_if_new(ctrl, sol);
526          if(newelement)
527          {
528            results = insert(results, sol);
529          }
530          kill sol;
531        }
532      }
533      kill G;
534      kill R_new;
535    }
536  }
537  return( results  );
538}
539example
540{
541  "EXAMPLE:";echo=2;
542  ring R = (2,a),x(1..3),lp;
543  minpoly=a2+a+1;
544  ideal I;
545  I[1]=x(1)^2*x(2)+(a)*x(1)*x(2)^2+(a+1);
546  I[2]=x(1)^2*x(2)*x(3)^2+(a)*x(1);
547  I[3]=(a+1)*x(1)*x(3)+(a+1)*x(1);
548  GBsolve(I);
549}
550////////////////////////////////////////////////////////////////////
551proc XLsolve(ideal I, list #)
552"USAGE:         XLsolve(I[, d]); I ideal, d optional integer
553                solve I (system of multivariate polynomials) with a
554                variant of the linearization technique, multiplying
555                the polynomials with monomials of degree at most d
556                (default is 2)
557RETURN:         list L of the common roots of I as ideals
558ASSUME:         basering is a finite field of type (p^n,a)"
559{
560  int i,j,k, D;
561  int SD = deg(I);
562  list solutions;
563  if(size(#))
564  {
565    if(typeof(#[1])=="int") { D = #[1]; }
566  }
567  else
568  { D = 2; }
569  list lMonomialsForMultiplying = monomialsOfDegreeAtMost(D+SD);
570
571  int m = ncols(I);
572  list extended_system;
573  list mm;
574  for(k=1; k<=size(lMonomialsForMultiplying)-SD; k++)
575  {
576    mm = lMonomialsForMultiplying[k];
577    for(i=1; i<=m; i++)
578    {
579      for(j=1; j<=size(mm); j++)
580      {
581        extended_system[size(extended_system)+1] = reduce(I[i]*mm[j], defaultIdeal());
582      }
583    }
584  }
585  ideal new_system = I;
586  for(i=1; i<=size(extended_system); i++)
587  {
588    new_system[m+i] = extended_system[i];
589  }
590  ideal reduced_system = linearReduce( new_system, lMonomialsForMultiplying);
591
592  solutions = simplesolver(reduced_system);
593
594  return(solutions);
595}
596example
597{
598  "EXAMPLE:";echo=2;
599  ring R = (2,a),x(1..3),lp;
600  minpoly=a2+a+1;
601  ideal I;
602  I[1]=(a)*x(1)^2+x(2)^2+(a+1);
603  I[2]=(a)*x(1)^2+(a)*x(1)*x(3)+(a)*x(2)^2+1;
604  I[3]=(a)*x(1)*x(3)+1;
605  I[4]=x(1)^2+x(1)*x(3)+(a);
606  XLsolve(I, 3);
607}
608
609////////////////////////////////////////////////////////////////////
610proc ZZsolve(ideal I)
611"USAGE:         ZZsolve(I); I ideal
612solve I (system of multivariate equations) over a
613finite field by mapping the polynomials to a single
614univariate polynomial over extension of the basering
615RETURN:         list, the common roots of I as ideal
616ASSUME:         basering is a finite field of type (p^n,a)
617"
618{
619  int i, j, nv, numeqs,r,l,e;
620  def original_ring = basering;
621  // check assumptions
622  if(npars(basering)>1)
623  {
624    ERROR("Basering must have at most one parameter");
625  }
626  if(char(basering)==0)
627  {
628    ERROR("Basering must have finite characteristic");
629  }
630  if(hasGFCoefficient(basering))
631  {
632    ERROR("not implemented for Galois fields");
633  }
634
635  nv = nvars(original_ring);
636  numeqs = ncols(I);
637  l = numeqs % nv;
638  if( l == 0)
639  {
640    r = numeqs div nv;
641  }
642  else
643  {
644    r = (numeqs div nv) +1;
645  }
646
647
648  list list_of_equations;
649  for(i=1; i<=r; i++)
650  {
651    list_of_equations[i] = ideal();
652  }
653  for(i=0; i<numeqs; i++)
654  {
655    list_of_equations[(i div nv)+1][(i % nv) +1] = I[i+1];
656  }
657
658  ring ring_for_matrix = (char(original_ring),@y),(x(1..nv),@X,@c(1..nv)(1..nv)),lp;
659  execute("minpoly="+Z_get_minpoly(size(original_ring)^nv, parstr(1))+";");
660
661  ideal IV;
662  for(i=1; i<=nv; i++)
663  {
664    IV[i] = var(i);
665  }
666
667  matrix M_C[nv][nv];
668  for(i=1;i<=nrows(M_C); i++)
669  {
670    for(j=1; j<=ncols(M_C); j++)
671    {
672      M_C[i,j] = @c(i)(j);
673    }
674  }
675
676  poly X = Z_phi(IV);
677  ideal IX_power_poly;
678  ideal IX_power_var;
679  for(i=1; i<=nv; i++)
680  {
681    e = (size(original_ring)^(i-1));
682    IX_power_poly[i] = X^e;
683    IX_power_var[i] = @X^e;
684  }
685  IX_power_poly = reduce(IX_power_poly, Z_default_ideal(nv, size(original_ring)));
686
687  def M = matrix(IX_power_poly,1,nv)*M_C - matrix(IV,1,nv);
688
689  ideal IC;
690  for(i=1; i<=ncols(M); i++)
691  {
692    for(j=1; j<=ncols(IV); j++)
693    {
694      IC[(i-1)*ncols(M)+j] = coeffs(M[1,i],IV[j])[2,1];
695    }
696  }
697
698  ideal IC_solultion = std(Presolve::solvelinearpart(IC));
699
700  matrix M_C_sol[nv][nv];
701  for(i=1;i<=nrows(M_C_sol); i++)
702  {
703    for(j=1; j<=ncols(M_C_sol); j++)
704    {
705      M_C_sol[i,j] = reduce(@c(i)(j), std(IC_solultion));
706    }
707  }
708  ideal I_subs;
709  I_subs = ideal(matrix(IX_power_var,1,nv)*M_C_sol);
710
711  setring original_ring;
712  string var_str = varstr(original_ring)+",@X,@y";
713  string minpoly_str = "minpoly="+string(minpoly)+";";
714  def ring_for_substitution = Ring::changevar(var_str, original_ring);
715
716  setring ring_for_substitution;
717  execute(minpoly_str);
718
719  ideal I_subs = imap(ring_for_matrix, I_subs);
720  ideal I = imap(original_ring, I);
721  list list_of_equations = imap(original_ring, list_of_equations);
722
723  list list_of_F;
724  for(i=1; i<=r; i++)
725  {
726    list_of_F[i] = Z_phi( list_of_equations[i] );
727  }
728
729  for(i=1; i<=nv; i++)
730  {
731    for(j=1; j<=r; j++)
732    {
733      list_of_F[j] = subst( list_of_F[j], var(i), I_subs[i] );
734    }
735  }
736  int s = size(original_ring);
737  if(npars(original_ring)==1)
738  {
739    for(j=1; j<=r; j++)
740    {
741      list_of_F[j] = subst(list_of_F[j], par(1), (@y^( (s^nv-1) div (s-1) )));
742    }
743  }
744
745  ring temp_ring = (char(original_ring),@y),@X,lp;
746  list list_of_F = imap(ring_for_substitution, list_of_F);
747
748  ring ring_for_factorization = (char(original_ring),@y),X,lp;
749  execute("minpoly="+Z_get_minpoly(size(original_ring)^nv, parstr(1))+";");
750  map rho = temp_ring,X;
751  list list_of_F = rho(list_of_F);
752  poly G = 0;
753  for(i=1; i<=r; i++)
754  {
755    G = gcd(G, list_of_F[i]);
756  }
757  if(G==1)
758  {
759    return(list());
760  }
761
762  list factors = Presolve::linearpart(factorize(G,1));
763
764  ideal check;
765  for(i=1; i<=nv; i++)
766  {
767    check[i] = X^(size(original_ring)^(i-1));
768  }
769  list fsols;
770
771  matrix sc;
772  list sl;
773  def sM;
774  matrix M_for_sol = fetch(ring_for_matrix, M_C_sol);
775  for(i=1; i<=size(factors[1]); i++)
776  {
777    sc = matrix(reduce(check, std(factors[1][i])), 1,nv  );
778
779    sl = list();
780    sM = sc*M_for_sol;
781    for(j=1; j<=ncols(sM); j++)
782    {
783      sl[j] = sM[1,j];
784    }
785    fsols[i] = sl;
786  }
787  if(size(fsols)==0)
788  {
789    return(list());
790  }
791  setring ring_for_substitution;
792  list ssols = imap(ring_for_factorization, fsols);
793  if(npars(original_ring)==1)
794  {
795    execute("poly P="+Z_get_minpoly(size(original_ring)^nv, "@y"));
796    poly RP = gcd(P,  (@y^( (s^nv-1) div (s-1) ))-a);
797    for(i=1; i<=size(ssols); i++)
798    {
799      for(j=1; j<=size(ssols[i]); j++)
800      {
801        ssols[i][j] = reduce( ssols[i][j], std(RP));
802      }
803    }
804  }
805  setring original_ring;
806  list solutions = imap(ring_for_substitution, ssols);
807  list final_solutions;
808  ideal ps;
809  for(i=1; i<=size(solutions); i++)
810  {
811    ps = ideal();
812    for(j=1; j<=nvars(original_ring); j++)
813    {
814      ps[j] = var(j)-solutions[i][j];
815    }
816    final_solutions = insert(final_solutions, std(ps));
817  }
818  return(final_solutions);
819}
820example
821{
822  "EXAMPLE:";echo=2;
823  ring R = (2,a),x(1..3),lp;
824  minpoly=a2+a+1;
825  ideal I;
826  I[1]=x(1)^2*x(2)+(a)*x(1)*x(2)^2+(a+1);
827  I[2]=x(1)^2*x(2)*x(3)^2+(a)*x(1);
828  I[3]=(a+1)*x(1)*x(3)+(a+1)*x(1);
829  ZZsolve(I);
830}
831////////////////////////////////////////////////////////////////////
832////////////////////////////////////////////////////////////////////
833static proc linearReduce(ideal I, list mons)
834{
835  //system("--no-warn", 1);
836  int LRtime = rtimer;
837  int i,j ;
838  int prime_field = 1;
839  list solutions, monomials;
840  for(i=1; i<=size(mons); i++)
841  {
842    monomials = reorderMonomials(mons[i])+monomials;
843  }
844  int number_of_monomials = size(monomials);
845
846  def original_ring = basering;
847  if(npars(basering)==1)
848  {
849    prime_field=0;
850    string minpolystr = "minpoly="
851    +get_minpoly_str(size(original_ring),parstr(original_ring,1))+";" ;
852  }
853  string old_vars = varstr(original_ring);
854  string new_vars = "@y(1.."+string( number_of_monomials )+")";
855
856  def ring_for_var_change = changevar( old_vars+","+new_vars, original_ring);
857  setring ring_for_var_change;
858  if( prime_field == 0)
859  {
860    execute(minpolystr);
861  }
862
863  list monomials = imap(original_ring, monomials);
864  ideal I = imap(original_ring, I);
865  ideal C;
866  intvec weights=1:nvars(original_ring);
867
868  for(i=1; i<= number_of_monomials; i++)
869  {
870    C[i] = monomials[i] - @y(i);
871    weights = weights,deg(monomials[i])+1;
872  }
873  ideal linear_eqs = I;
874  for(i=1; i<=ncols(C); i++)
875  {
876    linear_eqs = reduce(linear_eqs, C[i]);
877  }
878
879  def ring_for_elimination = changevar( new_vars, ring_for_var_change);
880  setring ring_for_elimination;
881  if( prime_field == 0)
882  {
883    execute(minpolystr);
884  }
885
886  ideal I = imap(ring_for_var_change, linear_eqs);
887  ideal lin_sol = solvelinearpart(I);
888  def ring_for_back_change = changeord( list(list("wp",weights),list("C",0:1)), ring_for_var_change);
889
890  setring ring_for_back_change;
891  if( prime_field == 0)
892  {
893    execute(minpolystr);
894  }
895
896  ideal lin_sol = imap(ring_for_elimination, lin_sol);
897  ideal C = imap(ring_for_var_change, C);
898  ideal J = lin_sol;
899  for(i=1; i<=ncols(C); i++)
900  {
901    J = reduce(J, C[i]);
902  }
903  setring original_ring;
904  ideal J = imap(ring_for_back_change, J);
905  return(J);
906}
907
908static proc monomialsOfDegreeAtMost(int k)
909{
910  int Mtime = rtimer;
911  list monomials, monoms, monoms_one, lower_monoms;
912  int n = nvars(basering);
913  int t,i,l,j,s;
914  for(i=1; i<=n; i++)
915  {
916    monoms_one[i] = var(i);
917  }
918  monomials = list(monoms_one);
919  if(1 < k)
920  {
921    for(t=2; t<=k; t++)
922    {
923      lower_monoms = monomials[t-1];
924      monoms = list();
925      s = 1;
926      for(i=1; i<=n; i++)
927      {
928        for(j=s; j<=size(lower_monoms); j++)
929        {
930          monoms = monoms+list(lower_monoms[j]*var(i));
931        }
932        s = s + int(binomial(n+t-2-i, t-2));
933      }
934      monomials[t] = monoms;
935    }
936  }
937  return(monomials);
938}
939
940static proc reorderMonomials(list monomials)
941{
942  list univar_monoms, mixed_monoms;
943
944  for(int j=1; j<=size(monomials); j++)
945  {
946    if( univariate(monomials[j])>0 )
947    {
948      univar_monoms = univar_monoms + list(monomials[j]);
949    }
950    else
951    {
952      mixed_monoms = mixed_monoms + list(monomials[j]);
953    }
954  }
955
956  return(univar_monoms + mixed_monoms);
957}
958
959static proc melyseg(poly g, list start)
960{
961  list gsub = g;
962  int i = 1;
963
964  while( start[1][1] <> char(basering) )
965  {
966    gsub[i+1] = subst( gsub[i], var(i), vec2elm(start[i]));
967    if( gsub[i+1] == 0 )
968    {
969      list new = increment(start,i);
970      for(int l=1; l<=size(start); l++)
971      {
972        if(start[l]<>new[l])
973        {
974          i = l;
975          break;
976        }
977      }
978      start = new;
979    }
980    else
981    {
982      if(i == nvars(basering))
983      {
984        return(start);
985      }else{
986        i++;
987      }
988    }
989  }
990  return(list());
991}
992
993static proc productOfEqs(ideal I)
994{
995  //system("--no-warn", 1);
996  ideal eqs = sort_ideal(I);
997  int i,q;
998  poly g = 1;
999  q = size(basering);
1000  ideal I = defaultIdeal();
1001
1002  for(i=1; i<=size(eqs); i++)
1003  {
1004    if(g==0){return(g);}
1005    g = reduce(g*(eqs[i]^(q-1)-1), I);
1006  }
1007  return( g );
1008}
1009
1010static proc clonering(list #)
1011{
1012  def original_ring = basering;
1013  int n = nvars(original_ring);
1014  int prime_field=npars(basering);
1015  if(prime_field)
1016  {
1017    string minpolystr = "minpoly="+
1018    get_minpoly_str(size(original_ring),parstr(original_ring,1))+";" ;
1019  }
1020
1021  if(size(#))
1022  {
1023    int newvars = #[1];
1024  }
1025  else
1026  {
1027    int newvars = nvars(original_ring);
1028  }
1029  string newvarstr = "v(1.."+string(newvars)+")";
1030  def newring = changevar(newvarstr, original_ring);
1031  setring newring;
1032  if( prime_field )
1033  {
1034    execute(minpolystr);
1035  }
1036  return(newring);
1037}
1038
1039static proc defaultIdeal()
1040{
1041  ideal I;
1042  for(int i=1; i<=nvars(basering); i++)
1043  {
1044    I[i] = var(i)^size(basering)-var(i);
1045  }
1046  return( std(I) );
1047}
1048
1049static proc order_of_extension()
1050{
1051  int oe=1;
1052  list rl = ringlist(basering);
1053  if( size(rl[1]) <> 1)
1054  {
1055    oe = deg( subst(minpoly,par(1),var(1)) );
1056  }
1057  return(oe);
1058}
1059
1060static proc vec2elm(intvec v)
1061{
1062  number g = 1;
1063  if(npars(basering) == 1) { g=par(1); }
1064  number e=0;
1065  int oe = size(v);
1066  for(int i=1; i<=oe; i++)
1067  {
1068    e = e+v[i]*g^(oe-i);
1069  }
1070  return(e);
1071}
1072
1073static proc increment(list l, list #)
1074{
1075  int c, i, j, oe;
1076  oe = order_of_extension();
1077  c = char(basering);
1078
1079  if( size(#) == 1 )
1080  {
1081    i = #[1];
1082  }
1083  else
1084  {
1085    i = size(l);
1086  }
1087
1088  l[i] = nextVec(l[i]);
1089  while( l[i][1] == c && i>1 )
1090  {
1091    l[i] = 0:oe;
1092    i--;
1093    l[i] = nextVec(l[i]);
1094  }
1095  if( i < size(l) )
1096  {
1097    for(j=i+1; j<=size(l); j++)
1098    {
1099      l[j] = 0:oe;
1100    }
1101  }
1102  return(l);
1103}
1104
1105static proc nextVec(intvec l)
1106{
1107  int c, i, j;
1108  i = size(l);
1109  c = char(basering);
1110  l[i] = l[i] + 1;
1111  while( l[i] == c && i>1 )
1112  {
1113    l[i] = 0;
1114    i--;
1115    l[i] = l[i] + 1;
1116  }
1117  return(l);
1118}
1119
1120static proc every_vector()
1121{
1122  list element, list_of_elements;
1123
1124  for(int i=1; i<=nvars(basering); i++)
1125  {
1126    element[i] = 0:order_of_extension();
1127  }
1128
1129  while(size(list_of_elements) < size(basering)^nvars(basering))
1130  {
1131    list_of_elements = list_of_elements + list(element);
1132    element = increment(element);
1133  }
1134  for(int i=1; i<=size(list_of_elements); i++)
1135  {
1136    for(int j=1; j<=size(list_of_elements[i]); j++)
1137    {
1138      list_of_elements[i][j] = vec2elm(list_of_elements[i][j]);
1139    }
1140  }
1141  return(list_of_elements);
1142}
1143
1144static proc num2int(number a)
1145{
1146  int N=0;
1147  if(order_of_extension() == 1)
1148  {
1149    N = int(a);
1150    if(N<0)
1151    {
1152      N = N + char(basering);
1153    }
1154  }
1155  else
1156  {
1157    ideal C = coeffs(subst(a,par(1),var(1)),var(1));
1158    for(int i=1; i<=ncols(C); i++)
1159    {
1160      int c = int(C[i]);
1161      if(c<0) { c = c + char(basering); }
1162      N = N + c*char(basering)^(i-1);
1163    }
1164  }
1165  return(N);
1166}
1167
1168static proc get_minpoly_str(int size_of_ring, string parname)
1169{
1170  def original_ring = basering;
1171  ring new_ring = (size_of_ring, A),x,lp;
1172  string S = string(minpoly);
1173  string SMP;
1174  if(S=="0")
1175  {
1176    SMP = SMP+parname;
1177  }
1178  else
1179  {
1180    for(int i=1; i<=size(S); i++)
1181    {
1182      if(S[i]=="A")
1183      {
1184        SMP = SMP+parname;
1185      }
1186      else
1187      {
1188        SMP = SMP+S[i];
1189      }
1190    }
1191  }
1192  return(SMP);
1193}
1194
1195static proc sort_ideal(ideal I)
1196{
1197  ideal OI;
1198  int i,j,M;
1199  poly P;
1200  M = ncols(I);
1201  OI = I;
1202  for(i=2; i<=M; i++)
1203  {
1204    j=i;
1205    while(size(OI[j-1])>size(OI[j]))
1206    {
1207      P = OI[j-1];
1208      OI[j-1] = OI[j];
1209      OI[j] = P;
1210      j--;
1211      if(j==1) break;
1212    }
1213  }
1214  return(OI);
1215}
1216
1217static proc add_if_new(list L, ideal I)
1218{
1219  int i, newelement;
1220  poly P;
1221
1222  I=std(I);
1223  for(i=1; i<=nvars(basering); i++)
1224  {
1225    P = P +  reduce(var(i),I)*var(1)^(i-1);
1226  }
1227  newelement=1;
1228  for(i=1; i<=size(L); i++)
1229  {
1230    if(L[i]==P)
1231    {
1232      newelement=0;
1233      break;
1234    }
1235  }
1236  if(newelement)
1237  {
1238    L = insert(L, P);
1239  }
1240  return(L,newelement);
1241}
1242
1243static proc Z_get_minpoly(int size_of_ring, string parname)
1244{
1245  def original_ring = basering;
1246  ring new_ring = (size_of_ring, A),x,lp;
1247  string S = string(minpoly);
1248  string SMP;
1249  if(S=="0")
1250  {
1251    SMP = SMP+parname;
1252  }
1253  else
1254  {
1255    for(int i=1; i<=size(S); i++)
1256    {
1257      if(S[i]=="A")
1258      {
1259        SMP = SMP+parname;
1260      }
1261      else
1262      {
1263        SMP = SMP+S[i];
1264      }
1265    }
1266  }
1267  return(SMP);
1268}
1269
1270static proc Z_phi(ideal I)
1271{
1272  poly f;
1273  for(int i=1; i<= ncols(I); i++)
1274{
1275    f = f+I[i]*@y^(i-1);
1276  }
1277  return(f);
1278}
1279
1280static proc Z_default_ideal(int number_of_variables, int q)
1281{
1282  ideal DI;
1283  for(int i=1; i<=number_of_variables; i++)
1284  {
1285    DI[i] = var(i)^q-var(i);
1286  }
1287  return(std(DI));
1288}
Note: See TracBrowser for help on using the repository browser.