source: git/kernel/linear_algebra/interpolation.cc @ fec0d17

spielwiese
Last change on this file since fec0d17 was fec0d17, checked in by Hans Schoenemann <hannes@…>, 3 years ago
fix: malloc -> omAlloc
  • Property mode set to 100644
File size: 49.2 KB
Line 
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
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
36static inline modp_number modp_mul (modp_number x,modp_number y)
37{
38//    return (x*y)%myp;
39   return (modp_number)(((unsigned long)x)*((unsigned long)y)%((unsigned long)myp));
40}
41static inline modp_number modp_sub (modp_number x,modp_number y)
42{
43//   if (x>=y) return x-y; else return x+myp-y;
44     return (x+myp-y)%myp;
45}
46
47typedef exponent *mono_type;
48typedef struct {mono_type mon; unsigned int point_ref;} condition_type;
49typedef modp_number *coordinate_products;
50typedef coordinate_products *coordinates;
51
52struct mon_list_entry_struct {mono_type mon;
53                              struct mon_list_entry_struct *next;};
54typedef struct mon_list_entry_struct mon_list_entry;
55
56struct row_list_entry_struct {modp_number *row_matrix;
57                              modp_number *row_solve;
58                              int first_col;
59                              struct row_list_entry_struct *next;};
60
61typedef struct row_list_entry_struct row_list_entry;
62
63struct generator_struct {modp_number *coef;
64                         mono_type lt;
65                         modp_number ltcoef;
66                         struct generator_struct *next;};
67
68typedef struct generator_struct generator_entry;
69
70struct modp_result_struct {modp_number p;
71                           generator_entry *generator;
72                           int n_generators;
73                           struct modp_result_struct *next;
74                           struct modp_result_struct *prev;};
75
76typedef struct modp_result_struct modp_result_entry;
77
78typedef modp_number *modp_coordinates;
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
116struct gen_list_struct {mpz_t *polycoef;
117                        mono_type *polyexp;
118                        struct gen_list_struct *next;};
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 comparizons by Singular
132STATIC_VAR poly comparizon_p2;
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;
166     modp_number *solve;
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("?");
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
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;
256    for (i=0;i<n_points;i++) max_coord=max_coord+multiplicity[i];
257    for (j=0;j<n_points;j++)
258    {
259        c=1;
260        for (i=1;i<=variables;i++)
261        {
262            c=(c*(multiplicity[j]+i-1))/i;
263        }
264        s=s+c;
265    }
266    return s;
267}
268
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     }
292     pSetm(comparizon_p1);
293     pSetm(comparizon_p2);
294     bool res=(pLmCmp(comparizon_p1,comparizon_p2)>0);
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;
334     m=(exponent*)omAlloc0(sizeof(exponent)*variables);
335     return m;
336}
337
338static void GeneralInit ()  // general initialization
339{
340     int i,j;
341     points=(coordinates*)omAlloc(sizeof(coordinates)*n_points);
342     for (i=0;i<n_points;i++)
343     {
344         points[i]=(coordinate_products*)omAlloc(sizeof(coordinate_products)*variables);
345         for (j=0;j<variables;j++) points[i][j]=(modp_number*)omAlloc0(sizeof(modp_number)*(max_coord));
346     }
347     condition_list=(condition_type*)omAlloc0(sizeof(condition_type)*final_base_dim);
348     for (i=0;i<final_base_dim;i++) condition_list[i].mon=(exponent*)omAlloc0(sizeof(exponent)*variables);
349     modp_points=(modp_coordinates*)omAlloc(sizeof(modp_coordinates)*n_points);
350     for (i=0;i<n_points;i++) modp_points[i]=(modp_number*)omAlloc0(sizeof(modp_number)*variables);
351     if (!only_modp)
352     {
353        q_points=(q_coordinates*)omAlloc0(sizeof(q_coordinates)*n_points);
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        }
359        int_points=(int_coordinates*)omAlloc0(sizeof(int_coordinates)*n_points);
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     }
366     coord_exist=(coord_exist_table*)omAlloc(sizeof(coord_exist_table)*n_points);
367     for (i=0;i<n_points;i++)
368     {
369         coord_exist[i]=(bool*)omAlloc0(sizeof(bool)*variables);
370     }
371     generic_column_name=(mono_type*)omAlloc(sizeof(mono_type)*final_base_dim);
372     for (i=0;i<final_base_dim;i++) generic_column_name[i]=ZeroMonomial ();
373     good_primes=0;
374     bad_primes=1;
375     generic_n_generators=0;
376     if (!only_modp)
377     {
378        polycoef=(mpz_t*)omAlloc(sizeof(mpz_t)*(final_base_dim+1));
379        polyexp=(mono_type*)omAlloc(sizeof(mono_type)*(final_base_dim+1));
380        for (i=0;i<=final_base_dim;i++)
381        {
382             mpz_init(polycoef[i]);
383             polyexp[i]=ZeroMonomial ();
384        }
385        mpz_init(common_denom);
386     }
387
388// set all globally used lists to be initially empty
389     modp_result=NULL;
390     cur_result=NULL;
391     gen_list=NULL;
392     n_results=0;
393
394// creates polynomials to compare by Singular
395     comparizon_p1=pOne();
396     comparizon_p2=pOne();
397}
398
399static void InitProcData ()  // initialization for procedures which do computations mod p
400{
401     int i;
402     check_list=MonListAdd (check_list,ZeroMonomial ());
403     my_row=(modp_number*)omAlloc0(sizeof(modp_number)*final_base_dim);
404     my_solve_row=(modp_number*)omAlloc0(sizeof(modp_number)*final_base_dim);
405     column_name=(mono_type*)omAlloc(sizeof(mono_type)*final_base_dim);
406     for (i=0;i<final_base_dim;i++) column_name[i]=ZeroMonomial ();
407     last_solve_column=0;
408     modp_number *gen_table;
409     modp_number pos_gen;
410     bool gen_ok;
411     modp_Reverse=(modp_number*)omAlloc0(sizeof(modp_number)*myp);
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);
462     omFree(condition_list);
463     for (i=0;i<n_points;i++) omFree(modp_points[i]);
464     omFree(modp_points);
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        }
478        omFree(int_points);
479        generic_lt=FreeMonList (generic_lt);
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     {
491         omFree(generic_column_name[i]);
492     }
493     omFree(generic_column_name);
494     for (i=0;i<n_points;i++) omFree(coord_exist[i]);
495     omFree(coord_exist);
496     pDelete(&comparizon_p1);
497     pDelete(&comparizon_p2);
498     omFree(multiplicity);
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;
506     check_list=FreeMonList (check_list);
507     lt_list=FreeMonList (lt_list);
508     base_list=FreeMonList (base_list);
509     omFree(my_row);
510     my_row=NULL;
511     omFree(my_solve_row);
512     my_solve_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]);
524     omFree(column_name);
525     omFree(modp_Reverse);
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);
830        omFree(check_list);
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]++;
841         check_list=MonListAdd (check_list,m);
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, shoul 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 debuging, 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--;
1041    return cf_getSmallPrime(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;
1050     modp_number prod;
1051     in_gamma=(modp_number*)omAlloc0(sizeof(modp_number)*n);
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];
1064         in_gamma[i]=OneInverse(prod,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{
1073     omFree(in_gamma);
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 occures
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        }
1199        modp_result=cur_result;
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;
1206        generic_lt=FreeMonList (generic_lt);
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     {
1266         if (!EqualMon(generic_column_name[i],column_name[i]))
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);
1339         omFree(gen_list);
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);
1444              ReduceCheckListByLTs ();
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 );
1513  variables=currRing->N;
1514  only_modp=rField_is_Zp(currRing);
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
1523  final_base_dim = CalcBaseDim ();
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;
1625  myp_index=cf_getNumSmallPrimes ();
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 ();
1639            if (only_modp) modp_PrepareProducts (); else int_PrepareProducts ();
1640
1641            modp_Main ();
1642
1643            if (!only_modp)
1644            {
1645               MultGenerators ();
1646               CheckColumnSequence ();
1647            }
1648            else
1649            {
1650               modp_SetColumnNames ();
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  {
1737    ret=idInit(generic_n_generators,1);
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
Note: See TracBrowser for help on using the repository browser.