source: git/Singular/LIB/ffsolve.lib @ a57b65

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