[0bc2dd] | 1 | /**************************************** |
---|
| 2 | * Computer Algebra System SINGULAR * |
---|
| 3 | ****************************************/ |
---|
[073295f] | 4 | /* $Id: interpolation.cc,v 1.7 2009-08-13 15:30:51 Singular Exp $ */ |
---|
[0bc2dd] | 5 | |
---|
[b5e57e2] | 6 | #include "mod2.h" |
---|
| 7 | #include "structs.h" |
---|
| 8 | #include "febase.h" |
---|
| 9 | #include "ideals.h" |
---|
| 10 | #include "intvec.h" |
---|
| 11 | #include "polys.h" |
---|
| 12 | #include "lists.h" |
---|
| 13 | #include "longrat.h" |
---|
| 14 | #include "ipid.h" |
---|
[0bc2dd] | 15 | #include "ring.h" |
---|
[b5e57e2] | 16 | #include <factory.h> |
---|
| 17 | |
---|
| 18 | //memory management |
---|
| 19 | #define mdmALLOC(x) omAlloc0(x) |
---|
| 20 | #define mdmFREE(x) omFree(x) |
---|
| 21 | |
---|
| 22 | // parameters to debug |
---|
| 23 | //#define debb |
---|
| 24 | //#define shmat |
---|
| 25 | //#define checksize |
---|
| 26 | //#define writemsg |
---|
| 27 | |
---|
| 28 | // possible strategies |
---|
| 29 | #define unsortedmatrix |
---|
| 30 | //#define integerstrategy |
---|
| 31 | |
---|
| 32 | #define modp_number int |
---|
| 33 | #define exponent int |
---|
| 34 | |
---|
| 35 | modp_number myp; // all modp computation done mod myp |
---|
| 36 | int myp_index; // index of small prime in Singular table of primes |
---|
| 37 | |
---|
| 38 | static inline modp_number modp_mul (modp_number x,modp_number y) |
---|
| 39 | { |
---|
| 40 | // return (x*y)%myp; |
---|
| 41 | return (modp_number)(((unsigned long)x)*((unsigned long)y)%((unsigned long)myp)); |
---|
| 42 | } |
---|
| 43 | static inline modp_number modp_sub (modp_number x,modp_number y) |
---|
| 44 | { |
---|
| 45 | // if (x>=y) return x-y; else return x+myp-y; |
---|
| 46 | return (x+myp-y)%myp; |
---|
| 47 | } |
---|
| 48 | |
---|
| 49 | typedef exponent *mono_type; |
---|
| 50 | typedef struct {mono_type mon; unsigned int point_ref;} condition_type; |
---|
| 51 | typedef modp_number *coordinate_products; |
---|
| 52 | typedef coordinate_products *coordinates; |
---|
| 53 | |
---|
| 54 | struct mon_list_entry_struct {mono_type mon; |
---|
| 55 | struct mon_list_entry_struct *next;}; |
---|
| 56 | typedef struct mon_list_entry_struct mon_list_entry; |
---|
| 57 | |
---|
| 58 | struct row_list_entry_struct {modp_number *row_matrix; |
---|
| 59 | modp_number *row_solve; |
---|
| 60 | int first_col; |
---|
| 61 | struct row_list_entry_struct *next;}; |
---|
| 62 | |
---|
| 63 | typedef struct row_list_entry_struct row_list_entry; |
---|
| 64 | |
---|
| 65 | struct generator_struct {modp_number *coef; |
---|
| 66 | mono_type lt; |
---|
| 67 | modp_number ltcoef; |
---|
| 68 | struct generator_struct *next;}; |
---|
| 69 | |
---|
| 70 | typedef struct generator_struct generator_entry; |
---|
| 71 | |
---|
| 72 | struct modp_result_struct {modp_number p; |
---|
| 73 | generator_entry *generator; |
---|
| 74 | int n_generators; |
---|
| 75 | struct modp_result_struct *next; |
---|
| 76 | struct modp_result_struct *prev;}; |
---|
| 77 | |
---|
| 78 | typedef struct modp_result_struct modp_result_entry; |
---|
| 79 | |
---|
| 80 | typedef modp_number *modp_coordinates; |
---|
| 81 | typedef mpq_t *q_coordinates; |
---|
| 82 | typedef mpz_t *int_coordinates; |
---|
| 83 | typedef bool *coord_exist_table; |
---|
| 84 | |
---|
| 85 | int final_base_dim; // dimension of the quotient space, known from the beginning |
---|
| 86 | int last_solve_column; // last non-zero column in "solve" part of matrix, used for speed up |
---|
| 87 | int n_points; // modp_number of ideals (points) |
---|
| 88 | int *multiplicity; // multiplicities of points |
---|
| 89 | int variables; // modp_number of variables |
---|
| 90 | int max_coord; // maximal possible coordinate product used during Evaluation |
---|
| 91 | bool only_modp; // perform only one modp computations |
---|
| 92 | |
---|
| 93 | modp_coordinates *modp_points; // coordinates of points for modp problem - used by Evaluate (this is also initial data for only modp) |
---|
| 94 | q_coordinates *q_points; // coordinates of points for rational data (not used for modp) |
---|
| 95 | int_coordinates *int_points; // coordinates of points for integer data - used to check generators (not used for modp) |
---|
| 96 | coord_exist_table *coord_exist; // checks whether all coordinates has been initialized |
---|
| 97 | mon_list_entry *check_list; // monomials to be checked in next stages |
---|
| 98 | coordinates *points; // power products of coordinates of points used in modp cycles |
---|
| 99 | condition_type *condition_list; // conditions stored in an array |
---|
| 100 | mon_list_entry *lt_list; // leading terms found so far |
---|
| 101 | mon_list_entry *base_list; // standard monomials found so far |
---|
| 102 | row_list_entry *row_list; // rows of the matrix (both parts) |
---|
| 103 | modp_number *my_row; // one special row to perform operations |
---|
| 104 | modp_number *my_solve_row; // one special row to find the linear dependence ("solve" part) |
---|
| 105 | mono_type *column_name; // monomials assigned to columns in solve_row |
---|
| 106 | |
---|
| 107 | int n_results; // modp_number of performed modp computations (not discarded) |
---|
| 108 | modp_number modp_denom; // denominator of mod p computations |
---|
| 109 | modp_result_entry *modp_result; // list of results for various mod p calculations (used for modp - first result is the desired one) |
---|
| 110 | modp_result_entry *cur_result; // pointer to current result (as before) |
---|
| 111 | modp_number *congr; // primes used in computations (chinese remainder theorem) (not used for modp) |
---|
[1f18bb] | 112 | modp_number *in_gamma; // inverts used in chinese remainder theorem (not used for modp) |
---|
[b5e57e2] | 113 | mpz_t bigcongr; // result, in fact, is given in mod bigcongr (not used for modp) |
---|
| 114 | |
---|
| 115 | mpz_t *polycoef; // polynomial integercoefficients (not used for modp) |
---|
| 116 | mono_type *polyexp; // polynomial exponents |
---|
| 117 | |
---|
| 118 | struct gen_list_struct {mpz_t *polycoef; |
---|
| 119 | mono_type *polyexp; |
---|
| 120 | struct gen_list_struct *next;}; |
---|
| 121 | typedef struct gen_list_struct gen_list_entry; |
---|
| 122 | |
---|
| 123 | gen_list_entry *gen_list=NULL; // list of resulting generators - output data (integer version) |
---|
| 124 | |
---|
| 125 | int generic_n_generators; // modp_number of generators - should be the same for all modp comp (not used for modp) |
---|
| 126 | mono_type *generic_column_name; // monomials assigned to columns in solve_row - should be the same for all modp comp (!!! used for modp) |
---|
| 127 | mon_list_entry *generic_lt=NULL; // leading terms for ordered generators - should be the same for all modp comp (not used for modp) |
---|
| 128 | int good_primes; // modp_number of good primes so far; |
---|
| 129 | int bad_primes; // modp_number of bad primes so far; |
---|
| 130 | mpz_t common_denom; // common denominator used to force points coordinates to Z (not used for modp) |
---|
| 131 | bool denom_divisible; // common denominator is divisible by p (not used for modp) |
---|
| 132 | |
---|
| 133 | poly comparizon_p1; //polynomials used to do comparizons by Singular |
---|
| 134 | poly comparizon_p2; |
---|
| 135 | |
---|
| 136 | modp_number *modp_Reverse; // reverses in mod p |
---|
| 137 | |
---|
| 138 | bool protocol; // true to show the protocol |
---|
| 139 | |
---|
| 140 | #ifdef checksize |
---|
| 141 | int maximal_size=0; |
---|
| 142 | #endif |
---|
| 143 | |
---|
| 144 | void WriteMono (mono_type m) // Writes a monomial on the screen - only for debug |
---|
| 145 | { |
---|
| 146 | int i; |
---|
| 147 | for (i=0;i<variables;i++) Print("_%d", m[i]); |
---|
[6b4ff12] | 148 | PrintS(" "); |
---|
[b5e57e2] | 149 | } |
---|
| 150 | |
---|
| 151 | #ifdef debb |
---|
| 152 | void WriteMonoList (mon_list_entry *list) |
---|
| 153 | { |
---|
| 154 | mon_list_entry *ptr; |
---|
| 155 | ptr=list; |
---|
| 156 | while (ptr!=NULL) |
---|
| 157 | { |
---|
| 158 | WriteMono(ptr->mon); |
---|
| 159 | ptr=ptr->next; |
---|
| 160 | } |
---|
| 161 | } |
---|
| 162 | |
---|
| 163 | void Info () |
---|
| 164 | { |
---|
| 165 | int i; |
---|
| 166 | row_list_entry *curptr; |
---|
| 167 | modp_number *row; |
---|
| 168 | modp_number *solve; |
---|
| 169 | char ch; |
---|
| 170 | cout << endl << "quotient basis: "; |
---|
| 171 | WriteMonoList (base_list); |
---|
| 172 | cout << endl << "leading terms: "; |
---|
| 173 | WriteMonoList (lt_list); |
---|
| 174 | cout << endl << "to be checked: "; |
---|
| 175 | WriteMonoList (check_list); |
---|
| 176 | #ifdef shmat |
---|
| 177 | cout << endl << "Matrix:" << endl; |
---|
| 178 | cout << " "; |
---|
| 179 | for (i=0;i<final_base_dim;i++) |
---|
| 180 | { |
---|
| 181 | WriteMono (column_name[i]); |
---|
| 182 | } |
---|
| 183 | cout << endl; |
---|
| 184 | curptr=row_list; |
---|
| 185 | while (curptr!=NULL) |
---|
| 186 | { |
---|
| 187 | row=curptr->row_matrix; |
---|
| 188 | solve=curptr->row_solve; |
---|
| 189 | for (i=0;i<final_base_dim;i++) |
---|
| 190 | { |
---|
| 191 | cout << *row << " , "; |
---|
| 192 | row++; |
---|
| 193 | } |
---|
| 194 | cout << " "; |
---|
| 195 | for (i=0;i<final_base_dim;i++) |
---|
| 196 | { |
---|
| 197 | cout << *solve << " , "; |
---|
| 198 | solve++; |
---|
| 199 | } |
---|
| 200 | curptr=curptr->next; |
---|
| 201 | cout << endl;} |
---|
| 202 | cout << "Special row: Solve row:" << endl; |
---|
| 203 | row=my_row; |
---|
| 204 | for (i=0;i<final_base_dim;i++) |
---|
| 205 | { |
---|
| 206 | cout << *row << " , "; |
---|
| 207 | row++; |
---|
| 208 | } |
---|
| 209 | cout << " "; |
---|
| 210 | row=my_solve_row; |
---|
| 211 | for (i=0;i<final_base_dim;i++) |
---|
| 212 | { |
---|
| 213 | cout << *row << " , "; |
---|
| 214 | row++; |
---|
| 215 | } |
---|
| 216 | #endif |
---|
| 217 | cout << endl; |
---|
| 218 | cin >> ch; |
---|
| 219 | cout << endl << endl; |
---|
| 220 | } |
---|
| 221 | #endif |
---|
| 222 | |
---|
| 223 | modp_number OneInverse(modp_number a,modp_number p) // computes inverse of d mod p without using tables |
---|
| 224 | { |
---|
| 225 | long u, v, u0, u1, u2, q, r; |
---|
| 226 | u1=1; u2=0; |
---|
| 227 | u = a; v = p; |
---|
| 228 | while (v != 0) { |
---|
| 229 | q = u / v; |
---|
| 230 | r = u % v; |
---|
| 231 | u = v; |
---|
| 232 | v = r; |
---|
| 233 | u0 = u2; |
---|
| 234 | u2 = u1 - q*u2; |
---|
| 235 | u1 = u0; |
---|
| 236 | } |
---|
| 237 | if (u1 < 0) u1=u1+p; |
---|
| 238 | // now it should be ok, but for any case... |
---|
| 239 | if ((u1<0)||((u1*a)%p!=1)) |
---|
| 240 | { |
---|
[6b4ff12] | 241 | PrintS("?"); |
---|
[b5e57e2] | 242 | modp_number i; |
---|
| 243 | for (i=1;i<p;i++) |
---|
| 244 | { |
---|
| 245 | if ((a*i)%p==1) return i; |
---|
| 246 | } |
---|
| 247 | } |
---|
| 248 | return (modp_number)u1; |
---|
| 249 | } |
---|
| 250 | |
---|
| 251 | int CalcBaseDim () // returns the maximal (expected) dimension of quotient space |
---|
| 252 | { |
---|
| 253 | int c; |
---|
| 254 | int i,j; |
---|
| 255 | int s=0; |
---|
| 256 | max_coord=1; |
---|
| 257 | for (i=0;i<n_points;i++) max_coord=max_coord+multiplicity[i]; |
---|
| 258 | for (j=0;j<n_points;j++) |
---|
| 259 | { |
---|
| 260 | c=1; |
---|
| 261 | for (i=1;i<=variables;i++) |
---|
| 262 | { |
---|
| 263 | c=(c*(multiplicity[j]+i-1))/i; |
---|
| 264 | } |
---|
| 265 | s=s+c; |
---|
| 266 | } |
---|
| 267 | return s; |
---|
| 268 | } |
---|
| 269 | |
---|
| 270 | bool EqualMon (mono_type m1, mono_type m2) // compares two monomials, true if equal, false otherwise |
---|
| 271 | { |
---|
| 272 | int i; |
---|
| 273 | for (i=0;i<variables;i++) |
---|
| 274 | if (m1[i]!=m2[i]) return false;; |
---|
| 275 | return true; |
---|
| 276 | } |
---|
| 277 | |
---|
| 278 | exponent MonDegree (mono_type mon) // computes the degree of a monomial |
---|
| 279 | { |
---|
| 280 | exponent v=0; |
---|
| 281 | int i; |
---|
| 282 | for (i=0;i<variables;i++) v=v+mon[i]; |
---|
| 283 | return v; |
---|
| 284 | } |
---|
| 285 | |
---|
| 286 | bool Greater (mono_type m1, mono_type m2) // true if m1 > m2, false otherwise. uses Singular comparing function |
---|
| 287 | { |
---|
| 288 | for (int j = variables; j; j--) |
---|
| 289 | { |
---|
| 290 | pSetExp(comparizon_p1, j, m1[j-1]); |
---|
| 291 | pSetExp(comparizon_p2, j, m2[j-1]); |
---|
| 292 | } |
---|
| 293 | pSetm(comparizon_p1); |
---|
| 294 | pSetm(comparizon_p2); |
---|
| 295 | bool res=(pLmCmp(comparizon_p1,comparizon_p2)>0); |
---|
| 296 | return res; |
---|
| 297 | } |
---|
| 298 | |
---|
| 299 | mon_list_entry* MonListAdd (mon_list_entry *list, mono_type mon) // adds one entry to the list of monomials structure |
---|
| 300 | { |
---|
| 301 | mon_list_entry *curptr=list; |
---|
| 302 | mon_list_entry *prevptr=NULL; |
---|
| 303 | mon_list_entry *temp; |
---|
| 304 | |
---|
| 305 | while (curptr!=NULL) |
---|
| 306 | { |
---|
| 307 | if (EqualMon (mon,(*curptr).mon)) return list; |
---|
| 308 | if (Greater ((*curptr).mon,mon)) break; |
---|
| 309 | prevptr=curptr; |
---|
| 310 | curptr=curptr->next;} |
---|
| 311 | temp=(mon_list_entry*)mdmALLOC(sizeof(mon_list_entry)); |
---|
| 312 | (*temp).next=curptr; |
---|
| 313 | (*temp).mon=(exponent*)mdmALLOC(sizeof(exponent)*variables); |
---|
| 314 | memcpy(temp->mon,mon,sizeof(exponent)*variables); |
---|
| 315 | if (prevptr==NULL) return temp; else { |
---|
| 316 | prevptr->next=temp; |
---|
| 317 | return list;} |
---|
| 318 | } |
---|
| 319 | |
---|
| 320 | mono_type MonListElement (mon_list_entry *list, int n) // returns ith element from list of monomials - no error checking! |
---|
| 321 | { |
---|
| 322 | mon_list_entry *cur=list; |
---|
| 323 | int i; |
---|
| 324 | for (i=0;i<n;i++) cur=cur->next; |
---|
| 325 | return cur->mon; |
---|
| 326 | } |
---|
| 327 | |
---|
| 328 | mono_type ZeroMonomial () // allocates memory for new monomial with all enries equal 0 |
---|
| 329 | { |
---|
| 330 | mono_type m; |
---|
| 331 | m=(exponent*)mdmALLOC(sizeof(exponent)*variables); |
---|
| 332 | int i; |
---|
| 333 | for (i=0;i<variables;i++) m[i]=0; |
---|
| 334 | return m; |
---|
| 335 | } |
---|
| 336 | |
---|
| 337 | void GeneralInit () // general initialization |
---|
| 338 | { |
---|
| 339 | int i,j,k; |
---|
| 340 | points=(coordinates*)mdmALLOC(sizeof(coordinates)*n_points); |
---|
| 341 | for (i=0;i<n_points;i++) |
---|
| 342 | { |
---|
| 343 | points[i]=(coordinate_products*)mdmALLOC(sizeof(coordinate_products)*variables); |
---|
| 344 | for (j=0;j<variables;j++) points[i][j]=(modp_number*)mdmALLOC(sizeof(modp_number)*(max_coord)); |
---|
| 345 | } |
---|
| 346 | condition_list=(condition_type*)mdmALLOC(sizeof(condition_type)*final_base_dim); |
---|
| 347 | for (i=0;i<final_base_dim;i++) condition_list[i].mon=(exponent*)mdmALLOC(sizeof(exponent)*variables); |
---|
| 348 | modp_points=(modp_coordinates*)mdmALLOC(sizeof(modp_coordinates)*n_points); |
---|
| 349 | for (i=0;i<n_points;i++) modp_points[i]=(modp_number*)mdmALLOC(sizeof(modp_number)*variables); |
---|
| 350 | if (!only_modp) |
---|
| 351 | { |
---|
| 352 | q_points=(q_coordinates*)mdmALLOC(sizeof(q_coordinates)*n_points); |
---|
| 353 | for (i=0;i<n_points;i++) |
---|
| 354 | { |
---|
| 355 | q_points[i]=(mpq_t*)mdmALLOC(sizeof(mpq_t)*variables); |
---|
| 356 | for (j=0;j<variables;j++) mpq_init(q_points[i][j]); |
---|
| 357 | } |
---|
| 358 | int_points=(int_coordinates*)mdmALLOC(sizeof(int_coordinates)*n_points); |
---|
| 359 | for (i=0;i<n_points;i++) |
---|
| 360 | { |
---|
| 361 | int_points[i]=(mpz_t*)mdmALLOC(sizeof(mpz_t)*variables); |
---|
| 362 | for (j=0;j<variables;j++) mpz_init(int_points[i][j]); |
---|
| 363 | } |
---|
| 364 | } |
---|
| 365 | coord_exist=(coord_exist_table*)mdmALLOC(sizeof(coord_exist_table)*n_points); |
---|
| 366 | for (i=0;i<n_points;i++) |
---|
| 367 | { |
---|
| 368 | coord_exist[i]=(bool*)mdmALLOC(sizeof(bool)*variables); |
---|
| 369 | for (j=0;j<variables;j++) coord_exist[i][j]=false; |
---|
| 370 | } |
---|
| 371 | generic_column_name=(mono_type*)mdmALLOC(sizeof(mono_type)*final_base_dim); |
---|
| 372 | for (i=0;i<final_base_dim;i++) generic_column_name[i]=ZeroMonomial (); |
---|
| 373 | good_primes=0; |
---|
| 374 | bad_primes=1; |
---|
| 375 | generic_n_generators=0; |
---|
| 376 | if (!only_modp) |
---|
| 377 | { |
---|
| 378 | polycoef=(mpz_t*)mdmALLOC(sizeof(mpz_t)*(final_base_dim+1)); |
---|
| 379 | polyexp=(mono_type*)mdmALLOC(sizeof(mono_type)*(final_base_dim+1)); |
---|
| 380 | for (i=0;i<=final_base_dim;i++) |
---|
| 381 | { |
---|
| 382 | mpz_init(polycoef[i]); |
---|
| 383 | polyexp[i]=ZeroMonomial (); |
---|
| 384 | } |
---|
| 385 | mpz_init(common_denom); |
---|
| 386 | } |
---|
| 387 | |
---|
| 388 | // set all globally used lists to be initially empty |
---|
| 389 | modp_result=NULL; |
---|
| 390 | cur_result=NULL; |
---|
| 391 | gen_list=NULL; |
---|
| 392 | n_results=0; |
---|
| 393 | |
---|
| 394 | // creates polynomials to compare by Singular |
---|
| 395 | comparizon_p1=pOne(); |
---|
| 396 | comparizon_p2=pOne(); |
---|
| 397 | } |
---|
| 398 | |
---|
| 399 | void InitProcData () // initialization for procedures which do computations mod p |
---|
| 400 | { |
---|
| 401 | int i,j; |
---|
| 402 | check_list=MonListAdd (check_list,ZeroMonomial ()); |
---|
| 403 | my_row=(modp_number*)mdmALLOC(sizeof(modp_number)*final_base_dim); |
---|
| 404 | my_solve_row=(modp_number*)mdmALLOC(sizeof(modp_number)*final_base_dim); |
---|
| 405 | column_name=(mono_type*)mdmALLOC(sizeof(mono_type)*final_base_dim); |
---|
| 406 | for (i=0;i<final_base_dim;i++) column_name[i]=ZeroMonomial (); |
---|
| 407 | last_solve_column=0; |
---|
| 408 | modp_number *gen_table; |
---|
| 409 | modp_number pos_gen; |
---|
| 410 | bool gen_ok; |
---|
| 411 | modp_Reverse=(modp_number*)mdmALLOC(sizeof(modp_number)*myp); |
---|
| 412 | |
---|
| 413 | // produces table of modp inverts by finding a generator of (Z_myp*,*) |
---|
| 414 | gen_table=(modp_number*)mdmALLOC(sizeof(modp_number)*myp); |
---|
| 415 | gen_table[1]=1; |
---|
| 416 | for (pos_gen=2;pos_gen<myp;pos_gen++) |
---|
| 417 | { |
---|
| 418 | gen_ok=true; |
---|
| 419 | for (i=2;i<myp;i++) |
---|
| 420 | { |
---|
| 421 | gen_table[i]=modp_mul(gen_table[i-1],pos_gen); |
---|
| 422 | if (gen_table[i]==1) |
---|
| 423 | { |
---|
| 424 | gen_ok=false; |
---|
| 425 | break; |
---|
| 426 | } |
---|
| 427 | } |
---|
| 428 | if (gen_ok) break; |
---|
| 429 | } |
---|
| 430 | for (i=2;i<myp;i++) modp_Reverse[gen_table[i]]=gen_table[myp-i+1]; |
---|
| 431 | modp_Reverse[1]=1; |
---|
| 432 | mdmFREE(gen_table); |
---|
| 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; |
---|
| 443 | mdmFREE(ptr->mon); |
---|
| 444 | mdmFREE(ptr); |
---|
| 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 | { |
---|
| 456 | mdmFREE(points[i][j]); |
---|
| 457 | } |
---|
| 458 | mdmFREE(points[i]); |
---|
| 459 | } |
---|
| 460 | mdmFREE(points); |
---|
| 461 | for (i=0;i<final_base_dim;i++) mdmFREE(condition_list[i].mon); |
---|
| 462 | mdmFREE(condition_list); |
---|
| 463 | for (i=0;i<n_points;i++) mdmFREE(modp_points[i]); |
---|
| 464 | mdmFREE(modp_points); |
---|
| 465 | if (!only_modp) |
---|
| 466 | { |
---|
| 467 | for (i=0;i<n_points;i++) |
---|
| 468 | { |
---|
| 469 | for (j=0;j<variables;j++) mpq_clear(q_points[i][j]); |
---|
| 470 | mdmFREE(q_points[i]); |
---|
| 471 | } |
---|
| 472 | mdmFREE(q_points); |
---|
| 473 | for (i=0;i<n_points;i++) |
---|
| 474 | { |
---|
| 475 | for (j=0;j<variables;j++) mpz_clear(int_points[i][j]); |
---|
| 476 | mdmFREE(int_points[i]); |
---|
| 477 | } |
---|
| 478 | mdmFREE(int_points); |
---|
| 479 | generic_lt=FreeMonList (generic_lt); |
---|
| 480 | for (i=0;i<=final_base_dim;i++) |
---|
| 481 | { |
---|
| 482 | mpz_clear(polycoef[i]); |
---|
| 483 | mdmFREE(polyexp[i]); |
---|
| 484 | } |
---|
| 485 | mdmFREE(polycoef); |
---|
| 486 | mdmFREE(polyexp); |
---|
| 487 | if (!only_modp) mpz_clear(common_denom); |
---|
| 488 | } |
---|
| 489 | for (i=0;i<final_base_dim;i++) |
---|
| 490 | { |
---|
| 491 | mdmFREE(generic_column_name[i]); |
---|
| 492 | } |
---|
| 493 | mdmFREE(generic_column_name); |
---|
| 494 | for (i=0;i<n_points;i++) mdmFREE(coord_exist[i]); |
---|
| 495 | mdmFREE(coord_exist); |
---|
| 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); |
---|
| 508 | mdmFREE(my_row); |
---|
| 509 | my_row=NULL; |
---|
| 510 | mdmFREE(my_solve_row); |
---|
| 511 | my_solve_row=NULL; |
---|
| 512 | ptr=row_list; |
---|
| 513 | while (ptr!=NULL) |
---|
| 514 | { |
---|
| 515 | pptr=ptr->next; |
---|
| 516 | mdmFREE(ptr->row_matrix); |
---|
| 517 | mdmFREE(ptr->row_solve); |
---|
| 518 | mdmFREE(ptr); |
---|
| 519 | ptr=pptr; |
---|
| 520 | } |
---|
| 521 | row_list=NULL; |
---|
| 522 | for (i=0;i<final_base_dim;i++) mdmFREE(column_name[i]); |
---|
| 523 | mdmFREE(column_name); |
---|
| 524 | mdmFREE(modp_Reverse); |
---|
| 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; |
---|
| 536 | mn=(exponent*)mdmALLOC(sizeof(exponent)*variables); |
---|
| 537 | memcpy(mn,mon,sizeof(exponent)*variables); |
---|
| 538 | for (k=0;k<variables;k++) |
---|
| 539 | { |
---|
| 540 | for (j=1;j<=con.mon[k];j++) // this loop computes the coefficient from derivation |
---|
| 541 | { |
---|
| 542 | *ev=modp_mul(*ev,mn[k]); |
---|
| 543 | mn[k]--; |
---|
| 544 | } |
---|
| 545 | *ev=modp_mul(*ev,points[con.point_ref][k][mn[k]]); |
---|
| 546 | } |
---|
| 547 | mdmFREE(mn); |
---|
| 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; |
---|
| 561 | mn=(exponent*)mdmALLOC(sizeof(exponent)*variables); |
---|
| 562 | memcpy(mn,mon,sizeof(exponent)*variables); |
---|
| 563 | for (k=0;k<variables;k++) |
---|
| 564 | { |
---|
| 565 | for (j=1;j<=con.mon[k];j++) // this loop computes the coefficient from derivation |
---|
| 566 | { |
---|
| 567 | mpz_set_si(mon_conv,mn[k]); // (***) |
---|
| 568 | mpz_mul(ev,ev,mon_conv); |
---|
| 569 | mn[k]--; |
---|
| 570 | } |
---|
| 571 | for (j=1;j<=mn[k];j++) mpz_mul(ev,ev,int_points[con.point_ref][k]); // this loop computes the product of coordinate |
---|
| 572 | } |
---|
| 573 | mdmFREE(mn); |
---|
| 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 | } |
---|
| 693 | mdmFREE(mon); |
---|
| 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; |
---|
| 765 | #ifdef debb |
---|
[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 | { |
---|
| 805 | if (p_ptr==NULL) |
---|
| 806 | check_list=c_ptr->next; |
---|
| 807 | else |
---|
| 808 | p_ptr->next=c_ptr->next; |
---|
| 809 | n_ptr=c_ptr->next; |
---|
| 810 | mdmFREE(c_ptr->mon); |
---|
| 811 | mdmFREE(c_ptr); |
---|
| 812 | c_ptr=n_ptr; |
---|
| 813 | } |
---|
| 814 | else |
---|
| 815 | { |
---|
| 816 | p_ptr=c_ptr; |
---|
| 817 | c_ptr=c_ptr->next; |
---|
| 818 | } |
---|
| 819 | } |
---|
| 820 | } |
---|
| 821 | |
---|
| 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; |
---|
| 829 | mdmFREE(check_list->mon); |
---|
| 830 | mdmFREE(check_list); |
---|
| 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 | } |
---|
| 872 | temp=(row_list_entry*)mdmALLOC(sizeof(row_list_entry)); |
---|
| 873 | (*temp).row_matrix=(modp_number*)mdmALLOC(sizeof(modp_number)*final_base_dim); |
---|
| 874 | memcpy((*temp).row_matrix,my_row,sizeof(modp_number)*(final_base_dim)); |
---|
| 875 | (*temp).row_solve=(modp_number*)mdmALLOC(sizeof(modp_number)*final_base_dim); |
---|
| 876 | memcpy((*temp).row_solve,my_solve_row,sizeof(modp_number)*(final_base_dim)); |
---|
| 877 | (*temp).first_col=first_col; |
---|
| 878 | temp->next=ptr; |
---|
| 879 | if (pptr==NULL) row_list=temp; else pptr->next=temp; |
---|
| 880 | memcpy(column_name[first_col],mon,sizeof(exponent)*variables); |
---|
| 881 | } |
---|
| 882 | |
---|
| 883 | |
---|
| 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; |
---|
| 941 | temp=(modp_result_entry*)mdmALLOC(sizeof(modp_result_entry)); |
---|
| 942 | if (cur_result==NULL) |
---|
| 943 | { |
---|
| 944 | modp_result=temp; |
---|
| 945 | temp->prev=NULL; |
---|
| 946 | } |
---|
| 947 | else |
---|
| 948 | { |
---|
| 949 | temp->prev=cur_result; |
---|
| 950 | cur_result->next=temp; |
---|
| 951 | } |
---|
| 952 | cur_result=temp; |
---|
| 953 | cur_result->next=NULL; |
---|
| 954 | cur_result->p=myp; |
---|
| 955 | cur_result->generator=NULL; |
---|
| 956 | cur_result->n_generators=0; |
---|
| 957 | n_results++; |
---|
| 958 | } |
---|
| 959 | |
---|
| 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; |
---|
| 968 | mdmFREE(cur_gen->coef); |
---|
| 969 | mdmFREE(cur_gen->lt); |
---|
| 970 | mdmFREE(cur_gen); |
---|
| 971 | cur_gen=next_gen; |
---|
| 972 | } |
---|
| 973 | mdmFREE(e); |
---|
| 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 | } |
---|
| 990 | temp=(generator_entry*)mdmALLOC(sizeof(generator_entry)); |
---|
| 991 | if (prev_ptr==NULL) cur_result->generator=temp; else prev_ptr->next=temp; |
---|
| 992 | temp->next=NULL; |
---|
| 993 | temp->coef=(modp_number*)mdmALLOC(sizeof(modp_number)*final_base_dim); |
---|
| 994 | memcpy(temp->coef,my_solve_row,sizeof(modp_number)*final_base_dim); |
---|
| 995 | temp->lt=ZeroMonomial (); |
---|
| 996 | memcpy(temp->lt,mon,sizeof(exponent)*variables); |
---|
| 997 | temp->ltcoef=1; |
---|
| 998 | cur_result->n_generators++; |
---|
| 999 | } |
---|
| 1000 | |
---|
| 1001 | void MultGenerators () // before reconstructing, all denominators must be cancelled |
---|
| 1002 | { |
---|
| 1003 | #ifndef integerstrategy |
---|
| 1004 | int i; |
---|
| 1005 | generator_entry *cur_ptr; |
---|
| 1006 | cur_ptr=cur_result->generator; |
---|
| 1007 | while (cur_ptr!=NULL) |
---|
| 1008 | { |
---|
| 1009 | for (i=0;i<final_base_dim;i++) |
---|
| 1010 | cur_ptr->coef[i]=modp_mul(cur_ptr->coef[i],modp_denom); |
---|
| 1011 | cur_ptr->ltcoef=modp_denom; |
---|
| 1012 | cur_ptr=cur_ptr->next; |
---|
| 1013 | } |
---|
| 1014 | #endif |
---|
| 1015 | } |
---|
| 1016 | |
---|
| 1017 | void PresentGenerator (int i) // only for debuging, writes a generator in its form in program |
---|
| 1018 | { |
---|
| 1019 | int j; |
---|
| 1020 | modp_result_entry *cur_ptr; |
---|
| 1021 | generator_entry *cur_gen; |
---|
| 1022 | cur_ptr=modp_result; |
---|
| 1023 | while (cur_ptr!=NULL) |
---|
| 1024 | { |
---|
| 1025 | cur_gen=cur_ptr->generator; |
---|
| 1026 | for (j=0;j<i;j++) cur_gen=cur_gen->next; |
---|
| 1027 | for (j=0;j<final_base_dim;j++) |
---|
| 1028 | { |
---|
| 1029 | Print("%d;", cur_gen->coef[j]); |
---|
| 1030 | } |
---|
[6b4ff12] | 1031 | PrintS(" and LT = "); |
---|
[b5e57e2] | 1032 | WriteMono (cur_gen->lt); |
---|
| 1033 | Print(" ( %d ) prime = %d\n", cur_gen->ltcoef, cur_ptr->p); |
---|
| 1034 | cur_ptr=cur_ptr->next; |
---|
| 1035 | } |
---|
| 1036 | } |
---|
| 1037 | |
---|
| 1038 | modp_number TakePrime (modp_number p) // takes "previous" (smaller) prime |
---|
| 1039 | { |
---|
| 1040 | modp_number d; |
---|
| 1041 | myp_index--; |
---|
| 1042 | return cf_getSmallPrime(myp_index); |
---|
| 1043 | } |
---|
| 1044 | |
---|
| 1045 | void PrepareChinese (int n) // initialization for CRA |
---|
| 1046 | { |
---|
| 1047 | int i,k; |
---|
| 1048 | modp_result_entry *cur_ptr; |
---|
| 1049 | cur_ptr=modp_result; |
---|
| 1050 | modp_number *congr_ptr; |
---|
| 1051 | modp_number prod; |
---|
[1f18bb] | 1052 | in_gamma=(modp_number*)mdmALLOC(sizeof(modp_number)*n); |
---|
[b5e57e2] | 1053 | congr=(modp_number*)mdmALLOC(sizeof(modp_number)*n); |
---|
| 1054 | congr_ptr=congr; |
---|
| 1055 | while (cur_ptr!=NULL) |
---|
| 1056 | { |
---|
| 1057 | *congr_ptr=cur_ptr->p; |
---|
| 1058 | cur_ptr=cur_ptr->next; |
---|
| 1059 | congr_ptr++; |
---|
| 1060 | } |
---|
| 1061 | for (k=1;k<n;k++) |
---|
| 1062 | { |
---|
| 1063 | prod=congr[0]%congr[k]; |
---|
| 1064 | for (i=1;i<=k-1;i++) prod=(prod*congr[i])%congr[k]; |
---|
[1f18bb] | 1065 | in_gamma[i]=OneInverse(prod,congr[k]); |
---|
[b5e57e2] | 1066 | } |
---|
| 1067 | mpz_init(bigcongr); |
---|
| 1068 | mpz_set_ui(bigcongr,congr[0]); |
---|
| 1069 | for (k=1;k<n;k++) mpz_mul_ui(bigcongr,bigcongr,congr[k]); |
---|
| 1070 | } |
---|
| 1071 | |
---|
| 1072 | void CloseChinese (int n) // after CRA |
---|
| 1073 | { |
---|
[1f18bb] | 1074 | mdmFREE(in_gamma); |
---|
[b5e57e2] | 1075 | mdmFREE(congr); |
---|
| 1076 | mpz_clear(bigcongr); |
---|
| 1077 | } |
---|
| 1078 | |
---|
| 1079 | void ClearGCD () // divides polynomials in basis by gcd of coefficients |
---|
| 1080 | { |
---|
| 1081 | bool first_gcd=true; |
---|
| 1082 | int i; |
---|
| 1083 | mpz_t g; |
---|
| 1084 | mpz_init(g); |
---|
| 1085 | for (i=0;i<=final_base_dim;i++) |
---|
| 1086 | { |
---|
| 1087 | if (mpz_sgn(polycoef[i])!=0) |
---|
| 1088 | { |
---|
| 1089 | if (first_gcd) |
---|
| 1090 | { |
---|
| 1091 | first_gcd=false; |
---|
| 1092 | mpz_set(g,polycoef[i]); |
---|
| 1093 | } |
---|
| 1094 | else |
---|
| 1095 | mpz_gcd(g,g,polycoef[i]); |
---|
| 1096 | } |
---|
| 1097 | } |
---|
| 1098 | for (i=0;i<=final_base_dim;i++) mpz_divexact(polycoef[i],polycoef[i],g); |
---|
| 1099 | mpz_clear(g); |
---|
| 1100 | } |
---|
| 1101 | |
---|
| 1102 | void ReconstructGenerator (int ngen,int n,bool show) // recostruction of generator from various modp comp |
---|
| 1103 | { |
---|
| 1104 | int size; |
---|
| 1105 | int i,j,k; |
---|
| 1106 | int coef; |
---|
| 1107 | char *str=NULL; |
---|
| 1108 | str=(char*)mdmALLOC(sizeof(char)*1000); |
---|
| 1109 | modp_result_entry *cur_ptr; |
---|
| 1110 | generator_entry *cur_gen; |
---|
| 1111 | modp_number *u; |
---|
| 1112 | modp_number *v; |
---|
| 1113 | modp_number temp; |
---|
| 1114 | mpz_t sol; |
---|
| 1115 | mpz_t nsol; |
---|
| 1116 | mpz_init(sol); |
---|
| 1117 | mpz_init(nsol); |
---|
| 1118 | u=(modp_number*)mdmALLOC(sizeof(modp_number)*n); |
---|
| 1119 | v=(modp_number*)mdmALLOC(sizeof(modp_number)*n); |
---|
| 1120 | for (coef=0;coef<=final_base_dim;coef++) |
---|
| 1121 | { |
---|
| 1122 | i=0; |
---|
| 1123 | cur_ptr=modp_result; |
---|
| 1124 | while (cur_ptr!=NULL) |
---|
| 1125 | { |
---|
| 1126 | cur_gen=cur_ptr->generator; |
---|
| 1127 | for (j=0;j<ngen;j++) cur_gen=cur_gen->next; // we have to take jth generator |
---|
| 1128 | if (coef<final_base_dim) u[i]=cur_gen->coef[coef]; else u[i]=cur_gen->ltcoef; |
---|
| 1129 | cur_ptr=cur_ptr->next; |
---|
| 1130 | i++; |
---|
| 1131 | } |
---|
| 1132 | v[0]=u[0]; // now CRA begins |
---|
| 1133 | for (k=1;k<n;k++) |
---|
| 1134 | { |
---|
| 1135 | temp=v[k-1]; |
---|
| 1136 | for (j=k-2;j>=0;j--) temp=(temp*congr[j]+v[j])%congr[k]; |
---|
| 1137 | temp=u[k]-temp; |
---|
| 1138 | if (temp<0) temp=temp+congr[k]; |
---|
[1f18bb] | 1139 | v[k]=(temp*in_gamma[k])%congr[k]; |
---|
[b5e57e2] | 1140 | } |
---|
| 1141 | mpz_set_si(sol,v[n-1]); |
---|
| 1142 | for (k=n-2;k>=0;k--) |
---|
| 1143 | { |
---|
| 1144 | mpz_mul_ui(sol,sol,congr[k]); |
---|
| 1145 | mpz_add_ui(sol,sol,v[k]); |
---|
| 1146 | } // now CRA ends |
---|
| 1147 | mpz_sub(nsol,sol,bigcongr); |
---|
| 1148 | int s=mpz_cmpabs(sol,nsol); |
---|
| 1149 | if (s>0) mpz_set(sol,nsol); // chooses representation closer to 0 |
---|
| 1150 | mpz_set(polycoef[coef],sol); |
---|
| 1151 | if (coef<final_base_dim) |
---|
| 1152 | memcpy(polyexp[coef],generic_column_name[coef],sizeof(exponent)*variables); |
---|
| 1153 | else |
---|
| 1154 | memcpy(polyexp[coef],MonListElement (generic_lt,ngen),sizeof(exponent)*variables); |
---|
| 1155 | #ifdef checksize |
---|
| 1156 | size=mpz_sizeinbase(sol,10); |
---|
| 1157 | if (size>maximal_size) maximal_size=size; |
---|
| 1158 | #endif |
---|
| 1159 | } |
---|
| 1160 | mdmFREE(u); |
---|
| 1161 | mdmFREE(v); |
---|
| 1162 | mdmFREE(str); |
---|
| 1163 | ClearGCD (); |
---|
| 1164 | mpz_clear(sol); |
---|
| 1165 | mpz_clear(nsol); |
---|
| 1166 | } |
---|
| 1167 | |
---|
| 1168 | void Discard () // some unlucky prime occures |
---|
| 1169 | { |
---|
| 1170 | modp_result_entry *temp; |
---|
| 1171 | #ifdef writemsg |
---|
| 1172 | Print(", prime=%d", cur_result->p); |
---|
| 1173 | #endif |
---|
| 1174 | bad_primes++; |
---|
| 1175 | if (good_primes>bad_primes) |
---|
| 1176 | { |
---|
| 1177 | #ifdef writemsg |
---|
| 1178 | Print("-discarding this comp (%dth)\n", n_results); |
---|
| 1179 | #endif |
---|
| 1180 | temp=cur_result; |
---|
| 1181 | cur_result=cur_result->prev; |
---|
| 1182 | cur_result->next=NULL; |
---|
| 1183 | n_results--; |
---|
| 1184 | FreeResultEntry (temp); |
---|
| 1185 | } |
---|
| 1186 | else |
---|
| 1187 | { |
---|
| 1188 | #ifdef writemsg |
---|
[6b4ff12] | 1189 | PrintS("-discarding ALL.\n"); |
---|
[b5e57e2] | 1190 | #endif |
---|
| 1191 | int i; |
---|
| 1192 | modp_result_entry *ntfree; |
---|
| 1193 | generator_entry *cur_gen; |
---|
| 1194 | temp=cur_result->prev; |
---|
| 1195 | while (temp!=NULL) |
---|
| 1196 | { |
---|
| 1197 | ntfree=temp->prev; |
---|
| 1198 | FreeResultEntry (temp); |
---|
| 1199 | temp=ntfree; |
---|
| 1200 | } |
---|
| 1201 | modp_result=cur_result; |
---|
| 1202 | cur_result->prev=NULL; |
---|
| 1203 | n_results=1; |
---|
| 1204 | good_primes=1; |
---|
| 1205 | bad_primes=0; |
---|
| 1206 | generic_n_generators=cur_result->n_generators; |
---|
| 1207 | cur_gen=cur_result->generator; |
---|
| 1208 | generic_lt=FreeMonList (generic_lt); |
---|
| 1209 | for (i=0;i<generic_n_generators;i++) |
---|
| 1210 | { |
---|
| 1211 | generic_lt=MonListAdd (generic_lt,cur_gen->lt); |
---|
| 1212 | cur_gen=cur_gen->next; |
---|
| 1213 | } |
---|
| 1214 | for (i=0;i<final_base_dim;i++) memcpy(generic_column_name[i],column_name[i],sizeof(exponent)*variables); |
---|
| 1215 | } |
---|
| 1216 | } |
---|
| 1217 | |
---|
| 1218 | void modp_SetColumnNames () // used by modp - sets column names, the old table will be destroyed |
---|
| 1219 | { |
---|
| 1220 | int i; |
---|
| 1221 | for (i=0;i<final_base_dim;i++) memcpy(generic_column_name[i],column_name[i],sizeof(exponent)*variables); |
---|
| 1222 | } |
---|
| 1223 | |
---|
| 1224 | void CheckColumnSequence () // checks if scheme of computations is as generic one |
---|
| 1225 | { |
---|
| 1226 | int i; |
---|
| 1227 | if (cur_result->n_generators!=generic_n_generators) |
---|
| 1228 | { |
---|
| 1229 | #ifdef writemsg |
---|
[6b4ff12] | 1230 | PrintS("wrong number of generators occured"); |
---|
[b5e57e2] | 1231 | #else |
---|
[6b4ff12] | 1232 | if (protocol) PrintS("ng"); |
---|
[b5e57e2] | 1233 | #endif |
---|
| 1234 | Discard (); |
---|
| 1235 | return; |
---|
| 1236 | } |
---|
| 1237 | if (denom_divisible) |
---|
| 1238 | { |
---|
| 1239 | #ifdef writemsg |
---|
[6b4ff12] | 1240 | PrintS("denom of coef divisible by p"); |
---|
[b5e57e2] | 1241 | #else |
---|
[6b4ff12] | 1242 | if (protocol) PrintS("dp"); |
---|
[b5e57e2] | 1243 | #endif |
---|
| 1244 | Discard (); |
---|
| 1245 | return; |
---|
| 1246 | } |
---|
| 1247 | generator_entry *cur_gen; |
---|
| 1248 | mon_list_entry *cur_mon; |
---|
| 1249 | cur_gen=cur_result->generator; |
---|
| 1250 | cur_mon=generic_lt; |
---|
| 1251 | for (i=0;i<generic_n_generators;i++) |
---|
| 1252 | { |
---|
| 1253 | if (!EqualMon(cur_mon->mon,cur_gen->lt)) |
---|
| 1254 | { |
---|
| 1255 | #ifdef writemsg |
---|
[6b4ff12] | 1256 | PrintS("wrong leading term occured"); |
---|
[b5e57e2] | 1257 | #else |
---|
[6b4ff12] | 1258 | if (protocol) PrintS("lt"); |
---|
[b5e57e2] | 1259 | #endif |
---|
| 1260 | Discard (); |
---|
| 1261 | return; |
---|
| 1262 | } |
---|
| 1263 | cur_gen=cur_gen->next; |
---|
| 1264 | cur_mon=cur_mon->next; |
---|
| 1265 | } |
---|
| 1266 | for (i=0;i<final_base_dim;i++) |
---|
| 1267 | { |
---|
| 1268 | if (!EqualMon(generic_column_name[i],column_name[i])) |
---|
| 1269 | { |
---|
| 1270 | #ifdef writemsg |
---|
[6b4ff12] | 1271 | PrintS("wrong seq of cols occured"); |
---|
[b5e57e2] | 1272 | #else |
---|
[6b4ff12] | 1273 | if (protocol) PrintS("sc"); |
---|
[b5e57e2] | 1274 | #endif |
---|
| 1275 | Discard (); |
---|
| 1276 | return; |
---|
| 1277 | } |
---|
| 1278 | } |
---|
| 1279 | good_primes++; |
---|
| 1280 | } |
---|
| 1281 | |
---|
| 1282 | void WriteGenerator () // writes generator (only for debugging) |
---|
| 1283 | { |
---|
| 1284 | char *str; |
---|
| 1285 | str=(char*)mdmALLOC(sizeof(char)*1000); |
---|
| 1286 | int i; |
---|
| 1287 | for (i=0;i<=final_base_dim;i++) |
---|
| 1288 | { |
---|
| 1289 | str=mpz_get_str(str,10,polycoef[i]); |
---|
[6b4ff12] | 1290 | PrintS(str); |
---|
| 1291 | PrintS("*"); |
---|
[b5e57e2] | 1292 | WriteMono(polyexp[i]); |
---|
[6b4ff12] | 1293 | PrintS(" "); |
---|
[b5e57e2] | 1294 | } |
---|
| 1295 | mdmFREE(str); |
---|
[6b4ff12] | 1296 | PrintLn(); |
---|
[b5e57e2] | 1297 | } |
---|
| 1298 | |
---|
| 1299 | bool CheckGenerator () // evaluates generator to check whether it is good |
---|
| 1300 | { |
---|
| 1301 | mpz_t val,sum,temp; |
---|
| 1302 | int con,i; |
---|
| 1303 | mpz_init(val); |
---|
| 1304 | mpz_init(sum); |
---|
| 1305 | for (con=0;con<final_base_dim;con++) |
---|
| 1306 | { |
---|
| 1307 | mpz_set_si(sum,0); |
---|
| 1308 | for (i=0;i<=final_base_dim;i++) |
---|
| 1309 | { |
---|
| 1310 | int_Evaluate(val, polyexp[i], condition_list[con]); |
---|
| 1311 | mpz_mul(val,val,polycoef[i]); |
---|
| 1312 | mpz_add(sum,sum,val); |
---|
| 1313 | } |
---|
| 1314 | if (mpz_sgn(sum)!=0) |
---|
| 1315 | { |
---|
| 1316 | mpz_clear(val); |
---|
| 1317 | mpz_clear(sum); |
---|
| 1318 | return false; |
---|
| 1319 | } |
---|
| 1320 | } |
---|
| 1321 | mpz_clear(val); |
---|
| 1322 | mpz_clear(sum); |
---|
| 1323 | return true; |
---|
| 1324 | } |
---|
| 1325 | |
---|
| 1326 | void ClearGenList () |
---|
| 1327 | { |
---|
| 1328 | gen_list_entry *temp; |
---|
| 1329 | int i; |
---|
| 1330 | while (gen_list!=NULL) |
---|
| 1331 | { |
---|
| 1332 | temp=gen_list->next; |
---|
| 1333 | for (i=0;i<=final_base_dim;i++) |
---|
| 1334 | { |
---|
| 1335 | mpz_clear(gen_list->polycoef[i]); |
---|
| 1336 | mdmFREE(gen_list->polyexp[i]); |
---|
| 1337 | } |
---|
| 1338 | mdmFREE(gen_list->polycoef); |
---|
| 1339 | mdmFREE(gen_list->polyexp); |
---|
| 1340 | mdmFREE(gen_list); |
---|
| 1341 | gen_list=temp; |
---|
| 1342 | } |
---|
| 1343 | } |
---|
| 1344 | |
---|
| 1345 | void UpdateGenList () |
---|
| 1346 | { |
---|
| 1347 | gen_list_entry *temp,*prev; |
---|
| 1348 | int i,j; |
---|
| 1349 | prev=NULL; |
---|
| 1350 | temp=gen_list; |
---|
| 1351 | exponent deg; |
---|
| 1352 | for (i=0;i<=final_base_dim;i++) |
---|
| 1353 | { |
---|
| 1354 | deg=MonDegree(polyexp[i]); |
---|
| 1355 | for (j=0;j<deg;j++) |
---|
| 1356 | { |
---|
| 1357 | mpz_mul(polycoef[i],polycoef[i],common_denom); |
---|
| 1358 | } |
---|
| 1359 | } |
---|
| 1360 | ClearGCD (); |
---|
| 1361 | while (temp!=NULL) |
---|
| 1362 | { |
---|
| 1363 | prev=temp; |
---|
| 1364 | temp=temp->next; |
---|
| 1365 | } |
---|
| 1366 | temp=(gen_list_entry*)mdmALLOC(sizeof(gen_list_entry)); |
---|
| 1367 | if (prev==NULL) gen_list=temp; else prev->next=temp; |
---|
| 1368 | temp->next=NULL; |
---|
| 1369 | temp->polycoef=(mpz_t*)mdmALLOC(sizeof(mpz_t)*(final_base_dim+1)); |
---|
| 1370 | temp->polyexp=(mono_type*)mdmALLOC(sizeof(mono_type)*(final_base_dim+1)); |
---|
| 1371 | for (i=0;i<=final_base_dim;i++) |
---|
| 1372 | { |
---|
| 1373 | mpz_init(temp->polycoef[i]); |
---|
| 1374 | mpz_set(temp->polycoef[i],polycoef[i]); |
---|
| 1375 | temp->polyexp[i]=ZeroMonomial (); |
---|
| 1376 | memcpy(temp->polyexp[i],polyexp[i],sizeof(exponent)*variables); |
---|
| 1377 | } |
---|
| 1378 | } |
---|
| 1379 | |
---|
| 1380 | void ShowGenList () |
---|
| 1381 | { |
---|
| 1382 | gen_list_entry *temp; |
---|
| 1383 | int i; |
---|
| 1384 | char *str; |
---|
| 1385 | str=(char*)mdmALLOC(sizeof(char)*1000); |
---|
| 1386 | temp=gen_list; |
---|
| 1387 | while (temp!=NULL) |
---|
| 1388 | { |
---|
[6b4ff12] | 1389 | PrintS("generator: "); |
---|
[b5e57e2] | 1390 | for (i=0;i<=final_base_dim;i++) |
---|
| 1391 | { |
---|
| 1392 | str=mpz_get_str(str,10,temp->polycoef[i]); |
---|
[6b4ff12] | 1393 | PrintS(str); |
---|
| 1394 | PrintS("*"); |
---|
[b5e57e2] | 1395 | WriteMono(temp->polyexp[i]); |
---|
| 1396 | } |
---|
[6b4ff12] | 1397 | PrintLn(); |
---|
[b5e57e2] | 1398 | temp=temp->next; |
---|
| 1399 | } |
---|
| 1400 | mdmFREE(str); |
---|
| 1401 | } |
---|
| 1402 | |
---|
| 1403 | |
---|
| 1404 | void modp_Main () |
---|
| 1405 | { |
---|
| 1406 | mono_type cur_mon; |
---|
| 1407 | cur_mon= ZeroMonomial (); |
---|
| 1408 | modp_denom=1; |
---|
| 1409 | bool row_is_zero; |
---|
| 1410 | |
---|
| 1411 | #ifdef debb |
---|
| 1412 | Info (); |
---|
| 1413 | #endif |
---|
| 1414 | |
---|
| 1415 | while (check_list!=NULL) |
---|
| 1416 | { |
---|
| 1417 | TakeNextMonomial (cur_mon); |
---|
| 1418 | ProduceRow (cur_mon); |
---|
| 1419 | #ifdef debb |
---|
| 1420 | cout << "row produced for monomial "; |
---|
| 1421 | WriteMono (cur_mon); |
---|
| 1422 | cout << endl; |
---|
| 1423 | Info (); |
---|
| 1424 | #endif |
---|
| 1425 | ReduceRow (); |
---|
| 1426 | row_is_zero = RowIsZero (); |
---|
| 1427 | if (row_is_zero) |
---|
| 1428 | { |
---|
| 1429 | lt_list=MonListAdd (lt_list,cur_mon); |
---|
| 1430 | ReduceCheckListByMon (cur_mon); |
---|
| 1431 | NewGenerator (cur_mon); |
---|
| 1432 | #ifdef debb |
---|
| 1433 | cout << "row is zero - linear dependence found (should be seen in my_solve_row)" << endl; |
---|
| 1434 | cout << "monomial added to leading terms list" << endl; |
---|
| 1435 | cout << "check list updated" << endl; |
---|
| 1436 | Info (); |
---|
| 1437 | #endif |
---|
| 1438 | } |
---|
| 1439 | else |
---|
| 1440 | { |
---|
| 1441 | base_list= MonListAdd (base_list,cur_mon); |
---|
| 1442 | UpdateCheckList (cur_mon); |
---|
| 1443 | ReduceCheckListByLTs (); |
---|
| 1444 | #ifdef debb |
---|
| 1445 | cout << "row is non-zero" << endl; |
---|
| 1446 | cout << "monomial added to quotient basis list" << endl; |
---|
| 1447 | cout << "new monomials added to check list" << endl; |
---|
| 1448 | cout << "check list reduced by monomials from leading term list" << endl; |
---|
| 1449 | Info (); |
---|
| 1450 | #endif |
---|
| 1451 | PrepareRow (cur_mon); |
---|
| 1452 | #ifdef debb |
---|
| 1453 | cout << "row prepared and put into matrix" << endl; |
---|
| 1454 | Info (); |
---|
| 1455 | #endif |
---|
| 1456 | } |
---|
| 1457 | } |
---|
| 1458 | mdmFREE(cur_mon); |
---|
| 1459 | } |
---|
| 1460 | |
---|
| 1461 | void ResolveCoeff (mpq_t c, number m) |
---|
| 1462 | { |
---|
| 1463 | if ((long)m & SR_INT) |
---|
| 1464 | { |
---|
| 1465 | long m_val=SR_TO_INT(m); |
---|
| 1466 | mpq_set_si(c,m_val,1); |
---|
| 1467 | } |
---|
| 1468 | else |
---|
| 1469 | { |
---|
| 1470 | if (m->s<2) |
---|
| 1471 | { |
---|
| 1472 | mpz_set(mpq_numref(c),&m->z); |
---|
| 1473 | mpz_set(mpq_denref(c),&m->n); |
---|
| 1474 | mpq_canonicalize(c); |
---|
| 1475 | } |
---|
| 1476 | else |
---|
| 1477 | { |
---|
| 1478 | mpq_set_z(c,&m->z); |
---|
| 1479 | } |
---|
| 1480 | } |
---|
| 1481 | } |
---|
| 1482 | |
---|
| 1483 | ideal interpolation(lists L, intvec *v) |
---|
| 1484 | { |
---|
| 1485 | protocol=TEST_OPT_PROT; // should be set if option(prot) is enabled |
---|
| 1486 | |
---|
| 1487 | bool data_ok=true; |
---|
| 1488 | |
---|
| 1489 | // reading the ring data *************************************************** |
---|
[dbbae8] | 1490 | if ((currRing==NULL) || ((!rField_is_Zp ())&&(!rField_is_Q ()))) |
---|
[b5e57e2] | 1491 | { |
---|
| 1492 | WerrorS("coefficient field should be Zp or Q!"); |
---|
| 1493 | return NULL; |
---|
| 1494 | } |
---|
| 1495 | if ((currRing->qideal)!=NULL) |
---|
| 1496 | { |
---|
| 1497 | WerrorS("quotient ring not supported!"); |
---|
| 1498 | return NULL; |
---|
| 1499 | } |
---|
| 1500 | if ((currRing->OrdSgn)!=1) |
---|
| 1501 | { |
---|
| 1502 | WerrorS("ordering must be global!"); |
---|
| 1503 | return NULL; |
---|
| 1504 | } |
---|
| 1505 | n_points=v->length (); |
---|
| 1506 | if (n_points!=(L->nr+1)) |
---|
| 1507 | { |
---|
| 1508 | WerrorS("list and intvec must have the same length!"); |
---|
| 1509 | return NULL; |
---|
| 1510 | } |
---|
| 1511 | variables=currRing->N; |
---|
| 1512 | only_modp=rField_is_Zp(); |
---|
| 1513 | if (only_modp) myp=rChar(); |
---|
| 1514 | // ring data read ********************************************************** |
---|
| 1515 | |
---|
| 1516 | |
---|
| 1517 | multiplicity=(int*)malloc(sizeof(int)*n_points); |
---|
| 1518 | int i; |
---|
| 1519 | for (i=0;i<n_points;i++) multiplicity[i]=(*v)[i]; |
---|
| 1520 | |
---|
| 1521 | final_base_dim = CalcBaseDim (); |
---|
| 1522 | |
---|
| 1523 | #ifdef writemsg |
---|
| 1524 | Print("number of variables: %d\n", variables); |
---|
| 1525 | Print("number of points: %d\n", n_points); |
---|
[6b4ff12] | 1526 | PrintS("multiplicities: "); |
---|
[b5e57e2] | 1527 | for (i=0;i<n_points;i++) Print("%d ", multiplicity[i]); |
---|
[6b4ff12] | 1528 | PrintLn(); |
---|
[b5e57e2] | 1529 | Print("general initialization for dimension %d ...\n", final_base_dim); |
---|
| 1530 | #endif |
---|
| 1531 | |
---|
| 1532 | GeneralInit (); |
---|
| 1533 | |
---|
| 1534 | // reading coordinates of points from ideals ********************************** |
---|
| 1535 | mpq_t divisor; |
---|
| 1536 | if (!only_modp) mpq_init(divisor); |
---|
| 1537 | int j; |
---|
| 1538 | for(i=0; i<=L->nr;i++) |
---|
| 1539 | { |
---|
| 1540 | ideal I=(ideal)L->m[i].Data(); |
---|
| 1541 | for(j=0;j<IDELEMS(I);j++) |
---|
| 1542 | { |
---|
| 1543 | poly p=I->m[j]; |
---|
| 1544 | if (p!=NULL) |
---|
| 1545 | { |
---|
| 1546 | poly ph=pHead(p); |
---|
| 1547 | int pcvar=pVar(ph); |
---|
| 1548 | if (pcvar!=0) |
---|
| 1549 | { |
---|
| 1550 | pcvar--; |
---|
| 1551 | if (coord_exist[i][pcvar]) |
---|
| 1552 | { |
---|
| 1553 | Print("coordinate %d for point %d initialized twice!\n",pcvar+1,i+1); |
---|
| 1554 | data_ok=false; |
---|
| 1555 | } |
---|
| 1556 | number m; |
---|
| 1557 | m=pGetCoeff(p); // possible coefficient standing by a leading monomial |
---|
| 1558 | if (!only_modp) ResolveCoeff (divisor,m); |
---|
| 1559 | number n; |
---|
| 1560 | if (pNext(p)!=NULL) n=pGetCoeff(pNext(p)); |
---|
| 1561 | else n=nInit(0); |
---|
| 1562 | if (only_modp) |
---|
| 1563 | { |
---|
| 1564 | n=nNeg(n); |
---|
| 1565 | n=nDiv(n,m); |
---|
| 1566 | modp_points[i][pcvar]=(int)((long)n); |
---|
| 1567 | } |
---|
| 1568 | else |
---|
| 1569 | { |
---|
| 1570 | ResolveCoeff (q_points[i][pcvar],n); |
---|
| 1571 | mpq_neg(q_points[i][pcvar],q_points[i][pcvar]); |
---|
| 1572 | mpq_div(q_points[i][pcvar],q_points[i][pcvar],divisor); |
---|
| 1573 | } |
---|
| 1574 | coord_exist[i][pcvar]=true; |
---|
| 1575 | } |
---|
| 1576 | else |
---|
| 1577 | { |
---|
[6b4ff12] | 1578 | PrintS("not a variable? "); |
---|
[b5e57e2] | 1579 | wrp(p); |
---|
| 1580 | PrintLn(); |
---|
| 1581 | data_ok=false; |
---|
| 1582 | } |
---|
| 1583 | pDelete(&ph); |
---|
| 1584 | } |
---|
| 1585 | } |
---|
| 1586 | } |
---|
| 1587 | if (!only_modp) mpq_clear(divisor); |
---|
| 1588 | // data from ideal read ******************************************************* |
---|
| 1589 | |
---|
| 1590 | // ckecking if all coordinates are initialized |
---|
| 1591 | for (i=0;i<n_points;i++) |
---|
| 1592 | { |
---|
| 1593 | for (j=0;j<variables;j++) |
---|
| 1594 | { |
---|
| 1595 | if (!coord_exist[i][j]) |
---|
| 1596 | { |
---|
| 1597 | Print("coordinate %d for point %d not known!\n",j+1,i+1); |
---|
| 1598 | data_ok=false; |
---|
| 1599 | } |
---|
| 1600 | } |
---|
| 1601 | } |
---|
| 1602 | |
---|
| 1603 | if (!data_ok) |
---|
| 1604 | { |
---|
[f4fea85] | 1605 | GeneralDone(); |
---|
[b5e57e2] | 1606 | WerrorS("data structure is invalid"); |
---|
| 1607 | return NULL; |
---|
| 1608 | } |
---|
| 1609 | |
---|
| 1610 | if (!only_modp) IntegerPoints (); |
---|
| 1611 | MakeConditions (); |
---|
| 1612 | #ifdef writemsg |
---|
[6b4ff12] | 1613 | PrintS("done.\n"); |
---|
[b5e57e2] | 1614 | #else |
---|
| 1615 | if (protocol) Print("[vdim %d]",final_base_dim); |
---|
| 1616 | #endif |
---|
| 1617 | |
---|
| 1618 | |
---|
| 1619 | // main procedure ********************************************************************* |
---|
| 1620 | int modp_cycles=10; |
---|
| 1621 | bool correct_gen=false; |
---|
| 1622 | if (only_modp) modp_cycles=1; |
---|
| 1623 | myp_index=cf_getNumSmallPrimes (); |
---|
| 1624 | |
---|
| 1625 | while ((!correct_gen)&&(myp_index>1)) |
---|
| 1626 | { |
---|
| 1627 | #ifdef writemsg |
---|
| 1628 | Print("trying %d cycles mod p...\n",modp_cycles); |
---|
| 1629 | #else |
---|
| 1630 | if (protocol) Print("(%d)",modp_cycles); |
---|
| 1631 | #endif |
---|
| 1632 | while ((n_results<modp_cycles)&&(myp_index>1)) // some computations mod p |
---|
| 1633 | { |
---|
| 1634 | if (!only_modp) myp=TakePrime (myp); |
---|
| 1635 | NewResultEntry (); |
---|
| 1636 | InitProcData (); |
---|
| 1637 | if (only_modp) modp_PrepareProducts (); else int_PrepareProducts (); |
---|
| 1638 | |
---|
| 1639 | modp_Main (); |
---|
| 1640 | |
---|
| 1641 | if (!only_modp) |
---|
| 1642 | { |
---|
| 1643 | MultGenerators (); |
---|
| 1644 | CheckColumnSequence (); |
---|
| 1645 | } |
---|
| 1646 | else |
---|
| 1647 | { |
---|
| 1648 | modp_SetColumnNames (); |
---|
| 1649 | } |
---|
| 1650 | FreeProcData (); |
---|
| 1651 | } |
---|
| 1652 | |
---|
| 1653 | if (!only_modp) |
---|
| 1654 | { |
---|
| 1655 | PrepareChinese (modp_cycles); |
---|
| 1656 | correct_gen=true; |
---|
| 1657 | for (i=0;i<generic_n_generators;i++) |
---|
| 1658 | { |
---|
| 1659 | ReconstructGenerator (i,modp_cycles,false); |
---|
| 1660 | correct_gen=CheckGenerator (); |
---|
| 1661 | if (!correct_gen) |
---|
| 1662 | { |
---|
| 1663 | #ifdef writemsg |
---|
[6b4ff12] | 1664 | PrintS("wrong generator!\n"); |
---|
[b5e57e2] | 1665 | #else |
---|
[6b4ff12] | 1666 | // if (protocol) PrintS("!g"); |
---|
[b5e57e2] | 1667 | #endif |
---|
| 1668 | ClearGenList (); |
---|
| 1669 | break; |
---|
| 1670 | } |
---|
| 1671 | else |
---|
| 1672 | { |
---|
| 1673 | UpdateGenList (); |
---|
| 1674 | } |
---|
| 1675 | } |
---|
| 1676 | #ifdef checksize |
---|
| 1677 | Print("maximal size of output: %d, precision bound: %d.\n",maximalsize,mpz_sizeinbase(bigcongr,10)); |
---|
| 1678 | #endif |
---|
| 1679 | CloseChinese (modp_cycles); |
---|
| 1680 | modp_cycles=modp_cycles+10; |
---|
| 1681 | } |
---|
| 1682 | else |
---|
| 1683 | { |
---|
| 1684 | correct_gen=true; |
---|
| 1685 | } |
---|
| 1686 | } |
---|
| 1687 | // end of main procedure ************************************************************************************ |
---|
| 1688 | |
---|
| 1689 | #ifdef writemsg |
---|
[6b4ff12] | 1690 | PrintS("computations finished.\n"); |
---|
[b5e57e2] | 1691 | #else |
---|
[6b4ff12] | 1692 | if (protocol) PrintLn(); |
---|
[b5e57e2] | 1693 | #endif |
---|
| 1694 | |
---|
| 1695 | if (!correct_gen) |
---|
| 1696 | { |
---|
| 1697 | GeneralDone (); |
---|
| 1698 | ClearGenList (); |
---|
| 1699 | WerrorS("internal error - coefficient too big!"); |
---|
| 1700 | return NULL; |
---|
| 1701 | } |
---|
| 1702 | |
---|
| 1703 | // passing data to ideal ************************************************************************************* |
---|
| 1704 | ideal ret; |
---|
| 1705 | |
---|
| 1706 | if (only_modp) |
---|
| 1707 | { |
---|
| 1708 | mono_type mon; |
---|
| 1709 | ret=idInit(modp_result->n_generators,1); |
---|
| 1710 | generator_entry *cur_gen=modp_result->generator; |
---|
| 1711 | for(i=0;i<IDELEMS(ret);i++) |
---|
| 1712 | { |
---|
| 1713 | poly p,sum; |
---|
| 1714 | sum=NULL; |
---|
| 1715 | int a; |
---|
| 1716 | int cf; |
---|
| 1717 | for (a=final_base_dim;a>=0;a--) |
---|
| 1718 | { |
---|
| 1719 | if (a==final_base_dim) cf=cur_gen->ltcoef; else cf=cur_gen->coef[a]; |
---|
| 1720 | if (cf!=0) |
---|
| 1721 | { |
---|
| 1722 | p=pISet(cf); |
---|
| 1723 | if (a==final_base_dim) mon=cur_gen->lt; else mon=generic_column_name[a]; |
---|
| 1724 | for (j=0;j<variables;j++) pSetExp(p,j+1,mon[j]); |
---|
| 1725 | pSetm(p); |
---|
| 1726 | sum=pAdd(sum,p); |
---|
| 1727 | } |
---|
| 1728 | } |
---|
| 1729 | ret->m[i]=sum; |
---|
| 1730 | cur_gen=cur_gen->next; |
---|
| 1731 | } |
---|
| 1732 | } |
---|
| 1733 | else |
---|
| 1734 | { |
---|
| 1735 | ret=idInit(generic_n_generators,1); |
---|
| 1736 | gen_list_entry *temp=gen_list; |
---|
| 1737 | for(i=0;i<IDELEMS(ret);i++) |
---|
| 1738 | { |
---|
| 1739 | poly p,sum; |
---|
| 1740 | sum=NULL; |
---|
| 1741 | int a; |
---|
| 1742 | for (a=final_base_dim;a>=0;a--) // build one polynomial |
---|
| 1743 | { |
---|
| 1744 | if (mpz_sgn(temp->polycoef[a])!=0) |
---|
| 1745 | { |
---|
| 1746 | p=pOne(); //a monomial |
---|
| 1747 | for (j=0;j<variables;j++) pSetExp(p,j+1,temp->polyexp[a][j]); |
---|
| 1748 | pSetm(p); // after all pSetExp |
---|
| 1749 | number n=(number)omAllocBin(rnumber_bin); |
---|
| 1750 | #ifdef LDEBUG |
---|
| 1751 | n->debug=123456; |
---|
| 1752 | #endif |
---|
| 1753 | mpz_init_set(&n->z,temp->polycoef[a]); |
---|
| 1754 | n->s=3; |
---|
| 1755 | nlNormalize(n); |
---|
| 1756 | pSetCoeff(p,n); |
---|
| 1757 | sum=pAdd(sum,p); |
---|
| 1758 | } |
---|
| 1759 | } |
---|
| 1760 | ret->m[i]=sum; |
---|
| 1761 | temp=temp->next; |
---|
| 1762 | } |
---|
| 1763 | } |
---|
| 1764 | // data transferred **************************************************************************** |
---|
| 1765 | |
---|
| 1766 | |
---|
| 1767 | GeneralDone (); |
---|
| 1768 | ClearGenList (); |
---|
| 1769 | return ret; |
---|
| 1770 | } |
---|
| 1771 | |
---|
| 1772 | |
---|
| 1773 | BOOLEAN jjINTERPOLATION (leftv res, leftv l, leftv v) |
---|
| 1774 | { |
---|
| 1775 | res->data=interpolation((lists)l->Data(),(intvec*)v->Data()); |
---|
| 1776 | setFlag(res,FLAG_STD); |
---|
| 1777 | return errorreported; |
---|
| 1778 | } |
---|