source: git/Singular/interpolation.cc @ 880433

spielwiese
Last change on this file since 880433 was 880433, checked in by Hans Schoenemann <hannes@…>, 11 years ago
fix: made local proc static, removed unsed variables (interpolation.cc)
  • Property mode set to 100644
File size: 49.3 KB
RevLine 
[0bc2dd]1/****************************************
2*  Computer Algebra System SINGULAR     *
3****************************************/
4
[762407]5#include "config.h"
[b1dfaf]6#include <kernel/mod2.h>
[0fb34ba]7#include <misc/options.h>
[cffea8f]8#ifdef HAVE_FACTORY
9#define SI_DONT_HAVE_GLOBAL_VARS
10#  include <factory/factory.h>
11#endif /* HAVE_FACTORY */
[5ff1d3]12
[0fb34ba]13#include <misc/intvec.h>
[5ff1d3]14
15#include <polys/monomials/ring.h>
[737a68]16#include <kernel/polys.h>
[5ff1d3]17
[0519cc]18#include <coeffs/longrat.h>
[5ff1d3]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"
[b5e57e2]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
[880433]39static modp_number myp;  // all modp computation done mod myp
40static int myp_index; // index of small prime in Singular table of primes
[b5e57e2]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
[880433]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
[b5e57e2]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
[880433]127static gen_list_entry *gen_list=NULL; // list of resulting generators - output data (integer version)
[b5e57e2]128
[880433]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)
[b5e57e2]136
[880433]137static poly comparizon_p1;  //polynomials used to do comparizons by Singular
138static poly comparizon_p2;
[b5e57e2]139
[880433]140static modp_number *modp_Reverse; // reverses in mod p
[b5e57e2]141
[880433]142static bool protocol; // true to show the protocol
[b5e57e2]143
144#ifdef checksize
[880433]145static int maximal_size=0;
[b5e57e2]146#endif
147
[30d574]148#if 0  /* only for debuggig*/
[b5e57e2]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]);
[6b4ff12]153     PrintS(" ");
[b5e57e2]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
[880433]227static modp_number OneInverse(modp_number a,modp_number p) // computes inverse of d mod p without using tables
[b5e57e2]228{
229   long  u, v, u0, u1, u2, q, r;
230   u1=1; u2=0;
231   u = a; v = p;
[30d574]232   while (v != 0)
233   {
[b5e57e2]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   {
[6b4ff12]246     PrintS("?");
[b5e57e2]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
[880433]256static int CalcBaseDim () // returns the maximal (expected) dimension of quotient space
[b5e57e2]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
[880433]275static bool EqualMon (mono_type m1, mono_type m2)  // compares two monomials, true if equal, false otherwise
[b5e57e2]276{
277     int i;
278     for (i=0;i<variables;i++)
279         if (m1[i]!=m2[i]) return false;;
280     return true;
281}
282
[880433]283static exponent MonDegree (mono_type mon)  // computes the degree of a monomial
[b5e57e2]284{
285     exponent v=0;
286     int i;
287     for (i=0;i<variables;i++) v=v+mon[i];
288     return v;
289}
290
[880433]291static bool Greater (mono_type m1, mono_type m2) // true if m1 > m2, false otherwise. uses Singular comparing function
[b5e57e2]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
[880433]304static mon_list_entry* MonListAdd (mon_list_entry *list, mono_type mon)  // adds one entry to the list of monomials structure
[b5e57e2]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;
[30d574]315           curptr=curptr->next;
316     }
317     temp=(mon_list_entry*)omAlloc0(sizeof(mon_list_entry));
[b5e57e2]318     (*temp).next=curptr;
[aad4e14]319     (*temp).mon=(exponent*)omAlloc(sizeof(exponent)*variables);
[b5e57e2]320     memcpy(temp->mon,mon,sizeof(exponent)*variables);
[30d574]321     if (prevptr==NULL) return temp;
322     else
323     {
[b5e57e2]324          prevptr->next=temp;
[30d574]325          return list;
326     }
[b5e57e2]327}
328
[880433]329static mono_type MonListElement (mon_list_entry *list, int n)  // returns ith element from list of monomials - no error checking!
[b5e57e2]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
[880433]337static mono_type ZeroMonomial ()  // allocates memory for new monomial with all enries equal 0
[b5e57e2]338{
339     mono_type m;
[30d574]340     m=(exponent*)omAlloc0(sizeof(exponent)*variables);
[b5e57e2]341     return m;
342}
343
[880433]344static void GeneralInit ()  // general initialization
[b5e57e2]345{
[880433]346     int i,j;
[aad4e14]347     points=(coordinates*)omAlloc(sizeof(coordinates)*n_points);
[b5e57e2]348     for (i=0;i<n_points;i++)
349     {
[aad4e14]350         points[i]=(coordinate_products*)omAlloc(sizeof(coordinate_products)*variables);
[30d574]351         for (j=0;j<variables;j++) points[i][j]=(modp_number*)omAlloc0(sizeof(modp_number)*(max_coord));
[b5e57e2]352     }
[30d574]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);
[aad4e14]355     modp_points=(modp_coordinates*)omAlloc(sizeof(modp_coordinates)*n_points);
[30d574]356     for (i=0;i<n_points;i++) modp_points[i]=(modp_number*)omAlloc0(sizeof(modp_number)*variables);
[b5e57e2]357     if (!only_modp)
358     {
[30d574]359        q_points=(q_coordinates*)omAlloc0(sizeof(q_coordinates)*n_points);
[b5e57e2]360        for (i=0;i<n_points;i++)
361        {
[aad4e14]362            q_points[i]=(mpq_t*)omAlloc(sizeof(mpq_t)*variables);
[b5e57e2]363            for (j=0;j<variables;j++) mpq_init(q_points[i][j]);
364        }
[30d574]365        int_points=(int_coordinates*)omAlloc0(sizeof(int_coordinates)*n_points);
[b5e57e2]366        for (i=0;i<n_points;i++)
367        {
[aad4e14]368            int_points[i]=(mpz_t*)omAlloc(sizeof(mpz_t)*variables);
[b5e57e2]369            for (j=0;j<variables;j++) mpz_init(int_points[i][j]);
370        }
371     }
[aad4e14]372     coord_exist=(coord_exist_table*)omAlloc(sizeof(coord_exist_table)*n_points);
[b5e57e2]373     for (i=0;i<n_points;i++)
374     {
[30d574]375         coord_exist[i]=(bool*)omAlloc0(sizeof(bool)*variables);
[b5e57e2]376     }
[aad4e14]377     generic_column_name=(mono_type*)omAlloc(sizeof(mono_type)*final_base_dim);
[b5e57e2]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     {
[aad4e14]384        polycoef=(mpz_t*)omAlloc(sizeof(mpz_t)*(final_base_dim+1));
385        polyexp=(mono_type*)omAlloc(sizeof(mono_type)*(final_base_dim+1));
[b5e57e2]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
[880433]405static void InitProcData ()  // initialization for procedures which do computations mod p
[b5e57e2]406{
[880433]407     int i;
[b5e57e2]408     check_list=MonListAdd (check_list,ZeroMonomial ());
[30d574]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);
[aad4e14]411     column_name=(mono_type*)omAlloc(sizeof(mono_type)*final_base_dim);
[b5e57e2]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;
[30d574]417     modp_Reverse=(modp_number*)omAlloc0(sizeof(modp_number)*myp);
[b5e57e2]418
419// produces table of modp inverts by finding a generator of (Z_myp*,*)
[aad4e14]420     gen_table=(modp_number*)omAlloc(sizeof(modp_number)*myp);
[b5e57e2]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;
[30d574]438     omFree(gen_table);
[b5e57e2]439}
440
[880433]441static mon_list_entry* FreeMonList (mon_list_entry *list)  // destroys a list of monomials, returns NULL pointer
[b5e57e2]442{
443     mon_list_entry *ptr;
444     mon_list_entry *pptr;
445     ptr=list;
446     while (ptr!=NULL)
447     {
448           pptr=ptr->next;
[30d574]449           omFree(ptr->mon);
450           omFree(ptr);
[b5e57e2]451           ptr=pptr;}
452     return NULL;
453}
454
[880433]455static void GeneralDone ()  // to be called before exit to free memory
[b5e57e2]456{
[880433]457     int i,j;
[b5e57e2]458     for (i=0;i<n_points;i++)
459     {
460         for (j=0;j<variables;j++)
461         {
[30d574]462             omFree(points[i][j]);
[b5e57e2]463         }
[30d574]464         omFree(points[i]);
[b5e57e2]465     }
[30d574]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);
[b5e57e2]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]);
[30d574]476            omFree(q_points[i]);
[b5e57e2]477        }
[30d574]478        omFree(q_points);
[b5e57e2]479        for (i=0;i<n_points;i++)
480        {
481            for (j=0;j<variables;j++) mpz_clear(int_points[i][j]);
[30d574]482            omFree(int_points[i]);
[b5e57e2]483        }
[30d574]484        omFree(int_points);
[b5e57e2]485        generic_lt=FreeMonList (generic_lt);
486        for (i=0;i<=final_base_dim;i++)
487        {
488            mpz_clear(polycoef[i]);
[30d574]489            omFree(polyexp[i]);
[b5e57e2]490        }
[30d574]491        omFree(polycoef);
492        omFree(polyexp);
[b5e57e2]493        if (!only_modp) mpz_clear(common_denom);
494     }
495     for (i=0;i<final_base_dim;i++)
496     {
[30d574]497         omFree(generic_column_name[i]);
[b5e57e2]498     }
[30d574]499     omFree(generic_column_name);
500     for (i=0;i<n_points;i++) omFree(coord_exist[i]);
501     omFree(coord_exist);
[b5e57e2]502     pDelete(&comparizon_p1);
503     pDelete(&comparizon_p2);
504}
505
[880433]506static void FreeProcData ()  // to be called after one modp computation to free memory
[b5e57e2]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);
[30d574]514     omFree(my_row);
[b5e57e2]515     my_row=NULL;
[30d574]516     omFree(my_solve_row);
[b5e57e2]517     my_solve_row=NULL;
518     ptr=row_list;
519     while (ptr!=NULL)
520     {
521           pptr=ptr->next;
[30d574]522           omFree(ptr->row_matrix);
523           omFree(ptr->row_solve);
524           omFree(ptr);
[b5e57e2]525           ptr=pptr;
526     }
527     row_list=NULL;
[30d574]528     for (i=0;i<final_base_dim;i++) omFree(column_name[i]);
529     omFree(column_name);
530     omFree(modp_Reverse);
[b5e57e2]531}
532
[880433]533static void modp_Evaluate(modp_number *ev, mono_type mon, condition_type con)  // evaluates monomial on condition (modp)
[b5e57e2]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;
[aad4e14]542    mn=(exponent*)omAlloc(sizeof(exponent)*variables);
[b5e57e2]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    }
[30d574]553    omFree(mn);
[b5e57e2]554}
555
[880433]556static void int_Evaluate(mpz_t ev, mono_type mon, condition_type con) // (***) evaluates monomial on condition for integer numbers
[b5e57e2]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;
[aad4e14]567    mn=(exponent*)omAlloc(sizeof(exponent)*variables);
[b5e57e2]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    }
[30d574]579    omFree(mn);
[b5e57e2]580    mpz_clear(mon_conv);
581}
582
[880433]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
[b5e57e2]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
[880433]604static void IntegerPoints ()  // produces integer points from rationals by scaling the coordinate system
[b5e57e2]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
[880433]632static void int_PrepareProducts ()  // prepares coordinates of points and products for modp case (from integer coefs)
[b5e57e2]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
[880433]656static void modp_PrepareProducts () // prepares products for modp computation from modp data
[b5e57e2]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
[880433]670static void MakeConditions ()  // prepares a list of conditions from list of multiplicities
[b5e57e2]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     }
[30d574]698     omFree(mon);
[b5e57e2]699}
700
[880433]701static void ReduceRow ()  // reduces my_row by previous rows, does the same for solve row
[b5e57e2]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;
[30d574]770#if 0 /* only debugging */
[6b4ff12]771           PrintS("reduction by row ");
[b5e57e2]772           Info ();
773#endif
774     }
775}
776
[880433]777static bool RowIsZero () // check whether a row is zero
[b5e57e2]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
[880433]791static bool DivisibleMon (mono_type m1, mono_type m2) // checks whether m1 is divisible by m2
[b5e57e2]792{
793     int i;
794     for (i=0;i<variables;i++)
795         if (m1[i]>m2[i]) return false;;
796     return true;
797}
798
[880433]799static void ReduceCheckListByMon (mono_type m)  // from check_list monomials divisible by m are thrown out
[b5e57e2]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          {
[30d574]810             if (p_ptr==NULL)
[b5e57e2]811                check_list=c_ptr->next;
812             else
813                p_ptr->next=c_ptr->next;
814             n_ptr=c_ptr->next;
[30d574]815             omFree(c_ptr->mon);
816             omFree(c_ptr);
[b5e57e2]817             c_ptr=n_ptr;
818          }
819          else
820          {
821              p_ptr=c_ptr;
822              c_ptr=c_ptr->next;
823          }
824     }
825}
826
[880433]827static void TakeNextMonomial (mono_type mon)  // reads first monomial from check_list, then it is deleted
[b5e57e2]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;
[30d574]834        omFree(check_list->mon);
835        omFree(check_list);
[b5e57e2]836        check_list=n_check_list;
837     }
838}
839
[880433]840static void UpdateCheckList (mono_type m) // adds all monomials which are "next" to m (divisible by m and degree +1)
[b5e57e2]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
[880433]851static void ReduceCheckListByLTs ()  // deletes all monomials from check_list which are divisible by one of the leading terms
[b5e57e2]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
[880433]862static void RowListAdd (int first_col, mono_type mon) // puts a row into matrix
[b5e57e2]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     }
[30d574]877     temp=(row_list_entry*)omAlloc0(sizeof(row_list_entry));
878     (*temp).row_matrix=(modp_number*)omAlloc0(sizeof(modp_number)*final_base_dim);
[b5e57e2]879     memcpy((*temp).row_matrix,my_row,sizeof(modp_number)*(final_base_dim));
[30d574]880     (*temp).row_solve=(modp_number*)omAlloc0(sizeof(modp_number)*final_base_dim);
[b5e57e2]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
[880433]889static void PrepareRow (mono_type mon) // prepares a row and puts it into matrix
[b5e57e2]890{
891     modp_number *row;
892     int first_col=-1;
893     int col;
[073295f]894     modp_number red_val=1;
[b5e57e2]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
[880433]942static void NewResultEntry ()  // creates an entry for new modp result
[b5e57e2]943{
944    modp_result_entry *temp;
[30d574]945    temp=(modp_result_entry*)omAlloc0(sizeof(modp_result_entry));
[b5e57e2]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
[880433]964static void FreeResultEntry (modp_result_entry *e) // destroys the result entry, without worrying about where it is
[b5e57e2]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;
[30d574]972         omFree(cur_gen->coef);
973         omFree(cur_gen->lt);
974         omFree(cur_gen);
[b5e57e2]975         cur_gen=next_gen;
976     }
[30d574]977     omFree(e);
[b5e57e2]978}
979
980
[880433]981static void NewGenerator (mono_type mon)  // new generator in modp comp found, shoul be stored on the list
[b5e57e2]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     }
[30d574]993     temp=(generator_entry*)omAlloc0(sizeof(generator_entry));
[aad4e14]994     if (prev_ptr==NULL) cur_result->generator=temp;
995     else prev_ptr->next=temp;
[b5e57e2]996     temp->next=NULL;
[30d574]997     temp->coef=(modp_number*)omAlloc0(sizeof(modp_number)*final_base_dim);
[b5e57e2]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
[880433]1005static void MultGenerators () // before reconstructing, all denominators must be cancelled
[b5e57e2]1006{
1007#ifndef integerstrategy
1008     int i;
1009     generator_entry *cur_ptr;
1010     cur_ptr=cur_result->generator;
1011     while (cur_ptr!=NULL)
1012     {
[30d574]1013         for (i=0;i<final_base_dim;i++)
[b5e57e2]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}
[30d574]1020#if 0 /* only debbuging */
[b5e57e2]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        }
[6b4ff12]1035        PrintS(" and LT = ");
[b5e57e2]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}
[30d574]1041#endif
[b5e57e2]1042
[880433]1043static modp_number TakePrime (modp_number p)  // takes "previous" (smaller) prime
[b5e57e2]1044{
[c16aa0]1045#ifdef HAVE_FACTORY
[b5e57e2]1046    myp_index--;
1047    return cf_getSmallPrime(myp_index);
[c16aa0]1048#else
1049    return IsPrime(p-1);
1050#endif
[b5e57e2]1051}
1052
[880433]1053static void PrepareChinese (int n) // initialization for CRA
[b5e57e2]1054{
1055     int i,k;
1056     modp_result_entry *cur_ptr;
1057     cur_ptr=modp_result;
1058     modp_number *congr_ptr;
1059     modp_number prod;
[30d574]1060     in_gamma=(modp_number*)omAlloc0(sizeof(modp_number)*n);
1061     congr=(modp_number*)omAlloc0(sizeof(modp_number)*n);
[b5e57e2]1062     congr_ptr=congr;
1063     while (cur_ptr!=NULL)
1064     {
1065        *congr_ptr=cur_ptr->p;
1066        cur_ptr=cur_ptr->next;
1067        congr_ptr++;
1068     }
1069     for (k=1;k<n;k++)
1070     {
1071         prod=congr[0]%congr[k];
1072         for (i=1;i<=k-1;i++) prod=(prod*congr[i])%congr[k];
[1f18bb]1073         in_gamma[i]=OneInverse(prod,congr[k]);
[b5e57e2]1074     }
1075     mpz_init(bigcongr);
1076     mpz_set_ui(bigcongr,congr[0]);
1077     for (k=1;k<n;k++) mpz_mul_ui(bigcongr,bigcongr,congr[k]);
1078}
1079
[880433]1080static void CloseChinese () // after CRA
[b5e57e2]1081{
[30d574]1082     omFree(in_gamma);
1083     omFree(congr);
[b5e57e2]1084     mpz_clear(bigcongr);
1085}
1086
[880433]1087static void ClearGCD () // divides polynomials in basis by gcd of coefficients
[b5e57e2]1088{
1089     bool first_gcd=true;
1090     int i;
1091     mpz_t g;
1092     mpz_init(g);
1093     for (i=0;i<=final_base_dim;i++)
1094     {
1095         if (mpz_sgn(polycoef[i])!=0)
1096         {
1097            if (first_gcd)
1098            {
1099               first_gcd=false;
1100               mpz_set(g,polycoef[i]);
1101            }
1102            else
1103               mpz_gcd(g,g,polycoef[i]);
1104         }
1105     }
1106     for (i=0;i<=final_base_dim;i++) mpz_divexact(polycoef[i],polycoef[i],g);
1107     mpz_clear(g);
1108}
1109
[880433]1110static void ReconstructGenerator (int ngen,int n) // recostruction of generator from various modp comp
[b5e57e2]1111{
1112     int i,j,k;
1113     int coef;
1114     char *str=NULL;
[30d574]1115     str=(char*)omAlloc0(sizeof(char)*1000);
[b5e57e2]1116     modp_result_entry *cur_ptr;
1117     generator_entry *cur_gen;
1118     modp_number *u;
1119     modp_number *v;
1120     modp_number temp;
1121     mpz_t sol;
1122     mpz_t nsol;
1123     mpz_init(sol);
1124     mpz_init(nsol);
[30d574]1125     u=(modp_number*)omAlloc0(sizeof(modp_number)*n);
1126     v=(modp_number*)omAlloc0(sizeof(modp_number)*n);
[b5e57e2]1127     for (coef=0;coef<=final_base_dim;coef++)
1128     {
1129         i=0;
1130         cur_ptr=modp_result;
1131         while (cur_ptr!=NULL)
1132         {
1133            cur_gen=cur_ptr->generator;
1134            for (j=0;j<ngen;j++) cur_gen=cur_gen->next; // we have to take jth generator
1135            if (coef<final_base_dim) u[i]=cur_gen->coef[coef]; else u[i]=cur_gen->ltcoef;
1136            cur_ptr=cur_ptr->next;
1137            i++;
1138         }
1139         v[0]=u[0]; // now CRA begins
1140         for (k=1;k<n;k++)
1141         {
1142             temp=v[k-1];
1143             for (j=k-2;j>=0;j--) temp=(temp*congr[j]+v[j])%congr[k];
1144             temp=u[k]-temp;
1145             if (temp<0) temp=temp+congr[k];
[1f18bb]1146             v[k]=(temp*in_gamma[k])%congr[k];
[b5e57e2]1147         }
1148         mpz_set_si(sol,v[n-1]);
1149         for (k=n-2;k>=0;k--)
1150         {
1151             mpz_mul_ui(sol,sol,congr[k]);
1152             mpz_add_ui(sol,sol,v[k]);
1153         }          // now CRA ends
1154         mpz_sub(nsol,sol,bigcongr);
1155         int s=mpz_cmpabs(sol,nsol);
1156         if (s>0) mpz_set(sol,nsol); // chooses representation closer to 0
1157         mpz_set(polycoef[coef],sol);
1158         if (coef<final_base_dim)
1159            memcpy(polyexp[coef],generic_column_name[coef],sizeof(exponent)*variables);
1160         else
1161            memcpy(polyexp[coef],MonListElement (generic_lt,ngen),sizeof(exponent)*variables);
1162#ifdef checksize
1163         size=mpz_sizeinbase(sol,10);
1164         if (size>maximal_size) maximal_size=size;
1165#endif
1166     }
[30d574]1167     omFree(u);
1168     omFree(v);
1169     omFree(str);
[b5e57e2]1170     ClearGCD ();
1171     mpz_clear(sol);
1172     mpz_clear(nsol);
1173}
1174
[880433]1175static void Discard ()  // some unlucky prime occures
[b5e57e2]1176{
1177     modp_result_entry *temp;
1178#ifdef writemsg
1179     Print(", prime=%d", cur_result->p);
1180#endif
1181     bad_primes++;
1182     if (good_primes>bad_primes)
1183     {
1184#ifdef writemsg
1185        Print("-discarding this comp (%dth)\n", n_results);
1186#endif
1187        temp=cur_result;
1188        cur_result=cur_result->prev;
1189        cur_result->next=NULL;
1190        n_results--;
1191        FreeResultEntry (temp);
1192     }
1193     else
1194     {
1195#ifdef writemsg
[6b4ff12]1196        PrintS("-discarding ALL.\n");
[b5e57e2]1197#endif
1198        int i;
1199        modp_result_entry *ntfree;
1200        generator_entry *cur_gen;
1201        temp=cur_result->prev;
1202        while (temp!=NULL)
1203        {
1204            ntfree=temp->prev;
1205            FreeResultEntry (temp);
1206            temp=ntfree;
1207        }
1208        modp_result=cur_result;
1209        cur_result->prev=NULL;
1210        n_results=1;
1211        good_primes=1;
1212        bad_primes=0;
1213        generic_n_generators=cur_result->n_generators;
1214        cur_gen=cur_result->generator;
1215        generic_lt=FreeMonList (generic_lt);
1216        for (i=0;i<generic_n_generators;i++)
1217        {
1218            generic_lt=MonListAdd (generic_lt,cur_gen->lt);
1219            cur_gen=cur_gen->next;
1220        }
1221        for (i=0;i<final_base_dim;i++) memcpy(generic_column_name[i],column_name[i],sizeof(exponent)*variables);
1222     }
1223}
1224
[880433]1225static void modp_SetColumnNames ()  // used by modp - sets column names, the old table will be destroyed
[b5e57e2]1226{
1227    int i;
1228    for (i=0;i<final_base_dim;i++) memcpy(generic_column_name[i],column_name[i],sizeof(exponent)*variables);
1229}
1230
[880433]1231static void CheckColumnSequence () // checks if scheme of computations is as generic one
[b5e57e2]1232{
1233     int i;
1234     if (cur_result->n_generators!=generic_n_generators)
1235     {
1236#ifdef writemsg
[6b4ff12]1237        PrintS("wrong number of generators occured");
[b5e57e2]1238#else
[6b4ff12]1239        if (protocol) PrintS("ng");
[b5e57e2]1240#endif
1241        Discard ();
1242        return;
1243     }
1244     if (denom_divisible)
1245     {
1246#ifdef writemsg
[6b4ff12]1247        PrintS("denom of coef divisible by p");
[b5e57e2]1248#else
[6b4ff12]1249        if (protocol) PrintS("dp");
[b5e57e2]1250#endif
1251        Discard ();
1252        return;
1253     }
1254     generator_entry *cur_gen;
1255     mon_list_entry *cur_mon;
1256     cur_gen=cur_result->generator;
1257     cur_mon=generic_lt;
1258     for (i=0;i<generic_n_generators;i++)
1259     {
1260         if (!EqualMon(cur_mon->mon,cur_gen->lt))
1261         {
1262#ifdef writemsg
[6b4ff12]1263            PrintS("wrong leading term occured");
[b5e57e2]1264#else
[6b4ff12]1265            if (protocol) PrintS("lt");
[b5e57e2]1266#endif
1267            Discard ();
1268            return;
1269         }
1270         cur_gen=cur_gen->next;
1271         cur_mon=cur_mon->next;
1272     }
1273     for (i=0;i<final_base_dim;i++)
1274     {
1275         if (!EqualMon(generic_column_name[i],column_name[i]))
1276         {
1277#ifdef writemsg
[6b4ff12]1278            PrintS("wrong seq of cols occured");
[b5e57e2]1279#else
[6b4ff12]1280            if (protocol) PrintS("sc");
[b5e57e2]1281#endif
1282            Discard ();
1283            return;
1284         }
1285     }
1286     good_primes++;
1287}
[30d574]1288#if 0 /* only debuggig */
[b5e57e2]1289void WriteGenerator () // writes generator (only for debugging)
1290{
1291     char *str;
[30d574]1292     str=(char*)omAlloc0(sizeof(char)*1000);
[b5e57e2]1293     int i;
1294     for (i=0;i<=final_base_dim;i++)
1295     {
1296         str=mpz_get_str(str,10,polycoef[i]);
[6b4ff12]1297         PrintS(str);
1298         PrintS("*");
[b5e57e2]1299         WriteMono(polyexp[i]);
[6b4ff12]1300         PrintS(" ");
[b5e57e2]1301     }
[30d574]1302     omFree(str);
[6b4ff12]1303     PrintLn();
[b5e57e2]1304}
[30d574]1305#endif
[b5e57e2]1306
[880433]1307static bool CheckGenerator () // evaluates generator to check whether it is good
[b5e57e2]1308{
[880433]1309     mpz_t val,sum;
[b5e57e2]1310     int con,i;
1311     mpz_init(val);
1312     mpz_init(sum);
1313     for (con=0;con<final_base_dim;con++)
1314     {
1315       mpz_set_si(sum,0);
1316       for (i=0;i<=final_base_dim;i++)
1317       {
1318         int_Evaluate(val, polyexp[i], condition_list[con]);
1319         mpz_mul(val,val,polycoef[i]);
1320         mpz_add(sum,sum,val);
1321       }
1322       if (mpz_sgn(sum)!=0)
1323       {
1324          mpz_clear(val);
1325          mpz_clear(sum);
1326          return false;
1327       }
1328    }
1329    mpz_clear(val);
1330    mpz_clear(sum);
1331    return true;
1332}
1333
[880433]1334static void ClearGenList ()
[b5e57e2]1335{
1336     gen_list_entry *temp;
1337     int i;
1338     while (gen_list!=NULL)
1339     {
1340         temp=gen_list->next;
1341         for (i=0;i<=final_base_dim;i++)
1342         {
1343             mpz_clear(gen_list->polycoef[i]);
[30d574]1344             omFree(gen_list->polyexp[i]);
[b5e57e2]1345         }
[30d574]1346         omFree(gen_list->polycoef);
1347         omFree(gen_list->polyexp);
1348         omFree(gen_list);
[b5e57e2]1349         gen_list=temp;
1350      }
1351}
1352
[880433]1353static void UpdateGenList ()
[b5e57e2]1354{
1355     gen_list_entry *temp,*prev;
1356     int i,j;
1357     prev=NULL;
1358     temp=gen_list;
1359     exponent deg;
1360     for (i=0;i<=final_base_dim;i++)
1361     {
1362         deg=MonDegree(polyexp[i]);
1363         for (j=0;j<deg;j++)
1364         {
1365             mpz_mul(polycoef[i],polycoef[i],common_denom);
1366         }
1367     }
1368     ClearGCD ();
1369     while (temp!=NULL)
1370     {
1371         prev=temp;
1372         temp=temp->next;
1373     }
[30d574]1374     temp=(gen_list_entry*)omAlloc0(sizeof(gen_list_entry));
[b5e57e2]1375     if (prev==NULL) gen_list=temp; else prev->next=temp;
1376     temp->next=NULL;
[aad4e14]1377     temp->polycoef=(mpz_t*)omAlloc(sizeof(mpz_t)*(final_base_dim+1));
1378     temp->polyexp=(mono_type*)omAlloc(sizeof(mono_type)*(final_base_dim+1));
[b5e57e2]1379     for (i=0;i<=final_base_dim;i++)
1380     {
1381         mpz_init(temp->polycoef[i]);
1382         mpz_set(temp->polycoef[i],polycoef[i]);
1383         temp->polyexp[i]=ZeroMonomial ();
1384         memcpy(temp->polyexp[i],polyexp[i],sizeof(exponent)*variables);
1385     }
1386}
1387
[30d574]1388#if 0 /* only debugging */
[b5e57e2]1389void ShowGenList ()
1390{
1391     gen_list_entry *temp;
1392     int i;
1393     char *str;
[30d574]1394     str=(char*)omAlloc0(sizeof(char)*1000);
[b5e57e2]1395     temp=gen_list;
1396     while (temp!=NULL)
1397     {
[6b4ff12]1398         PrintS("generator: ");
[b5e57e2]1399         for (i=0;i<=final_base_dim;i++)
1400         {
1401             str=mpz_get_str(str,10,temp->polycoef[i]);
[6b4ff12]1402             PrintS(str);
1403             PrintS("*");
[b5e57e2]1404             WriteMono(temp->polyexp[i]);
1405         }
[6b4ff12]1406         PrintLn();
[b5e57e2]1407         temp=temp->next;
1408     }
[30d574]1409     omFree(str);
[b5e57e2]1410}
[30d574]1411#endif
[b5e57e2]1412
1413
[880433]1414static void modp_Main ()
[b5e57e2]1415{
1416     mono_type cur_mon;
1417     cur_mon= ZeroMonomial ();
1418     modp_denom=1;
1419     bool row_is_zero;
1420
[30d574]1421#if 0 /* only debugging */
[b5e57e2]1422     Info ();
1423#endif
1424
1425     while (check_list!=NULL)
1426     {
1427           TakeNextMonomial (cur_mon);
1428           ProduceRow (cur_mon);
[30d574]1429#if 0 /* only debugging */
[b5e57e2]1430           cout << "row produced for monomial ";
1431           WriteMono (cur_mon);
1432           cout << endl;
1433           Info ();
1434#endif
1435           ReduceRow ();
1436           row_is_zero = RowIsZero ();
1437           if (row_is_zero)
1438           {
1439              lt_list=MonListAdd (lt_list,cur_mon);
1440              ReduceCheckListByMon (cur_mon);
1441              NewGenerator (cur_mon);
[30d574]1442#if 0 /* only debugging */
[b5e57e2]1443              cout << "row is zero - linear dependence found (should be seen in my_solve_row)" << endl;
1444              cout << "monomial added to leading terms list" << endl;
1445              cout << "check list updated" << endl;
1446              Info ();
1447#endif
1448           }
1449           else
1450           {
1451              base_list= MonListAdd (base_list,cur_mon);
1452              UpdateCheckList (cur_mon);
1453              ReduceCheckListByLTs ();
[30d574]1454#if 0 /* only debugging */
[b5e57e2]1455              cout << "row is non-zero" << endl;
1456              cout << "monomial added to quotient basis list" << endl;
1457              cout << "new monomials added to check list" << endl;
1458              cout << "check list reduced by monomials from leading term list" << endl;
1459              Info ();
1460#endif
1461              PrepareRow (cur_mon);
[30d574]1462#if 0 /* only debugging */
[b5e57e2]1463              cout << "row prepared and put into matrix" << endl;
1464              Info ();
1465#endif
1466           }
1467        }
[30d574]1468        omFree(cur_mon);
[b5e57e2]1469}
1470
[880433]1471static void ResolveCoeff (mpq_t c, number m)
[b5e57e2]1472{
[880433]1473  if ((long)m & SR_INT)
1474  {
1475    long m_val=SR_TO_INT(m);
1476    mpq_set_si(c,m_val,1);
1477  }
1478  else
1479  {
1480    if (m->s<2)
1481    {
1482      mpz_set(mpq_numref(c),m->z);
1483      mpz_set(mpq_denref(c),m->n);
1484      mpq_canonicalize(c);
1485    }
1486    else
1487    {
1488      mpq_set_z(c,m->z);
1489    }
1490  }
[b5e57e2]1491}
1492
1493ideal interpolation(lists L, intvec *v)
1494{
1495  protocol=TEST_OPT_PROT;  // should be set if option(prot) is enabled
1496
1497  bool data_ok=true;
1498
1499  // reading the ring data ***************************************************
[0519cc]1500  if ((currRing==NULL) || ((!rField_is_Zp (currRing))&&(!rField_is_Q (currRing))))
[b5e57e2]1501  {
1502     WerrorS("coefficient field should be Zp or Q!");
1503     return NULL;
1504  }
1505  if ((currRing->qideal)!=NULL)
1506  {
1507     WerrorS("quotient ring not supported!");
1508     return NULL;
1509  }
1510  if ((currRing->OrdSgn)!=1)
1511  {
1512     WerrorS("ordering must be global!");
1513     return NULL;
1514  }
1515  n_points=v->length ();
1516  if (n_points!=(L->nr+1))
1517  {
1518     WerrorS("list and intvec must have the same length!");
1519     return NULL;
1520  }
1521  variables=currRing->N;
[0519cc]1522  only_modp=rField_is_Zp(currRing);
1523  if (only_modp) myp=rChar(currRing);
[b5e57e2]1524  // ring data read **********************************************************
1525
1526
1527  multiplicity=(int*)malloc(sizeof(int)*n_points);
1528  int i;
1529  for (i=0;i<n_points;i++) multiplicity[i]=(*v)[i];
1530
1531  final_base_dim = CalcBaseDim ();
1532
1533#ifdef writemsg
1534  Print("number of variables: %d\n", variables);
1535  Print("number of points: %d\n", n_points);
[6b4ff12]1536  PrintS("multiplicities: ");
[b5e57e2]1537  for (i=0;i<n_points;i++) Print("%d ", multiplicity[i]);
[6b4ff12]1538  PrintLn();
[b5e57e2]1539  Print("general initialization for dimension %d ...\n", final_base_dim);
1540#endif
1541
1542  GeneralInit ();
1543
1544// reading coordinates of points from ideals **********************************
1545  mpq_t divisor;
1546  if (!only_modp) mpq_init(divisor);
1547  int j;
1548  for(i=0; i<=L->nr;i++)
1549  {
1550    ideal I=(ideal)L->m[i].Data();
1551    for(j=0;j<IDELEMS(I);j++)
1552    {
1553      poly p=I->m[j];
1554      if (p!=NULL)
1555      {
1556        poly ph=pHead(p);
1557        int pcvar=pVar(ph);
1558        if (pcvar!=0)
1559        {
1560          pcvar--;
1561          if (coord_exist[i][pcvar])
1562          {
1563             Print("coordinate %d for point %d initialized twice!\n",pcvar+1,i+1);
1564             data_ok=false;
1565          }
1566          number m;
1567          m=pGetCoeff(p); // possible coefficient standing by a leading monomial
1568          if (!only_modp) ResolveCoeff (divisor,m);
1569          number n;
1570          if (pNext(p)!=NULL) n=pGetCoeff(pNext(p));
1571          else n=nInit(0);
1572          if (only_modp)
1573          {
1574            n=nNeg(n);
1575            n=nDiv(n,m);
1576            modp_points[i][pcvar]=(int)((long)n);
1577          }
1578          else
1579          {
1580            ResolveCoeff (q_points[i][pcvar],n);
1581            mpq_neg(q_points[i][pcvar],q_points[i][pcvar]);
1582            mpq_div(q_points[i][pcvar],q_points[i][pcvar],divisor);
1583          }
1584          coord_exist[i][pcvar]=true;
1585        }
1586        else
1587        {
[6b4ff12]1588          PrintS("not a variable? ");
[b5e57e2]1589          wrp(p);
1590          PrintLn();
1591          data_ok=false;
1592        }
1593        pDelete(&ph);
1594      }
1595    }
1596  }
1597  if (!only_modp) mpq_clear(divisor);
1598  // data from ideal read *******************************************************
1599
1600  // ckecking if all coordinates are initialized
1601  for (i=0;i<n_points;i++)
1602  {
1603      for (j=0;j<variables;j++)
1604      {
1605          if (!coord_exist[i][j])
1606          {
1607             Print("coordinate %d for point %d not known!\n",j+1,i+1);
1608             data_ok=false;
1609          }
1610      }
1611  }
1612
1613  if (!data_ok)
1614  {
[f4fea85]1615     GeneralDone();
[b5e57e2]1616     WerrorS("data structure is invalid");
1617     return NULL;
1618  }
1619
1620  if (!only_modp) IntegerPoints ();
1621  MakeConditions ();
1622#ifdef writemsg
[6b4ff12]1623  PrintS("done.\n");
[b5e57e2]1624#else
1625  if (protocol) Print("[vdim %d]",final_base_dim);
1626#endif
1627
1628
1629// main procedure *********************************************************************
1630  int modp_cycles=10;
1631  bool correct_gen=false;
1632  if (only_modp) modp_cycles=1;
[c16aa0]1633  #ifdef HAVE_FACTORY
[b5e57e2]1634  myp_index=cf_getNumSmallPrimes ();
[c16aa0]1635  #endif
[b5e57e2]1636
1637  while ((!correct_gen)&&(myp_index>1))
1638  {
1639#ifdef writemsg
1640        Print("trying %d cycles mod p...\n",modp_cycles);
1641#else
1642        if (protocol) Print("(%d)",modp_cycles);
1643#endif
1644        while ((n_results<modp_cycles)&&(myp_index>1))  // some computations mod p
1645        {
1646            if (!only_modp) myp=TakePrime (myp);
1647            NewResultEntry ();
1648            InitProcData ();
1649            if (only_modp) modp_PrepareProducts (); else int_PrepareProducts ();
1650
1651            modp_Main ();
1652
1653            if (!only_modp)
1654            {
1655               MultGenerators ();
1656               CheckColumnSequence ();
1657            }
1658            else
1659            {
1660               modp_SetColumnNames ();
1661            }
1662            FreeProcData ();
1663        }
1664
1665        if (!only_modp)
1666        {
1667           PrepareChinese (modp_cycles);
1668           correct_gen=true;
1669           for (i=0;i<generic_n_generators;i++)
1670           {
[880433]1671               ReconstructGenerator (i,modp_cycles);
[b5e57e2]1672               correct_gen=CheckGenerator ();
1673               if (!correct_gen)
1674               {
1675#ifdef writemsg
[6b4ff12]1676                  PrintS("wrong generator!\n");
[b5e57e2]1677#else
[6b4ff12]1678//                  if (protocol) PrintS("!g");
[b5e57e2]1679#endif
1680                  ClearGenList ();
1681                  break;
1682               }
1683               else
1684               {
1685                  UpdateGenList ();
1686               }
1687           }
1688#ifdef checksize
1689           Print("maximal size of output: %d, precision bound: %d.\n",maximalsize,mpz_sizeinbase(bigcongr,10));
1690#endif
[880433]1691           CloseChinese ();
[b5e57e2]1692           modp_cycles=modp_cycles+10;
1693        }
1694        else
1695        {
1696           correct_gen=true;
1697        }
1698  }
1699// end of main procedure ************************************************************************************
1700
1701#ifdef writemsg
[6b4ff12]1702  PrintS("computations finished.\n");
[b5e57e2]1703#else
[6b4ff12]1704  if (protocol) PrintLn();
[b5e57e2]1705#endif
1706
1707  if (!correct_gen)
1708  {
1709     GeneralDone ();
1710     ClearGenList ();
1711     WerrorS("internal error - coefficient too big!");
1712     return NULL;
1713  }
1714
1715// passing data to ideal *************************************************************************************
1716  ideal ret;
1717
1718  if (only_modp)
1719  {
1720    mono_type mon;
1721    ret=idInit(modp_result->n_generators,1);
1722    generator_entry *cur_gen=modp_result->generator;
1723    for(i=0;i<IDELEMS(ret);i++)
1724    {
1725      poly p,sum;
1726      sum=NULL;
1727      int a;
1728      int cf;
1729      for (a=final_base_dim;a>=0;a--)
1730      {
1731        if (a==final_base_dim) cf=cur_gen->ltcoef; else cf=cur_gen->coef[a];
1732        if (cf!=0)
1733        {
1734            p=pISet(cf);
1735            if (a==final_base_dim) mon=cur_gen->lt; else mon=generic_column_name[a];
1736            for (j=0;j<variables;j++) pSetExp(p,j+1,mon[j]);
1737            pSetm(p);
1738            sum=pAdd(sum,p);
1739        }
1740      }
1741      ret->m[i]=sum;
1742      cur_gen=cur_gen->next;
1743    }
1744  }
1745  else
1746  {
1747    ret=idInit(generic_n_generators,1);
1748    gen_list_entry *temp=gen_list;
1749    for(i=0;i<IDELEMS(ret);i++)
1750    {
1751      poly p,sum;
1752      sum=NULL;
1753      int a;
1754      for (a=final_base_dim;a>=0;a--) // build one polynomial
1755      {
1756          if (mpz_sgn(temp->polycoef[a])!=0)
1757          {
[896561]1758             number n=ALLOC_RNUMBER();
[b5e57e2]1759#ifdef LDEBUG
1760             n->debug=123456;
1761#endif
[0f7301]1762             mpz_init_set(n->z,temp->polycoef[a]);
[b5e57e2]1763             n->s=3;
[0519cc]1764             nlNormalize(n, currRing->cf);
[3eabdb]1765             p=pNSet(n); //a monomial
1766             for (j=0;j<variables;j++) pSetExp(p,j+1,temp->polyexp[a][j]);
1767             pSetm(p); // after all pSetExp
[b5e57e2]1768             sum=pAdd(sum,p);
1769          }
1770      }
1771      ret->m[i]=sum;
1772      temp=temp->next;
1773    }
1774  }
1775// data transferred ****************************************************************************
1776
1777
1778  GeneralDone ();
1779  ClearGenList ();
1780  return ret;
1781}
1782
1783
1784BOOLEAN jjINTERPOLATION (leftv res, leftv l, leftv v)
1785{
1786  res->data=interpolation((lists)l->Data(),(intvec*)v->Data());
1787  setFlag(res,FLAG_STD);
1788  return errorreported;
1789}
Note: See TracBrowser for help on using the repository browser.