My Project
Loading...
Searching...
No Matches
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
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
33STATIC_VAR modp_number myp; // all modp computation done mod myp
34STATIC_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
48typedef struct {mono_type mon; unsigned int point_ref;} condition_type;
51
54typedef struct mon_list_entry_struct mon_list_entry;
55
60
61typedef struct row_list_entry_struct row_list_entry;
62
67
68typedef struct generator_struct generator_entry;
69
71 generator_entry *generator;
75
76typedef struct modp_result_struct modp_result_entry;
77
79typedef mpq_t *q_coordinates;
80typedef mpz_t *int_coordinates;
81typedef bool *coord_exist_table;
82
83STATIC_VAR int final_base_dim; // dimension of the quotient space, known from the beginning
84STATIC_VAR int last_solve_column; // last non-zero column in "solve" part of matrix, used for speed up
85STATIC_VAR int n_points; // modp_number of ideals (points)
86STATIC_VAR int *multiplicity; // multiplicities of points
87STATIC_VAR int variables; // modp_number of variables
88STATIC_VAR int max_coord; // maximal possible coordinate product used during Evaluation
89STATIC_VAR bool only_modp; // perform only one modp computations
90
91STATIC_VAR modp_coordinates *modp_points; // coordinates of points for modp problem - used by Evaluate (this is also initial data for only modp)
92STATIC_VAR q_coordinates *q_points; // coordinates of points for rational data (not used for modp)
93STATIC_VAR int_coordinates *int_points; // coordinates of points for integer data - used to check generators (not used for modp)
94STATIC_VAR coord_exist_table *coord_exist; // checks whether all coordinates has been initialized
95STATIC_VAR mon_list_entry *check_list; // monomials to be checked in next stages
96STATIC_VAR coordinates *points; // power products of coordinates of points used in modp cycles
97STATIC_VAR condition_type *condition_list; // conditions stored in an array
98STATIC_VAR mon_list_entry *lt_list; // leading terms found so far
99STATIC_VAR mon_list_entry *base_list; // standard monomials found so far
100STATIC_VAR row_list_entry *row_list; // rows of the matrix (both parts)
101STATIC_VAR modp_number *my_row; // one special row to perform operations
102STATIC_VAR modp_number *my_solve_row; // one special row to find the linear dependence ("solve" part)
103STATIC_VAR mono_type *column_name; // monomials assigned to columns in solve_row
104
105STATIC_VAR int n_results; // modp_number of performed modp computations (not discarded)
106STATIC_VAR modp_number modp_denom; // denominator of mod p computations
107STATIC_VAR modp_result_entry *modp_result; // list of results for various mod p calculations (used for modp - first result is the desired one)
108STATIC_VAR modp_result_entry *cur_result; // pointer to current result (as before)
109STATIC_VAR modp_number *congr; // primes used in computations (chinese remainder theorem) (not used for modp)
110STATIC_VAR modp_number *in_gamma; // inverts used in chinese remainder theorem (not used for modp)
111STATIC_VAR mpz_t bigcongr; // result, in fact, is given in mod bigcongr (not used for modp)
112
113STATIC_VAR mpz_t *polycoef; // polynomial integercoefficients (not used for modp)
114STATIC_VAR mono_type *polyexp; // polynomial exponents
115
119typedef struct gen_list_struct gen_list_entry;
120
121STATIC_VAR gen_list_entry *gen_list=NULL; // list of resulting generators - output data (integer version)
122
123STATIC_VAR int generic_n_generators; // modp_number of generators - should be the same for all modp comp (not used for modp)
124STATIC_VAR mono_type *generic_column_name; // monomials assigned to columns in solve_row - should be the same for all modp comp (!!! used for modp)
125STATIC_VAR mon_list_entry *generic_lt=NULL; // leading terms for ordered generators - should be the same for all modp comp (not used for modp)
126STATIC_VAR int good_primes; // modp_number of good primes so far;
127STATIC_VAR int bad_primes; // modp_number of bad primes so far;
128STATIC_VAR mpz_t common_denom; // common denominator used to force points coordinates to Z (not used for modp)
129STATIC_VAR bool denom_divisible; // common denominator is divisible by p (not used for modp)
130
131STATIC_VAR poly comparizon_p1; //polynomials used to do comparisons by Singular
133
134STATIC_VAR modp_number *modp_Reverse; // reverses in mod p
135
136STATIC_VAR bool protocol; // true to show the protocol
137
138#ifdef checksize
139STATIC_VAR int maximal_size=0;
140#endif
141
142#if 0 /* only for debuggig*/
143void 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
150void 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
161void 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
221static 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("?");
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
250static 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
269static 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
277static 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
285static 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
298static 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
323static 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
331static mono_type ZeroMonomial () // allocates memory for new monomial with all enries equal 0
332{
333 mono_type m;
335 return m;
336}
337
338static 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]);
384 }
385 mpz_init(common_denom);
386 }
387
388// set all globally used lists to be initially empty
392 n_results=0;
393
394// creates polynomials to compare by Singular
397}
398
399static void InitProcData () // initialization for procedures which do computations mod p
400{
401 int i;
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
435static 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
449static 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 }
473 for (i=0;i<n_points;i++)
474 {
475 for (j=0;j<variables;j++) mpz_clear(int_points[i][j]);
477 }
480 for (i=0;i<=final_base_dim;i++)
481 {
482 mpz_clear(polycoef[i]);
483 omFree(polyexp[i]);
484 }
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
501static 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 }
523 for (i=0;i<final_base_dim;i++) omFree(column_name[i]);
526}
527
528static 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
551static 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
578static 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
599static 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
627static 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
651static 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
665static 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
696static 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
772static 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
786static 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
794static 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
822static 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
835static 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
846static 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
857static 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
884static 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
937static 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
959static 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
976static 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
1000static 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 */
1016void 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
1038static modp_number TakePrime (modp_number /*p*/) // takes "previous" (smaller) prime
1039{
1040 myp_index--;
1042}
1043
1044static 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;
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
1071static void CloseChinese () // after CRA
1072{
1074 omFree(congr);
1075 mpz_clear(bigcongr);
1076}
1077
1078static 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
1101static 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
1166static 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++;
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
1216static 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
1222static 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 */
1280void 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
1298static 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
1325static 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);
1340 gen_list=temp;
1341 }
1342}
1343
1344static 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 */
1380void 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
1405static 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
1462static 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
1484ideal 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 );
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
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
for(j=0;j< factors.length();j++)
Definition: facHensel.cc:129
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 mon_list_entry * FreeMonList(mon_list_entry *list)
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_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 mon_list_entry * MonListAdd(mon_list_entry *list, mono_type mon)
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_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