source: git/Singular/interpolation.cc @ e4e36c

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