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