source: git/Singular/interpolation.cc @ 0519cc

spielwiese
Last change on this file since 0519cc was 0519cc, checked in by Martin Lee <martinlee84@…>, 13 years ago
fixed interpolation.cc
  • Property mode set to 100644
File size: 48.9 KB
Line 
1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4/* $Id$ */
5
6#include <kernel/mod2.h>
7#include <misc/options.h>
8#include <kernel/febase.h>
9#include <kernel/ideals.h>
10#include <misc/intvec.h>
11#include <polys/polys.h>
12#include <Singular/lists.h>
13#include <coeffs/longrat.h>
14#include <Singular/ipid.h>
15#include <polys/monomials/ring.h>
16#ifdef HAVE_FACTORY
17#define SI_DONT_HAVE_GLOBAL_VARS
18#  include <factory/factory.h>
19#endif /* HAVE_FACTORY */
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
33modp_number myp;  // all modp computation done mod myp
34int 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
83int final_base_dim;    // dimension of the quotient space, known from the beginning
84int last_solve_column;  // last non-zero column in "solve" part of matrix, used for speed up
85int n_points;  // modp_number of ideals (points)
86int *multiplicity;  // multiplicities of points
87int variables;  // modp_number of variables
88int max_coord;  // maximal possible coordinate product used during Evaluation
89bool only_modp;  // perform only one modp computations
90
91modp_coordinates *modp_points; // coordinates of points for modp problem - used by Evaluate (this is also initial data for only modp)
92q_coordinates *q_points; // coordinates of points for rational data (not used for modp)
93int_coordinates *int_points; // coordinates of points for integer data - used to check generators (not used for modp)
94coord_exist_table *coord_exist; // checks whether all coordinates has been initialized
95mon_list_entry *check_list; // monomials to be checked in next stages
96coordinates *points; // power products of coordinates of points used in modp cycles
97condition_type *condition_list; // conditions stored in an array
98mon_list_entry *lt_list; // leading terms found so far
99mon_list_entry *base_list; // standard monomials found so far
100row_list_entry *row_list; // rows of the matrix (both parts)
101modp_number *my_row; // one special row to perform operations
102modp_number *my_solve_row; // one special row to find the linear dependence ("solve" part)
103mono_type *column_name; // monomials assigned to columns in solve_row
104
105int n_results;  // modp_number of performed modp computations (not discarded)
106modp_number modp_denom; // denominator of mod p computations
107modp_result_entry *modp_result; // list of results for various mod p calculations (used for modp - first result is the desired one)
108modp_result_entry *cur_result; // pointer to current result (as before)
109modp_number *congr; // primes used in computations (chinese remainder theorem) (not used for modp)
110modp_number *in_gamma; // inverts used in chinese remainder theorem (not used for modp)
111mpz_t bigcongr; // result, in fact, is given in mod bigcongr (not used for modp)
112
113mpz_t *polycoef; // polynomial integercoefficients (not used for modp)
114mono_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
121gen_list_entry *gen_list=NULL; // list of resulting generators - output data (integer version)
122
123int generic_n_generators; // modp_number of generators - should be the same for all modp comp (not used for modp)
124mono_type *generic_column_name; // monomials assigned to columns in solve_row - should be the same for all modp comp (!!! used for modp)
125mon_list_entry *generic_lt=NULL; // leading terms for ordered generators - should be the same for all modp comp (not used for modp)
126int good_primes; // modp_number of good primes so far;
127int bad_primes; // modp_number of bad primes so far;
128mpz_t common_denom; // common denominator used to force points coordinates to Z (not used for modp)
129bool denom_divisible; // common denominator is divisible by p (not used for modp)
130
131poly comparizon_p1;  //polynomials used to do comparizons by Singular
132poly comparizon_p2;
133
134modp_number *modp_Reverse; // reverses in mod p
135
136bool protocol; // true to show the protocol
137
138#ifdef checksize
139int 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
221modp_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
250int 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
269bool 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
277exponent 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
285bool 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
298mon_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
323mono_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
331mono_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
338void GeneralInit ()  // general initialization
339{
340     int i,j,k;
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
399void InitProcData ()  // initialization for procedures which do computations mod p
400{
401     int i,j;
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
435mon_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
449void GeneralDone ()  // to be called before exit to free memory
450{
451     int i,j,k;
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}
499
500void FreeProcData ()  // to be called after one modp computation to free memory
501{
502     int i;
503     row_list_entry *ptr;
504     row_list_entry *pptr;
505     check_list=FreeMonList (check_list);
506     lt_list=FreeMonList (lt_list);
507     base_list=FreeMonList (base_list);
508     omFree(my_row);
509     my_row=NULL;
510     omFree(my_solve_row);
511     my_solve_row=NULL;
512     ptr=row_list;
513     while (ptr!=NULL)
514     {
515           pptr=ptr->next;
516           omFree(ptr->row_matrix);
517           omFree(ptr->row_solve);
518           omFree(ptr);
519           ptr=pptr;
520     }
521     row_list=NULL;
522     for (i=0;i<final_base_dim;i++) omFree(column_name[i]);
523     omFree(column_name);
524     omFree(modp_Reverse);
525}
526
527void modp_Evaluate(modp_number *ev, mono_type mon, condition_type con)  // evaluates monomial on condition (modp)
528{
529    int i;
530    *ev=0;
531    for (i=0;i<variables;i++)
532        if (con.mon[i] > mon[i]) return ;;
533    *ev=1;
534    int j,k;
535    mono_type mn;
536    mn=(exponent*)omAlloc(sizeof(exponent)*variables);
537    memcpy(mn,mon,sizeof(exponent)*variables);
538    for (k=0;k<variables;k++)
539    {
540        for (j=1;j<=con.mon[k];j++) // this loop computes the coefficient from derivation
541        {
542            *ev=modp_mul(*ev,mn[k]);
543            mn[k]--;
544        }
545        *ev=modp_mul(*ev,points[con.point_ref][k][mn[k]]);
546    }
547    omFree(mn);
548}
549
550void int_Evaluate(mpz_t ev, mono_type mon, condition_type con) // (***) evaluates monomial on condition for integer numbers
551{
552    int i;
553    mpz_set_si(ev,0);
554    for (i=0;i<variables;i++)
555        if (con.mon[i] > mon[i]) return ;;
556    mpz_set_si(ev,1);
557    mpz_t mon_conv;
558    mpz_init(mon_conv);
559    int j,k;
560    mono_type mn;
561    mn=(exponent*)omAlloc(sizeof(exponent)*variables);
562    memcpy(mn,mon,sizeof(exponent)*variables);
563    for (k=0;k<variables;k++)
564    {
565        for (j=1;j<=con.mon[k];j++) // this loop computes the coefficient from derivation
566        {
567            mpz_set_si(mon_conv,mn[k]); // (***)
568            mpz_mul(ev,ev,mon_conv);
569            mn[k]--;
570        }
571        for (j=1;j<=mn[k];j++) mpz_mul(ev,ev,int_points[con.point_ref][k]);  // this loop computes the product of coordinate
572    }
573    omFree(mn);
574    mpz_clear(mon_conv);
575}
576
577void ProduceRow(mono_type mon)  // produces a row for monomial - first part is an evaluated part, second 0 to obtain the coefs of dependence
578{
579     modp_number *row;
580     condition_type *con;
581     int i;
582     row=my_row;
583     con=condition_list;
584     for (i=0;i<final_base_dim;i++)
585     {
586         modp_Evaluate (row, mon, *con);
587         row++;
588         con++;
589     }
590     row=my_solve_row;
591     for (i=0;i<final_base_dim;i++)
592     {
593         *row=0;
594         row++;
595     }
596}
597
598void IntegerPoints ()  // produces integer points from rationals by scaling the coordinate system
599{
600     int i,j;
601     mpz_set_si(common_denom,1); // this is common scaling factor
602     for (i=0;i<n_points;i++)
603     {
604         for (j=0;j<variables;j++)
605         {
606             mpz_lcm(common_denom,common_denom,mpq_denref(q_points[i][j]));
607         }
608     }
609     mpq_t temp;
610     mpq_init(temp);
611     mpq_t denom_q;
612     mpq_init(denom_q);
613     mpq_set_z(denom_q,common_denom);
614     for (i=0;i<n_points;i++)
615     {
616         for (j=0;j<variables;j++)
617         {
618             mpq_mul(temp,q_points[i][j],denom_q);  // multiplies by common denominator
619             mpz_set(int_points[i][j],mpq_numref(temp));  // and changes into integer
620         }
621     }
622     mpq_clear(temp);
623     mpq_clear(denom_q);
624}
625
626void int_PrepareProducts ()  // prepares coordinates of points and products for modp case (from integer coefs)
627{
628     bool ok=true;
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
651void 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
665void 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
696void 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
772bool 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
786bool 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
794void 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
822void 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
835void 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
846void 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
857void 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
884void 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
937void NewResultEntry ()  // creates an entry for new modp result
938{
939    modp_result_entry *temp;
940    int i;
941    temp=(modp_result_entry*)omAlloc0(sizeof(modp_result_entry));
942    if (cur_result==NULL)
943    {
944       modp_result=temp;
945       temp->prev=NULL;
946    }
947    else
948    {
949       temp->prev=cur_result;
950       cur_result->next=temp;
951    }
952    cur_result=temp;
953    cur_result->next=NULL;
954    cur_result->p=myp;
955    cur_result->generator=NULL;
956    cur_result->n_generators=0;
957    n_results++;
958}
959
960void FreeResultEntry (modp_result_entry *e) // destroys the result entry, without worrying about where it is
961{
962     generator_entry *cur_gen;
963     generator_entry *next_gen;
964     cur_gen=e->generator;
965     while (cur_gen!=NULL)
966     {
967         next_gen=cur_gen->next;
968         omFree(cur_gen->coef);
969         omFree(cur_gen->lt);
970         omFree(cur_gen);
971         cur_gen=next_gen;
972     }
973     omFree(e);
974}
975
976
977void NewGenerator (mono_type mon)  // new generator in modp comp found, shoul be stored on the list
978{
979     int i;
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
1002void 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
1040modp_number TakePrime (modp_number p)  // takes "previous" (smaller) prime
1041{
1042#ifdef HAVE_FACTORY
1043    myp_index--;
1044    return cf_getSmallPrime(myp_index);
1045#else
1046    return IsPrime(p-1);
1047#endif
1048}
1049
1050void PrepareChinese (int n) // initialization for CRA
1051{
1052     int i,k;
1053     modp_result_entry *cur_ptr;
1054     cur_ptr=modp_result;
1055     modp_number *congr_ptr;
1056     modp_number prod;
1057     in_gamma=(modp_number*)omAlloc0(sizeof(modp_number)*n);
1058     congr=(modp_number*)omAlloc0(sizeof(modp_number)*n);
1059     congr_ptr=congr;
1060     while (cur_ptr!=NULL)
1061     {
1062        *congr_ptr=cur_ptr->p;
1063        cur_ptr=cur_ptr->next;
1064        congr_ptr++;
1065     }
1066     for (k=1;k<n;k++)
1067     {
1068         prod=congr[0]%congr[k];
1069         for (i=1;i<=k-1;i++) prod=(prod*congr[i])%congr[k];
1070         in_gamma[i]=OneInverse(prod,congr[k]);
1071     }
1072     mpz_init(bigcongr);
1073     mpz_set_ui(bigcongr,congr[0]);
1074     for (k=1;k<n;k++) mpz_mul_ui(bigcongr,bigcongr,congr[k]);
1075}
1076
1077void CloseChinese (int n) // after CRA
1078{
1079     omFree(in_gamma);
1080     omFree(congr);
1081     mpz_clear(bigcongr);
1082}
1083
1084void ClearGCD () // divides polynomials in basis by gcd of coefficients
1085{
1086     bool first_gcd=true;
1087     int i;
1088     mpz_t g;
1089     mpz_init(g);
1090     for (i=0;i<=final_base_dim;i++)
1091     {
1092         if (mpz_sgn(polycoef[i])!=0)
1093         {
1094            if (first_gcd)
1095            {
1096               first_gcd=false;
1097               mpz_set(g,polycoef[i]);
1098            }
1099            else
1100               mpz_gcd(g,g,polycoef[i]);
1101         }
1102     }
1103     for (i=0;i<=final_base_dim;i++) mpz_divexact(polycoef[i],polycoef[i],g);
1104     mpz_clear(g);
1105}
1106
1107void ReconstructGenerator (int ngen,int n,bool show) // recostruction of generator from various modp comp
1108{
1109     int size;
1110     int i,j,k;
1111     int coef;
1112     char *str=NULL;
1113     str=(char*)omAlloc0(sizeof(char)*1000);
1114     modp_result_entry *cur_ptr;
1115     generator_entry *cur_gen;
1116     modp_number *u;
1117     modp_number *v;
1118     modp_number temp;
1119     mpz_t sol;
1120     mpz_t nsol;
1121     mpz_init(sol);
1122     mpz_init(nsol);
1123     u=(modp_number*)omAlloc0(sizeof(modp_number)*n);
1124     v=(modp_number*)omAlloc0(sizeof(modp_number)*n);
1125     for (coef=0;coef<=final_base_dim;coef++)
1126     {
1127         i=0;
1128         cur_ptr=modp_result;
1129         while (cur_ptr!=NULL)
1130         {
1131            cur_gen=cur_ptr->generator;
1132            for (j=0;j<ngen;j++) cur_gen=cur_gen->next; // we have to take jth generator
1133            if (coef<final_base_dim) u[i]=cur_gen->coef[coef]; else u[i]=cur_gen->ltcoef;
1134            cur_ptr=cur_ptr->next;
1135            i++;
1136         }
1137         v[0]=u[0]; // now CRA begins
1138         for (k=1;k<n;k++)
1139         {
1140             temp=v[k-1];
1141             for (j=k-2;j>=0;j--) temp=(temp*congr[j]+v[j])%congr[k];
1142             temp=u[k]-temp;
1143             if (temp<0) temp=temp+congr[k];
1144             v[k]=(temp*in_gamma[k])%congr[k];
1145         }
1146         mpz_set_si(sol,v[n-1]);
1147         for (k=n-2;k>=0;k--)
1148         {
1149             mpz_mul_ui(sol,sol,congr[k]);
1150             mpz_add_ui(sol,sol,v[k]);
1151         }          // now CRA ends
1152         mpz_sub(nsol,sol,bigcongr);
1153         int s=mpz_cmpabs(sol,nsol);
1154         if (s>0) mpz_set(sol,nsol); // chooses representation closer to 0
1155         mpz_set(polycoef[coef],sol);
1156         if (coef<final_base_dim)
1157            memcpy(polyexp[coef],generic_column_name[coef],sizeof(exponent)*variables);
1158         else
1159            memcpy(polyexp[coef],MonListElement (generic_lt,ngen),sizeof(exponent)*variables);
1160#ifdef checksize
1161         size=mpz_sizeinbase(sol,10);
1162         if (size>maximal_size) maximal_size=size;
1163#endif
1164     }
1165     omFree(u);
1166     omFree(v);
1167     omFree(str);
1168     ClearGCD ();
1169     mpz_clear(sol);
1170     mpz_clear(nsol);
1171}
1172
1173void Discard ()  // some unlucky prime occures
1174{
1175     modp_result_entry *temp;
1176#ifdef writemsg
1177     Print(", prime=%d", cur_result->p);
1178#endif
1179     bad_primes++;
1180     if (good_primes>bad_primes)
1181     {
1182#ifdef writemsg
1183        Print("-discarding this comp (%dth)\n", n_results);
1184#endif
1185        temp=cur_result;
1186        cur_result=cur_result->prev;
1187        cur_result->next=NULL;
1188        n_results--;
1189        FreeResultEntry (temp);
1190     }
1191     else
1192     {
1193#ifdef writemsg
1194        PrintS("-discarding ALL.\n");
1195#endif
1196        int i;
1197        modp_result_entry *ntfree;
1198        generator_entry *cur_gen;
1199        temp=cur_result->prev;
1200        while (temp!=NULL)
1201        {
1202            ntfree=temp->prev;
1203            FreeResultEntry (temp);
1204            temp=ntfree;
1205        }
1206        modp_result=cur_result;
1207        cur_result->prev=NULL;
1208        n_results=1;
1209        good_primes=1;
1210        bad_primes=0;
1211        generic_n_generators=cur_result->n_generators;
1212        cur_gen=cur_result->generator;
1213        generic_lt=FreeMonList (generic_lt);
1214        for (i=0;i<generic_n_generators;i++)
1215        {
1216            generic_lt=MonListAdd (generic_lt,cur_gen->lt);
1217            cur_gen=cur_gen->next;
1218        }
1219        for (i=0;i<final_base_dim;i++) memcpy(generic_column_name[i],column_name[i],sizeof(exponent)*variables);
1220     }
1221}
1222
1223void modp_SetColumnNames ()  // used by modp - sets column names, the old table will be destroyed
1224{
1225    int i;
1226    for (i=0;i<final_base_dim;i++) memcpy(generic_column_name[i],column_name[i],sizeof(exponent)*variables);
1227}
1228
1229void CheckColumnSequence () // checks if scheme of computations is as generic one
1230{
1231     int i;
1232     if (cur_result->n_generators!=generic_n_generators)
1233     {
1234#ifdef writemsg
1235        PrintS("wrong number of generators occured");
1236#else
1237        if (protocol) PrintS("ng");
1238#endif
1239        Discard ();
1240        return;
1241     }
1242     if (denom_divisible)
1243     {
1244#ifdef writemsg
1245        PrintS("denom of coef divisible by p");
1246#else
1247        if (protocol) PrintS("dp");
1248#endif
1249        Discard ();
1250        return;
1251     }
1252     generator_entry *cur_gen;
1253     mon_list_entry *cur_mon;
1254     cur_gen=cur_result->generator;
1255     cur_mon=generic_lt;
1256     for (i=0;i<generic_n_generators;i++)
1257     {
1258         if (!EqualMon(cur_mon->mon,cur_gen->lt))
1259         {
1260#ifdef writemsg
1261            PrintS("wrong leading term occured");
1262#else
1263            if (protocol) PrintS("lt");
1264#endif
1265            Discard ();
1266            return;
1267         }
1268         cur_gen=cur_gen->next;
1269         cur_mon=cur_mon->next;
1270     }
1271     for (i=0;i<final_base_dim;i++)
1272     {
1273         if (!EqualMon(generic_column_name[i],column_name[i]))
1274         {
1275#ifdef writemsg
1276            PrintS("wrong seq of cols occured");
1277#else
1278            if (protocol) PrintS("sc");
1279#endif
1280            Discard ();
1281            return;
1282         }
1283     }
1284     good_primes++;
1285}
1286#if 0 /* only debuggig */
1287void WriteGenerator () // writes generator (only for debugging)
1288{
1289     char *str;
1290     str=(char*)omAlloc0(sizeof(char)*1000);
1291     int i;
1292     for (i=0;i<=final_base_dim;i++)
1293     {
1294         str=mpz_get_str(str,10,polycoef[i]);
1295         PrintS(str);
1296         PrintS("*");
1297         WriteMono(polyexp[i]);
1298         PrintS(" ");
1299     }
1300     omFree(str);
1301     PrintLn();
1302}
1303#endif
1304
1305bool CheckGenerator () // evaluates generator to check whether it is good
1306{
1307     mpz_t val,sum,temp;
1308     int con,i;
1309     mpz_init(val);
1310     mpz_init(sum);
1311     for (con=0;con<final_base_dim;con++)
1312     {
1313       mpz_set_si(sum,0);
1314       for (i=0;i<=final_base_dim;i++)
1315       {
1316         int_Evaluate(val, polyexp[i], condition_list[con]);
1317         mpz_mul(val,val,polycoef[i]);
1318         mpz_add(sum,sum,val);
1319       }
1320       if (mpz_sgn(sum)!=0)
1321       {
1322          mpz_clear(val);
1323          mpz_clear(sum);
1324          return false;
1325       }
1326    }
1327    mpz_clear(val);
1328    mpz_clear(sum);
1329    return true;
1330}
1331
1332void ClearGenList ()
1333{
1334     gen_list_entry *temp;
1335     int i;
1336     while (gen_list!=NULL)
1337     {
1338         temp=gen_list->next;
1339         for (i=0;i<=final_base_dim;i++)
1340         {
1341             mpz_clear(gen_list->polycoef[i]);
1342             omFree(gen_list->polyexp[i]);
1343         }
1344         omFree(gen_list->polycoef);
1345         omFree(gen_list->polyexp);
1346         omFree(gen_list);
1347         gen_list=temp;
1348      }
1349}
1350
1351void UpdateGenList ()
1352{
1353     gen_list_entry *temp,*prev;
1354     int i,j;
1355     prev=NULL;
1356     temp=gen_list;
1357     exponent deg;
1358     for (i=0;i<=final_base_dim;i++)
1359     {
1360         deg=MonDegree(polyexp[i]);
1361         for (j=0;j<deg;j++)
1362         {
1363             mpz_mul(polycoef[i],polycoef[i],common_denom);
1364         }
1365     }
1366     ClearGCD ();
1367     while (temp!=NULL)
1368     {
1369         prev=temp;
1370         temp=temp->next;
1371     }
1372     temp=(gen_list_entry*)omAlloc0(sizeof(gen_list_entry));
1373     if (prev==NULL) gen_list=temp; else prev->next=temp;
1374     temp->next=NULL;
1375     temp->polycoef=(mpz_t*)omAlloc(sizeof(mpz_t)*(final_base_dim+1));
1376     temp->polyexp=(mono_type*)omAlloc(sizeof(mono_type)*(final_base_dim+1));
1377     for (i=0;i<=final_base_dim;i++)
1378     {
1379         mpz_init(temp->polycoef[i]);
1380         mpz_set(temp->polycoef[i],polycoef[i]);
1381         temp->polyexp[i]=ZeroMonomial ();
1382         memcpy(temp->polyexp[i],polyexp[i],sizeof(exponent)*variables);
1383     }
1384}
1385
1386#if 0 /* only debugging */
1387void ShowGenList ()
1388{
1389     gen_list_entry *temp;
1390     int i;
1391     char *str;
1392     str=(char*)omAlloc0(sizeof(char)*1000);
1393     temp=gen_list;
1394     while (temp!=NULL)
1395     {
1396         PrintS("generator: ");
1397         for (i=0;i<=final_base_dim;i++)
1398         {
1399             str=mpz_get_str(str,10,temp->polycoef[i]);
1400             PrintS(str);
1401             PrintS("*");
1402             WriteMono(temp->polyexp[i]);
1403         }
1404         PrintLn();
1405         temp=temp->next;
1406     }
1407     omFree(str);
1408}
1409#endif
1410
1411
1412void modp_Main ()
1413{
1414     mono_type cur_mon;
1415     cur_mon= ZeroMonomial ();
1416     modp_denom=1;
1417     bool row_is_zero;
1418
1419#if 0 /* only debugging */
1420     Info ();
1421#endif
1422
1423     while (check_list!=NULL)
1424     {
1425           TakeNextMonomial (cur_mon);
1426           ProduceRow (cur_mon);
1427#if 0 /* only debugging */
1428           cout << "row produced for monomial ";
1429           WriteMono (cur_mon);
1430           cout << endl;
1431           Info ();
1432#endif
1433           ReduceRow ();
1434           row_is_zero = RowIsZero ();
1435           if (row_is_zero)
1436           {
1437              lt_list=MonListAdd (lt_list,cur_mon);
1438              ReduceCheckListByMon (cur_mon);
1439              NewGenerator (cur_mon);
1440#if 0 /* only debugging */
1441              cout << "row is zero - linear dependence found (should be seen in my_solve_row)" << endl;
1442              cout << "monomial added to leading terms list" << endl;
1443              cout << "check list updated" << endl;
1444              Info ();
1445#endif
1446           }
1447           else
1448           {
1449              base_list= MonListAdd (base_list,cur_mon);
1450              UpdateCheckList (cur_mon);
1451              ReduceCheckListByLTs ();
1452#if 0 /* only debugging */
1453              cout << "row is non-zero" << endl;
1454              cout << "monomial added to quotient basis list" << endl;
1455              cout << "new monomials added to check list" << endl;
1456              cout << "check list reduced by monomials from leading term list" << endl;
1457              Info ();
1458#endif
1459              PrepareRow (cur_mon);
1460#if 0 /* only debugging */
1461              cout << "row prepared and put into matrix" << endl;
1462              Info ();
1463#endif
1464           }
1465        }
1466        omFree(cur_mon);
1467}
1468
1469void ResolveCoeff (mpq_t c, number m)
1470{
1471     if ((long)m & SR_INT)
1472     {
1473        long m_val=SR_TO_INT(m);
1474        mpq_set_si(c,m_val,1);
1475     }
1476     else
1477     {
1478        if (m->s<2)
1479        {
1480           mpz_set(mpq_numref(c),m->z);
1481           mpz_set(mpq_denref(c),m->n);
1482           mpq_canonicalize(c);
1483        }
1484        else
1485        {
1486           mpq_set_z(c,m->z);
1487        }
1488     }
1489}
1490
1491ideal interpolation(lists L, intvec *v)
1492{
1493  protocol=TEST_OPT_PROT;  // should be set if option(prot) is enabled
1494
1495  bool data_ok=true;
1496
1497  // reading the ring data ***************************************************
1498  if ((currRing==NULL) || ((!rField_is_Zp (currRing))&&(!rField_is_Q (currRing))))
1499  {
1500     WerrorS("coefficient field should be Zp or Q!");
1501     return NULL;
1502  }
1503  if ((currRing->qideal)!=NULL)
1504  {
1505     WerrorS("quotient ring not supported!");
1506     return NULL;
1507  }
1508  if ((currRing->OrdSgn)!=1)
1509  {
1510     WerrorS("ordering must be global!");
1511     return NULL;
1512  }
1513  n_points=v->length ();
1514  if (n_points!=(L->nr+1))
1515  {
1516     WerrorS("list and intvec must have the same length!");
1517     return NULL;
1518  }
1519  variables=currRing->N;
1520  only_modp=rField_is_Zp(currRing);
1521  if (only_modp) myp=rChar(currRing);
1522  // ring data read **********************************************************
1523
1524
1525  multiplicity=(int*)malloc(sizeof(int)*n_points);
1526  int i;
1527  for (i=0;i<n_points;i++) multiplicity[i]=(*v)[i];
1528
1529  final_base_dim = CalcBaseDim ();
1530
1531#ifdef writemsg
1532  Print("number of variables: %d\n", variables);
1533  Print("number of points: %d\n", n_points);
1534  PrintS("multiplicities: ");
1535  for (i=0;i<n_points;i++) Print("%d ", multiplicity[i]);
1536  PrintLn();
1537  Print("general initialization for dimension %d ...\n", final_base_dim);
1538#endif
1539
1540  GeneralInit ();
1541
1542// reading coordinates of points from ideals **********************************
1543  mpq_t divisor;
1544  if (!only_modp) mpq_init(divisor);
1545  int j;
1546  for(i=0; i<=L->nr;i++)
1547  {
1548    ideal I=(ideal)L->m[i].Data();
1549    for(j=0;j<IDELEMS(I);j++)
1550    {
1551      poly p=I->m[j];
1552      if (p!=NULL)
1553      {
1554        poly ph=pHead(p);
1555        int pcvar=pVar(ph);
1556        if (pcvar!=0)
1557        {
1558          pcvar--;
1559          if (coord_exist[i][pcvar])
1560          {
1561             Print("coordinate %d for point %d initialized twice!\n",pcvar+1,i+1);
1562             data_ok=false;
1563          }
1564          number m;
1565          m=pGetCoeff(p); // possible coefficient standing by a leading monomial
1566          if (!only_modp) ResolveCoeff (divisor,m);
1567          number n;
1568          if (pNext(p)!=NULL) n=pGetCoeff(pNext(p));
1569          else n=nInit(0);
1570          if (only_modp)
1571          {
1572            n=nNeg(n);
1573            n=nDiv(n,m);
1574            modp_points[i][pcvar]=(int)((long)n);
1575          }
1576          else
1577          {
1578            ResolveCoeff (q_points[i][pcvar],n);
1579            mpq_neg(q_points[i][pcvar],q_points[i][pcvar]);
1580            mpq_div(q_points[i][pcvar],q_points[i][pcvar],divisor);
1581          }
1582          coord_exist[i][pcvar]=true;
1583        }
1584        else
1585        {
1586          PrintS("not a variable? ");
1587          wrp(p);
1588          PrintLn();
1589          data_ok=false;
1590        }
1591        pDelete(&ph);
1592      }
1593    }
1594  }
1595  if (!only_modp) mpq_clear(divisor);
1596  // data from ideal read *******************************************************
1597
1598  // ckecking if all coordinates are initialized
1599  for (i=0;i<n_points;i++)
1600  {
1601      for (j=0;j<variables;j++)
1602      {
1603          if (!coord_exist[i][j])
1604          {
1605             Print("coordinate %d for point %d not known!\n",j+1,i+1);
1606             data_ok=false;
1607          }
1608      }
1609  }
1610
1611  if (!data_ok)
1612  {
1613     GeneralDone();
1614     WerrorS("data structure is invalid");
1615     return NULL;
1616  }
1617
1618  if (!only_modp) IntegerPoints ();
1619  MakeConditions ();
1620#ifdef writemsg
1621  PrintS("done.\n");
1622#else
1623  if (protocol) Print("[vdim %d]",final_base_dim);
1624#endif
1625
1626
1627// main procedure *********************************************************************
1628  int modp_cycles=10;
1629  bool correct_gen=false;
1630  if (only_modp) modp_cycles=1;
1631  #ifdef HAVE_FACTORY
1632  myp_index=cf_getNumSmallPrimes ();
1633  #endif
1634
1635  while ((!correct_gen)&&(myp_index>1))
1636  {
1637#ifdef writemsg
1638        Print("trying %d cycles mod p...\n",modp_cycles);
1639#else
1640        if (protocol) Print("(%d)",modp_cycles);
1641#endif
1642        while ((n_results<modp_cycles)&&(myp_index>1))  // some computations mod p
1643        {
1644            if (!only_modp) myp=TakePrime (myp);
1645            NewResultEntry ();
1646            InitProcData ();
1647            if (only_modp) modp_PrepareProducts (); else int_PrepareProducts ();
1648
1649            modp_Main ();
1650
1651            if (!only_modp)
1652            {
1653               MultGenerators ();
1654               CheckColumnSequence ();
1655            }
1656            else
1657            {
1658               modp_SetColumnNames ();
1659            }
1660            FreeProcData ();
1661        }
1662
1663        if (!only_modp)
1664        {
1665           PrepareChinese (modp_cycles);
1666           correct_gen=true;
1667           for (i=0;i<generic_n_generators;i++)
1668           {
1669               ReconstructGenerator (i,modp_cycles,false);
1670               correct_gen=CheckGenerator ();
1671               if (!correct_gen)
1672               {
1673#ifdef writemsg
1674                  PrintS("wrong generator!\n");
1675#else
1676//                  if (protocol) PrintS("!g");
1677#endif
1678                  ClearGenList ();
1679                  break;
1680               }
1681               else
1682               {
1683                  UpdateGenList ();
1684               }
1685           }
1686#ifdef checksize
1687           Print("maximal size of output: %d, precision bound: %d.\n",maximalsize,mpz_sizeinbase(bigcongr,10));
1688#endif
1689           CloseChinese (modp_cycles);
1690           modp_cycles=modp_cycles+10;
1691        }
1692        else
1693        {
1694           correct_gen=true;
1695        }
1696  }
1697// end of main procedure ************************************************************************************
1698
1699#ifdef writemsg
1700  PrintS("computations finished.\n");
1701#else
1702  if (protocol) PrintLn();
1703#endif
1704
1705  if (!correct_gen)
1706  {
1707     GeneralDone ();
1708     ClearGenList ();
1709     WerrorS("internal error - coefficient too big!");
1710     return NULL;
1711  }
1712
1713// passing data to ideal *************************************************************************************
1714  ideal ret;
1715
1716  if (only_modp)
1717  {
1718    mono_type mon;
1719    ret=idInit(modp_result->n_generators,1);
1720    generator_entry *cur_gen=modp_result->generator;
1721    for(i=0;i<IDELEMS(ret);i++)
1722    {
1723      poly p,sum;
1724      sum=NULL;
1725      int a;
1726      int cf;
1727      for (a=final_base_dim;a>=0;a--)
1728      {
1729        if (a==final_base_dim) cf=cur_gen->ltcoef; else cf=cur_gen->coef[a];
1730        if (cf!=0)
1731        {
1732            p=pISet(cf);
1733            if (a==final_base_dim) mon=cur_gen->lt; else mon=generic_column_name[a];
1734            for (j=0;j<variables;j++) pSetExp(p,j+1,mon[j]);
1735            pSetm(p);
1736            sum=pAdd(sum,p);
1737        }
1738      }
1739      ret->m[i]=sum;
1740      cur_gen=cur_gen->next;
1741    }
1742  }
1743  else
1744  {
1745    ret=idInit(generic_n_generators,1);
1746    gen_list_entry *temp=gen_list;
1747    for(i=0;i<IDELEMS(ret);i++)
1748    {
1749      poly p,sum;
1750      sum=NULL;
1751      int a;
1752      for (a=final_base_dim;a>=0;a--) // build one polynomial
1753      {
1754          if (mpz_sgn(temp->polycoef[a])!=0)
1755          {
1756             number n=ALLOC_RNUMBER();
1757#ifdef LDEBUG
1758             n->debug=123456;
1759#endif
1760             mpz_init_set(n->z,temp->polycoef[a]);
1761             n->s=3;
1762             nlNormalize(n, currRing->cf);
1763             p=pNSet(n); //a monomial
1764             for (j=0;j<variables;j++) pSetExp(p,j+1,temp->polyexp[a][j]);
1765             pSetm(p); // after all pSetExp
1766             sum=pAdd(sum,p);
1767          }
1768      }
1769      ret->m[i]=sum;
1770      temp=temp->next;
1771    }
1772  }
1773// data transferred ****************************************************************************
1774
1775
1776  GeneralDone ();
1777  ClearGenList ();
1778  return ret;
1779}
1780
1781
1782BOOLEAN jjINTERPOLATION (leftv res, leftv l, leftv v)
1783{
1784  res->data=interpolation((lists)l->Data(),(intvec*)v->Data());
1785  setFlag(res,FLAG_STD);
1786  return errorreported;
1787}
Note: See TracBrowser for help on using the repository browser.