interpolation.cc
Go to the documentation of this file.
1 /****************************************
2 * Computer Algebra System SINGULAR *
3 ****************************************/
4 
5 #include "kernel/mod2.h"
6 
7 #include "factory/factory.h"
8 
9 #include "misc/options.h"
10 #include "misc/intvec.h"
11 
12 #include "coeffs/longrat.h" // snumber ...
13 
14 #include "polys/monomials/ring.h"
15 
16 #include "kernel/polys.h"
17 #include "kernel/ideals.h"
18 
19 #include "interpolation.h"
20 
21 // parameters to debug
22 //#define shmat
23 //#define checksize
24 //#define writemsg
25 
26 // possible strategies
27 #define unsortedmatrix
28 //#define integerstrategy
29 
30 #define modp_number int
31 #define exponent int
32 
33 static modp_number myp; // all modp computation done mod myp
34 static int myp_index; // index of small prime in Singular table of primes
35 
37 {
38 // return (x*y)%myp;
39  return (modp_number)(((unsigned long)x)*((unsigned long)y)%((unsigned long)myp));
40 }
42 {
43 // if (x>=y) return x-y; else return x+myp-y;
44  return (x+myp-y)%myp;
45 }
46 
48 typedef struct {mono_type mon; unsigned int point_ref;} condition_type;
51 
54 typedef struct mon_list_entry_struct mon_list_entry;
55 
58  int first_col;
60 
61 typedef struct row_list_entry_struct row_list_entry;
62 
67 
68 typedef struct generator_struct generator_entry;
69 
71  generator_entry *generator;
75 
76 typedef struct modp_result_struct modp_result_entry;
77 
79 typedef mpq_t *q_coordinates;
80 typedef mpz_t *int_coordinates;
81 typedef bool *coord_exist_table;
82 
83 static int final_base_dim; // dimension of the quotient space, known from the beginning
84 static int last_solve_column; // last non-zero column in "solve" part of matrix, used for speed up
85 static int n_points; // modp_number of ideals (points)
86 static int *multiplicity; // multiplicities of points
87 static int variables; // modp_number of variables
88 static int max_coord; // maximal possible coordinate product used during Evaluation
89 static bool only_modp; // perform only one modp computations
90 
91 static modp_coordinates *modp_points; // coordinates of points for modp problem - used by Evaluate (this is also initial data for only modp)
92 static q_coordinates *q_points; // coordinates of points for rational data (not used for modp)
93 static int_coordinates *int_points; // coordinates of points for integer data - used to check generators (not used for modp)
94 static coord_exist_table *coord_exist; // checks whether all coordinates has been initialized
95 static mon_list_entry *check_list; // monomials to be checked in next stages
96 static coordinates *points; // power products of coordinates of points used in modp cycles
97 static condition_type *condition_list; // conditions stored in an array
98 static mon_list_entry *lt_list; // leading terms found so far
99 static mon_list_entry *base_list; // standard monomials found so far
100 static row_list_entry *row_list; // rows of the matrix (both parts)
101 static modp_number *my_row; // one special row to perform operations
102 static modp_number *my_solve_row; // one special row to find the linear dependence ("solve" part)
103 static mono_type *column_name; // monomials assigned to columns in solve_row
104 
105 static int n_results; // modp_number of performed modp computations (not discarded)
106 static modp_number modp_denom; // denominator of mod p computations
107 static modp_result_entry *modp_result; // list of results for various mod p calculations (used for modp - first result is the desired one)
108 static modp_result_entry *cur_result; // pointer to current result (as before)
109 static modp_number *congr; // primes used in computations (chinese remainder theorem) (not used for modp)
110 static modp_number *in_gamma; // inverts used in chinese remainder theorem (not used for modp)
111 static mpz_t bigcongr; // result, in fact, is given in mod bigcongr (not used for modp)
112 
113 static mpz_t *polycoef; // polynomial integercoefficients (not used for modp)
114 static mono_type *polyexp; // polynomial exponents
115 
116 struct gen_list_struct {mpz_t *polycoef;
119 typedef struct gen_list_struct gen_list_entry;
120 
121 static gen_list_entry *gen_list=NULL; // list of resulting generators - output data (integer version)
122 
123 static int generic_n_generators; // modp_number of generators - should be the same for all modp comp (not used for modp)
124 static mono_type *generic_column_name; // monomials assigned to columns in solve_row - should be the same for all modp comp (!!! used for modp)
125 static mon_list_entry *generic_lt=NULL; // leading terms for ordered generators - should be the same for all modp comp (not used for modp)
126 static int good_primes; // modp_number of good primes so far;
127 static int bad_primes; // modp_number of bad primes so far;
128 static mpz_t common_denom; // common denominator used to force points coordinates to Z (not used for modp)
129 static bool denom_divisible; // common denominator is divisible by p (not used for modp)
130 
131 static poly comparizon_p1; //polynomials used to do comparizons by Singular
132 static poly comparizon_p2;
133 
134 static modp_number *modp_Reverse; // reverses in mod p
135 
136 static bool protocol; // true to show the protocol
137 
138 #ifdef checksize
139 static int maximal_size=0;
140 #endif
141 
142 #if 0 /* only for debuggig*/
143 void WriteMono (mono_type m) // Writes a monomial on the screen - only for debug
144 {
145  int i;
146  for (i=0;i<variables;i++) Print("_%d", m[i]);
147  PrintS(" ");
148 }
149 
150 void WriteMonoList (mon_list_entry *list)
151 {
152  mon_list_entry *ptr;
153  ptr=list;
154  while (ptr!=NULL)
155  {
156  WriteMono(ptr->mon);
157  ptr=ptr->next;
158  }
159 }
160 
161 void Info ()
162 {
163  int i;
164  row_list_entry *curptr;
165  modp_number *row;
167  char ch;
168  cout << endl << "quotient basis: ";
169  WriteMonoList (base_list);
170  cout << endl << "leading terms: ";
171  WriteMonoList (lt_list);
172  cout << endl << "to be checked: ";
173  WriteMonoList (check_list);
174 #ifdef shmat
175  cout << endl << "Matrix:" << endl;
176  cout << " ";
177  for (i=0;i<final_base_dim;i++)
178  {
179  WriteMono (column_name[i]);
180  }
181  cout << endl;
182  curptr=row_list;
183  while (curptr!=NULL)
184  {
185  row=curptr->row_matrix;
186  solve=curptr->row_solve;
187  for (i=0;i<final_base_dim;i++)
188  {
189  cout << *row << " , ";
190  row++;
191  }
192  cout << " ";
193  for (i=0;i<final_base_dim;i++)
194  {
195  cout << *solve << " , ";
196  solve++;
197  }
198  curptr=curptr->next;
199  cout << endl;}
200  cout << "Special row: Solve row:" << endl;
201  row=my_row;
202  for (i=0;i<final_base_dim;i++)
203  {
204  cout << *row << " , ";
205  row++;
206  }
207  cout << " ";
208  row=my_solve_row;
209  for (i=0;i<final_base_dim;i++)
210  {
211  cout << *row << " , ";
212  row++;
213  }
214 #endif
215  cout << endl;
216  cin >> ch;
217  cout << endl << endl;
218 }
219 #endif
220 
221 static modp_number OneInverse(modp_number a,modp_number p) // computes inverse of d mod p without using tables
222 {
223  long u, v, u0, u1, u2, q, r;
224  u1=1; u2=0;
225  u = a; v = p;
226  while (v != 0)
227  {
228  q = u / v;
229  r = u % v;
230  u = v;
231  v = r;
232  u0 = u2;
233  u2 = u1 - q*u2;
234  u1 = u0;
235  }
236  if (u1 < 0) u1=u1+p;
237 // now it should be ok, but for any case...
238  if ((u1<0)||((u1*a)%p!=1))
239  {
240  PrintS("?");
241  modp_number i;
242  for (i=1;i<p;i++)
243  {
244  if ((a*i)%p==1) return i;
245  }
246  }
247  return (modp_number)u1;
248 }
249 
250 static int CalcBaseDim () // returns the maximal (expected) dimension of quotient space
251 {
252  int c;
253  int i,j;
254  int s=0;
255  max_coord=1;
256  for (i=0;i<n_points;i++) max_coord=max_coord+multiplicity[i];
257  for (j=0;j<n_points;j++)
258  {
259  c=1;
260  for (i=1;i<=variables;i++)
261  {
262  c=(c*(multiplicity[j]+i-1))/i;
263  }
264  s=s+c;
265  }
266  return s;
267 }
268 
269 static bool EqualMon (mono_type m1, mono_type m2) // compares two monomials, true if equal, false otherwise
270 {
271  int i;
272  for (i=0;i<variables;i++)
273  if (m1[i]!=m2[i]) return false;;
274  return true;
275 }
276 
277 static exponent MonDegree (mono_type mon) // computes the degree of a monomial
278 {
279  exponent v=0;
280  int i;
281  for (i=0;i<variables;i++) v=v+mon[i];
282  return v;
283 }
284 
285 static bool Greater (mono_type m1, mono_type m2) // true if m1 > m2, false otherwise. uses Singular comparing function
286 {
287  for (int j = variables; j; j--)
288  {
289  pSetExp(comparizon_p1, j, m1[j-1]);
290  pSetExp(comparizon_p2, j, m2[j-1]);
291  }
295  return res;
296 }
297 
298 static mon_list_entry* MonListAdd (mon_list_entry *list, mono_type mon) // adds one entry to the list of monomials structure
299 {
300  mon_list_entry *curptr=list;
301  mon_list_entry *prevptr=NULL;
302  mon_list_entry *temp;
303 
304  while (curptr!=NULL)
305  {
306  if (EqualMon (mon,(*curptr).mon)) return list;
307  if (Greater ((*curptr).mon,mon)) break;
308  prevptr=curptr;
309  curptr=curptr->next;
310  }
311  temp=(mon_list_entry*)omAlloc0(sizeof(mon_list_entry));
312  (*temp).next=curptr;
313  (*temp).mon=(exponent*)omAlloc(sizeof(exponent)*variables);
314  memcpy(temp->mon,mon,sizeof(exponent)*variables);
315  if (prevptr==NULL) return temp;
316  else
317  {
318  prevptr->next=temp;
319  return list;
320  }
321 }
322 
323 static mono_type MonListElement (mon_list_entry *list, int n) // returns ith element from list of monomials - no error checking!
324 {
325  mon_list_entry *cur=list;
326  int i;
327  for (i=0;i<n;i++) cur=cur->next;
328  return cur->mon;
329 }
330 
331 static mono_type ZeroMonomial () // allocates memory for new monomial with all enries equal 0
332 {
333  mono_type m;
334  m=(exponent*)omAlloc0(sizeof(exponent)*variables);
335  return m;
336 }
337 
338 static void GeneralInit () // general initialization
339 {
340  int i,j;
342  for (i=0;i<n_points;i++)
343  {
345  for (j=0;j<variables;j++) points[i][j]=(modp_number*)omAlloc0(sizeof(modp_number)*(max_coord));
346  }
347  condition_list=(condition_type*)omAlloc0(sizeof(condition_type)*final_base_dim);
348  for (i=0;i<final_base_dim;i++) condition_list[i].mon=(exponent*)omAlloc0(sizeof(exponent)*variables);
350  for (i=0;i<n_points;i++) modp_points[i]=(modp_number*)omAlloc0(sizeof(modp_number)*variables);
351  if (!only_modp)
352  {
354  for (i=0;i<n_points;i++)
355  {
356  q_points[i]=(mpq_t*)omAlloc(sizeof(mpq_t)*variables);
357  for (j=0;j<variables;j++) mpq_init(q_points[i][j]);
358  }
360  for (i=0;i<n_points;i++)
361  {
362  int_points[i]=(mpz_t*)omAlloc(sizeof(mpz_t)*variables);
363  for (j=0;j<variables;j++) mpz_init(int_points[i][j]);
364  }
365  }
367  for (i=0;i<n_points;i++)
368  {
369  coord_exist[i]=(bool*)omAlloc0(sizeof(bool)*variables);
370  }
372  for (i=0;i<final_base_dim;i++) generic_column_name[i]=ZeroMonomial ();
373  good_primes=0;
374  bad_primes=1;
376  if (!only_modp)
377  {
378  polycoef=(mpz_t*)omAlloc(sizeof(mpz_t)*(final_base_dim+1));
379  polyexp=(mono_type*)omAlloc(sizeof(mono_type)*(final_base_dim+1));
380  for (i=0;i<=final_base_dim;i++)
381  {
382  mpz_init(polycoef[i]);
383  polyexp[i]=ZeroMonomial ();
384  }
385  mpz_init(common_denom);
386  }
387 
388 // set all globally used lists to be initially empty
391  gen_list=NULL;
392  n_results=0;
393 
394 // creates polynomials to compare by Singular
397 }
398 
399 static void InitProcData () // initialization for procedures which do computations mod p
400 {
401  int i;
406  for (i=0;i<final_base_dim;i++) column_name[i]=ZeroMonomial ();
408  modp_number *gen_table;
409  modp_number pos_gen;
410  bool gen_ok;
412 
413 // produces table of modp inverts by finding a generator of (Z_myp*,*)
414  gen_table=(modp_number*)omAlloc(sizeof(modp_number)*myp);
415  gen_table[1]=1;
416  for (pos_gen=2;pos_gen<myp;pos_gen++)
417  {
418  gen_ok=true;
419  for (i=2;i<myp;i++)
420  {
421  gen_table[i]=modp_mul(gen_table[i-1],pos_gen);
422  if (gen_table[i]==1)
423  {
424  gen_ok=false;
425  break;
426  }
427  }
428  if (gen_ok) break;
429  }
430  for (i=2;i<myp;i++) modp_Reverse[gen_table[i]]=gen_table[myp-i+1];
431  modp_Reverse[1]=1;
432  omFree(gen_table);
433 }
434 
435 static mon_list_entry* FreeMonList (mon_list_entry *list) // destroys a list of monomials, returns NULL pointer
436 {
437  mon_list_entry *ptr;
438  mon_list_entry *pptr;
439  ptr=list;
440  while (ptr!=NULL)
441  {
442  pptr=ptr->next;
443  omFree(ptr->mon);
444  omFree(ptr);
445  ptr=pptr;}
446  return NULL;
447 }
448 
449 static void GeneralDone () // to be called before exit to free memory
450 {
451  int i,j;
452  for (i=0;i<n_points;i++)
453  {
454  for (j=0;j<variables;j++)
455  {
456  omFree(points[i][j]);
457  }
458  omFree(points[i]);
459  }
460  omFree(points);
461  for (i=0;i<final_base_dim;i++) omFree(condition_list[i].mon);
462  omFree(condition_list);
463  for (i=0;i<n_points;i++) omFree(modp_points[i]);
465  if (!only_modp)
466  {
467  for (i=0;i<n_points;i++)
468  {
469  for (j=0;j<variables;j++) mpq_clear(q_points[i][j]);
470  omFree(q_points[i]);
471  }
472  omFree(q_points);
473  for (i=0;i<n_points;i++)
474  {
475  for (j=0;j<variables;j++) mpz_clear(int_points[i][j]);
476  omFree(int_points[i]);
477  }
480  for (i=0;i<=final_base_dim;i++)
481  {
482  mpz_clear(polycoef[i]);
483  omFree(polyexp[i]);
484  }
485  omFree(polycoef);
486  omFree(polyexp);
487  if (!only_modp) mpz_clear(common_denom);
488  }
489  for (i=0;i<final_base_dim;i++)
490  {
492  }
494  for (i=0;i<n_points;i++) omFree(coord_exist[i]);
498 }
499 
500 static void FreeProcData () // to be called after one modp computation to free memory
501 {
502  int i;
503  row_list_entry *ptr;
504  row_list_entry *pptr;
508  omFree(my_row);
509  my_row=NULL;
512  ptr=row_list;
513  while (ptr!=NULL)
514  {
515  pptr=ptr->next;
516  omFree(ptr->row_matrix);
517  omFree(ptr->row_solve);
518  omFree(ptr);
519  ptr=pptr;
520  }
521  row_list=NULL;
522  for (i=0;i<final_base_dim;i++) omFree(column_name[i]);
525 }
526 
527 static void modp_Evaluate(modp_number *ev, mono_type mon, condition_type con) // evaluates monomial on condition (modp)
528 {
529  int i;
530  *ev=0;
531  for (i=0;i<variables;i++)
532  if (con.mon[i] > mon[i]) return ;;
533  *ev=1;
534  int j,k;
535  mono_type mn;
536  mn=(exponent*)omAlloc(sizeof(exponent)*variables);
537  memcpy(mn,mon,sizeof(exponent)*variables);
538  for (k=0;k<variables;k++)
539  {
540  for (j=1;j<=con.mon[k];j++) // this loop computes the coefficient from derivation
541  {
542  *ev=modp_mul(*ev,mn[k]);
543  mn[k]--;
544  }
545  *ev=modp_mul(*ev,points[con.point_ref][k][mn[k]]);
546  }
547  omFree(mn);
548 }
549 
550 static void int_Evaluate(mpz_t ev, mono_type mon, condition_type con) // (***) evaluates monomial on condition for integer numbers
551 {
552  int i;
553  mpz_set_ui(ev,0);
554  for (i=0;i<variables;i++)
555  if (con.mon[i] > mon[i]) return ;;
556  mpz_set_ui(ev,1);
557  mpz_t mon_conv;
558  mpz_init(mon_conv);
559  int j,k;
560  mono_type mn;
561  mn=(exponent*)omAlloc(sizeof(exponent)*variables);
562  memcpy(mn,mon,sizeof(exponent)*variables);
563  for (k=0;k<variables;k++)
564  {
565  for (j=1;j<=con.mon[k];j++) // this loop computes the coefficient from derivation
566  {
567  mpz_set_si(mon_conv,mn[k]); // (***)
568  mpz_mul(ev,ev,mon_conv);
569  mn[k]--;
570  }
571  for (j=1;j<=mn[k];j++) mpz_mul(ev,ev,int_points[con.point_ref][k]); // this loop computes the product of coordinate
572  }
573  omFree(mn);
574  mpz_clear(mon_conv);
575 }
576 
577 static void ProduceRow(mono_type mon) // produces a row for monomial - first part is an evaluated part, second 0 to obtain the coefs of dependence
578 {
579  modp_number *row;
580  condition_type *con;
581  int i;
582  row=my_row;
583  con=condition_list;
584  for (i=0;i<final_base_dim;i++)
585  {
586  modp_Evaluate (row, mon, *con);
587  row++;
588  con++;
589  }
590  row=my_solve_row;
591  for (i=0;i<final_base_dim;i++)
592  {
593  *row=0;
594  row++;
595  }
596 }
597 
598 static void IntegerPoints () // produces integer points from rationals by scaling the coordinate system
599 {
600  int i,j;
601  mpz_set_ui(common_denom,1); // this is common scaling factor
602  for (i=0;i<n_points;i++)
603  {
604  for (j=0;j<variables;j++)
605  {
606  mpz_lcm(common_denom,common_denom,mpq_denref(q_points[i][j]));
607  }
608  }
609  mpq_t temp;
610  mpq_init(temp);
611  mpq_t denom_q;
612  mpq_init(denom_q);
613  mpq_set_z(denom_q,common_denom);
614  for (i=0;i<n_points;i++)
615  {
616  for (j=0;j<variables;j++)
617  {
618  mpq_mul(temp,q_points[i][j],denom_q); // multiplies by common denominator
619  mpz_set(int_points[i][j],mpq_numref(temp)); // and changes into integer
620  }
621  }
622  mpq_clear(temp);
623  mpq_clear(denom_q);
624 }
625 
626 static void int_PrepareProducts () // prepares coordinates of points and products for modp case (from integer coefs)
627 {
628  int i,j,k;
629  mpz_t big_myp;
630  mpz_init(big_myp);
631  mpz_set_si(big_myp,myp);
632  mpz_t temp;
633  mpz_init(temp);
634  for (i=0;i<n_points;i++)
635  {
636  for (j=0;j<variables;j++)
637  {
638  mpz_mod(temp,int_points[i][j],big_myp); // coordinate is now modulo myp
639  points[i][j][1]=mpz_get_ui(temp);
640  points[i][j][0]=1;
641  for (k=2;k<max_coord;k++) points[i][j][k]=modp_mul(points[i][j][k-1],points[i][j][1]);
642  }
643  }
644  mpz_mod(temp,common_denom,big_myp); // computes the common denominator (from rational data) modulo myp
645  denom_divisible=(mpz_sgn(temp)==0); // checks whether it is divisible by modp
646  mpz_clear(temp);
647  mpz_clear(big_myp);
648 }
649 
650 static void modp_PrepareProducts () // prepares products for modp computation from modp data
651 {
652  int i,j,k;
653  for (i=0;i<n_points;i++)
654  {
655  for (j=0;j<variables;j++)
656  {
657  points[i][j][1]=modp_points[i][j];
658  points[i][j][0]=1;
659  for (k=2;k<max_coord;k++) points[i][j][k]=modp_mul(points[i][j][k-1],points[i][j][1]);
660  }
661  }
662 }
663 
664 static void MakeConditions () // prepares a list of conditions from list of multiplicities
665 {
666  condition_type *con;
667  int n,i,k;
668  mono_type mon;
669  mon=ZeroMonomial ();
670  con=condition_list;
671  for (n=0;n<n_points;n++)
672  {
673  for (i=0;i<variables;i++) mon[i]=0;
674  while (mon[0]<multiplicity[n])
675  {
676  if (MonDegree (mon) < multiplicity[n])
677  {
678  memcpy(con->mon,mon,sizeof(exponent)*variables);
679  con->point_ref=n;
680  con++;
681  }
682  k=variables-1;
683  mon[k]++;
684  while ((k>0)&&(mon[k]>=multiplicity[n]))
685  {
686  mon[k]=0;
687  k--;
688  mon[k]++;
689  }
690  }
691  }
692  omFree(mon);
693 }
694 
695 static void ReduceRow () // reduces my_row by previous rows, does the same for solve row
696 {
697  if (row_list==NULL) return ;
698  row_list_entry *row_ptr;
699  modp_number *cur_row_ptr;
700  modp_number *solve_row_ptr;
701  modp_number *my_row_ptr;
702  modp_number *my_solve_row_ptr;
703  int first_col;
704  int i;
705  modp_number red_val;
706  modp_number mul_val;
707 #ifdef integerstrategy
708  modp_number *m_row_ptr;
709  modp_number prep_val;
710 #endif
711  row_ptr=row_list;
712  while (row_ptr!=NULL)
713  {
714  cur_row_ptr=row_ptr->row_matrix;
715  solve_row_ptr=row_ptr->row_solve;
716  my_row_ptr=my_row;
717  my_solve_row_ptr=my_solve_row;
718  first_col=row_ptr->first_col;
719  cur_row_ptr=cur_row_ptr+first_col; // reduction begins at first_col position
720  my_row_ptr=my_row_ptr+first_col;
721  red_val=*my_row_ptr;
722  if (red_val!=0)
723  {
724 #ifdef integerstrategy
725  prep_val=*cur_row_ptr;
726  if (prep_val!=1)
727  {
728  m_row_ptr=my_row;
729  for (i=0;i<final_base_dim;i++)
730  {
731  if (*m_row_ptr!=0) *m_row_ptr=modp_mul(*m_row_ptr,prep_val);
732  m_row_ptr++;
733  }
734  m_row_ptr=my_solve_row;
735  for (i=0;i<last_solve_column;i++)
736  {
737  if (*m_row_ptr!=0) *m_row_ptr=modp_mul(*m_row_ptr,prep_val);
738  m_row_ptr++;
739  }
740  }
741 #endif
742  for (i=first_col;i<final_base_dim;i++)
743  {
744  if (*cur_row_ptr!=0)
745  {
746  mul_val=modp_mul(*cur_row_ptr,red_val);
747  *my_row_ptr=modp_sub(*my_row_ptr,mul_val);
748  }
749  my_row_ptr++;
750  cur_row_ptr++;
751  }
752  for (i=0;i<=last_solve_column;i++) // last_solve_column stores the last non-enmpty entry in solve matrix
753  {
754  if (*solve_row_ptr!=0)
755  {
756  mul_val=modp_mul(*solve_row_ptr,red_val);
757  *my_solve_row_ptr=modp_sub(*my_solve_row_ptr,mul_val);
758  }
759  solve_row_ptr++;
760  my_solve_row_ptr++;
761  }
762  }
763  row_ptr=row_ptr->next;
764 #if 0 /* only debugging */
765  PrintS("reduction by row ");
766  Info ();
767 #endif
768  }
769 }
770 
771 static bool RowIsZero () // check whether a row is zero
772 {
773  bool zero=1;
774  int i;
775  modp_number *row;
776  row=my_row;
777  for (i=0;i<final_base_dim;i++)
778  {
779  if (*row!=0) {zero=0; break;}
780  row++;
781  }
782  return zero;
783 }
784 
785 static bool DivisibleMon (mono_type m1, mono_type m2) // checks whether m1 is divisible by m2
786 {
787  int i;
788  for (i=0;i<variables;i++)
789  if (m1[i]>m2[i]) return false;;
790  return true;
791 }
792 
793 static void ReduceCheckListByMon (mono_type m) // from check_list monomials divisible by m are thrown out
794 {
795  mon_list_entry *c_ptr;
796  mon_list_entry *p_ptr;
797  mon_list_entry *n_ptr;
798  c_ptr=check_list;
799  p_ptr=NULL;
800  while (c_ptr!=NULL)
801  {
802  if (DivisibleMon (m,c_ptr->mon))
803  {
804  if (p_ptr==NULL)
805  check_list=c_ptr->next;
806  else
807  p_ptr->next=c_ptr->next;
808  n_ptr=c_ptr->next;
809  omFree(c_ptr->mon);
810  omFree(c_ptr);
811  c_ptr=n_ptr;
812  }
813  else
814  {
815  p_ptr=c_ptr;
816  c_ptr=c_ptr->next;
817  }
818  }
819 }
820 
821 static void TakeNextMonomial (mono_type mon) // reads first monomial from check_list, then it is deleted
822 {
823  mon_list_entry *n_check_list;
824  if (check_list!=NULL)
825  {
826  memcpy(mon,check_list->mon,sizeof(exponent)*variables);
827  n_check_list=check_list->next;
828  omFree(check_list->mon);
830  check_list=n_check_list;
831  }
832 }
833 
834 static void UpdateCheckList (mono_type m) // adds all monomials which are "next" to m (divisible by m and degree +1)
835 {
836  int i;
837  for (i=0;i<variables;i++)
838  {
839  m[i]++;
841  m[i]--;
842  }
843 }
844 
845 static void ReduceCheckListByLTs () // deletes all monomials from check_list which are divisible by one of the leading terms
846 {
847  mon_list_entry *ptr;
848  ptr=lt_list;
849  while (ptr!=NULL)
850  {
851  ReduceCheckListByMon (ptr->mon);
852  ptr=ptr->next;
853  }
854 }
855 
856 static void RowListAdd (int first_col, mono_type mon) // puts a row into matrix
857 {
858  row_list_entry *ptr;
859  row_list_entry *pptr;
860  row_list_entry *temp;
861  ptr=row_list;
862  pptr=NULL;
863  while (ptr!=NULL)
864  {
865 #ifndef unsortedmatrix
866  if ( first_col <= ptr->first_col ) break;
867 #endif
868  pptr=ptr;
869  ptr=ptr->next;
870  }
871  temp=(row_list_entry*)omAlloc0(sizeof(row_list_entry));
872  (*temp).row_matrix=(modp_number*)omAlloc0(sizeof(modp_number)*final_base_dim);
873  memcpy((*temp).row_matrix,my_row,sizeof(modp_number)*(final_base_dim));
874  (*temp).row_solve=(modp_number*)omAlloc0(sizeof(modp_number)*final_base_dim);
875  memcpy((*temp).row_solve,my_solve_row,sizeof(modp_number)*(final_base_dim));
876  (*temp).first_col=first_col;
877  temp->next=ptr;
878  if (pptr==NULL) row_list=temp; else pptr->next=temp;
879  memcpy(column_name[first_col],mon,sizeof(exponent)*variables);
880 }
881 
882 
883 static void PrepareRow (mono_type mon) // prepares a row and puts it into matrix
884 {
885  modp_number *row;
886  int first_col=-1;
887  int col;
888  modp_number red_val=1;
889  row=my_row;
890 #ifdef integerstrategy
891  for (col=0;col<final_base_dim;col++)
892  {
893  if (*row!=0)
894  {
895  first_col=col;
896  break;
897  }
898  row++;
899  }
900  my_solve_row[first_col]=1;
901  if (first_col > last_solve_column) last_solve_column=first_col;
902 #else
903  for (col=0;col<final_base_dim;col++)
904  {
905  if (*row!=0)
906  {
907  first_col=col;
908  red_val=modp_Reverse[*row]; // first non-zero entry, should multiply all row by inverse to obtain 1
909  modp_denom=modp_mul(modp_denom,*row); // remembers the divisor
910  *row=1;
911  break;
912  }
913  row++;
914  }
915  my_solve_row[first_col]=1;
916  if (first_col > last_solve_column) last_solve_column=first_col;
917  if (red_val!=1)
918  {
919  row++;
920  for (col=first_col+1;col<final_base_dim;col++)
921  {
922  if (*row!=0) *row=modp_mul(*row,red_val);
923  row++;
924  }
925  row=my_solve_row;
926  for (col=0;col<=last_solve_column;col++)
927  {
928  if (*row!=0) *row=modp_mul(*row,red_val);
929  row++;
930  }
931  }
932 #endif
933  RowListAdd (first_col, mon);
934 }
935 
936 static void NewResultEntry () // creates an entry for new modp result
937 {
938  modp_result_entry *temp;
939  temp=(modp_result_entry*)omAlloc0(sizeof(modp_result_entry));
940  if (cur_result==NULL)
941  {
942  modp_result=temp;
943  temp->prev=NULL;
944  }
945  else
946  {
947  temp->prev=cur_result;
948  cur_result->next=temp;
949  }
950  cur_result=temp;
951  cur_result->next=NULL;
952  cur_result->p=myp;
953  cur_result->generator=NULL;
954  cur_result->n_generators=0;
955  n_results++;
956 }
957 
958 static void FreeResultEntry (modp_result_entry *e) // destroys the result entry, without worrying about where it is
959 {
960  generator_entry *cur_gen;
961  generator_entry *next_gen;
962  cur_gen=e->generator;
963  while (cur_gen!=NULL)
964  {
965  next_gen=cur_gen->next;
966  omFree(cur_gen->coef);
967  omFree(cur_gen->lt);
968  omFree(cur_gen);
969  cur_gen=next_gen;
970  }
971  omFree(e);
972 }
973 
974 
975 static void NewGenerator (mono_type mon) // new generator in modp comp found, shoul be stored on the list
976 {
977  generator_entry *cur_ptr;
978  generator_entry *prev_ptr;
979  generator_entry *temp;
980  cur_ptr=cur_result->generator;
981  prev_ptr=NULL;
982  while (cur_ptr!=NULL)
983  {
984  prev_ptr=cur_ptr;
985  cur_ptr=cur_ptr->next;
986  }
987  temp=(generator_entry*)omAlloc0(sizeof(generator_entry));
988  if (prev_ptr==NULL) cur_result->generator=temp;
989  else prev_ptr->next=temp;
990  temp->next=NULL;
991  temp->coef=(modp_number*)omAlloc0(sizeof(modp_number)*final_base_dim);
992  memcpy(temp->coef,my_solve_row,sizeof(modp_number)*final_base_dim);
993  temp->lt=ZeroMonomial ();
994  memcpy(temp->lt,mon,sizeof(exponent)*variables);
995  temp->ltcoef=1;
996  cur_result->n_generators++;
997 }
998 
999 static void MultGenerators () // before reconstructing, all denominators must be cancelled
1000 {
1001 #ifndef integerstrategy
1002  int i;
1003  generator_entry *cur_ptr;
1004  cur_ptr=cur_result->generator;
1005  while (cur_ptr!=NULL)
1006  {
1007  for (i=0;i<final_base_dim;i++)
1008  cur_ptr->coef[i]=modp_mul(cur_ptr->coef[i],modp_denom);
1009  cur_ptr->ltcoef=modp_denom;
1010  cur_ptr=cur_ptr->next;
1011  }
1012 #endif
1013 }
1014 #if 0 /* only debbuging */
1015 void PresentGenerator (int i) // only for debuging, writes a generator in its form in program
1016 {
1017  int j;
1018  modp_result_entry *cur_ptr;
1019  generator_entry *cur_gen;
1020  cur_ptr=modp_result;
1021  while (cur_ptr!=NULL)
1022  {
1023  cur_gen=cur_ptr->generator;
1024  for (j=0;j<i;j++) cur_gen=cur_gen->next;
1025  for (j=0;j<final_base_dim;j++)
1026  {
1027  Print("%d;", cur_gen->coef[j]);
1028  }
1029  PrintS(" and LT = ");
1030  WriteMono (cur_gen->lt);
1031  Print(" ( %d ) prime = %d\n", cur_gen->ltcoef, cur_ptr->p);
1032  cur_ptr=cur_ptr->next;
1033  }
1034 }
1035 #endif
1036 
1037 static modp_number TakePrime (modp_number /*p*/) // takes "previous" (smaller) prime
1038 {
1039  myp_index--;
1040  return cf_getSmallPrime(myp_index);
1041 }
1042 
1043 static void PrepareChinese (int n) // initialization for CRA
1044 {
1045  int i,k;
1046  modp_result_entry *cur_ptr;
1047  cur_ptr=modp_result;
1048  modp_number *congr_ptr;
1049  modp_number prod;
1051  congr=(modp_number*)omAlloc0(sizeof(modp_number)*n);
1052  congr_ptr=congr;
1053  while (cur_ptr!=NULL)
1054  {
1055  *congr_ptr=cur_ptr->p;
1056  cur_ptr=cur_ptr->next;
1057  congr_ptr++;
1058  }
1059  for (k=1;k<n;k++)
1060  {
1061  prod=congr[0]%congr[k];
1062  for (i=1;i<=k-1;i++) prod=(prod*congr[i])%congr[k];
1063  in_gamma[i]=OneInverse(prod,congr[k]);
1064  }
1065  mpz_init(bigcongr);
1066  mpz_set_si(bigcongr,congr[0]);
1067  for (k=1;k<n;k++) mpz_mul_ui(bigcongr,bigcongr,congr[k]);
1068 }
1069 
1070 static void CloseChinese () // after CRA
1071 {
1072  omFree(in_gamma);
1073  omFree(congr);
1074  mpz_clear(bigcongr);
1075 }
1076 
1077 static void ClearGCD () // divides polynomials in basis by gcd of coefficients
1078 {
1079  bool first_gcd=true;
1080  int i;
1081  mpz_t g;
1082  mpz_init(g);
1083  for (i=0;i<=final_base_dim;i++)
1084  {
1085  if (mpz_sgn(polycoef[i])!=0)
1086  {
1087  if (first_gcd)
1088  {
1089  first_gcd=false;
1090  mpz_set(g,polycoef[i]);
1091  }
1092  else
1093  mpz_gcd(g,g,polycoef[i]);
1094  }
1095  }
1096  for (i=0;i<=final_base_dim;i++) mpz_divexact(polycoef[i],polycoef[i],g);
1097  mpz_clear(g);
1098 }
1099 
1100 static void ReconstructGenerator (int ngen,int n) // recostruction of generator from various modp comp
1101 {
1102  int i,j,k;
1103  int coef;
1104  char *str=NULL;
1105  str=(char*)omAlloc0(sizeof(char)*1000);
1106  modp_result_entry *cur_ptr;
1107  generator_entry *cur_gen;
1108  modp_number *u;
1109  modp_number *v;
1110  modp_number temp;
1111  mpz_t sol;
1112  mpz_t nsol;
1113  mpz_init(sol);
1114  mpz_init(nsol);
1115  u=(modp_number*)omAlloc0(sizeof(modp_number)*n);
1116  v=(modp_number*)omAlloc0(sizeof(modp_number)*n);
1117  for (coef=0;coef<=final_base_dim;coef++)
1118  {
1119  i=0;
1120  cur_ptr=modp_result;
1121  while (cur_ptr!=NULL)
1122  {
1123  cur_gen=cur_ptr->generator;
1124  for (j=0;j<ngen;j++) cur_gen=cur_gen->next; // we have to take jth generator
1125  if (coef<final_base_dim) u[i]=cur_gen->coef[coef]; else u[i]=cur_gen->ltcoef;
1126  cur_ptr=cur_ptr->next;
1127  i++;
1128  }
1129  v[0]=u[0]; // now CRA begins
1130  for (k=1;k<n;k++)
1131  {
1132  temp=v[k-1];
1133  for (j=k-2;j>=0;j--) temp=(temp*congr[j]+v[j])%congr[k];
1134  temp=u[k]-temp;
1135  if (temp<0) temp=temp+congr[k];
1136  v[k]=(temp*in_gamma[k])%congr[k];
1137  }
1138  mpz_set_si(sol,v[n-1]);
1139  for (k=n-2;k>=0;k--)
1140  {
1141  mpz_mul_ui(sol,sol,congr[k]);
1142  mpz_add_ui(sol,sol,v[k]);
1143  } // now CRA ends
1144  mpz_sub(nsol,sol,bigcongr);
1145  int s=mpz_cmpabs(sol,nsol);
1146  if (s>0) mpz_set(sol,nsol); // chooses representation closer to 0
1147  mpz_set(polycoef[coef],sol);
1148  if (coef<final_base_dim)
1149  memcpy(polyexp[coef],generic_column_name[coef],sizeof(exponent)*variables);
1150  else
1151  memcpy(polyexp[coef],MonListElement (generic_lt,ngen),sizeof(exponent)*variables);
1152 #ifdef checksize
1153  size=mpz_sizeinbase(sol,10);
1154  if (size>maximal_size) maximal_size=size;
1155 #endif
1156  }
1157  omFree(u);
1158  omFree(v);
1159  omFree(str);
1160  ClearGCD ();
1161  mpz_clear(sol);
1162  mpz_clear(nsol);
1163 }
1164 
1165 static void Discard () // some unlucky prime occures
1166 {
1167  modp_result_entry *temp;
1168 #ifdef writemsg
1169  Print(", prime=%d", cur_result->p);
1170 #endif
1171  bad_primes++;
1172  if (good_primes>bad_primes)
1173  {
1174 #ifdef writemsg
1175  Print("-discarding this comp (%dth)\n", n_results);
1176 #endif
1177  temp=cur_result;
1178  cur_result=cur_result->prev;
1179  cur_result->next=NULL;
1180  n_results--;
1181  FreeResultEntry (temp);
1182  }
1183  else
1184  {
1185 #ifdef writemsg
1186  PrintS("-discarding ALL.\n");
1187 #endif
1188  int i;
1189  modp_result_entry *ntfree;
1190  generator_entry *cur_gen;
1191  temp=cur_result->prev;
1192  while (temp!=NULL)
1193  {
1194  ntfree=temp->prev;
1195  FreeResultEntry (temp);
1196  temp=ntfree;
1197  }
1199  cur_result->prev=NULL;
1200  n_results=1;
1201  good_primes=1;
1202  bad_primes=0;
1203  generic_n_generators=cur_result->n_generators;
1204  cur_gen=cur_result->generator;
1206  for (i=0;i<generic_n_generators;i++)
1207  {
1208  generic_lt=MonListAdd (generic_lt,cur_gen->lt);
1209  cur_gen=cur_gen->next;
1210  }
1211  for (i=0;i<final_base_dim;i++) memcpy(generic_column_name[i],column_name[i],sizeof(exponent)*variables);
1212  }
1213 }
1214 
1215 static void modp_SetColumnNames () // used by modp - sets column names, the old table will be destroyed
1216 {
1217  int i;
1218  for (i=0;i<final_base_dim;i++) memcpy(generic_column_name[i],column_name[i],sizeof(exponent)*variables);
1219 }
1220 
1221 static void CheckColumnSequence () // checks if scheme of computations is as generic one
1222 {
1223  int i;
1224  if (cur_result->n_generators!=generic_n_generators)
1225  {
1226 #ifdef writemsg
1227  PrintS("wrong number of generators occurred");
1228 #else
1229  if (protocol) PrintS("ng");
1230 #endif
1231  Discard ();
1232  return;
1233  }
1234  if (denom_divisible)
1235  {
1236 #ifdef writemsg
1237  PrintS("denom of coef divisible by p");
1238 #else
1239  if (protocol) PrintS("dp");
1240 #endif
1241  Discard ();
1242  return;
1243  }
1244  generator_entry *cur_gen;
1245  mon_list_entry *cur_mon;
1246  cur_gen=cur_result->generator;
1247  cur_mon=generic_lt;
1248  for (i=0;i<generic_n_generators;i++)
1249  {
1250  if (!EqualMon(cur_mon->mon,cur_gen->lt))
1251  {
1252 #ifdef writemsg
1253  PrintS("wrong leading term occurred");
1254 #else
1255  if (protocol) PrintS("lt");
1256 #endif
1257  Discard ();
1258  return;
1259  }
1260  cur_gen=cur_gen->next;
1261  cur_mon=cur_mon->next;
1262  }
1263  for (i=0;i<final_base_dim;i++)
1264  {
1266  {
1267 #ifdef writemsg
1268  PrintS("wrong seq of cols occurred");
1269 #else
1270  if (protocol) PrintS("sc");
1271 #endif
1272  Discard ();
1273  return;
1274  }
1275  }
1276  good_primes++;
1277 }
1278 #if 0 /* only debuggig */
1279 void WriteGenerator () // writes generator (only for debugging)
1280 {
1281  char *str;
1282  str=(char*)omAlloc0(sizeof(char)*1000);
1283  int i;
1284  for (i=0;i<=final_base_dim;i++)
1285  {
1286  str=mpz_get_str(str,10,polycoef[i]);
1287  PrintS(str);
1288  PrintS("*");
1289  WriteMono(polyexp[i]);
1290  PrintS(" ");
1291  }
1292  omFree(str);
1293  PrintLn();
1294 }
1295 #endif
1296 
1297 static bool CheckGenerator () // evaluates generator to check whether it is good
1298 {
1299  mpz_t val,sum;
1300  int con,i;
1301  mpz_init(val);
1302  mpz_init(sum);
1303  for (con=0;con<final_base_dim;con++)
1304  {
1305  mpz_set_ui(sum,0);
1306  for (i=0;i<=final_base_dim;i++)
1307  {
1308  int_Evaluate(val, polyexp[i], condition_list[con]);
1309  mpz_mul(val,val,polycoef[i]);
1310  mpz_add(sum,sum,val);
1311  }
1312  if (mpz_sgn(sum)!=0)
1313  {
1314  mpz_clear(val);
1315  mpz_clear(sum);
1316  return false;
1317  }
1318  }
1319  mpz_clear(val);
1320  mpz_clear(sum);
1321  return true;
1322 }
1323 
1324 static void ClearGenList ()
1325 {
1326  gen_list_entry *temp;
1327  int i;
1328  while (gen_list!=NULL)
1329  {
1330  temp=gen_list->next;
1331  for (i=0;i<=final_base_dim;i++)
1332  {
1333  mpz_clear(gen_list->polycoef[i]);
1334  omFree(gen_list->polyexp[i]);
1335  }
1336  omFree(gen_list->polycoef);
1337  omFree(gen_list->polyexp);
1338  omFree(gen_list);
1339  gen_list=temp;
1340  }
1341 }
1342 
1343 static void UpdateGenList ()
1344 {
1345  gen_list_entry *temp,*prev;
1346  int i,j;
1347  prev=NULL;
1348  temp=gen_list;
1349  exponent deg;
1350  for (i=0;i<=final_base_dim;i++)
1351  {
1352  deg=MonDegree(polyexp[i]);
1353  for (j=0;j<deg;j++)
1354  {
1355  mpz_mul(polycoef[i],polycoef[i],common_denom);
1356  }
1357  }
1358  ClearGCD ();
1359  while (temp!=NULL)
1360  {
1361  prev=temp;
1362  temp=temp->next;
1363  }
1364  temp=(gen_list_entry*)omAlloc0(sizeof(gen_list_entry));
1365  if (prev==NULL) gen_list=temp; else prev->next=temp;
1366  temp->next=NULL;
1367  temp->polycoef=(mpz_t*)omAlloc(sizeof(mpz_t)*(final_base_dim+1));
1368  temp->polyexp=(mono_type*)omAlloc(sizeof(mono_type)*(final_base_dim+1));
1369  for (i=0;i<=final_base_dim;i++)
1370  {
1371  mpz_init(temp->polycoef[i]);
1372  mpz_set(temp->polycoef[i],polycoef[i]);
1373  temp->polyexp[i]=ZeroMonomial ();
1374  memcpy(temp->polyexp[i],polyexp[i],sizeof(exponent)*variables);
1375  }
1376 }
1377 
1378 #if 0 /* only debugging */
1379 void ShowGenList ()
1380 {
1381  gen_list_entry *temp;
1382  int i;
1383  char *str;
1384  str=(char*)omAlloc0(sizeof(char)*1000);
1385  temp=gen_list;
1386  while (temp!=NULL)
1387  {
1388  PrintS("generator: ");
1389  for (i=0;i<=final_base_dim;i++)
1390  {
1391  str=mpz_get_str(str,10,temp->polycoef[i]);
1392  PrintS(str);
1393  PrintS("*");
1394  WriteMono(temp->polyexp[i]);
1395  }
1396  PrintLn();
1397  temp=temp->next;
1398  }
1399  omFree(str);
1400 }
1401 #endif
1402 
1403 
1404 static void modp_Main ()
1405 {
1406  mono_type cur_mon;
1407  cur_mon= ZeroMonomial ();
1408  modp_denom=1;
1409  bool row_is_zero;
1410 
1411 #if 0 /* only debugging */
1412  Info ();
1413 #endif
1414 
1415  while (check_list!=NULL)
1416  {
1417  TakeNextMonomial (cur_mon);
1418  ProduceRow (cur_mon);
1419 #if 0 /* only debugging */
1420  cout << "row produced for monomial ";
1421  WriteMono (cur_mon);
1422  cout << endl;
1423  Info ();
1424 #endif
1425  ReduceRow ();
1426  row_is_zero = RowIsZero ();
1427  if (row_is_zero)
1428  {
1429  lt_list=MonListAdd (lt_list,cur_mon);
1430  ReduceCheckListByMon (cur_mon);
1431  NewGenerator (cur_mon);
1432 #if 0 /* only debugging */
1433  cout << "row is zero - linear dependence found (should be seen in my_solve_row)" << endl;
1434  cout << "monomial added to leading terms list" << endl;
1435  cout << "check list updated" << endl;
1436  Info ();
1437 #endif
1438  }
1439  else
1440  {
1441  base_list= MonListAdd (base_list,cur_mon);
1442  UpdateCheckList (cur_mon);
1444 #if 0 /* only debugging */
1445  cout << "row is non-zero" << endl;
1446  cout << "monomial added to quotient basis list" << endl;
1447  cout << "new monomials added to check list" << endl;
1448  cout << "check list reduced by monomials from leading term list" << endl;
1449  Info ();
1450 #endif
1451  PrepareRow (cur_mon);
1452 #if 0 /* only debugging */
1453  cout << "row prepared and put into matrix" << endl;
1454  Info ();
1455 #endif
1456  }
1457  }
1458  omFree(cur_mon);
1459 }
1460 
1461 static void ResolveCoeff (mpq_t c, number m)
1462 {
1463  if ((long)m & SR_INT)
1464  {
1465  long m_val=SR_TO_INT(m);
1466  mpq_set_si(c,m_val,1);
1467  }
1468  else
1469  {
1470  if (m->s<2)
1471  {
1472  mpz_set(mpq_numref(c),m->z);
1473  mpz_set(mpq_denref(c),m->n);
1474  mpq_canonicalize(c);
1475  }
1476  else
1477  {
1478  mpq_set_z(c,m->z);
1479  }
1480  }
1481 }
1482 
1483 ideal interpolation(const std::vector<ideal>& L, intvec *v)
1484 {
1485  protocol=TEST_OPT_PROT; // should be set if option(prot) is enabled
1486 
1487  bool data_ok=true;
1488 
1489  // reading the ring data ***************************************************
1490  if ((currRing==NULL) || ((!rField_is_Zp (currRing))&&(!rField_is_Q (currRing))))
1491  {
1492  WerrorS("coefficient field should be Zp or Q!");
1493  return NULL;
1494  }
1495  if ((currRing->qideal)!=NULL)
1496  {
1497  WerrorS("quotient ring not supported!");
1498  return NULL;
1499  }
1500  if ((currRing->OrdSgn)!=1)
1501  {
1502  WerrorS("ordering must be global!");
1503  return NULL;
1504  }
1505  n_points=v->length ();
1506  if (n_points!=L.size())
1507  {
1508  WerrorS("list and intvec must have the same length!");
1509  return NULL;
1510  }
1511  assume( n_points > 0 );
1512  variables=currRing->N;
1514  if (only_modp) myp=rChar(currRing);
1515  // ring data read **********************************************************
1516 
1517 
1518  multiplicity=(int*)malloc(sizeof(int)*n_points); // TODO: use omalloc!
1519  int i;
1520  for (i=0;i<n_points;i++) multiplicity[i]=(*v)[i];
1521 
1523 
1524 #ifdef writemsg
1525  Print("number of variables: %d\n", variables);
1526  Print("number of points: %d\n", n_points);
1527  PrintS("multiplicities: ");
1528  for (i=0;i<n_points;i++) Print("%d ", multiplicity[i]);
1529  PrintLn();
1530  Print("general initialization for dimension %d ...\n", final_base_dim);
1531 #endif
1532 
1533  GeneralInit ();
1534 
1535 // reading coordinates of points from ideals **********************************
1536  mpq_t divisor;
1537  if (!only_modp) mpq_init(divisor);
1538  int j;
1539  for(i=0; i<L.size();i++)
1540  {
1541  ideal I = L[i];
1542  for(j=0;j<IDELEMS(I);j++)
1543  {
1544  poly p=I->m[j];
1545  if (p!=NULL)
1546  {
1547  poly ph=pHead(p);
1548  int pcvar=pVar(ph);
1549  if (pcvar!=0)
1550  {
1551  pcvar--;
1552  if (coord_exist[i][pcvar])
1553  {
1554  Print("coordinate %d for point %d initialized twice!\n",pcvar+1,i+1);
1555  data_ok=false;
1556  }
1557  number m;
1558  m=pGetCoeff(p); // possible coefficient standing by a leading monomial
1559  if (!only_modp) ResolveCoeff (divisor,m);
1560  number n;
1561  if (pNext(p)!=NULL) n=pGetCoeff(pNext(p));
1562  else n=nInit(0);
1563  if (only_modp)
1564  {
1565  n=nInpNeg(n);
1566  n=nDiv(n,m);
1567  modp_points[i][pcvar]=(int)((long)n);
1568  }
1569  else
1570  {
1571  ResolveCoeff (q_points[i][pcvar],n);
1572  mpq_neg(q_points[i][pcvar],q_points[i][pcvar]);
1573  mpq_div(q_points[i][pcvar],q_points[i][pcvar],divisor);
1574  }
1575  coord_exist[i][pcvar]=true;
1576  }
1577  else
1578  {
1579  PrintS("not a variable? ");
1580  wrp(p);
1581  PrintLn();
1582  data_ok=false;
1583  }
1584  pDelete(&ph);
1585  }
1586  }
1587  }
1588  if (!only_modp) mpq_clear(divisor);
1589  // data from ideal read *******************************************************
1590 
1591  // ckecking if all coordinates are initialized
1592  for (i=0;i<n_points;i++)
1593  {
1594  for (j=0;j<variables;j++)
1595  {
1596  if (!coord_exist[i][j])
1597  {
1598  Print("coordinate %d for point %d not known!\n",j+1,i+1);
1599  data_ok=false;
1600  }
1601  }
1602  }
1603 
1604  if (!data_ok)
1605  {
1606  GeneralDone();
1607  WerrorS("data structure is invalid");
1608  return NULL;
1609  }
1610 
1611  if (!only_modp) IntegerPoints ();
1612  MakeConditions ();
1613 #ifdef writemsg
1614  PrintS("done.\n");
1615 #else
1616  if (protocol) Print("[vdim %d]",final_base_dim);
1617 #endif
1618 
1619 
1620 // main procedure *********************************************************************
1621  int modp_cycles=10;
1622  bool correct_gen=false;
1623  if (only_modp) modp_cycles=1;
1625 
1626  while ((!correct_gen)&&(myp_index>1))
1627  {
1628 #ifdef writemsg
1629  Print("trying %d cycles mod p...\n",modp_cycles);
1630 #else
1631  if (protocol) Print("(%d)",modp_cycles);
1632 #endif
1633  while ((n_results<modp_cycles)&&(myp_index>1)) // some computations mod p
1634  {
1635  if (!only_modp) myp=TakePrime (myp);
1636  NewResultEntry ();
1637  InitProcData ();
1639 
1640  modp_Main ();
1641 
1642  if (!only_modp)
1643  {
1644  MultGenerators ();
1646  }
1647  else
1648  {
1650  }
1651  FreeProcData ();
1652  }
1653 
1654  if (!only_modp)
1655  {
1656  PrepareChinese (modp_cycles);
1657  correct_gen=true;
1658  for (i=0;i<generic_n_generators;i++)
1659  {
1660  ReconstructGenerator (i,modp_cycles);
1661  correct_gen=CheckGenerator ();
1662  if (!correct_gen)
1663  {
1664 #ifdef writemsg
1665  PrintS("wrong generator!\n");
1666 #else
1667 // if (protocol) PrintS("!g");
1668 #endif
1669  ClearGenList ();
1670  break;
1671  }
1672  else
1673  {
1674  UpdateGenList ();
1675  }
1676  }
1677 #ifdef checksize
1678  Print("maximal size of output: %d, precision bound: %d.\n",maximalsize,mpz_sizeinbase(bigcongr,10));
1679 #endif
1680  CloseChinese ();
1681  modp_cycles=modp_cycles+10;
1682  }
1683  else
1684  {
1685  correct_gen=true;
1686  }
1687  }
1688 // end of main procedure ************************************************************************************
1689 
1690 #ifdef writemsg
1691  PrintS("computations finished.\n");
1692 #else
1693  if (protocol) PrintLn();
1694 #endif
1695 
1696  if (!correct_gen)
1697  {
1698  GeneralDone ();
1699  ClearGenList ();
1700  WerrorS("internal error - coefficient too big!");
1701  return NULL;
1702  }
1703 
1704 // passing data to ideal *************************************************************************************
1705  ideal ret;
1706 
1707  if (only_modp)
1708  {
1709  mono_type mon;
1710  ret=idInit(modp_result->n_generators,1);
1711  generator_entry *cur_gen=modp_result->generator;
1712  for(i=0;i<IDELEMS(ret);i++)
1713  {
1714  poly p,sum;
1715  sum=NULL;
1716  int a;
1717  int cf;
1718  for (a=final_base_dim;a>=0;a--)
1719  {
1720  if (a==final_base_dim) cf=cur_gen->ltcoef; else cf=cur_gen->coef[a];
1721  if (cf!=0)
1722  {
1723  p=pISet(cf);
1724  if (a==final_base_dim) mon=cur_gen->lt; else mon=generic_column_name[a];
1725  for (j=0;j<variables;j++) pSetExp(p,j+1,mon[j]);
1726  pSetm(p);
1727  sum=pAdd(sum,p);
1728  }
1729  }
1730  ret->m[i]=sum;
1731  cur_gen=cur_gen->next;
1732  }
1733  }
1734  else
1735  {
1737  gen_list_entry *temp=gen_list;
1738  for(i=0;i<IDELEMS(ret);i++)
1739  {
1740  poly p,sum;
1741  sum=NULL;
1742  int a;
1743  for (a=final_base_dim;a>=0;a--) // build one polynomial
1744  {
1745  if (mpz_sgn(temp->polycoef[a])!=0)
1746  {
1747  number n=ALLOC_RNUMBER();
1748 #ifdef LDEBUG
1749  n->debug=123456;
1750 #endif
1751  mpz_init_set(n->z,temp->polycoef[a]);
1752  n->s=3;
1753  n_Normalize(n, currRing->cf);
1754  p=pNSet(n); //a monomial
1755  for (j=0;j<variables;j++) pSetExp(p,j+1,temp->polyexp[a][j]);
1756  pSetm(p); // after all pSetExp
1757  sum=pAdd(sum,p);
1758  }
1759  }
1760  ret->m[i]=sum;
1761  temp=temp->next;
1762  }
1763  }
1764 // data transferred ****************************************************************************
1765 
1766 
1767  GeneralDone ();
1768  ClearGenList ();
1769  return ret;
1770 }
1771 
1772 
static int final_base_dim
static void CheckColumnSequence()
static void modp_Evaluate(modp_number *ev, mono_type mon, condition_type con)
static int variables
int cf_getSmallPrime(int i)
Definition: cf_primes.cc:28
static int myp_index
#define modp_number
const CanonicalForm int s
Definition: facAbsFact.cc:55
static void GeneralInit()
static modp_number modp_denom
#define pVar(m)
Definition: polys.h:376
static void FreeProcData()
const CanonicalForm int const CFList const Variable & y
Definition: facAbsFact.cc:57
int j
Definition: facHensel.cc:105
static int good_primes
#define pSetm(p)
Definition: polys.h:266
static void UpdateGenList()
static modp_number * in_gamma
void PrintLn()
Definition: reporter.cc:310
#define Print
Definition: emacs.cc:80
static void CloseChinese()
#define pAdd(p, q)
Definition: polys.h:198
static void GeneralDone()
static void ReduceRow()
static modp_result_entry * cur_result
struct modp_result_struct * next
static bool EqualMon(mono_type m1, mono_type m2)
#define pNSet(n)
Definition: polys.h:308
struct mon_list_entry_struct * next
ideal interpolation(const std::vector< ideal > &L, intvec *v)
#define TEST_OPT_PROT
Definition: options.h:102
static modp_result_entry * modp_result
#define pSetExp(p, i, v)
Definition: polys.h:42
static coord_exist_table * coord_exist
Compatiblity layer for legacy polynomial operations (over currRing)
static mono_type MonListElement(mon_list_entry *list, int n)
static void MultGenerators()
#define pLmCmp(p, q)
returns 0|1|-1 if p=q|p>q|p<q w.r.t monomial ordering
Definition: polys.h:105
struct row_list_entry_struct * next
struct gen_list_struct * next
bool * coord_exist_table
static void ReduceCheckListByLTs()
static mpz_t common_denom
static void MakeConditions()
int rChar(ring r)
Definition: ring.cc:713
modp_number * coef
static modp_number * congr
static void ClearGCD()
static int * multiplicity
static bool only_modp
static int n_results
static poly comparizon_p2
static void modp_SetColumnNames()
static FORCE_INLINE void n_Normalize(number &n, const coeffs r)
inplace-normalization of n; produces some canonical representation of n;
Definition: coeffs.h:578
static coordinates * points
bool solve(int **extmat, int nrows, int ncols)
Definition: cf_linsys.cc:504
static mon_list_entry * MonListAdd(mon_list_entry *list, mono_type mon)
static mon_list_entry * generic_lt
static void modp_PrepareProducts()
g
Definition: cfModGcd.cc:4031
void WerrorS(const char *s)
Definition: feFopen.cc:24
int k
Definition: cfEzgcd.cc:92
static number & pGetCoeff(poly p)
return an alias to the leading coefficient of p assumes that p != NULL NOTE: not copy ...
Definition: monomials.h:44
static bool protocol
#define omAlloc(size)
Definition: omAllocDecl.h:210
static int_coordinates * int_points
mono_type * polyexp
static int generic_n_generators
static poly comparizon_p1
static modp_number * modp_Reverse
static exponent MonDegree(mono_type mon)
static modp_number myp
static void NewGenerator(mono_type mon)
if(yy_init)
Definition: libparse.cc:1418
static bool DivisibleMon(mono_type m1, mono_type m2)
modp_number * modp_coordinates
static mpz_t * polycoef
Definition: intvec.h:19
CanonicalForm res
Definition: facAbsFact.cc:64
#define omFree(addr)
Definition: omAllocDecl.h:261
struct modp_result_struct * prev
static void Discard()
#define assume(x)
Definition: mod2.h:390
static gen_list_entry * gen_list
static void FreeResultEntry(modp_result_entry *e)
struct generator_struct * next
modp_number * coordinate_products
#define nInpNeg(n)
Definition: numbers.h:21
for(int i=0;i<=n;i++) degsf[i]
Definition: cfEzgcd.cc:65
static mpz_t bigcongr
static void TakeNextMonomial(mono_type mon)
void * malloc(size_t size)
Definition: omalloc.c:92
modp_number * row_matrix
static mono_type ZeroMonomial()
static int bad_primes
#define exponent
int m
Definition: cfEzgcd.cc:121
unsigned int point_ref
static modp_number TakePrime(modp_number)
int cf_getNumSmallPrimes()
Definition: cf_primes.cc:34
int i
Definition: cfEzgcd.cc:125
static bool denom_divisible
void PrintS(const char *s)
Definition: reporter.cc:284
static BOOLEAN rField_is_Q(const ring r)
Definition: ring.h:501
#define pOne()
Definition: polys.h:310
static bool CheckGenerator()
modp_number ltcoef
#define pHead(p)
returns newly allocated copy of Lm(p), coef is copied, next=NULL, p might be NULL ...
Definition: polys.h:67
#define IDELEMS(i)
Definition: simpleideals.h:23
static modp_number OneInverse(modp_number a, modp_number p)
static void ClearGenList()
int size(const CanonicalForm &f, const Variable &v)
int size ( const CanonicalForm & f, const Variable & v )
Definition: cf_ops.cc:600
static BOOLEAN rField_is_Zp(const ring r)
Definition: ring.h:495
static mon_list_entry * check_list
static bool Greater(mono_type m1, mono_type m2)
static void ProduceRow(mono_type mon)
static condition_type * condition_list
#define SR_TO_INT(SR)
Definition: longrat.h:68
static void RowListAdd(int first_col, mono_type mon)
static mono_type * polyexp
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:35
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:37
static void InitProcData()
static modp_number * my_solve_row
#define nDiv(a, b)
Definition: numbers.h:32
static modp_number modp_sub(modp_number x, modp_number y)
static mon_list_entry * lt_list
CanonicalForm cf
Definition: cfModGcd.cc:4024
static row_list_entry * row_list
static void modp_Main()
static q_coordinates * q_points
exponent * mono_type
static void UpdateCheckList(mono_type m)
#define NULL
Definition: omList.c:12
static void int_Evaluate(mpz_t ev, mono_type mon, condition_type con)
int length() const
Definition: intvec.h:94
static void int_PrepareProducts()
static void ReduceCheckListByMon(mono_type m)
fq_nmod_poly_t prod
Definition: facHensel.cc:95
mpz_t * int_coordinates
static mono_type * column_name
#define SR_INT
Definition: longrat.h:66
mpq_t * q_coordinates
#define ALLOC_RNUMBER()
Definition: coeffs.h:87
#define pDelete(p_ptr)
Definition: polys.h:181
static int max_coord
Variable x
Definition: cfModGcd.cc:4023
static mono_type * generic_column_name
#define pNext(p)
Definition: monomials.h:36
ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:12
static modp_number modp_mul(modp_number x, modp_number y)
static void PrepareChinese(int n)
static int last_solve_column
static void PrepareRow(mono_type mon)
static void IntegerPoints()
int p
Definition: cfModGcd.cc:4019
void wrp(poly p)
Definition: polys.h:305
static void ResolveCoeff(mpq_t c, number m)
static void ReconstructGenerator(int ngen, int n)
#define pISet(i)
Definition: polys.h:307
#define nInit(i)
Definition: numbers.h:24
static modp_number * my_row
modp_number * row_solve
static bool RowIsZero()
generator_entry * generator
static void NewResultEntry()
coordinate_products * coordinates
static int n_points
#define omAlloc0(size)
Definition: omAllocDecl.h:211
static mon_list_entry * base_list
static int CalcBaseDim()
static modp_coordinates * modp_points
return
Definition: cfGcdAlgExt.cc:218
static mon_list_entry * FreeMonList(mon_list_entry *list)