[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] | 39 | static modp_number myp; // all modp computation done mod myp |
---|
| 40 | static int myp_index; // index of small prime in Singular table of primes |
---|
[b5e57e2] | 41 | |
---|
| 42 | static 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 | } |
---|
| 47 | static 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 | |
---|
| 53 | typedef exponent *mono_type; |
---|
| 54 | typedef struct {mono_type mon; unsigned int point_ref;} condition_type; |
---|
| 55 | typedef modp_number *coordinate_products; |
---|
| 56 | typedef coordinate_products *coordinates; |
---|
| 57 | |
---|
| 58 | struct mon_list_entry_struct {mono_type mon; |
---|
| 59 | struct mon_list_entry_struct *next;}; |
---|
| 60 | typedef struct mon_list_entry_struct mon_list_entry; |
---|
| 61 | |
---|
| 62 | struct 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 | |
---|
| 67 | typedef struct row_list_entry_struct row_list_entry; |
---|
| 68 | |
---|
| 69 | struct generator_struct {modp_number *coef; |
---|
| 70 | mono_type lt; |
---|
| 71 | modp_number ltcoef; |
---|
| 72 | struct generator_struct *next;}; |
---|
| 73 | |
---|
| 74 | typedef struct generator_struct generator_entry; |
---|
| 75 | |
---|
| 76 | struct 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 | |
---|
| 82 | typedef struct modp_result_struct modp_result_entry; |
---|
| 83 | |
---|
| 84 | typedef modp_number *modp_coordinates; |
---|
| 85 | typedef mpq_t *q_coordinates; |
---|
| 86 | typedef mpz_t *int_coordinates; |
---|
| 87 | typedef bool *coord_exist_table; |
---|
| 88 | |
---|
[880433] | 89 | static int final_base_dim; // dimension of the quotient space, known from the beginning |
---|
| 90 | static int last_solve_column; // last non-zero column in "solve" part of matrix, used for speed up |
---|
| 91 | static int n_points; // modp_number of ideals (points) |
---|
| 92 | static int *multiplicity; // multiplicities of points |
---|
| 93 | static int variables; // modp_number of variables |
---|
| 94 | static int max_coord; // maximal possible coordinate product used during Evaluation |
---|
| 95 | static bool only_modp; // perform only one modp computations |
---|
| 96 | |
---|
| 97 | static modp_coordinates *modp_points; // coordinates of points for modp problem - used by Evaluate (this is also initial data for only modp) |
---|
| 98 | static q_coordinates *q_points; // coordinates of points for rational data (not used for modp) |
---|
| 99 | static int_coordinates *int_points; // coordinates of points for integer data - used to check generators (not used for modp) |
---|
| 100 | static coord_exist_table *coord_exist; // checks whether all coordinates has been initialized |
---|
| 101 | static mon_list_entry *check_list; // monomials to be checked in next stages |
---|
| 102 | static coordinates *points; // power products of coordinates of points used in modp cycles |
---|
| 103 | static condition_type *condition_list; // conditions stored in an array |
---|
| 104 | static mon_list_entry *lt_list; // leading terms found so far |
---|
| 105 | static mon_list_entry *base_list; // standard monomials found so far |
---|
| 106 | static row_list_entry *row_list; // rows of the matrix (both parts) |
---|
| 107 | static modp_number *my_row; // one special row to perform operations |
---|
| 108 | static modp_number *my_solve_row; // one special row to find the linear dependence ("solve" part) |
---|
| 109 | static mono_type *column_name; // monomials assigned to columns in solve_row |
---|
| 110 | |
---|
| 111 | static int n_results; // modp_number of performed modp computations (not discarded) |
---|
| 112 | static modp_number modp_denom; // denominator of mod p computations |
---|
| 113 | static modp_result_entry *modp_result; // list of results for various mod p calculations (used for modp - first result is the desired one) |
---|
| 114 | static modp_result_entry *cur_result; // pointer to current result (as before) |
---|
| 115 | static modp_number *congr; // primes used in computations (chinese remainder theorem) (not used for modp) |
---|
| 116 | static modp_number *in_gamma; // inverts used in chinese remainder theorem (not used for modp) |
---|
| 117 | static mpz_t bigcongr; // result, in fact, is given in mod bigcongr (not used for modp) |
---|
| 118 | |
---|
| 119 | static mpz_t *polycoef; // polynomial integercoefficients (not used for modp) |
---|
| 120 | static mono_type *polyexp; // polynomial exponents |
---|
[b5e57e2] | 121 | |
---|
| 122 | struct gen_list_struct {mpz_t *polycoef; |
---|
| 123 | mono_type *polyexp; |
---|
| 124 | struct gen_list_struct *next;}; |
---|
| 125 | typedef struct gen_list_struct gen_list_entry; |
---|
| 126 | |
---|
[880433] | 127 | static gen_list_entry *gen_list=NULL; // list of resulting generators - output data (integer version) |
---|
[b5e57e2] | 128 | |
---|
[880433] | 129 | static int generic_n_generators; // modp_number of generators - should be the same for all modp comp (not used for modp) |
---|
| 130 | static mono_type *generic_column_name; // monomials assigned to columns in solve_row - should be the same for all modp comp (!!! used for modp) |
---|
| 131 | static mon_list_entry *generic_lt=NULL; // leading terms for ordered generators - should be the same for all modp comp (not used for modp) |
---|
| 132 | static int good_primes; // modp_number of good primes so far; |
---|
| 133 | static int bad_primes; // modp_number of bad primes so far; |
---|
| 134 | static mpz_t common_denom; // common denominator used to force points coordinates to Z (not used for modp) |
---|
| 135 | static bool denom_divisible; // common denominator is divisible by p (not used for modp) |
---|
[b5e57e2] | 136 | |
---|
[880433] | 137 | static poly comparizon_p1; //polynomials used to do comparizons by Singular |
---|
| 138 | static poly comparizon_p2; |
---|
[b5e57e2] | 139 | |
---|
[880433] | 140 | static modp_number *modp_Reverse; // reverses in mod p |
---|
[b5e57e2] | 141 | |
---|
[880433] | 142 | static bool protocol; // true to show the protocol |
---|
[b5e57e2] | 143 | |
---|
| 144 | #ifdef checksize |
---|
[880433] | 145 | static int maximal_size=0; |
---|
[b5e57e2] | 146 | #endif |
---|
| 147 | |
---|
[30d574] | 148 | #if 0 /* only for debuggig*/ |
---|
[b5e57e2] | 149 | void 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 | |
---|
| 156 | void 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 | |
---|
| 167 | void 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] | 227 | static 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] | 256 | static 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] | 275 | static 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] | 283 | static 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] | 291 | static 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] | 304 | static 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] | 329 | static 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] | 337 | static 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] | 344 | static 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] | 405 | static 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] | 441 | static 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] | 455 | static 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] | 506 | static 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] | 533 | static 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] | 556 | static 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] | 583 | static 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] | 604 | static 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] | 632 | static 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] | 656 | static 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] | 670 | static 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] | 701 | static 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] | 777 | static 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] | 791 | static 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] | 799 | static 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] | 827 | static 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] | 840 | static 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] | 851 | static 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] | 862 | static 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] | 889 | static 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] | 942 | static 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] | 964 | static 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] | 981 | static 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] | 1005 | static 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] | 1021 | void 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] | 1043 | static 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] | 1053 | static 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] | 1080 | static void CloseChinese () // after CRA |
---|
[b5e57e2] | 1081 | { |
---|
[30d574] | 1082 | omFree(in_gamma); |
---|
| 1083 | omFree(congr); |
---|
[b5e57e2] | 1084 | mpz_clear(bigcongr); |
---|
| 1085 | } |
---|
| 1086 | |
---|
[880433] | 1087 | static 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] | 1110 | static 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] | 1175 | static 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] | 1225 | static 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] | 1231 | static 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] | 1289 | void 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] | 1307 | static 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] | 1334 | static 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] | 1353 | static 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] | 1389 | void 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] | 1414 | static 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] | 1471 | static 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 | |
---|
| 1493 | ideal 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 | |
---|
| 1784 | BOOLEAN 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 | } |
---|