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

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