Changeset 30d574 in git
- Timestamp:
- Nov 11, 2010, 10:23:22 AM (13 years ago)
- Branches:
- (u'spielwiese', '0d6b7fcd9813a1ca1ed4220cfa2b104b97a0a003')
- Children:
- a80025e8d85c97a2cc2a6b86ef219a392df4f2b5
- Parents:
- d882518b93f286d7db8ab355fea6f1781b41f6d1
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
Singular/interpolation.cc
rd88251 r30d574 15 15 #include <kernel/ring.h> 16 16 #ifdef HAVE_FACTORY 17 #define SI_DONT_HAVE_GLOBAL_VARS 17 18 # include <factory/factory.h> 18 19 #endif /* HAVE_FACTORY */ 19 20 20 //memory management21 #define mdmALLOC(x) omAlloc0(x)22 #define mdmFREE(x) omFree(x)23 24 21 // parameters to debug 25 //#define debb26 22 //#define shmat 27 23 //#define checksize … … 144 140 #endif 145 141 142 #if 0 /* only for debuggig*/ 146 143 void WriteMono (mono_type m) // Writes a monomial on the screen - only for debug 147 144 { … … 151 148 } 152 149 153 #ifdef debb154 150 void WriteMonoList (mon_list_entry *list) 155 151 { … … 228 224 u1=1; u2=0; 229 225 u = a; v = p; 230 while (v != 0) { 226 while (v != 0) 227 { 231 228 q = u / v; 232 229 r = u % v; … … 310 307 if (Greater ((*curptr).mon,mon)) break; 311 308 prevptr=curptr; 312 curptr=curptr->next;} 313 temp=(mon_list_entry*)mdmALLOC(sizeof(mon_list_entry)); 309 curptr=curptr->next; 310 } 311 temp=(mon_list_entry*)omAlloc0(sizeof(mon_list_entry)); 314 312 (*temp).next=curptr; 315 (*temp).mon=(exponent*) mdmALLOC(sizeof(exponent)*variables);313 (*temp).mon=(exponent*)omAlloc0(sizeof(exponent)*variables); 316 314 memcpy(temp->mon,mon,sizeof(exponent)*variables); 317 if (prevptr==NULL) return temp; else { 315 if (prevptr==NULL) return temp; 316 else 317 { 318 318 prevptr->next=temp; 319 return list;} 319 return list; 320 } 320 321 } 321 322 … … 331 332 { 332 333 mono_type m; 333 m=(exponent*) mdmALLOC(sizeof(exponent)*variables);334 m=(exponent*)omAlloc0(sizeof(exponent)*variables); 334 335 int i; 335 336 for (i=0;i<variables;i++) m[i]=0; … … 340 341 { 341 342 int i,j,k; 342 points=(coordinates*) mdmALLOC(sizeof(coordinates)*n_points);343 points=(coordinates*)omAlloc0(sizeof(coordinates)*n_points); 343 344 for (i=0;i<n_points;i++) 344 345 { 345 points[i]=(coordinate_products*) mdmALLOC(sizeof(coordinate_products)*variables);346 for (j=0;j<variables;j++) points[i][j]=(modp_number*) mdmALLOC(sizeof(modp_number)*(max_coord));347 } 348 condition_list=(condition_type*) mdmALLOC(sizeof(condition_type)*final_base_dim);349 for (i=0;i<final_base_dim;i++) condition_list[i].mon=(exponent*) mdmALLOC(sizeof(exponent)*variables);350 modp_points=(modp_coordinates*) mdmALLOC(sizeof(modp_coordinates)*n_points);351 for (i=0;i<n_points;i++) modp_points[i]=(modp_number*) mdmALLOC(sizeof(modp_number)*variables);346 points[i]=(coordinate_products*)omAlloc0(sizeof(coordinate_products)*variables); 347 for (j=0;j<variables;j++) points[i][j]=(modp_number*)omAlloc0(sizeof(modp_number)*(max_coord)); 348 } 349 condition_list=(condition_type*)omAlloc0(sizeof(condition_type)*final_base_dim); 350 for (i=0;i<final_base_dim;i++) condition_list[i].mon=(exponent*)omAlloc0(sizeof(exponent)*variables); 351 modp_points=(modp_coordinates*)omAlloc0(sizeof(modp_coordinates)*n_points); 352 for (i=0;i<n_points;i++) modp_points[i]=(modp_number*)omAlloc0(sizeof(modp_number)*variables); 352 353 if (!only_modp) 353 354 { 354 q_points=(q_coordinates*) mdmALLOC(sizeof(q_coordinates)*n_points);355 q_points=(q_coordinates*)omAlloc0(sizeof(q_coordinates)*n_points); 355 356 for (i=0;i<n_points;i++) 356 357 { 357 q_points[i]=(mpq_t*) mdmALLOC(sizeof(mpq_t)*variables);358 q_points[i]=(mpq_t*)omAlloc0(sizeof(mpq_t)*variables); 358 359 for (j=0;j<variables;j++) mpq_init(q_points[i][j]); 359 360 } 360 int_points=(int_coordinates*) mdmALLOC(sizeof(int_coordinates)*n_points);361 int_points=(int_coordinates*)omAlloc0(sizeof(int_coordinates)*n_points); 361 362 for (i=0;i<n_points;i++) 362 363 { 363 int_points[i]=(mpz_t*) mdmALLOC(sizeof(mpz_t)*variables);364 int_points[i]=(mpz_t*)omAlloc0(sizeof(mpz_t)*variables); 364 365 for (j=0;j<variables;j++) mpz_init(int_points[i][j]); 365 366 } 366 367 } 367 coord_exist=(coord_exist_table*) mdmALLOC(sizeof(coord_exist_table)*n_points);368 coord_exist=(coord_exist_table*)omAlloc0(sizeof(coord_exist_table)*n_points); 368 369 for (i=0;i<n_points;i++) 369 370 { 370 coord_exist[i]=(bool*) mdmALLOC(sizeof(bool)*variables);371 coord_exist[i]=(bool*)omAlloc0(sizeof(bool)*variables); 371 372 for (j=0;j<variables;j++) coord_exist[i][j]=false; 372 373 } 373 generic_column_name=(mono_type*) mdmALLOC(sizeof(mono_type)*final_base_dim);374 generic_column_name=(mono_type*)omAlloc0(sizeof(mono_type)*final_base_dim); 374 375 for (i=0;i<final_base_dim;i++) generic_column_name[i]=ZeroMonomial (); 375 376 good_primes=0; … … 378 379 if (!only_modp) 379 380 { 380 polycoef=(mpz_t*) mdmALLOC(sizeof(mpz_t)*(final_base_dim+1));381 polyexp=(mono_type*) mdmALLOC(sizeof(mono_type)*(final_base_dim+1));381 polycoef=(mpz_t*)omAlloc0(sizeof(mpz_t)*(final_base_dim+1)); 382 polyexp=(mono_type*)omAlloc0(sizeof(mono_type)*(final_base_dim+1)); 382 383 for (i=0;i<=final_base_dim;i++) 383 384 { … … 403 404 int i,j; 404 405 check_list=MonListAdd (check_list,ZeroMonomial ()); 405 my_row=(modp_number*) mdmALLOC(sizeof(modp_number)*final_base_dim);406 my_solve_row=(modp_number*) mdmALLOC(sizeof(modp_number)*final_base_dim);407 column_name=(mono_type*) mdmALLOC(sizeof(mono_type)*final_base_dim);406 my_row=(modp_number*)omAlloc0(sizeof(modp_number)*final_base_dim); 407 my_solve_row=(modp_number*)omAlloc0(sizeof(modp_number)*final_base_dim); 408 column_name=(mono_type*)omAlloc0(sizeof(mono_type)*final_base_dim); 408 409 for (i=0;i<final_base_dim;i++) column_name[i]=ZeroMonomial (); 409 410 last_solve_column=0; … … 411 412 modp_number pos_gen; 412 413 bool gen_ok; 413 modp_Reverse=(modp_number*) mdmALLOC(sizeof(modp_number)*myp);414 modp_Reverse=(modp_number*)omAlloc0(sizeof(modp_number)*myp); 414 415 415 416 // produces table of modp inverts by finding a generator of (Z_myp*,*) 416 gen_table=(modp_number*) mdmALLOC(sizeof(modp_number)*myp);417 gen_table=(modp_number*)omAlloc0(sizeof(modp_number)*myp); 417 418 gen_table[1]=1; 418 419 for (pos_gen=2;pos_gen<myp;pos_gen++) … … 432 433 for (i=2;i<myp;i++) modp_Reverse[gen_table[i]]=gen_table[myp-i+1]; 433 434 modp_Reverse[1]=1; 434 mdmFREE(gen_table);435 omFree(gen_table); 435 436 } 436 437 … … 443 444 { 444 445 pptr=ptr->next; 445 mdmFREE(ptr->mon);446 mdmFREE(ptr);446 omFree(ptr->mon); 447 omFree(ptr); 447 448 ptr=pptr;} 448 449 return NULL; … … 456 457 for (j=0;j<variables;j++) 457 458 { 458 mdmFREE(points[i][j]);459 omFree(points[i][j]); 459 460 } 460 mdmFREE(points[i]);461 } 462 mdmFREE(points);463 for (i=0;i<final_base_dim;i++) mdmFREE(condition_list[i].mon);464 mdmFREE(condition_list);465 for (i=0;i<n_points;i++) mdmFREE(modp_points[i]);466 mdmFREE(modp_points);461 omFree(points[i]); 462 } 463 omFree(points); 464 for (i=0;i<final_base_dim;i++) omFree(condition_list[i].mon); 465 omFree(condition_list); 466 for (i=0;i<n_points;i++) omFree(modp_points[i]); 467 omFree(modp_points); 467 468 if (!only_modp) 468 469 { … … 470 471 { 471 472 for (j=0;j<variables;j++) mpq_clear(q_points[i][j]); 472 mdmFREE(q_points[i]);473 } 474 mdmFREE(q_points);473 omFree(q_points[i]); 474 } 475 omFree(q_points); 475 476 for (i=0;i<n_points;i++) 476 477 { 477 478 for (j=0;j<variables;j++) mpz_clear(int_points[i][j]); 478 mdmFREE(int_points[i]);479 } 480 mdmFREE(int_points);479 omFree(int_points[i]); 480 } 481 omFree(int_points); 481 482 generic_lt=FreeMonList (generic_lt); 482 483 for (i=0;i<=final_base_dim;i++) 483 484 { 484 485 mpz_clear(polycoef[i]); 485 mdmFREE(polyexp[i]);486 } 487 mdmFREE(polycoef);488 mdmFREE(polyexp);486 omFree(polyexp[i]); 487 } 488 omFree(polycoef); 489 omFree(polyexp); 489 490 if (!only_modp) mpz_clear(common_denom); 490 491 } 491 492 for (i=0;i<final_base_dim;i++) 492 493 { 493 mdmFREE(generic_column_name[i]);494 } 495 mdmFREE(generic_column_name);496 for (i=0;i<n_points;i++) mdmFREE(coord_exist[i]);497 mdmFREE(coord_exist);494 omFree(generic_column_name[i]); 495 } 496 omFree(generic_column_name); 497 for (i=0;i<n_points;i++) omFree(coord_exist[i]); 498 omFree(coord_exist); 498 499 pDelete(&comparizon_p1); 499 500 pDelete(&comparizon_p2); … … 508 509 lt_list=FreeMonList (lt_list); 509 510 base_list=FreeMonList (base_list); 510 mdmFREE(my_row);511 omFree(my_row); 511 512 my_row=NULL; 512 mdmFREE(my_solve_row);513 omFree(my_solve_row); 513 514 my_solve_row=NULL; 514 515 ptr=row_list; … … 516 517 { 517 518 pptr=ptr->next; 518 mdmFREE(ptr->row_matrix);519 mdmFREE(ptr->row_solve);520 mdmFREE(ptr);519 omFree(ptr->row_matrix); 520 omFree(ptr->row_solve); 521 omFree(ptr); 521 522 ptr=pptr; 522 523 } 523 524 row_list=NULL; 524 for (i=0;i<final_base_dim;i++) mdmFREE(column_name[i]);525 mdmFREE(column_name);526 mdmFREE(modp_Reverse);525 for (i=0;i<final_base_dim;i++) omFree(column_name[i]); 526 omFree(column_name); 527 omFree(modp_Reverse); 527 528 } 528 529 … … 536 537 int j,k; 537 538 mono_type mn; 538 mn=(exponent*) mdmALLOC(sizeof(exponent)*variables);539 mn=(exponent*)omAlloc0(sizeof(exponent)*variables); 539 540 memcpy(mn,mon,sizeof(exponent)*variables); 540 541 for (k=0;k<variables;k++) … … 547 548 *ev=modp_mul(*ev,points[con.point_ref][k][mn[k]]); 548 549 } 549 mdmFREE(mn);550 omFree(mn); 550 551 } 551 552 … … 561 562 int j,k; 562 563 mono_type mn; 563 mn=(exponent*) mdmALLOC(sizeof(exponent)*variables);564 mn=(exponent*)omAlloc0(sizeof(exponent)*variables); 564 565 memcpy(mn,mon,sizeof(exponent)*variables); 565 566 for (k=0;k<variables;k++) … … 573 574 for (j=1;j<=mn[k];j++) mpz_mul(ev,ev,int_points[con.point_ref][k]); // this loop computes the product of coordinate 574 575 } 575 mdmFREE(mn);576 omFree(mn); 576 577 mpz_clear(mon_conv); 577 578 } … … 693 694 } 694 695 } 695 mdmFREE(mon);696 omFree(mon); 696 697 } 697 698 … … 765 766 } 766 767 row_ptr=row_ptr->next; 767 #if def debb768 #if 0 /* only debugging */ 768 769 PrintS("reduction by row "); 769 770 Info (); … … 805 806 if (DivisibleMon (m,c_ptr->mon)) 806 807 { 807 if (p_ptr==NULL) 808 if (p_ptr==NULL) 808 809 check_list=c_ptr->next; 809 810 else 810 811 p_ptr->next=c_ptr->next; 811 812 n_ptr=c_ptr->next; 812 mdmFREE(c_ptr->mon);813 mdmFREE(c_ptr);813 omFree(c_ptr->mon); 814 omFree(c_ptr); 814 815 c_ptr=n_ptr; 815 816 } … … 829 830 memcpy(mon,check_list->mon,sizeof(exponent)*variables); 830 831 n_check_list=check_list->next; 831 mdmFREE(check_list->mon);832 mdmFREE(check_list);832 omFree(check_list->mon); 833 omFree(check_list); 833 834 check_list=n_check_list; 834 835 } … … 872 873 ptr=ptr->next; 873 874 } 874 temp=(row_list_entry*) mdmALLOC(sizeof(row_list_entry));875 (*temp).row_matrix=(modp_number*) mdmALLOC(sizeof(modp_number)*final_base_dim);875 temp=(row_list_entry*)omAlloc0(sizeof(row_list_entry)); 876 (*temp).row_matrix=(modp_number*)omAlloc0(sizeof(modp_number)*final_base_dim); 876 877 memcpy((*temp).row_matrix,my_row,sizeof(modp_number)*(final_base_dim)); 877 (*temp).row_solve=(modp_number*) mdmALLOC(sizeof(modp_number)*final_base_dim);878 (*temp).row_solve=(modp_number*)omAlloc0(sizeof(modp_number)*final_base_dim); 878 879 memcpy((*temp).row_solve,my_solve_row,sizeof(modp_number)*(final_base_dim)); 879 880 (*temp).first_col=first_col; … … 941 942 modp_result_entry *temp; 942 943 int i; 943 temp=(modp_result_entry*) mdmALLOC(sizeof(modp_result_entry));944 temp=(modp_result_entry*)omAlloc0(sizeof(modp_result_entry)); 944 945 if (cur_result==NULL) 945 946 { … … 968 969 { 969 970 next_gen=cur_gen->next; 970 mdmFREE(cur_gen->coef);971 mdmFREE(cur_gen->lt);972 mdmFREE(cur_gen);971 omFree(cur_gen->coef); 972 omFree(cur_gen->lt); 973 omFree(cur_gen); 973 974 cur_gen=next_gen; 974 975 } 975 mdmFREE(e);976 omFree(e); 976 977 } 977 978 … … 990 991 cur_ptr=cur_ptr->next; 991 992 } 992 temp=(generator_entry*) mdmALLOC(sizeof(generator_entry));993 temp=(generator_entry*)omAlloc0(sizeof(generator_entry)); 993 994 if (prev_ptr==NULL) cur_result->generator=temp; else prev_ptr->next=temp; 994 995 temp->next=NULL; 995 temp->coef=(modp_number*) mdmALLOC(sizeof(modp_number)*final_base_dim);996 temp->coef=(modp_number*)omAlloc0(sizeof(modp_number)*final_base_dim); 996 997 memcpy(temp->coef,my_solve_row,sizeof(modp_number)*final_base_dim); 997 998 temp->lt=ZeroMonomial (); … … 1009 1010 while (cur_ptr!=NULL) 1010 1011 { 1011 for (i=0;i<final_base_dim;i++) 1012 for (i=0;i<final_base_dim;i++) 1012 1013 cur_ptr->coef[i]=modp_mul(cur_ptr->coef[i],modp_denom); 1013 1014 cur_ptr->ltcoef=modp_denom; … … 1016 1017 #endif 1017 1018 } 1018 1019 #if 0 /* only debbuging */ 1019 1020 void PresentGenerator (int i) // only for debuging, writes a generator in its form in program 1020 1021 { … … 1037 1038 } 1038 1039 } 1040 #endif 1039 1041 1040 1042 modp_number TakePrime (modp_number p) // takes "previous" (smaller) prime … … 1055 1057 modp_number *congr_ptr; 1056 1058 modp_number prod; 1057 in_gamma=(modp_number*) mdmALLOC(sizeof(modp_number)*n);1058 congr=(modp_number*) mdmALLOC(sizeof(modp_number)*n);1059 in_gamma=(modp_number*)omAlloc0(sizeof(modp_number)*n); 1060 congr=(modp_number*)omAlloc0(sizeof(modp_number)*n); 1059 1061 congr_ptr=congr; 1060 1062 while (cur_ptr!=NULL) … … 1077 1079 void CloseChinese (int n) // after CRA 1078 1080 { 1079 mdmFREE(in_gamma);1080 mdmFREE(congr);1081 omFree(in_gamma); 1082 omFree(congr); 1081 1083 mpz_clear(bigcongr); 1082 1084 } … … 1111 1113 int coef; 1112 1114 char *str=NULL; 1113 str=(char*) mdmALLOC(sizeof(char)*1000);1115 str=(char*)omAlloc0(sizeof(char)*1000); 1114 1116 modp_result_entry *cur_ptr; 1115 1117 generator_entry *cur_gen; … … 1121 1123 mpz_init(sol); 1122 1124 mpz_init(nsol); 1123 u=(modp_number*) mdmALLOC(sizeof(modp_number)*n);1124 v=(modp_number*) mdmALLOC(sizeof(modp_number)*n);1125 u=(modp_number*)omAlloc0(sizeof(modp_number)*n); 1126 v=(modp_number*)omAlloc0(sizeof(modp_number)*n); 1125 1127 for (coef=0;coef<=final_base_dim;coef++) 1126 1128 { … … 1163 1165 #endif 1164 1166 } 1165 mdmFREE(u);1166 mdmFREE(v);1167 mdmFREE(str);1167 omFree(u); 1168 omFree(v); 1169 omFree(str); 1168 1170 ClearGCD (); 1169 1171 mpz_clear(sol); … … 1284 1286 good_primes++; 1285 1287 } 1286 1288 #if 0 /* only debuggig */ 1287 1289 void WriteGenerator () // writes generator (only for debugging) 1288 1290 { 1289 1291 char *str; 1290 str=(char*) mdmALLOC(sizeof(char)*1000);1292 str=(char*)omAlloc0(sizeof(char)*1000); 1291 1293 int i; 1292 1294 for (i=0;i<=final_base_dim;i++) … … 1298 1300 PrintS(" "); 1299 1301 } 1300 mdmFREE(str);1302 omFree(str); 1301 1303 PrintLn(); 1302 1304 } 1305 #endif 1303 1306 1304 1307 bool CheckGenerator () // evaluates generator to check whether it is good … … 1339 1342 { 1340 1343 mpz_clear(gen_list->polycoef[i]); 1341 mdmFREE(gen_list->polyexp[i]);1344 omFree(gen_list->polyexp[i]); 1342 1345 } 1343 mdmFREE(gen_list->polycoef);1344 mdmFREE(gen_list->polyexp);1345 mdmFREE(gen_list);1346 omFree(gen_list->polycoef); 1347 omFree(gen_list->polyexp); 1348 omFree(gen_list); 1346 1349 gen_list=temp; 1347 1350 } … … 1369 1372 temp=temp->next; 1370 1373 } 1371 temp=(gen_list_entry*) mdmALLOC(sizeof(gen_list_entry));1374 temp=(gen_list_entry*)omAlloc0(sizeof(gen_list_entry)); 1372 1375 if (prev==NULL) gen_list=temp; else prev->next=temp; 1373 1376 temp->next=NULL; 1374 temp->polycoef=(mpz_t*) mdmALLOC(sizeof(mpz_t)*(final_base_dim+1));1375 temp->polyexp=(mono_type*) mdmALLOC(sizeof(mono_type)*(final_base_dim+1));1377 temp->polycoef=(mpz_t*)omAlloc0(sizeof(mpz_t)*(final_base_dim+1)); 1378 temp->polyexp=(mono_type*)omAlloc0(sizeof(mono_type)*(final_base_dim+1)); 1376 1379 for (i=0;i<=final_base_dim;i++) 1377 1380 { … … 1383 1386 } 1384 1387 1388 #if 0 /* only debugging */ 1385 1389 void ShowGenList () 1386 1390 { … … 1388 1392 int i; 1389 1393 char *str; 1390 str=(char*) mdmALLOC(sizeof(char)*1000);1394 str=(char*)omAlloc0(sizeof(char)*1000); 1391 1395 temp=gen_list; 1392 1396 while (temp!=NULL) … … 1403 1407 temp=temp->next; 1404 1408 } 1405 mdmFREE(str); 1406 } 1409 omFree(str); 1410 } 1411 #endif 1407 1412 1408 1413 … … 1414 1419 bool row_is_zero; 1415 1420 1416 #if def debb1421 #if 0 /* only debugging */ 1417 1422 Info (); 1418 1423 #endif … … 1422 1427 TakeNextMonomial (cur_mon); 1423 1428 ProduceRow (cur_mon); 1424 #if def debb1429 #if 0 /* only debugging */ 1425 1430 cout << "row produced for monomial "; 1426 1431 WriteMono (cur_mon); … … 1435 1440 ReduceCheckListByMon (cur_mon); 1436 1441 NewGenerator (cur_mon); 1437 #if def debb1442 #if 0 /* only debugging */ 1438 1443 cout << "row is zero - linear dependence found (should be seen in my_solve_row)" << endl; 1439 1444 cout << "monomial added to leading terms list" << endl; … … 1447 1452 UpdateCheckList (cur_mon); 1448 1453 ReduceCheckListByLTs (); 1449 #if def debb1454 #if 0 /* only debugging */ 1450 1455 cout << "row is non-zero" << endl; 1451 1456 cout << "monomial added to quotient basis list" << endl; … … 1455 1460 #endif 1456 1461 PrepareRow (cur_mon); 1457 #if def debb1462 #if 0 /* only debugging */ 1458 1463 cout << "row prepared and put into matrix" << endl; 1459 1464 Info (); … … 1461 1466 } 1462 1467 } 1463 mdmFREE(cur_mon);1468 omFree(cur_mon); 1464 1469 } 1465 1470
Note: See TracChangeset
for help on using the changeset viewer.