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

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