My Project
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_VAR modp_number myp; // all modp computation done mod myp
34 STATIC_VAR 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_VAR int final_base_dim; // dimension of the quotient space, known from the beginning
84 STATIC_VAR int last_solve_column; // last non-zero column in "solve" part of matrix, used for speed up
85 STATIC_VAR int n_points; // modp_number of ideals (points)
86 STATIC_VAR int *multiplicity; // multiplicities of points
87 STATIC_VAR int variables; // modp_number of variables
88 STATIC_VAR int max_coord; // maximal possible coordinate product used during Evaluation
89 STATIC_VAR bool only_modp; // perform only one modp computations
90 
91 STATIC_VAR modp_coordinates *modp_points; // coordinates of points for modp problem - used by Evaluate (this is also initial data for only modp)
92 STATIC_VAR q_coordinates *q_points; // coordinates of points for rational data (not used for modp)
93 STATIC_VAR int_coordinates *int_points; // coordinates of points for integer data - used to check generators (not used for modp)
94 STATIC_VAR coord_exist_table *coord_exist; // checks whether all coordinates has been initialized
95 STATIC_VAR mon_list_entry *check_list; // monomials to be checked in next stages
96 STATIC_VAR coordinates *points; // power products of coordinates of points used in modp cycles
97 STATIC_VAR condition_type *condition_list; // conditions stored in an array
98 STATIC_VAR mon_list_entry *lt_list; // leading terms found so far
99 STATIC_VAR mon_list_entry *base_list; // standard monomials found so far
100 STATIC_VAR row_list_entry *row_list; // rows of the matrix (both parts)
101 STATIC_VAR modp_number *my_row; // one special row to perform operations
102 STATIC_VAR modp_number *my_solve_row; // one special row to find the linear dependence ("solve" part)
103 STATIC_VAR mono_type *column_name; // monomials assigned to columns in solve_row
104 
105 STATIC_VAR int n_results; // modp_number of performed modp computations (not discarded)
106 STATIC_VAR modp_number modp_denom; // denominator of mod p computations
107 STATIC_VAR modp_result_entry *modp_result; // list of results for various mod p calculations (used for modp - first result is the desired one)
108 STATIC_VAR modp_result_entry *cur_result; // pointer to current result (as before)
109 STATIC_VAR modp_number *congr; // primes used in computations (chinese remainder theorem) (not used for modp)
110 STATIC_VAR modp_number *in_gamma; // inverts used in chinese remainder theorem (not used for modp)
111 STATIC_VAR mpz_t bigcongr; // result, in fact, is given in mod bigcongr (not used for modp)
112 
113 STATIC_VAR mpz_t *polycoef; // polynomial integercoefficients (not used for modp)
114 STATIC_VAR 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_VAR gen_list_entry *gen_list=NULL; // list of resulting generators - output data (integer version)
122 
123 STATIC_VAR int generic_n_generators; // modp_number of generators - should be the same for all modp comp (not used for modp)
124 STATIC_VAR 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_VAR 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_VAR int good_primes; // modp_number of good primes so far;
127 STATIC_VAR int bad_primes; // modp_number of bad primes so far;
128 STATIC_VAR mpz_t common_denom; // common denominator used to force points coordinates to Z (not used for modp)
129 STATIC_VAR bool denom_divisible; // common denominator is divisible by p (not used for modp)
130 
131 STATIC_VAR poly comparizon_p1; //polynomials used to do comparisons by Singular
133 
134 STATIC_VAR modp_number *modp_Reverse; // reverses in mod p
135 
136 STATIC_VAR bool protocol; // true to show the protocol
137 
138 #ifdef checksize
139 STATIC_VAR 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;
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  }
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  }
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));
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);
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]);
499 }
500 
501 static void FreeProcData () // to be called after one modp computation to free memory
502 {
503  int i;
504  row_list_entry *ptr;
505  row_list_entry *pptr;
509  omFree(my_row);
510  my_row=NULL;
513  ptr=row_list;
514  while (ptr!=NULL)
515  {
516  pptr=ptr->next;
517  omFree(ptr->row_matrix);
518  omFree(ptr->row_solve);
519  omFree(ptr);
520  ptr=pptr;
521  }
522  row_list=NULL;
523  for (i=0;i<final_base_dim;i++) omFree(column_name[i]);
526 }
527 
528 static void modp_Evaluate(modp_number *ev, mono_type mon, condition_type con) // evaluates monomial on condition (modp)
529 {
530  int i;
531  *ev=0;
532  for (i=0;i<variables;i++)
533  if (con.mon[i] > mon[i]) return ;;
534  *ev=1;
535  int j,k;
536  mono_type mn;
537  mn=(exponent*)omAlloc(sizeof(exponent)*variables);
538  memcpy(mn,mon,sizeof(exponent)*variables);
539  for (k=0;k<variables;k++)
540  {
541  for (j=1;j<=con.mon[k];j++) // this loop computes the coefficient from derivation
542  {
543  *ev=modp_mul(*ev,mn[k]);
544  mn[k]--;
545  }
546  *ev=modp_mul(*ev,points[con.point_ref][k][mn[k]]);
547  }
548  omFree(mn);
549 }
550 
551 static void int_Evaluate(mpz_t ev, mono_type mon, condition_type con) // (***) evaluates monomial on condition for integer numbers
552 {
553  int i;
554  mpz_set_ui(ev,0);
555  for (i=0;i<variables;i++)
556  if (con.mon[i] > mon[i]) return ;;
557  mpz_set_ui(ev,1);
558  mpz_t mon_conv;
559  mpz_init(mon_conv);
560  int j,k;
561  mono_type mn;
562  mn=(exponent*)omAlloc(sizeof(exponent)*variables);
563  memcpy(mn,mon,sizeof(exponent)*variables);
564  for (k=0;k<variables;k++)
565  {
566  for (j=1;j<=con.mon[k];j++) // this loop computes the coefficient from derivation
567  {
568  mpz_set_si(mon_conv,mn[k]); // (***)
569  mpz_mul(ev,ev,mon_conv);
570  mn[k]--;
571  }
572  for (j=1;j<=mn[k];j++) mpz_mul(ev,ev,int_points[con.point_ref][k]); // this loop computes the product of coordinate
573  }
574  omFree(mn);
575  mpz_clear(mon_conv);
576 }
577 
578 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
579 {
580  modp_number *row;
581  condition_type *con;
582  int i;
583  row=my_row;
584  con=condition_list;
585  for (i=0;i<final_base_dim;i++)
586  {
587  modp_Evaluate (row, mon, *con);
588  row++;
589  con++;
590  }
591  row=my_solve_row;
592  for (i=0;i<final_base_dim;i++)
593  {
594  *row=0;
595  row++;
596  }
597 }
598 
599 static void IntegerPoints () // produces integer points from rationals by scaling the coordinate system
600 {
601  int i,j;
602  mpz_set_ui(common_denom,1); // this is common scaling factor
603  for (i=0;i<n_points;i++)
604  {
605  for (j=0;j<variables;j++)
606  {
607  mpz_lcm(common_denom,common_denom,mpq_denref(q_points[i][j]));
608  }
609  }
610  mpq_t temp;
611  mpq_init(temp);
612  mpq_t denom_q;
613  mpq_init(denom_q);
614  mpq_set_z(denom_q,common_denom);
615  for (i=0;i<n_points;i++)
616  {
617  for (j=0;j<variables;j++)
618  {
619  mpq_mul(temp,q_points[i][j],denom_q); // multiplies by common denominator
620  mpz_set(int_points[i][j],mpq_numref(temp)); // and changes into integer
621  }
622  }
623  mpq_clear(temp);
624  mpq_clear(denom_q);
625 }
626 
627 static void int_PrepareProducts () // prepares coordinates of points and products for modp case (from integer coefs)
628 {
629  int i,j,k;
630  mpz_t big_myp;
631  mpz_init(big_myp);
632  mpz_set_si(big_myp,myp);
633  mpz_t temp;
634  mpz_init(temp);
635  for (i=0;i<n_points;i++)
636  {
637  for (j=0;j<variables;j++)
638  {
639  mpz_mod(temp,int_points[i][j],big_myp); // coordinate is now modulo myp
640  points[i][j][1]=mpz_get_ui(temp);
641  points[i][j][0]=1;
642  for (k=2;k<max_coord;k++) points[i][j][k]=modp_mul(points[i][j][k-1],points[i][j][1]);
643  }
644  }
645  mpz_mod(temp,common_denom,big_myp); // computes the common denominator (from rational data) modulo myp
646  denom_divisible=(mpz_sgn(temp)==0); // checks whether it is divisible by modp
647  mpz_clear(temp);
648  mpz_clear(big_myp);
649 }
650 
651 static void modp_PrepareProducts () // prepares products for modp computation from modp data
652 {
653  int i,j,k;
654  for (i=0;i<n_points;i++)
655  {
656  for (j=0;j<variables;j++)
657  {
658  points[i][j][1]=modp_points[i][j];
659  points[i][j][0]=1;
660  for (k=2;k<max_coord;k++) points[i][j][k]=modp_mul(points[i][j][k-1],points[i][j][1]);
661  }
662  }
663 }
664 
665 static void MakeConditions () // prepares a list of conditions from list of multiplicities
666 {
667  condition_type *con;
668  int n,i,k;
669  mono_type mon;
670  mon=ZeroMonomial ();
671  con=condition_list;
672  for (n=0;n<n_points;n++)
673  {
674  for (i=0;i<variables;i++) mon[i]=0;
675  while (mon[0]<multiplicity[n])
676  {
677  if (MonDegree (mon) < multiplicity[n])
678  {
679  memcpy(con->mon,mon,sizeof(exponent)*variables);
680  con->point_ref=n;
681  con++;
682  }
683  k=variables-1;
684  mon[k]++;
685  while ((k>0)&&(mon[k]>=multiplicity[n]))
686  {
687  mon[k]=0;
688  k--;
689  mon[k]++;
690  }
691  }
692  }
693  omFree(mon);
694 }
695 
696 static void ReduceRow () // reduces my_row by previous rows, does the same for solve row
697 {
698  if (row_list==NULL) return ;
699  row_list_entry *row_ptr;
700  modp_number *cur_row_ptr;
701  modp_number *solve_row_ptr;
702  modp_number *my_row_ptr;
703  modp_number *my_solve_row_ptr;
704  int first_col;
705  int i;
706  modp_number red_val;
707  modp_number mul_val;
708 #ifdef integerstrategy
709  modp_number *m_row_ptr;
710  modp_number prep_val;
711 #endif
712  row_ptr=row_list;
713  while (row_ptr!=NULL)
714  {
715  cur_row_ptr=row_ptr->row_matrix;
716  solve_row_ptr=row_ptr->row_solve;
717  my_row_ptr=my_row;
718  my_solve_row_ptr=my_solve_row;
719  first_col=row_ptr->first_col;
720  cur_row_ptr=cur_row_ptr+first_col; // reduction begins at first_col position
721  my_row_ptr=my_row_ptr+first_col;
722  red_val=*my_row_ptr;
723  if (red_val!=0)
724  {
725 #ifdef integerstrategy
726  prep_val=*cur_row_ptr;
727  if (prep_val!=1)
728  {
729  m_row_ptr=my_row;
730  for (i=0;i<final_base_dim;i++)
731  {
732  if (*m_row_ptr!=0) *m_row_ptr=modp_mul(*m_row_ptr,prep_val);
733  m_row_ptr++;
734  }
735  m_row_ptr=my_solve_row;
736  for (i=0;i<last_solve_column;i++)
737  {
738  if (*m_row_ptr!=0) *m_row_ptr=modp_mul(*m_row_ptr,prep_val);
739  m_row_ptr++;
740  }
741  }
742 #endif
743  for (i=first_col;i<final_base_dim;i++)
744  {
745  if (*cur_row_ptr!=0)
746  {
747  mul_val=modp_mul(*cur_row_ptr,red_val);
748  *my_row_ptr=modp_sub(*my_row_ptr,mul_val);
749  }
750  my_row_ptr++;
751  cur_row_ptr++;
752  }
753  for (i=0;i<=last_solve_column;i++) // last_solve_column stores the last non-enmpty entry in solve matrix
754  {
755  if (*solve_row_ptr!=0)
756  {
757  mul_val=modp_mul(*solve_row_ptr,red_val);
758  *my_solve_row_ptr=modp_sub(*my_solve_row_ptr,mul_val);
759  }
760  solve_row_ptr++;
761  my_solve_row_ptr++;
762  }
763  }
764  row_ptr=row_ptr->next;
765 #if 0 /* only debugging */
766  PrintS("reduction by row ");
767  Info ();
768 #endif
769  }
770 }
771 
772 static bool RowIsZero () // check whether a row is zero
773 {
774  bool zero=1;
775  int i;
776  modp_number *row;
777  row=my_row;
778  for (i=0;i<final_base_dim;i++)
779  {
780  if (*row!=0) {zero=0; break;}
781  row++;
782  }
783  return zero;
784 }
785 
786 static bool DivisibleMon (mono_type m1, mono_type m2) // checks whether m1 is divisible by m2
787 {
788  int i;
789  for (i=0;i<variables;i++)
790  if (m1[i]>m2[i]) return false;;
791  return true;
792 }
793 
794 static void ReduceCheckListByMon (mono_type m) // from check_list monomials divisible by m are thrown out
795 {
796  mon_list_entry *c_ptr;
797  mon_list_entry *p_ptr;
798  mon_list_entry *n_ptr;
799  c_ptr=check_list;
800  p_ptr=NULL;
801  while (c_ptr!=NULL)
802  {
803  if (DivisibleMon (m,c_ptr->mon))
804  {
805  if (p_ptr==NULL)
806  check_list=c_ptr->next;
807  else
808  p_ptr->next=c_ptr->next;
809  n_ptr=c_ptr->next;
810  omFree(c_ptr->mon);
811  omFree(c_ptr);
812  c_ptr=n_ptr;
813  }
814  else
815  {
816  p_ptr=c_ptr;
817  c_ptr=c_ptr->next;
818  }
819  }
820 }
821 
822 static void TakeNextMonomial (mono_type mon) // reads first monomial from check_list, then it is deleted
823 {
824  mon_list_entry *n_check_list;
825  if (check_list!=NULL)
826  {
827  memcpy(mon,check_list->mon,sizeof(exponent)*variables);
828  n_check_list=check_list->next;
829  omFree(check_list->mon);
831  check_list=n_check_list;
832  }
833 }
834 
835 static void UpdateCheckList (mono_type m) // adds all monomials which are "next" to m (divisible by m and degree +1)
836 {
837  int i;
838  for (i=0;i<variables;i++)
839  {
840  m[i]++;
842  m[i]--;
843  }
844 }
845 
846 static void ReduceCheckListByLTs () // deletes all monomials from check_list which are divisible by one of the leading terms
847 {
848  mon_list_entry *ptr;
849  ptr=lt_list;
850  while (ptr!=NULL)
851  {
852  ReduceCheckListByMon (ptr->mon);
853  ptr=ptr->next;
854  }
855 }
856 
857 static void RowListAdd (int first_col, mono_type mon) // puts a row into matrix
858 {
859  row_list_entry *ptr;
860  row_list_entry *pptr;
861  row_list_entry *temp;
862  ptr=row_list;
863  pptr=NULL;
864  while (ptr!=NULL)
865  {
866 #ifndef unsortedmatrix
867  if ( first_col <= ptr->first_col ) break;
868 #endif
869  pptr=ptr;
870  ptr=ptr->next;
871  }
872  temp=(row_list_entry*)omAlloc0(sizeof(row_list_entry));
873  (*temp).row_matrix=(modp_number*)omAlloc0(sizeof(modp_number)*final_base_dim);
874  memcpy((*temp).row_matrix,my_row,sizeof(modp_number)*(final_base_dim));
875  (*temp).row_solve=(modp_number*)omAlloc0(sizeof(modp_number)*final_base_dim);
876  memcpy((*temp).row_solve,my_solve_row,sizeof(modp_number)*(final_base_dim));
877  (*temp).first_col=first_col;
878  temp->next=ptr;
879  if (pptr==NULL) row_list=temp; else pptr->next=temp;
880  memcpy(column_name[first_col],mon,sizeof(exponent)*variables);
881 }
882 
883 
884 static void PrepareRow (mono_type mon) // prepares a row and puts it into matrix
885 {
886  modp_number *row;
887  int first_col=-1;
888  int col;
889  modp_number red_val=1;
890  row=my_row;
891 #ifdef integerstrategy
892  for (col=0;col<final_base_dim;col++)
893  {
894  if (*row!=0)
895  {
896  first_col=col;
897  break;
898  }
899  row++;
900  }
901  my_solve_row[first_col]=1;
902  if (first_col > last_solve_column) last_solve_column=first_col;
903 #else
904  for (col=0;col<final_base_dim;col++)
905  {
906  if (*row!=0)
907  {
908  first_col=col;
909  red_val=modp_Reverse[*row]; // first non-zero entry, should multiply all row by inverse to obtain 1
910  modp_denom=modp_mul(modp_denom,*row); // remembers the divisor
911  *row=1;
912  break;
913  }
914  row++;
915  }
916  my_solve_row[first_col]=1;
917  if (first_col > last_solve_column) last_solve_column=first_col;
918  if (red_val!=1)
919  {
920  row++;
921  for (col=first_col+1;col<final_base_dim;col++)
922  {
923  if (*row!=0) *row=modp_mul(*row,red_val);
924  row++;
925  }
926  row=my_solve_row;
927  for (col=0;col<=last_solve_column;col++)
928  {
929  if (*row!=0) *row=modp_mul(*row,red_val);
930  row++;
931  }
932  }
933 #endif
934  RowListAdd (first_col, mon);
935 }
936 
937 static void NewResultEntry () // creates an entry for new modp result
938 {
939  modp_result_entry *temp;
940  temp=(modp_result_entry*)omAlloc0(sizeof(modp_result_entry));
941  if (cur_result==NULL)
942  {
943  modp_result=temp;
944  temp->prev=NULL;
945  }
946  else
947  {
948  temp->prev=cur_result;
949  cur_result->next=temp;
950  }
951  cur_result=temp;
952  cur_result->next=NULL;
953  cur_result->p=myp;
954  cur_result->generator=NULL;
955  cur_result->n_generators=0;
956  n_results++;
957 }
958 
959 static void FreeResultEntry (modp_result_entry *e) // destroys the result entry, without worrying about where it is
960 {
961  generator_entry *cur_gen;
962  generator_entry *next_gen;
963  cur_gen=e->generator;
964  while (cur_gen!=NULL)
965  {
966  next_gen=cur_gen->next;
967  omFree(cur_gen->coef);
968  omFree(cur_gen->lt);
969  omFree(cur_gen);
970  cur_gen=next_gen;
971  }
972  omFree(e);
973 }
974 
975 
976 static void NewGenerator (mono_type mon) // new generator in modp comp found, should be stored on the list
977 {
978  generator_entry *cur_ptr;
979  generator_entry *prev_ptr;
980  generator_entry *temp;
981  cur_ptr=cur_result->generator;
982  prev_ptr=NULL;
983  while (cur_ptr!=NULL)
984  {
985  prev_ptr=cur_ptr;
986  cur_ptr=cur_ptr->next;
987  }
988  temp=(generator_entry*)omAlloc0(sizeof(generator_entry));
989  if (prev_ptr==NULL) cur_result->generator=temp;
990  else prev_ptr->next=temp;
991  temp->next=NULL;
992  temp->coef=(modp_number*)omAlloc0(sizeof(modp_number)*final_base_dim);
993  memcpy(temp->coef,my_solve_row,sizeof(modp_number)*final_base_dim);
994  temp->lt=ZeroMonomial ();
995  memcpy(temp->lt,mon,sizeof(exponent)*variables);
996  temp->ltcoef=1;
997  cur_result->n_generators++;
998 }
999 
1000 static void MultGenerators () // before reconstructing, all denominators must be cancelled
1001 {
1002 #ifndef integerstrategy
1003  int i;
1004  generator_entry *cur_ptr;
1005  cur_ptr=cur_result->generator;
1006  while (cur_ptr!=NULL)
1007  {
1008  for (i=0;i<final_base_dim;i++)
1009  cur_ptr->coef[i]=modp_mul(cur_ptr->coef[i],modp_denom);
1010  cur_ptr->ltcoef=modp_denom;
1011  cur_ptr=cur_ptr->next;
1012  }
1013 #endif
1014 }
1015 #if 0 /* only debbuging */
1016 void PresentGenerator (int i) // only for debugging, writes a generator in its form in program
1017 {
1018  int j;
1019  modp_result_entry *cur_ptr;
1020  generator_entry *cur_gen;
1021  cur_ptr=modp_result;
1022  while (cur_ptr!=NULL)
1023  {
1024  cur_gen=cur_ptr->generator;
1025  for (j=0;j<i;j++) cur_gen=cur_gen->next;
1026  for (j=0;j<final_base_dim;j++)
1027  {
1028  Print("%d;", cur_gen->coef[j]);
1029  }
1030  PrintS(" and LT = ");
1031  WriteMono (cur_gen->lt);
1032  Print(" ( %d ) prime = %d\n", cur_gen->ltcoef, cur_ptr->p);
1033  cur_ptr=cur_ptr->next;
1034  }
1035 }
1036 #endif
1037 
1038 static modp_number TakePrime (modp_number /*p*/) // takes "previous" (smaller) prime
1039 {
1040  myp_index--;
1041  return cf_getSmallPrime(myp_index);
1042 }
1043 
1044 static void PrepareChinese (int n) // initialization for CRA
1045 {
1046  int i,k;
1047  modp_result_entry *cur_ptr;
1048  cur_ptr=modp_result;
1049  modp_number *congr_ptr;
1050  modp_number prod;
1052  congr=(modp_number*)omAlloc0(sizeof(modp_number)*n);
1053  congr_ptr=congr;
1054  while (cur_ptr!=NULL)
1055  {
1056  *congr_ptr=cur_ptr->p;
1057  cur_ptr=cur_ptr->next;
1058  congr_ptr++;
1059  }
1060  for (k=1;k<n;k++)
1061  {
1062  prod=congr[0]%congr[k];
1063  for (i=1;i<=k-1;i++) prod=(prod*congr[i])%congr[k];
1065  }
1066  mpz_init(bigcongr);
1067  mpz_set_si(bigcongr,congr[0]);
1068  for (k=1;k<n;k++) mpz_mul_ui(bigcongr,bigcongr,congr[k]);
1069 }
1070 
1071 static void CloseChinese () // after CRA
1072 {
1073  omFree(in_gamma);
1074  omFree(congr);
1075  mpz_clear(bigcongr);
1076 }
1077 
1078 static void ClearGCD () // divides polynomials in basis by gcd of coefficients
1079 {
1080  bool first_gcd=true;
1081  int i;
1082  mpz_t g;
1083  mpz_init(g);
1084  for (i=0;i<=final_base_dim;i++)
1085  {
1086  if (mpz_sgn(polycoef[i])!=0)
1087  {
1088  if (first_gcd)
1089  {
1090  first_gcd=false;
1091  mpz_set(g,polycoef[i]);
1092  }
1093  else
1094  mpz_gcd(g,g,polycoef[i]);
1095  }
1096  }
1097  for (i=0;i<=final_base_dim;i++) mpz_divexact(polycoef[i],polycoef[i],g);
1098  mpz_clear(g);
1099 }
1100 
1101 static void ReconstructGenerator (int ngen,int n) // recostruction of generator from various modp comp
1102 {
1103  int i,j,k;
1104  int coef;
1105  char *str=NULL;
1106  str=(char*)omAlloc0(sizeof(char)*1000);
1107  modp_result_entry *cur_ptr;
1108  generator_entry *cur_gen;
1109  modp_number *u;
1110  modp_number *v;
1111  modp_number temp;
1112  mpz_t sol;
1113  mpz_t nsol;
1114  mpz_init(sol);
1115  mpz_init(nsol);
1116  u=(modp_number*)omAlloc0(sizeof(modp_number)*n);
1117  v=(modp_number*)omAlloc0(sizeof(modp_number)*n);
1118  for (coef=0;coef<=final_base_dim;coef++)
1119  {
1120  i=0;
1121  cur_ptr=modp_result;
1122  while (cur_ptr!=NULL)
1123  {
1124  cur_gen=cur_ptr->generator;
1125  for (j=0;j<ngen;j++) cur_gen=cur_gen->next; // we have to take jth generator
1126  if (coef<final_base_dim) u[i]=cur_gen->coef[coef]; else u[i]=cur_gen->ltcoef;
1127  cur_ptr=cur_ptr->next;
1128  i++;
1129  }
1130  v[0]=u[0]; // now CRA begins
1131  for (k=1;k<n;k++)
1132  {
1133  temp=v[k-1];
1134  for (j=k-2;j>=0;j--) temp=(temp*congr[j]+v[j])%congr[k];
1135  temp=u[k]-temp;
1136  if (temp<0) temp=temp+congr[k];
1137  v[k]=(temp*in_gamma[k])%congr[k];
1138  }
1139  mpz_set_si(sol,v[n-1]);
1140  for (k=n-2;k>=0;k--)
1141  {
1142  mpz_mul_ui(sol,sol,congr[k]);
1143  mpz_add_ui(sol,sol,v[k]);
1144  } // now CRA ends
1145  mpz_sub(nsol,sol,bigcongr);
1146  int s=mpz_cmpabs(sol,nsol);
1147  if (s>0) mpz_set(sol,nsol); // chooses representation closer to 0
1148  mpz_set(polycoef[coef],sol);
1149  if (coef<final_base_dim)
1150  memcpy(polyexp[coef],generic_column_name[coef],sizeof(exponent)*variables);
1151  else
1152  memcpy(polyexp[coef],MonListElement (generic_lt,ngen),sizeof(exponent)*variables);
1153 #ifdef checksize
1154  size=mpz_sizeinbase(sol,10);
1155  if (size>maximal_size) maximal_size=size;
1156 #endif
1157  }
1158  omFree(u);
1159  omFree(v);
1160  omFree(str);
1161  ClearGCD ();
1162  mpz_clear(sol);
1163  mpz_clear(nsol);
1164 }
1165 
1166 static void Discard () // some unlucky prime occurs
1167 {
1168  modp_result_entry *temp;
1169 #ifdef writemsg
1170  Print(", prime=%d", cur_result->p);
1171 #endif
1172  bad_primes++;
1173  if (good_primes>bad_primes)
1174  {
1175 #ifdef writemsg
1176  Print("-discarding this comp (%dth)\n", n_results);
1177 #endif
1178  temp=cur_result;
1179  cur_result=cur_result->prev;
1180  cur_result->next=NULL;
1181  n_results--;
1182  FreeResultEntry (temp);
1183  }
1184  else
1185  {
1186 #ifdef writemsg
1187  PrintS("-discarding ALL.\n");
1188 #endif
1189  int i;
1190  modp_result_entry *ntfree;
1191  generator_entry *cur_gen;
1192  temp=cur_result->prev;
1193  while (temp!=NULL)
1194  {
1195  ntfree=temp->prev;
1196  FreeResultEntry (temp);
1197  temp=ntfree;
1198  }
1200  cur_result->prev=NULL;
1201  n_results=1;
1202  good_primes=1;
1203  bad_primes=0;
1204  generic_n_generators=cur_result->n_generators;
1205  cur_gen=cur_result->generator;
1207  for (i=0;i<generic_n_generators;i++)
1208  {
1209  generic_lt=MonListAdd (generic_lt,cur_gen->lt);
1210  cur_gen=cur_gen->next;
1211  }
1212  for (i=0;i<final_base_dim;i++) memcpy(generic_column_name[i],column_name[i],sizeof(exponent)*variables);
1213  }
1214 }
1215 
1216 static void modp_SetColumnNames () // used by modp - sets column names, the old table will be destroyed
1217 {
1218  int i;
1219  for (i=0;i<final_base_dim;i++) memcpy(generic_column_name[i],column_name[i],sizeof(exponent)*variables);
1220 }
1221 
1222 static void CheckColumnSequence () // checks if scheme of computations is as generic one
1223 {
1224  int i;
1225  if (cur_result->n_generators!=generic_n_generators)
1226  {
1227 #ifdef writemsg
1228  PrintS("wrong number of generators occurred");
1229 #else
1230  if (protocol) PrintS("ng");
1231 #endif
1232  Discard ();
1233  return;
1234  }
1235  if (denom_divisible)
1236  {
1237 #ifdef writemsg
1238  PrintS("denom of coef divisible by p");
1239 #else
1240  if (protocol) PrintS("dp");
1241 #endif
1242  Discard ();
1243  return;
1244  }
1245  generator_entry *cur_gen;
1246  mon_list_entry *cur_mon;
1247  cur_gen=cur_result->generator;
1248  cur_mon=generic_lt;
1249  for (i=0;i<generic_n_generators;i++)
1250  {
1251  if (!EqualMon(cur_mon->mon,cur_gen->lt))
1252  {
1253 #ifdef writemsg
1254  PrintS("wrong leading term occurred");
1255 #else
1256  if (protocol) PrintS("lt");
1257 #endif
1258  Discard ();
1259  return;
1260  }
1261  cur_gen=cur_gen->next;
1262  cur_mon=cur_mon->next;
1263  }
1264  for (i=0;i<final_base_dim;i++)
1265  {
1267  {
1268 #ifdef writemsg
1269  PrintS("wrong seq of cols occurred");
1270 #else
1271  if (protocol) PrintS("sc");
1272 #endif
1273  Discard ();
1274  return;
1275  }
1276  }
1277  good_primes++;
1278 }
1279 #if 0 /* only debuggig */
1280 void WriteGenerator () // writes generator (only for debugging)
1281 {
1282  char *str;
1283  str=(char*)omAlloc0(sizeof(char)*1000);
1284  int i;
1285  for (i=0;i<=final_base_dim;i++)
1286  {
1287  str=mpz_get_str(str,10,polycoef[i]);
1288  PrintS(str);
1289  PrintS("*");
1290  WriteMono(polyexp[i]);
1291  PrintS(" ");
1292  }
1293  omFree(str);
1294  PrintLn();
1295 }
1296 #endif
1297 
1298 static bool CheckGenerator () // evaluates generator to check whether it is good
1299 {
1300  mpz_t val,sum;
1301  int con,i;
1302  mpz_init(val);
1303  mpz_init(sum);
1304  for (con=0;con<final_base_dim;con++)
1305  {
1306  mpz_set_ui(sum,0);
1307  for (i=0;i<=final_base_dim;i++)
1308  {
1309  int_Evaluate(val, polyexp[i], condition_list[con]);
1310  mpz_mul(val,val,polycoef[i]);
1311  mpz_add(sum,sum,val);
1312  }
1313  if (mpz_sgn(sum)!=0)
1314  {
1315  mpz_clear(val);
1316  mpz_clear(sum);
1317  return false;
1318  }
1319  }
1320  mpz_clear(val);
1321  mpz_clear(sum);
1322  return true;
1323 }
1324 
1325 static void ClearGenList ()
1326 {
1327  gen_list_entry *temp;
1328  int i;
1329  while (gen_list!=NULL)
1330  {
1331  temp=gen_list->next;
1332  for (i=0;i<=final_base_dim;i++)
1333  {
1334  mpz_clear(gen_list->polycoef[i]);
1335  omFree(gen_list->polyexp[i]);
1336  }
1337  omFree(gen_list->polycoef);
1338  omFree(gen_list->polyexp);
1339  omFree(gen_list);
1340  gen_list=temp;
1341  }
1342 }
1343 
1344 static void UpdateGenList ()
1345 {
1346  gen_list_entry *temp,*prev;
1347  int i,j;
1348  prev=NULL;
1349  temp=gen_list;
1350  exponent deg;
1351  for (i=0;i<=final_base_dim;i++)
1352  {
1353  deg=MonDegree(polyexp[i]);
1354  for (j=0;j<deg;j++)
1355  {
1356  mpz_mul(polycoef[i],polycoef[i],common_denom);
1357  }
1358  }
1359  ClearGCD ();
1360  while (temp!=NULL)
1361  {
1362  prev=temp;
1363  temp=temp->next;
1364  }
1365  temp=(gen_list_entry*)omAlloc0(sizeof(gen_list_entry));
1366  if (prev==NULL) gen_list=temp; else prev->next=temp;
1367  temp->next=NULL;
1368  temp->polycoef=(mpz_t*)omAlloc(sizeof(mpz_t)*(final_base_dim+1));
1369  temp->polyexp=(mono_type*)omAlloc(sizeof(mono_type)*(final_base_dim+1));
1370  for (i=0;i<=final_base_dim;i++)
1371  {
1372  mpz_init(temp->polycoef[i]);
1373  mpz_set(temp->polycoef[i],polycoef[i]);
1374  temp->polyexp[i]=ZeroMonomial ();
1375  memcpy(temp->polyexp[i],polyexp[i],sizeof(exponent)*variables);
1376  }
1377 }
1378 
1379 #if 0 /* only debugging */
1380 void ShowGenList ()
1381 {
1382  gen_list_entry *temp;
1383  int i;
1384  char *str;
1385  str=(char*)omAlloc0(sizeof(char)*1000);
1386  temp=gen_list;
1387  while (temp!=NULL)
1388  {
1389  PrintS("generator: ");
1390  for (i=0;i<=final_base_dim;i++)
1391  {
1392  str=mpz_get_str(str,10,temp->polycoef[i]);
1393  PrintS(str);
1394  PrintS("*");
1395  WriteMono(temp->polyexp[i]);
1396  }
1397  PrintLn();
1398  temp=temp->next;
1399  }
1400  omFree(str);
1401 }
1402 #endif
1403 
1404 
1405 static void modp_Main ()
1406 {
1407  mono_type cur_mon;
1408  cur_mon= ZeroMonomial ();
1409  modp_denom=1;
1410  bool row_is_zero;
1411 
1412 #if 0 /* only debugging */
1413  Info ();
1414 #endif
1415 
1416  while (check_list!=NULL)
1417  {
1418  TakeNextMonomial (cur_mon);
1419  ProduceRow (cur_mon);
1420 #if 0 /* only debugging */
1421  cout << "row produced for monomial ";
1422  WriteMono (cur_mon);
1423  cout << endl;
1424  Info ();
1425 #endif
1426  ReduceRow ();
1427  row_is_zero = RowIsZero ();
1428  if (row_is_zero)
1429  {
1430  lt_list=MonListAdd (lt_list,cur_mon);
1431  ReduceCheckListByMon (cur_mon);
1432  NewGenerator (cur_mon);
1433 #if 0 /* only debugging */
1434  cout << "row is zero - linear dependence found (should be seen in my_solve_row)" << endl;
1435  cout << "monomial added to leading terms list" << endl;
1436  cout << "check list updated" << endl;
1437  Info ();
1438 #endif
1439  }
1440  else
1441  {
1442  base_list= MonListAdd (base_list,cur_mon);
1443  UpdateCheckList (cur_mon);
1445 #if 0 /* only debugging */
1446  cout << "row is non-zero" << endl;
1447  cout << "monomial added to quotient basis list" << endl;
1448  cout << "new monomials added to check list" << endl;
1449  cout << "check list reduced by monomials from leading term list" << endl;
1450  Info ();
1451 #endif
1452  PrepareRow (cur_mon);
1453 #if 0 /* only debugging */
1454  cout << "row prepared and put into matrix" << endl;
1455  Info ();
1456 #endif
1457  }
1458  }
1459  omFree(cur_mon);
1460 }
1461 
1462 static void ResolveCoeff (mpq_t c, number m)
1463 {
1464  if ((long)m & SR_INT)
1465  {
1466  long m_val=SR_TO_INT(m);
1467  mpq_set_si(c,m_val,1);
1468  }
1469  else
1470  {
1471  if (m->s<2)
1472  {
1473  mpz_set(mpq_numref(c),m->z);
1474  mpz_set(mpq_denref(c),m->n);
1475  mpq_canonicalize(c);
1476  }
1477  else
1478  {
1479  mpq_set_z(c,m->z);
1480  }
1481  }
1482 }
1483 
1484 ideal interpolation(const std::vector<ideal>& L, intvec *v)
1485 {
1486  protocol=TEST_OPT_PROT; // should be set if option(prot) is enabled
1487 
1488  bool data_ok=true;
1489 
1490  // reading the ring data ***************************************************
1491  if ((currRing==NULL) || ((!rField_is_Zp (currRing))&&(!rField_is_Q (currRing))))
1492  {
1493  WerrorS("coefficient field should be Zp or Q!");
1494  return NULL;
1495  }
1496  if ((currRing->qideal)!=NULL)
1497  {
1498  WerrorS("quotient ring not supported!");
1499  return NULL;
1500  }
1501  if ((currRing->OrdSgn)!=1)
1502  {
1503  WerrorS("ordering must be global!");
1504  return NULL;
1505  }
1506  n_points=v->length ();
1507  if (n_points!=L.size())
1508  {
1509  WerrorS("list and intvec must have the same length!");
1510  return NULL;
1511  }
1512  assume( n_points > 0 );
1513  variables=currRing->N;
1515  if (only_modp) myp=rChar(currRing);
1516  // ring data read **********************************************************
1517 
1518 
1519  multiplicity=(int*)omAlloc(sizeof(int)*n_points);
1520  int i;
1521  for (i=0;i<n_points;i++) multiplicity[i]=(*v)[i];
1522 
1524 
1525 #ifdef writemsg
1526  Print("number of variables: %d\n", variables);
1527  Print("number of points: %d\n", n_points);
1528  PrintS("multiplicities: ");
1529  for (i=0;i<n_points;i++) Print("%d ", multiplicity[i]);
1530  PrintLn();
1531  Print("general initialization for dimension %d ...\n", final_base_dim);
1532 #endif
1533 
1534  GeneralInit ();
1535 
1536 // reading coordinates of points from ideals **********************************
1537  mpq_t divisor;
1538  if (!only_modp) mpq_init(divisor);
1539  int j;
1540  for(i=0; i<L.size();i++)
1541  {
1542  ideal I = L[i];
1543  for(j=0;j<IDELEMS(I);j++)
1544  {
1545  poly p=I->m[j];
1546  if (p!=NULL)
1547  {
1548  poly ph=pHead(p);
1549  int pcvar=pVar(ph);
1550  if (pcvar!=0)
1551  {
1552  pcvar--;
1553  if (coord_exist[i][pcvar])
1554  {
1555  Print("coordinate %d for point %d initialized twice!\n",pcvar+1,i+1);
1556  data_ok=false;
1557  }
1558  number m;
1559  m=pGetCoeff(p); // possible coefficient standing by a leading monomial
1560  if (!only_modp) ResolveCoeff (divisor,m);
1561  number n;
1562  if (pNext(p)!=NULL) n=pGetCoeff(pNext(p));
1563  else n=nInit(0);
1564  if (only_modp)
1565  {
1566  n=nInpNeg(n);
1567  n=nDiv(n,m);
1568  modp_points[i][pcvar]=(int)((long)n);
1569  }
1570  else
1571  {
1572  ResolveCoeff (q_points[i][pcvar],n);
1573  mpq_neg(q_points[i][pcvar],q_points[i][pcvar]);
1574  mpq_div(q_points[i][pcvar],q_points[i][pcvar],divisor);
1575  }
1576  coord_exist[i][pcvar]=true;
1577  }
1578  else
1579  {
1580  PrintS("not a variable? ");
1581  wrp(p);
1582  PrintLn();
1583  data_ok=false;
1584  }
1585  pDelete(&ph);
1586  }
1587  }
1588  }
1589  if (!only_modp) mpq_clear(divisor);
1590  // data from ideal read *******************************************************
1591 
1592  // ckecking if all coordinates are initialized
1593  for (i=0;i<n_points;i++)
1594  {
1595  for (j=0;j<variables;j++)
1596  {
1597  if (!coord_exist[i][j])
1598  {
1599  Print("coordinate %d for point %d not known!\n",j+1,i+1);
1600  data_ok=false;
1601  }
1602  }
1603  }
1604 
1605  if (!data_ok)
1606  {
1607  GeneralDone();
1608  WerrorS("data structure is invalid");
1609  return NULL;
1610  }
1611 
1612  if (!only_modp) IntegerPoints ();
1613  MakeConditions ();
1614 #ifdef writemsg
1615  PrintS("done.\n");
1616 #else
1617  if (protocol) Print("[vdim %d]",final_base_dim);
1618 #endif
1619 
1620 
1621 // main procedure *********************************************************************
1622  int modp_cycles=10;
1623  bool correct_gen=false;
1624  if (only_modp) modp_cycles=1;
1626 
1627  while ((!correct_gen)&&(myp_index>1))
1628  {
1629 #ifdef writemsg
1630  Print("trying %d cycles mod p...\n",modp_cycles);
1631 #else
1632  if (protocol) Print("(%d)",modp_cycles);
1633 #endif
1634  while ((n_results<modp_cycles)&&(myp_index>1)) // some computations mod p
1635  {
1636  if (!only_modp) myp=TakePrime (myp);
1637  NewResultEntry ();
1638  InitProcData ();
1640 
1641  modp_Main ();
1642 
1643  if (!only_modp)
1644  {
1645  MultGenerators ();
1647  }
1648  else
1649  {
1651  }
1652  FreeProcData ();
1653  }
1654 
1655  if (!only_modp)
1656  {
1657  PrepareChinese (modp_cycles);
1658  correct_gen=true;
1659  for (i=0;i<generic_n_generators;i++)
1660  {
1661  ReconstructGenerator (i,modp_cycles);
1662  correct_gen=CheckGenerator ();
1663  if (!correct_gen)
1664  {
1665 #ifdef writemsg
1666  PrintS("wrong generator!\n");
1667 #else
1668 // if (protocol) PrintS("!g");
1669 #endif
1670  ClearGenList ();
1671  break;
1672  }
1673  else
1674  {
1675  UpdateGenList ();
1676  }
1677  }
1678 #ifdef checksize
1679  Print("maximal size of output: %d, precision bound: %d.\n",maximalsize,mpz_sizeinbase(bigcongr,10));
1680 #endif
1681  CloseChinese ();
1682  modp_cycles=modp_cycles+10;
1683  }
1684  else
1685  {
1686  correct_gen=true;
1687  }
1688  }
1689 // end of main procedure ************************************************************************************
1690 
1691 #ifdef writemsg
1692  PrintS("computations finished.\n");
1693 #else
1694  if (protocol) PrintLn();
1695 #endif
1696 
1697  if (!correct_gen)
1698  {
1699  GeneralDone ();
1700  ClearGenList ();
1701  WerrorS("internal error - coefficient too big!");
1702  return NULL;
1703  }
1704 
1705 // passing data to ideal *************************************************************************************
1706  ideal ret;
1707 
1708  if (only_modp)
1709  {
1710  mono_type mon;
1711  ret=idInit(modp_result->n_generators,1);
1712  generator_entry *cur_gen=modp_result->generator;
1713  for(i=0;i<IDELEMS(ret);i++)
1714  {
1715  poly p,sum;
1716  sum=NULL;
1717  int a;
1718  int cf;
1719  for (a=final_base_dim;a>=0;a--)
1720  {
1721  if (a==final_base_dim) cf=cur_gen->ltcoef; else cf=cur_gen->coef[a];
1722  if (cf!=0)
1723  {
1724  p=pISet(cf);
1725  if (a==final_base_dim) mon=cur_gen->lt; else mon=generic_column_name[a];
1726  for (j=0;j<variables;j++) pSetExp(p,j+1,mon[j]);
1727  pSetm(p);
1728  sum=pAdd(sum,p);
1729  }
1730  }
1731  ret->m[i]=sum;
1732  cur_gen=cur_gen->next;
1733  }
1734  }
1735  else
1736  {
1738  gen_list_entry *temp=gen_list;
1739  for(i=0;i<IDELEMS(ret);i++)
1740  {
1741  poly p,sum;
1742  sum=NULL;
1743  int a;
1744  for (a=final_base_dim;a>=0;a--) // build one polynomial
1745  {
1746  if (mpz_sgn(temp->polycoef[a])!=0)
1747  {
1748  number n=ALLOC_RNUMBER();
1749 #ifdef LDEBUG
1750  n->debug=123456;
1751 #endif
1752  mpz_init_set(n->z,temp->polycoef[a]);
1753  n->s=3;
1754  n_Normalize(n, currRing->cf);
1755  p=pNSet(n); //a monomial
1756  for (j=0;j<variables;j++) pSetExp(p,j+1,temp->polyexp[a][j]);
1757  pSetm(p); // after all pSetExp
1758  sum=pAdd(sum,p);
1759  }
1760  }
1761  ret->m[i]=sum;
1762  temp=temp->next;
1763  }
1764  }
1765 // data transferred ****************************************************************************
1766 
1767 
1768  GeneralDone ();
1769  ClearGenList ();
1770  return ret;
1771 }
1772 
1773 
int size(const CanonicalForm &f, const Variable &v)
int size ( const CanonicalForm & f, const Variable & v )
Definition: cf_ops.cc:600
int m
Definition: cfEzgcd.cc:128
for(int i=0;i<=n;i++) degsf[i]
Definition: cfEzgcd.cc:72
int i
Definition: cfEzgcd.cc:132
int k
Definition: cfEzgcd.cc:99
return
Definition: cfGcdAlgExt.cc:218
Variable x
Definition: cfModGcd.cc:4082
int p
Definition: cfModGcd.cc:4078
g
Definition: cfModGcd.cc:4090
CanonicalForm cf
Definition: cfModGcd.cc:4083
bool solve(int **extmat, int nrows, int ncols)
Definition: cf_linsys.cc:504
int cf_getNumSmallPrimes()
Definition: cf_primes.cc:34
int cf_getSmallPrime(int i)
Definition: cf_primes.cc:28
Definition: intvec.h:23
#define ALLOC_RNUMBER()
Definition: coeffs.h:87
static FORCE_INLINE void n_Normalize(number &n, const coeffs r)
inplace-normalization of n; produces some canonical representation of n;
Definition: coeffs.h:575
#define Print
Definition: emacs.cc:80
const CanonicalForm int s
Definition: facAbsFact.cc:51
const CanonicalForm int const CFList const Variable & y
Definition: facAbsFact.cc:53
CanonicalForm res
Definition: facAbsFact.cc:60
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:39
int j
Definition: facHensel.cc:110
fq_nmod_poly_t prod
Definition: facHensel.cc:100
void WerrorS(const char *s)
Definition: feFopen.cc:24
if(!FE_OPT_NO_SHELL_FLAG)(void) system(sys)
#define STATIC_VAR
Definition: globaldefs.h:7
static void NewGenerator(mono_type mon)
static void TakeNextMonomial(mono_type mon)
static modp_number TakePrime(modp_number)
STATIC_VAR int * multiplicity
static void ClearGenList()
static int CalcBaseDim()
static void Discard()
STATIC_VAR modp_number * my_row
static modp_number modp_sub(modp_number x, modp_number y)
STATIC_VAR int myp_index
static void ReduceCheckListByMon(mono_type m)
STATIC_VAR int last_solve_column
static void modp_Evaluate(modp_number *ev, mono_type mon, condition_type con)
unsigned int point_ref
static mono_type ZeroMonomial()
static void UpdateGenList()
static mono_type MonListElement(mon_list_entry *list, int n)
static void FreeProcData()
STATIC_VAR mpz_t bigcongr
static void ResolveCoeff(mpq_t c, number m)
STATIC_VAR mon_list_entry * check_list
static modp_number modp_mul(modp_number x, modp_number y)
STATIC_VAR mon_list_entry * generic_lt
ideal interpolation(const std::vector< ideal > &L, intvec *v)
#define exponent
static void IntegerPoints()
static bool EqualMon(mono_type m1, mono_type m2)
STATIC_VAR coordinates * points
struct modp_result_struct * prev
static void MakeConditions()
static void CloseChinese()
static void ClearGCD()
static void int_Evaluate(mpz_t ev, mono_type mon, condition_type con)
struct modp_result_struct * next
coordinate_products * coordinates
static void NewResultEntry()
static void PrepareChinese(int n)
STATIC_VAR modp_number myp
static void ReduceCheckListByLTs()
static mon_list_entry * MonListAdd(mon_list_entry *list, mono_type mon)
STATIC_VAR q_coordinates * q_points
static modp_number OneInverse(modp_number a, modp_number p)
static void ReconstructGenerator(int ngen, int n)
STATIC_VAR modp_result_entry * modp_result
STATIC_VAR int bad_primes
static void modp_PrepareProducts()
STATIC_VAR int n_points
static void PrepareRow(mono_type mon)
mpq_t * q_coordinates
modp_number * coef
STATIC_VAR modp_coordinates * modp_points
STATIC_VAR bool only_modp
static bool CheckGenerator()
#define modp_number
STATIC_VAR modp_number * my_solve_row
STATIC_VAR mono_type * polyexp
STATIC_VAR condition_type * condition_list
STATIC_VAR poly comparizon_p1
static exponent MonDegree(mono_type mon)
static void InitProcData()
static mon_list_entry * FreeMonList(mon_list_entry *list)
STATIC_VAR modp_number * congr
STATIC_VAR modp_number modp_denom
STATIC_VAR mpz_t common_denom
struct mon_list_entry_struct * next
STATIC_VAR modp_number * in_gamma
STATIC_VAR int good_primes
STATIC_VAR mon_list_entry * base_list
STATIC_VAR mon_list_entry * lt_list
exponent * mono_type
static void ProduceRow(mono_type mon)
struct row_list_entry_struct * next
STATIC_VAR int final_base_dim
static void MultGenerators()
STATIC_VAR mono_type * generic_column_name
STATIC_VAR modp_number * modp_Reverse
struct generator_struct * next
static void UpdateCheckList(mono_type m)
static void GeneralInit()
STATIC_VAR mpz_t * polycoef
modp_number * modp_coordinates
bool * coord_exist_table
STATIC_VAR mono_type * column_name
STATIC_VAR bool denom_divisible
modp_number * row_matrix
STATIC_VAR int_coordinates * int_points
static bool RowIsZero()
STATIC_VAR gen_list_entry * gen_list
static void FreeResultEntry(modp_result_entry *e)
STATIC_VAR int variables
static void modp_Main()
mono_type * polyexp
STATIC_VAR int max_coord
static void int_PrepareProducts()
mpz_t * int_coordinates
struct gen_list_struct * next
modp_number ltcoef
modp_number * row_solve
static void GeneralDone()
static void CheckColumnSequence()
static bool DivisibleMon(mono_type m1, mono_type m2)
static void ReduceRow()
generator_entry * generator
STATIC_VAR int generic_n_generators
STATIC_VAR bool protocol
STATIC_VAR int n_results
static void RowListAdd(int first_col, mono_type mon)
STATIC_VAR row_list_entry * row_list
modp_number * coordinate_products
STATIC_VAR coord_exist_table * coord_exist
static void modp_SetColumnNames()
STATIC_VAR poly comparizon_p2
STATIC_VAR modp_result_entry * cur_result
static bool Greater(mono_type m1, mono_type m2)
#define SR_INT
Definition: longrat.h:67
#define SR_TO_INT(SR)
Definition: longrat.h:69
#define assume(x)
Definition: mod2.h:389
#define pNext(p)
Definition: monomials.h:36
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
char * str(leftv arg)
Definition: shared.cc:704
#define nDiv(a, b)
Definition: numbers.h:32
#define nInpNeg(n)
Definition: numbers.h:21
#define nInit(i)
Definition: numbers.h:24
#define omAlloc(size)
Definition: omAllocDecl.h:210
#define omFree(addr)
Definition: omAllocDecl.h:261
#define omAlloc0(size)
Definition: omAllocDecl.h:211
#define NULL
Definition: omList.c:12
#define TEST_OPT_PROT
Definition: options.h:104
VAR ring currRing
Widely used global variable which specifies the current polynomial ring for Singular interpreter and ...
Definition: polys.cc:13
Compatibility layer for legacy polynomial operations (over currRing)
#define pAdd(p, q)
Definition: polys.h:203
#define pDelete(p_ptr)
Definition: polys.h:186
#define pHead(p)
returns newly allocated copy of Lm(p), coef is copied, next=NULL, p might be NULL
Definition: polys.h:67
#define pSetm(p)
Definition: polys.h:271
#define pNSet(n)
Definition: polys.h:313
#define pVar(m)
Definition: polys.h:380
void wrp(poly p)
Definition: polys.h:310
#define pSetExp(p, i, v)
Definition: polys.h:42
#define pLmCmp(p, q)
returns 0|1|-1 if p=q|p>q|p<q w.r.t monomial ordering
Definition: polys.h:105
#define pOne()
Definition: polys.h:315
#define pISet(i)
Definition: polys.h:312
void PrintS(const char *s)
Definition: reporter.cc:284
void PrintLn()
Definition: reporter.cc:310
int rChar(ring r)
Definition: ring.cc:713
static BOOLEAN rField_is_Zp(const ring r)
Definition: ring.h:500
static BOOLEAN rField_is_Q(const ring r)
Definition: ring.h:506
ideal idInit(int idsize, int rank)
initialise an ideal / module
Definition: simpleideals.cc:35
#define IDELEMS(i)
Definition: simpleideals.h:23