Changeset cc07066 in git
- Timestamp:
- Aug 31, 2020, 11:57:27 AM (4 years ago)
- Branches:
- (u'spielwiese', '5b153614cbc72bfa198d75b1e9e33dab2645d9fe')
- Children:
- fdea544a15d2eb16f952ba088de11e5efa5494d9
- Parents:
- 6924d452049dc8bfbb851672c4b3ceac05081b6c
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
Singular/LIB/pfd.lib
r6924d4 rcc07066 5 5 LIBRARY: pfd.lib Multivariate Partial Fraction Decomposition 6 6 7 AUTHOR: Marcel Wittmann, e-mail: mwittman@ rhrk.uni-kl.de7 AUTHOR: Marcel Wittmann, e-mail: mwittman@mathematik.uni-kl.de 8 8 9 9 OVERVIEW: … … 12 12 \"smaller\" numerators and denominators. 13 13 This can be used to shorten the IBP reduction coeffcients of multi-loop Feynman 14 integrals as shown in [J. B öhm, M. Wittmann, Z. Wu, Y. Xu, Y. Zhang: 'IBP14 integrals as shown in [J. Boehm, M. Wittmann, Z. Wu, Y. Xu, Y. Zhang: 'IBP 15 15 reduction coefficients made simple'] (preprint 2020). For this application, 16 16 we also provide a procedure that applies the algorithm to all entries of a … … 30 30 getStringpfd(); turn a decomposition gotten from @code{pfd} into one 31 31 string 32 getStringpfd_indexed(); like @code{getStringpfd}, but writ ingthe denominator33 factors just as @code{q1}, @code{q2},...32 getStringpfd_indexed(); like @code{getStringpfd}, but writes the denominator 33 factors just as @code{q1}, @code{q2}, ... 34 34 readInputTXT(); read a matrix of rational functions from a txt-file 35 35 pfdMat(); apply @code{pfd} to a matrix of rational functions … … 50 50 static proc mod_init() 51 51 { 52 printlevel = 2; 53 // export some static functions so they can be run in parallel using parallelWaitAll: 52 54 if (!defined(Tasks)) {LIB "tasks.lib";} 53 // export some static functions so they can be run in parallel using parallelWaitAll54 55 exportto(Tasks,pfdWrap); 55 56 importfrom(Tasks,pfdWrap); 56 57 exportto(Pfd,pfdWrap); 57 58 58 exportto(Tasks,testEntry); 59 59 importfrom(Tasks,testEntry); … … 66 66 pfd(arguments[, parallelize]); arguments list, parallelize int 67 67 RETURN: a partial fraction decomposition of f/g as a list @code{l} where 68 @code{l[1]} is an ideal generate by irreducible polynomials and68 @code{l[1]} is an ideal generated by irreducible polynomials and 69 69 @code{l[2]} is a list of fractions. 70 70 Each fraction is represented by a list of 71 71 @* 1) the numerator polynomial 72 @* 2) an intvec giving the indices @code{i} for which @code{l[1][i]}73 occursas a factor in the denominator72 @* 2) an intvec of indices @code{i} for which @code{l[1][i]} occurs 73 as a factor in the denominator 74 74 @* 3) an intvec containing the exponents of those irreducible factors. 75 @* Setting debugto a positive integer measures runtimes and76 creates a log file (default: debug=0).75 @* Setting @code{debug} to a positive integer measures runtimes and 76 creates a log file (default: @code{debug=0}). 77 77 @* The denominator g can also be given in factorized form as a list of 78 78 an ideal of irreducible non constant polynomials and an intvec of 79 79 exponents. This can save time since the first step in the algorithm is 80 to factorize g. ( g=list(ideal(0),intvec(0:0)) represents a denominator81 of 1.)80 to factorize g. (A list of the zero-ideal and an empty intvec 81 represents a denominator of 1.) 82 82 @* If instead of f and g, the input is a single list (or even a list of 83 lists) containing elements of the form list(f,g[,debug]) (f,g[,debug] 84 as above), the algorithm is applied to all entries in parallel (using 85 @ref{parallel_lib}), if @code{parallelize=1} (default) and in sequence 86 if @code{parallelize=0}. A list of the results is returned. 83 lists) containing elements of the form @code{list(f,g[,debug])} 84 (@code{f,g,debug} as above), the algorithm is applied to all entries 85 in parallel (using @ref{parallel_lib}), if @code{parallelize=1} 86 (default) and in sequence if @code{parallelize=0}. A list (or list of 87 lists) of the results is returned. 87 88 NOTE: The result depends on the monomial ordering. For \"small\" results 88 use dp.89 use @code{dp}. 89 90 EXAMPLE: example pfd; shows an example 90 91 SEE ALSO: checkpfd, evaluatepfd, displaypfd, displaypfd_long, pfdMat" … … 141 142 142 143 int debug=0; link l=":w "; 143 if(size(#)>2) 144 if(size(#)>3) {ERROR("wrong number of arguments, expected 2 or 3");} 145 if(size(#)==3) 144 146 { 145 147 debug=#[3]; … … 160 162 if(debug) 161 163 {fprintf(l,"%ntotal: 0 ms (numerator was 0)",0); close(l);} 162 if(voice<3) 163 {displaypfd(dec);} 164 if(voice<=printlevel) {displaypfd(dec);} 164 165 return(dec); 165 166 } … … 171 172 if(debug) 172 173 {fprintf(l,"%ntotal: 0 ms (denominator was constant)",0); close(l);} 173 if(voice<3) 174 {displaypfd(dec);} 174 if(voice<=printlevel) {displaypfd(dec);} 175 175 return(dec); 176 176 } 177 177 178 // (1) factorization of the denominator :178 // (1) factorization of the denominator //////////////////////////////////// 179 179 if(debug) {int t1 = rtimer; write(l,"factorizing ");} 180 180 list factor = factorize(g); 181 181 number lcoeff; 182 for(i=2; i<=size(factor[1]); i++) // making polynomials unique182 for(i=2; i<=size(factor[1]); i++) 183 183 { 184 184 lcoeff = leadcoef(factor[1][i]); … … 203 203 if(debug) 204 204 {fprintf(l,"%ntotal: 0 ms (denominator was constant)",0); close(l);} 205 if(voice<3) 206 {displaypfd(dec);} 205 if(voice<=printlevel) {displaypfd(dec);} 207 206 return(dec); 208 207 } … … 221 220 } 222 221 else 223 224 225 // (2) Nullstellensatz decomposition 222 {ERROR("wrong type for second argument, expected poly or list(ideal,intvec)");}} 223 224 // (2) Nullstellensatz decomposition ///////////////////////////////////////// 226 225 if(debug) 227 226 {int t2 = rtimer; write(l,"Nullstellensatz decomposition ");} … … 248 247 {t2 = rtimer-t2; fprintf(l,"done! (%s ms, %s terms)", t2, size(dec));} 249 248 250 // (3) short numerator decomposition 249 // (3) short numerator decomposition ///////////////////////////////////////// 251 250 if(debug) 252 251 {int t3 = rtimer; write(l,"short numerator decompositions "); counter=0;} … … 277 276 {t3 = rtimer-t3; fprintf(l,"done! (%s ms, %s terms)", t3, size(dec));} 278 277 279 // (4) algebraic dependence decomposition 278 // (4) algebraic dependence decomposition //////////////////////////////////// 280 279 if(debug) 281 280 {int t4 = rtimer; write(l,"algebraic dependence decomposition "); counter=0;} … … 305 304 {t4 = rtimer-t4; fprintf(l,"done! (%s ms, %s terms)", t4, size(dec));} 306 305 307 // (5) numerator decomposition 306 // (5) numerator decomposition /////////////////////////////////////////////// 308 307 if(debug) 309 308 {int t5 = rtimer; write(l,"numerator decompositions "); counter=0;} … … 335 334 dec = list(q,dec); 336 335 if(debug) {fprintf(l,"%ntotal: %s ms", t1+t2+t3+t4+t5); close(l);} 337 if(voice< 3) {displaypfd(dec);}336 if(voice<=printlevel) {displaypfd(dec);} 338 337 return(dec); 339 338 } … … 366 365 // a more complicated example 367 366 ring S = 0,(s12,s15,s23,s34,s45),dp; 368 poly f = (7*s12^4*s15^2 + 11*s12^3*s15^3 + 4*s12^2*s15^4 - 10*s12^4*s15*s23367 poly f = 7*s12^4*s15^2 + 11*s12^3*s15^3 + 4*s12^2*s15^4 - 10*s12^4*s15*s23 369 368 - 14*s12^3*s15^2*s23 - 4*s12^2*s15^3*s23 + 3*s12^4*s23^2 + 3*s12^3*s15*s23^2 370 369 + 13*s12^4*s15*s34 + 12*s12^3*s15^2*s34 + 2*s12^2*s15^3*s34 … … 402 401 - 8*s12*s15*s34*s45^3 - 10*s15^2*s34*s45^3 + 9*s12*s23*s34*s45^3 403 402 + s12*s34^2*s45^3 + 8*s15*s34^2*s45^3 - 9*s23*s34^2*s45^3 - s34^3*s45^3 404 - s15^2*s45^4 + s15*s34*s45^4 );403 - s15^2*s45^4 + s15*s34*s45^4; 405 404 poly g = 4*s12*s15*(s12 + s15 - s34)*(s15 - s23 - s34)*(s12 + s23 - s45) 406 405 *(s12 - s34 - s45)*(s12 + s15 - s34 - s45)*s45; … … 424 423 425 424 if(m==0) // denominator is 1 426 { 427 return(list(l)); //do nothing, return input 428 } 425 {return(list(l));} // do nothing, return input 429 426 430 427 ideal qe = q[indices]; 431 428 for(int i=1; i<=m; i++) 432 { 433 qe[i] = qe[i]^e[i]; 434 } 429 {qe[i] = qe[i]^e[i];} 435 430 matrix T; 436 431 ideal qe_std = liftstd(qe,T); … … 438 433 if(deg(qe_std) == 0) 439 434 { 440 T = T/qe_std[1]; // now 1 = T[1,1]*qe[1] + ... + T[m,1]*qe[m] is a Nullstellensatz certificate 435 T = T/qe_std[1]; 436 // now 1 = T[1,1]*qe[1] + ... + T[m,1]*qe[m] is a Nullstellensatz certificate 441 437 list dec; 442 438 poly h; … … 445 441 h = T[i,1]; 446 442 if(h != 0) 447 { 448 dec[size(dec)+1] = list(f*h,delete(indices,i),delete(e,i)); 449 } 443 {dec[size(dec)+1] = list(f*h,delete(indices,i),delete(e,i));} 450 444 } 451 445 return(dec); … … 453 447 else 454 448 { 455 return(list(l)); // do nothing, return input449 return(list(l)); // do nothing, return input 456 450 } 457 451 } … … 472 466 if(debug) 473 467 {fprintf(ll," shortNumeratorDecompStep: %s ms (m=%s, e=%s) " 474 +"--> constant denominator", rtimer-tt, m, e);}475 return(list(list(),list(l))); // do nothing, return input476 } 477 478 ideal q_denom = q[indices]; // factorsin the denominator468 +"--> constant denominator", rtimer-tt, m, e);} 469 return(list(list(),list(l))); // do nothing, return input 470 } 471 472 ideal q_denom = q[indices]; // factors occuring in the denominator 479 473 matrix T; 480 474 ideal q_std = liftstd(q_denom,T); … … 488 482 +"--> remainder is nonzero", rtimer-tt, m, e);} 489 483 return(list(list(),list(l))); // if there is a rest, decomposing further would 490 // not help in the next step (alg. depend. decomposition) 491 } 484 } // not help in the next step (alg. depend. decomposition) 492 485 493 486 matrix a = divrem[1]/divrem[3][1,1]; // now f = r + a[1,1]*q_std[1] + ... +a[m,1]*q_std[m] 494 a = T*a; // lift coefficients ==>now f = r + a[1,1]*q[1] + ... +a[m,1]*q[m]487 a = T*a; // lift coefficients ==> now f = r + a[1,1]*q[1] + ... +a[m,1]*q[m] 495 488 496 489 // reduce w.r.t. groebner basis of syz(q) to make the numerators "smaller": … … 520 513 } 521 514 dec[size(dec)+1] = fraction; 522 515 } 523 516 524 517 if(debug) … … 544 537 if(debug) 545 538 {fprintf(ll," algDependDecompStep: %s ms (m=%s, e=%s)", rtimer-tt, m, e);} 546 return(list(l)); // do nothing, return input539 return(list(l)); // do nothing, return input 547 540 } 548 541 549 542 if(m<=d) 550 543 { 551 if(size(syz(module(transpose(jacob(ideal(q[indices]))))))==0) // jacobian criterion544 if(size(syz(module(transpose(jacob(ideal(q[indices]))))))==0) // jacobian criterion 552 545 { 553 546 if(debug) 554 547 {fprintf(ll," algDependDecompStep: %s ms (m=%s, e=%s) " 555 548 +"--> alg. indep.", rtimer-tt, m, e);} 556 return(list(l)); // do nothing, return input549 return(list(l)); // do nothing, return input 557 550 } 558 551 } … … 567 560 ideal I; 568 561 for(i=1; i<=m; i++) 569 { 570 I[i] = y(i)-q[indices[i]];//^e[i]; 571 } 562 {I[i] = y(i)-q[indices[i]];} 572 563 573 564 ideal annihilatingPolys = eliminate(I,intvec(1..d)); … … 575 566 poly g = annihilatingPolys[1]; 576 567 577 poly tail = g[size(g)]; // term of lowest dp-order (thus lowest degree)568 poly tail = g[size(g)]; // term of lowest dp-order (thus lowest degree) 578 569 number tcoeff = leadcoef(tail); 579 570 intvec texpon = leadexp(tail); … … 597 588 for(i=1; i<=m; i++) 598 589 { 599 //pow = e[i]*(expon[i]-texpon[i]-1);600 590 pow = expon[i]-texpon[i]-e[i]; 601 591 if(pow>=0) … … 634 624 {fprintf(ll," numeratorDecompStep: %s ms (m=%s, e=%s) " 635 625 +"--> constant denominator", rtimer-tt, m, e);} 636 return(list(list(),list(l))); // do nothing, return input637 } 638 639 ideal q_denom = q[indices]; // factors in the denominator626 return(list(list(),list(l))); // do nothing, return input 627 } 628 629 ideal q_denom = q[indices]; // factors in the denominator 640 630 matrix T; 641 631 ideal q_std = liftstd(q_denom,T); … … 643 633 matrix a = divrem[1]/divrem[3][1,1]; 644 634 poly r = divrem[2][1]/divrem[3][1,1]; // now f = r + a[1,1]*q_std[1] + ... +a[m,1]*q_std[m] 645 a = T*a; // lift coefficients ==>now f = r + a[1,1]*q[1] + ... +a[m,1]*q[m]635 a = T*a; // lift coefficients ==> now f = r + a[1,1]*q[1] + ... +a[m,1]*q[m] 646 636 647 637 // reduce w.r.t. groebner basis of syz(q) to make the numerators "smaller" … … 652 642 list fraction,dec,rest; 653 643 if(r!=0) 654 { 655 rest[1] = list(r,indices,e); 656 } 644 {rest[1] = list(r,indices,e);} 657 645 for(i=1; i<=m; i++) 658 646 { … … 684 672 685 673 static proc mergepfd(list dec1, list dec2) 686 // assumes dec1 is already sorted w.r.t. dp in the denominator exponents 687 { 674 { 675 // Note: assumes dec1 is already sorted w.r.t. dp in the denominator exponents 688 676 int n1=size(dec1); 689 677 int n2=size(dec2); … … 717 705 } 718 706 if(is_dp_smaller(dec1[a][2], dec1[a][3], entry[2], entry[3])) 719 { // numerator exponent vector of dec1[a] is dp-smaller than entry 720 dec1=insert(dec1,entry,a); n1++; 721 } 707 {dec1=insert(dec1,entry,a); n1++;} 722 708 else{if(entry[2]==dec1[a][2] && entry[3]==dec1[a][3]) 723 709 { 724 dec1[a][1] = dec1[a][1] + entry[1]; //same denominator: add numerators710 dec1[a][1] = dec1[a][1] + entry[1]; //same denominator: add numerators 725 711 if(dec1[a][1]==0) 726 { 727 dec1 = delete(dec1,a); n1--; 728 } 712 {dec1 = delete(dec1,a); n1--;} 729 713 } 730 714 else 731 { 732 dec1=insert(dec1,entry,a-1); n1++; 733 }} 715 {dec1=insert(dec1,entry,a-1); n1++;}} 734 716 } 735 717 return(dec1); … … 758 740 "USAGE: checkpfd(list(f,g),dec[,N,C]); f,g poly, dec list, N,C int 759 741 RETURN: 0 or 1 760 PURPOSE: test sfor (mathematical) equality of f/g and a partial fraction742 PURPOSE: test for (mathematical) equality of f/g and a partial fraction 761 743 decomposition dec. The list dec has to have the same structure as the 762 744 output of @ref{pfd}. 763 745 @* The denominator g can also be given in factorized form as a list of 764 746 an ideal of irreducible non constant polynomials and an intvec of 765 exponents. (g=list(ideal(0),intvec(0:0)) represents a denominator 766 of 1.) 747 exponents. This can save time since the first step in the algorithm is 748 to factorize g. (a list of the zero-ideal and an empty intvec 749 represents a denominator of 1.) 767 750 @* By default the test is done (exactly) by bringing all terms of the 768 751 decomposition on the same denominator and comparing to f/g. 769 @* If additional parameters N [, C] are given and if N>0, a probabilistic770 method is chosen: evaluation at N random points with coordinates771 between -C and C. This may be faster for big polynomials.752 @* If additional parameters N [, C] are given and if @code{N>0}, a 753 probabilistic method is chosen: evaluation at N random points with 754 coordinates between -C and C. This may be faster for big polynomials. 772 755 SEE ALSO: pfd 773 756 EXAMPLE: example checkpfd; shows an example" 774 757 { 775 poly f =fraction[1];776 def g =fraction[2];758 poly f = fraction[1]; 759 def g = fraction[2]; 777 760 if(size(#)>0) 778 761 { … … 793 776 794 777 if(typeof(g)=="poly") 795 { 796 val1 = number(substitute(g,vars,values)); 797 } 798 else{if(typeof(g)=="list") //denominator given in factorized form 778 {val1 = number(substitute(g,vars,values));} 779 else{if(typeof(g)=="list") // denominator given in factorized form 799 780 { 800 781 val1 = number(1); 801 782 for(j=1; j<=size(g[1]); j++) 802 { 803 val1 = val1 * number(substitute(g[1][j]^g[2][j],vars,values)); 804 } 783 {val1 = val1 * number(substitute(g[1][j]^g[2][j],vars,values));} 805 784 } 806 785 else {ERROR("wrong type for second argument, expected poly or list");}} … … 811 790 if(div_by_0) {continue;} 812 791 if(val1 != val2) 813 { 814 return(0); 815 } 792 {return(0);} 816 793 } 817 794 return(1); … … 848 825 num = dec[i][1]; 849 826 for(j=1; j<=m; j++) 850 { 851 num = num * q[j]^(e[j]); 852 } 827 {num = num * q[j]^(e[j]);} 853 828 sum_of_numerators = sum_of_numerators + num; 854 829 } 855 856 830 // the decomposition is now equal to sum_of_numerators/product(q[i]^e_max[i]) (i from 1 to imax) 857 831 // now: check if this equals f/g: 858 859 832 if(typeof(g)=="poly") 860 833 { … … 864 837 num = f/fact[1][1]; 865 838 } 866 else{if(typeof(g)=="list") //denominator given in factorized form839 else{if(typeof(g)=="list") //denominator given in factorized form 867 840 { 868 841 ideal q_g = g[1]; … … 877 850 for(i=1; i<=m_g; i++) 878 851 { 879 880 852 j=0; 881 853 for(k=1; k<=m; k++) … … 884 856 if(c*q_g[i]==q[k]) {j=k; break;} 885 857 } 886 887 858 if(j==0) 888 { 889 sum_of_numerators = sum_of_numerators*q_g[i]^e_g[i]; 890 } 859 {sum_of_numerators = sum_of_numerators*q_g[i]^e_g[i];} 891 860 else 892 861 { 893 num = num*(c^e_g[i]); //fix lead coefficient 894 862 num = num*(c^e_g[i]); //fix lead coefficient 895 863 expon = e_g[i]-e_max[j]; 896 864 if(expon>0) 897 { 898 sum_of_numerators = sum_of_numerators*q[j]^expon; 899 } 865 {sum_of_numerators = sum_of_numerators*q[j]^expon;} 900 866 else{if(expon<0) 901 { 902 num = num*q[j]^(-expon); 903 }} 867 {num = num*q[j]^(-expon);}} 904 868 } 905 869 } … … 933 897 "USAGE: evaluatepfd(dec, values[, mode]); dec list, values ideal, mode int 934 898 RETURN: the number gotten by substituting the numbers generating the ideal 935 values for the variables in the partial fraction decomposition dec. 936 The list dec has to have the same structure as the output of @ref{pfd}. 937 @* mode = 1: raise Error in case of division by 0 (default) 938 @* mode = 2: return a second integer which is 1 if the denominator 899 @code{values} for the variables in the partial fraction decomposition 900 @code{dec}. The list @code{dec} has to have the same structure as the 901 output of @ref{pfd}. 902 @* @code{mode=1}: raise Error in case of division by 0 (default) 903 @* @code{mode=2}: return a second integer which is 1 if the denominator 939 904 becomes 0, and 0 otherwise. 940 905 SEE ALSO: pfd … … 944 909 if(size(#)>0) {mode = #[1];} 945 910 946 ideal vars = maxideal(1); // ideal generate by ring variables911 ideal vars = maxideal(1); // ideal generated by ring variables 947 912 number val=0; 948 913 number denom; … … 978 943 979 944 displaypfd_long(dec); 945 // evaluation at x=2, y=1: 980 946 ideal values = 2,1; 981 evaluate(dec,values); 982 // compare: f(2,1)/g(2,1) = (2+2*1)/(2^2-1^1) = 4/7 947 evaluatepfd(dec,values); 948 949 // compare: f(2,1)/g(2,1) = (2+2*1)/(2^2-1^1) = 4/3 983 950 } 984 951 … … 987 954 int n=size(l); 988 955 for(int i=1; i<=n; i++) 989 { 990 if(entry==l[i]) {return(i);} 991 } 956 {if(entry==l[i]) {return(i);}} 992 957 return(0); 993 958 } … … 995 960 proc displaypfd(list dec) 996 961 "USAGE: displaypfd(dec); dec list 997 PURPOSE: print a partial fraction decomposition dec in a readable way. 998 The list dec has to have the same structure as the output of @ref{pfd}. 962 PURPOSE: print a partial fraction decomposition @code{dec} in a readable way. 963 The list @code{dec} has to have the same structure as the output of 964 @ref{pfd}. 999 965 SEE ALSO: pfd, displaypfd_long, getStringpfd, getStringpfd_indexed 1000 966 EXAMPLE: example displaypfd; shows an example" … … 1032 998 printf("q%s = %s",i,q[i]); 1033 999 } 1034 printf("(%s terms)%n", size(dec), 0); 1000 if(size(dec)==1) {printf("(%s terms)%n", 1, 0);} 1001 else {printf("(%s terms)%n", size(dec), 0);} 1035 1002 } 1036 1003 example … … 1057 1024 print(" "+getStringFraction(dec[1],q)); 1058 1025 for(int i=2; i<=size(dec); i++) 1059 { 1060 print("+ "+getStringFraction(dec[i],q)); 1061 } 1062 printf("(%s terms)%n", size(dec), 0); 1026 {print("+ "+getStringFraction(dec[i],q));} 1027 if(size(dec)==1) {printf("(%s terms)%n", 1, 0);} 1028 else {printf("(%s terms)%n", size(dec), 0);} 1063 1029 } 1064 1030 example … … 1077 1043 proc getStringpfd(list dec) 1078 1044 "USAGE: getStringpfd(dec); dec list 1079 PURPOSE: turn a partial fraction decomposition dec into one string. The list 1080 dec has to have the same structure as the output of @ref{pfd}. 1045 PURPOSE: turn a partial fraction decomposition @code{dec} into one string. The 1046 list @code{dec} has to have the same structure as the output of 1047 @ref{pfd}. 1081 1048 SEE ALSO: pfd, getStringpfd_indexed, displaypfd, displaypfd_long 1082 1049 EXAMPLE: example getStringpfd; shows an example" … … 1086 1053 string s = getStringFraction(dec[1],q); 1087 1054 for(int i=2; i<=size(dec); i++) 1088 { 1089 s = s+" + "+getStringFraction(dec[i],q); 1090 } 1055 {s = s+" + "+getStringFraction(dec[i],q);} 1091 1056 return(s); 1092 1057 } … … 1107 1072 proc getStringpfd_indexed(list dec) 1108 1073 "USAGE: getStringpfd_indexed(dec); dec list 1109 PURPOSE: turn a partial fraction decomposition dec into one string, writing the 1110 denominator factors just as q1,q2,... . The list dec has to have the 1111 same structure as the output of @ref{pfd}. 1074 PURPOSE: turn a partial fraction decomposition @code{dec} into one string, 1075 writing the denominator factors just as @code{q1},@code{q2},... . The 1076 list @code{dec} has to have the same structure as the output of 1077 @ref{pfd}. 1112 1078 SEE ALSO: pfd, getStringpfd, displaypfd, displaypfd_long 1113 1079 EXAMPLE: example getStringpfd_indexed; shows an example" … … 1116 1082 string s = getStringFraction_indexed(dec[1]); 1117 1083 for(int i=2; i<=size(dec); i++) 1118 { 1119 s = s+" + "+getStringFraction_indexed(dec[i]); 1120 } 1084 {s = s+" + "+getStringFraction_indexed(dec[i]);} 1121 1085 return(s); 1122 1086 } … … 1178 1142 PURPOSE: read matrix of rational functions from a txt-file and turn it into a 1179 1143 matrix (i.e. a list of lists) of pairs of polynomials (numerators and 1180 denominators). The string file is the [directory +] name of the file, 1181 r is the ring chosen for the numerator/denominator polynomials. 1182 @* The input file should be a list of lists using the brackets \"{\", 1183 \"}\" and \",\" as seperation, e.g. 1184 @* \"{{(x+y)/(x^2-x*y), -(x^2*y+1)/(y), x^2}, {(x+1)/y, y/x, 0}}\". 1144 denominators). The string @code{file} should be the [directory +] name 1145 of the file in the form \"@code{<path-to-file>/<filename>.txt}\". 1146 @* The input file should be a list of lists separated by the characters 1147 \"{\", \"}\" and \",\". Example: 1148 @* \"{{(x+y)/(x^2-x*y), -(x^2*y+1)/(y), x^2}, {(x+1)/y, y/x, 0}}\" 1149 @* Each rational function has to be an expression of the form \"a\", 1150 \"(a)/(b)\", \"(b)^(-n)\" or \"(a)*(b)^(-n)\" where \"a\",\"b\" stand 1151 for polynomials (i.e. strings, that can be parsed as a polynomial with 1152 the @code{execute} command) and \"n\" stands for a positive integer. A 1153 minus sign \"-\" followed by such an expression is also allowed. 1154 @* IMPORTANT: The strings \"a\",\"b\" must NOT contain the symbol \"/\". 1155 (So in case the coefficient field is the rationals, all denominators 1156 in the coefficients of numerator and denominator polynomials should be 1157 cleared.) 1185 1158 @* The file should contain less than 2^31 characters (filesize < 2 GB). 1186 For bigger files the matrix should be split row-wise in multiple 1187 matrices and saved in different files. A list of the filenames 1188 (in the right order) can then be given as first argument instead. 1159 For bigger files the matrix should be split row-wise into multiple 1160 matrices and saved in different files (each smaller than 2 GB). A list 1161 of the filenames (in the right order) can then be given as first 1162 argument instead. 1189 1163 @* Also the basering has to match the variable names used in the 1190 1164 input file(s). 1191 @* mode = 1: save result to an ssi file of the same name (default)1192 @* mode = 2: return result1193 @* mode = 3: both save to ssi file andreturn result1165 @* @code{mode=1} (default): save result to an ssi-file of the same name 1166 @* @code{mode=2}: return result 1167 @* @code{mode=3}: save to ssi-file AND return result 1194 1168 SEE ALSO: pfdMat" 1195 1169 { … … 1198 1172 if(!defined(basering)) 1199 1173 {ERROR("no basering defined!");} 1200 int left__,right__,pos1__, pos2__,tmp__,i__,j__,t__,tt__;1174 int left__,right__,pos1__,mid__,pos2__,tmp__,i__,j__,k__,t__,tt__,depth__; 1201 1175 int mode__=1; 1202 1176 if(size(#)>0) {mode__=#[1];} … … 1207 1181 for(i__=1;i__<=n;i__++) 1208 1182 { 1209 printf(" file %s of %s:",i__,n);1183 dbprint(sprintf(" file %s of %s:",i__,n)); 1210 1184 if(find(file__[i__],".txt")==0) {ERROR("wrong file type, expected txt");} 1185 printlevel = printlevel+1; 1211 1186 mat__ = mat__ + readInputTXT(file__[i__],2); 1187 printlevel = printlevel-1; 1212 1188 } 1213 1189 … … 1215 1191 1216 1192 string filename__ = file__[1][1,find(file__[1],".txt")-1]; 1217 print(" saving to file "+filename__+".ssi "); t__ = rtimer;1193 dbprint(" saving to file "+filename__+".ssi "); t__ = rtimer; 1218 1194 write("ssi:w "+filename__+".ssi", mat__); 1219 printf(" done! (%s ms)", rtimer-t__);1195 dbprint(sprintf(" done! (%s ms)", rtimer-t__)); 1220 1196 1221 1197 if(mode__==3) {return(mat__);} … … 1225 1201 if(find(file__,".txt")==0) {ERROR("wrong file type, expected txt");} 1226 1202 1227 print(" reading file "); t__=rtimer;1203 dbprint(" reading file "); t__=rtimer; 1228 1204 string data__ = read(":r "+file__); 1229 printf(" done! (%s ms)", rtimer-t__); 1230 1231 1232 print(" processing input "); t__ = rtimer; 1205 dbprint(sprintf(" done! (%s ms)", rtimer-t__)); 1206 1207 dbprint(" processing input "); t__ = rtimer; 1233 1208 left__ = find(data__,"{"); 1234 1209 right__ = find(data__,"}"); … … 1236 1211 if(left__<tmp__ && tmp__<right__) {left__ = tmp__;} 1237 1212 1238 int finished__ ;1213 int finished__,n__; 1239 1214 poly p__,q__; 1240 1215 list row__,mat__; 1216 string s__,ss__; 1241 1217 i__=0; 1242 1218 while(1) … … 1251 1227 { 1252 1228 j__++; 1229 s__ = ""; 1253 1230 1254 1231 pos1__ = pos2__+1; 1255 pos2__ = find(data__,"/",pos1__); 1256 tmp__ = find(data__,",",pos1__); 1257 if(pos2__==0 || pos2__>right__) //no denominator 1258 { 1259 if(tmp__==0 || tmp__>right__) {pos2__ = right__; finished__ = 1;} 1260 else {pos2__ = tmp__;} 1261 execute("p__=" + fixBrackets(data__[pos1__,pos2__-pos1__])); 1262 q__=1; 1263 } 1264 else{if(tmp__>0 && tmp__<pos2__) // no denominator 1265 { 1266 pos2__ = tmp__; 1267 execute("p__=" + fixBrackets(data__[pos1__,pos2__-pos1__])); 1268 q__=1; 1269 } 1270 else 1271 { 1272 execute("p__=" + fixBrackets(data__[pos1__,pos2__-pos1__])); 1273 1274 pos1__ = pos2__+1; 1275 pos2__ = tmp__; 1276 if(pos2__==0 || pos2__>right__) 1232 pos2__ = find(data__,",",pos1__); 1233 if(pos2__==0 || pos2__>right__) // end of row 1234 { 1235 finished__ = 1; 1236 pos2__ = right__; 1237 } 1238 s__ = data__[pos1__,pos2__-pos1__]; 1239 mid__ = find(s__,"/"); 1240 if(mid__==0) 1241 { 1242 tmp__ = find(s__,"^(-"); 1243 if(tmp__==0) //no denominator 1277 1244 { 1278 pos2__=right__; 1279 finished__=1; 1245 execute("p__=" + s__); 1246 q__=1; 1247 row__[j__] = list(p__,q__); 1248 continue; 1280 1249 } 1281 execute("q__=" + fixBrackets(data__[pos1__,pos2__-pos1__])); 1282 }} 1250 else // denominator is given by a negative exponent 1251 { 1252 if(find(s__,"^(-",tmp__+1)>0) 1253 {ERROR(sprintf("invalid syntax in (%s,%s)-th entry:" 1254 +" more than one negative exponent",i__,j__));} 1255 for(k__=tmp__+3; s__[k__]!=")"; k__++) {} 1256 execute("n__=" + s__[tmp__+3,k__-tmp__-3]); 1257 while(k__<size(s__)) 1258 { 1259 k__++; 1260 if(s__[k__]!=" ") 1261 {ERROR(sprintf("invalid syntax in (%s,%s)-th entry",i__,j__));} 1262 } 1263 s__ = s__[1,tmp__-1]; 1264 depth__=0; 1265 for(k__=tmp__-1; 1; k__--) 1266 { 1267 ss__ = s__[k__]; 1268 if(ss__ == ")") {depth__++;} 1269 else{if(ss__ == "(") {depth__--;}} 1270 if(depth__==0) {break;} 1271 } 1272 if(k__>1) 1273 { 1274 while(1) 1275 { 1276 if(k__==0) 1277 {ERROR(sprintf("invalid syntax in (%s,%s)-th entry",i__,j__));} 1278 k__--; 1279 ss__ = s__[k__]; 1280 if(ss__=="*" || ss__==" " || ss__=="-") {break;} 1281 } 1282 } 1283 s__ = s__[1,k__] + "1/" + s__[k__+1,size(s__)-k__] + "^" + string(n__); 1284 mid__ = k__+2; // position of the character "/" 1285 } 1286 } 1287 if(find(s__,"/",mid__+1)>0) 1288 {ERROR(sprintf("invalid syntax in (%s,%s)-th entry:%n" 1289 +"no '/' allowed in the string representing the polynomials",i__,j__,0));} 1290 1291 execute("p__=" + fixBrackets(s__[1,mid__-1])); 1292 execute("q__=" + fixBrackets(s__[mid__+1,size(s__)-mid__])); 1283 1293 1284 1294 row__[j__] = list(p__,q__); 1285 1295 } 1286 1287 mat__[i__] = row__; // append row to matrix 1288 printf(" row %s done! (%s ms)",i__,rtimer-tt__); 1296 mat__[i__] = row__; // append row to matrix 1297 dbprint(sprintf(" row %s done! (%s ms)",i__,rtimer-tt__)); 1289 1298 1290 1299 left__ = find(data__,"{",right__); … … 1292 1301 right__ = find(data__,"}",left__); 1293 1302 } 1294 printf(" done! (%s ms)", rtimer-t__);1303 dbprint(sprintf(" done! (%s ms)", rtimer-t__)); 1295 1304 1296 1305 if(mode__==2) {return(mat__);} 1297 1306 1298 1307 string filename__ = file__[1,find(file__,".txt")-1]; 1299 print(" saving to file "+filename__+".ssi "); t__ = rtimer;1308 dbprint(" saving to file "+filename__+".ssi "); t__ = rtimer; 1300 1309 write("ssi:w "+filename__+".ssi", mat__); 1301 printf(" done! (%s ms)", rtimer-t__);1310 dbprint(sprintf(" done! (%s ms)", rtimer-t__)); 1302 1311 1303 1312 if(mode__==3) {return(mat__);} 1304 1313 } 1305 /* new version:1306 {1307 system("--ticks-per-sec",1000);1308 if(!defined(basering))1309 {ERROR("no basering defined!");}1310 int left__,right__,pos1__,pos2__,tmp__,i__,j__,k__,d__,t__,tt__;1311 int mode__=1;1312 if(size(#)>0) {mode__=#[1];}1313 if(typeof(file__)=="list") // list of filenames given --> apply function to each1314 { // file and concatenate the resulting matrices1315 int n = size(file__);1316 list mat__ = list();1317 for(i__=1;i__<=n;i__++)1318 {1319 printf(" file %s of %s:",i__,n);1320 if(find(file__[i__],".txt")==0) {ERROR("wrong file type, expected txt");}1321 mat__ = mat__ + readInputTXT(file__[i__],2);1322 }1323 1324 if(mode__==2) {return(mat__);}1325 1326 string filename__ = file__[1][1,find(file__[1],".txt")-1];1327 print(" saving to "+filename__+".ssi "); t__ = rtimer;1328 write("ssi:w "+filename__+".ssi", mat__);1329 printf(" done! (%s ms)", rtimer-t__);1330 1331 if(mode__==3) {return(mat__);}1332 }1333 if(typeof(file__)!="string")1334 {ERROR("wrong type for first argument (expected string or list)");}1335 if(find(file__,".txt")==0) {ERROR("wrong file type, expected txt");}1336 1337 print(" reading input matrix as string from "+file__); t__=rtimer;1338 string data__ = read(":r "+file__);1339 printf(" done! (%s ms)", rtimer-t__);1340 1341 1342 print(" processing input "); t__ = rtimer;1343 left__ = find(data__,"{");1344 right__ = find(data__,"}");1345 tmp__ = find(data__,"{",left__+1);1346 if(left__<tmp__ && tmp__<right__) {left__ = tmp__;}1347 1348 int finished__;1349 poly p__,q__;1350 int depth__;1351 list mat__;1352 string s__,s1__,s2__;1353 for(i__=1;1;i__++)1354 {1355 tt__ = rtimer;1356 mat__[i__] = list();1357 pos2__ = left__;1358 finished__ = 0;1359 for(j__=1; not finished__; j__++)1360 {1361 pos1__ = pos2__+1;1362 pos2__ = find(data__,",",pos1__);1363 if(pos2__==0 || pos2__>right__) {finished__=1; pos2__=right__;}1364 s__ = data__[pos1__,pos2__-pos1__];1365 for(k__=1; s__[k__]==" "; k__++) {}1366 s__ = s__[k__,size(s__)-k__+1]; // remove spaces1367 1368 tmp__=find(s__,"/");1369 //printf("s: %s",s__);1370 if(tmp__>0)1371 {1372 depth__=0;1373 d__=0;1374 for(k__=tmp__+1;k__<=size(s__);k__++)1375 {1376 if(s__[k__]=="(") {depth__++; k__++; continue;}1377 if(s__[k__]==")") {depth__--; k__++; continue;}1378 if(s__[k__]=="/" && depth__<=d__) {tmp__=k__; d__=depth__;}1379 }1380 s1__ = s__[1,tmp__-1];1381 s2__ = s__[tmp__+1,size(s__)-tmp__];1382 depth__=0;1383 for(k__=size(s1__); k__>1; k__--)1384 {1385 if(s1__[k__]==")") {depth__++; k__--; continue;}1386 if(s1__[k__]=="(") {depth__--; k__--; continue;}1387 if(depth__==0 && (s1__[k__]=="+" || s1__[k__]=="-")) {tmp__=0; break}1388 }1389 depth__=0;1390 for(k__=1; k__<=size(s2__); k__++)1391 {1392 if(s2__[k__]=="(") {depth__++; k__++; continue;}1393 if(s2__[k__]==")") {depth__--; k__++; continue;}1394 if(depth__==0 && (s2__[k__]=="+" || s2__[k__]=="-" || s2__[k__]=="*")) {tmp__=0; break}1395 }1396 //printf("s1: %s,%ns2: %s",s1__,s2__);1397 }1398 if(tmp__==0) // no denominator1399 {1400 execute("p__=" + s__);1401 q__=1;1402 }1403 else // s1__ = numerator, s2__ = denominator1404 {1405 execute("p__=" + fixBrackets(s1__));1406 execute("q__=" + fixBrackets(s2__));1407 if(deg(q__)==0) {p__=p__/q__; q__=1;}1408 }1409 ;1410 1411 mat__[i__][j__] = list(p__,q__);1412 }1413 printf(" row %s done! (%s ms)",i__,rtimer-tt__);1414 1415 left__ = find(data__,"{",right__);1416 if(left__==0) {break;}1417 right__ = find(data__,"}",left__);1418 }1419 printf(" done! (%s ms)", rtimer-t__);1420 1421 if(mode__==2) {return(mat__);}1422 1423 string filename__ = file__[1,find(file__,".txt")-1];1424 print(" saving to "+filename__+".ssi "); t__ = rtimer;1425 write("ssi:w "+filename__+".ssi", mat__);1426 printf(" done! (%s ms)", rtimer-t__);1427 1428 if(mode__==3) {return(mat__);}1429 } */1430 1431 1314 1432 1315 static proc fixBrackets(string data) … … 1489 1372 file string, dotest,ignore_nonlin,output_mode,parallelize int 1490 1373 PURPOSE: apply @code{pfd} to all entries of a matrix of rational functions 1491 saved in a file. The string file is the [directory +] name ofthe1492 file.1374 saved in a txt-file. The string @code{file} should be the 1375 [directory +] name of the file. 1493 1376 @* The input file can either be a txt-file or an ssi-file created with 1494 1377 @code{readInputTXT}. In case of a txt-file, the base ring has to match 1495 1378 and the matrix has to be in the same format specified in 1496 @ref{readInputTXT}. 1497 Also, txt-files that are bigger than 2 GB should be split as described 1498 for @code{readInputTXT} and a list of the filenames can be given as 1499 first argument instead. 1500 @* The result is written in two txt-files: once with the denominators 1501 written out and once with indexed denominator factors q1,q2,... . 1502 The polynomials q1,q2,... are saved in a separate txt-file. 1379 @ref{readInputTXT}. Also, txt-files that are bigger than 2 GB should 1380 be split as described for @code{readInputTXT} and a list of the 1381 filenames can be given as first argument instead. 1382 @* The result is saved in multiple txt- (and ssi-) files (see below) 1383 within the directory of the input file. 1503 1384 @* Also a logfile is created, which protocols the memory used and the 1504 runtimes of for each matrix entry in real-time.1385 runtimes of @code{pfd} for each matrix entry in real-time. 1505 1386 1506 1387 @* There are also 4 optional arguments: 1507 1388 @* If @code{dotest} is nonzero, test the results with checkpfd: 1508 dotest<0: exact test (may be slow), 1509 dotest>0: do this amount of probabilistic tests for each entry 1510 (see @ref{checkpfd}). 1511 @* (default: @code{dotest=-1}) 1512 @* If @code{ignore_nonlin} is nonzero, for each denominator, the 1513 nonlinear denominator factors are saved into a seperate file and 1514 removed before applying @code{pfd}. The nonlinear factors thus do not 1515 appear in the output files. (So the input matrix is equal to the 1516 output matrix divided element-wise by the matrix containing the 1517 nonlinear factors.) 1518 (default: @code{ignore_nonlin=1}) 1519 @* If @code{output_mode} is nonzero, the input and result are also 1520 saved to ssi-files containing matrices (as list of lists) of the input 1521 rational functions (as lists of numerator and denominator poly) and 1522 decompositions respectively. For the ssi-file containing the result, 1523 the first entry is an ideal @code{q} of denominator factors and the 1524 second entry is a matrix (as list of lists) containing the 1525 decompositions, each of which is a list of terms, where a term is 1526 represented as in the result of @ref{pfd} by a list @code{l} containing 1527 @* 1) the numerator polynomial 1528 @* 2) an intvec giving the indices @code{i} for which @code{q[i]} 1529 occurs as a factor in the denominator 1530 @* 3) an intvec containing the exponents of those irreducible factors. 1531 @* Also, the results of @code{pfd} will be saved directly in ssi-files 1532 named pfd_results_i_j where i,j are the matrix indices in the current 1533 directory. (default: @code{output_mode=0}) 1534 @* If @code{parallelize} is nonzero, the decompositions are calculated in 1535 parallel using @ref{parallel_lib}. (default: @code{parallelize=1}) 1389 @* @code{dotest<0} (default): exact test (may be slow), 1390 @* @code{dotest>0}: do this amount of probabilistic tests for each entry 1391 (see @ref{checkpfd}). 1392 1393 @* If @code{ignore_nonlin} is nonzero (default), for each denominator, 1394 the nonlinear factors in the factorization are removed before applying 1395 @code{pfd} (and added back in in the output files). 1396 1397 @* If @code{parallelize} is nonzero (default), the decompositions are 1398 calculated in parallel using @ref{parallel_lib}. 1399 1400 @* The parameter @code{output_mode} controls the output files created: 1401 @* @code{output_mode=1} (default): The result consists of two files: 1402 @code{<filename>_pfd_indexed.txt} contains the matrix of all 1403 decompositions (as list of lists separated by the characters \"{\", 1404 \"}\" and \",\") where all the denominators are written in factorized 1405 form depending on irreducible factors @code{q1}, @code{q2}, ... . 1406 The file @code{<filename>_denominator_factors.txt} lists all the 1407 polynomials @code{q1}, @code{q2}, ... . 1408 @* @code{output_mode=2}: Additionally to mode 1, the file 1409 @code{<filename>_pfd.txt} is created, which also contains the matrix 1410 of decompositions but the factors in the denominators are written out. 1411 @* @code{output_mode=3}: Additionally to mode 2, the result and some 1412 intermediate results are saved as SINGULAR objects in ssi-files: 1413 @* @code{<filename>.ssi}: contains the result of @code{readInputTXT} in 1414 case a txt-file was given as input. 1415 @* @code{<filename>_factorized_denominators.ssi}: like the first file, 1416 but the denominators are saved in factorized form, that is as a list 1417 of an ideal of irreducible non constant polynomials and an intvec of 1418 exponents. 1419 @* @code{<filename>_linear_part.ssi} (only if @code{ignore_nonlin} is 1420 nonzero): like the previous file, but all the irreducible denominator 1421 factors are removed 1422 @* @code{<filename>_non_linear_factors.ssi} (only if @code{ignore_nonlin} 1423 is nonzero): a list of an ideal @code{p} generated by irreducible 1424 polynomials and a matrix (list of lists) of the nonlinear denominator 1425 factors of each entry of the input matrix. These are represented as 1426 lists of an intvec of indices @code{i} for which @code{p[i]} occurs 1427 as a (nonlinear) factor in the denominator and an intvec containing 1428 the exponents of those factors. 1429 @* @code{<filename>_pfd.ssi}: a list, where the first entry is an ideal 1430 @code{q} of denominator factors and the second entry is a matrix (as 1431 list of lists) containing the decompositions, each of which is a list 1432 of terms, where a term is represented as in the result of @ref{pfd} 1433 by a list containing 1434 @* 1) the numerator polynomial 1435 @* 2) an intvec of indices @code{i} for which @code{q[i]} occurs 1436 as a factor in the denominator 1437 @* 3) an intvec containing the exponents of those irreducible factors. 1438 @* IMPORTANT: If @code{ignore_nonlin} is nonzero, this file contains the 1439 decompositions of the entries of the matrix in 1440 @code{<filename>_linear_part.ssi}. Thus the nonlinear factors, are 1441 NOT contained in this file. 1442 @* @code{output_mode=4}: Additionally to mode 3, the direct output of 1443 each call of @code{pfd} is saved in separate ssi-files called 1444 @code{pfd_results_i_j.ssi} where i,j are the matrix indices. This 1445 creates a lot of files, but may be useful in case the algorithm does 1446 not terminate in time for every matrix entry. Other than the files 1447 created in mode 1-3, these files are saved in the current directory, 1448 rather than the directory of the input file. 1536 1449 SEE ALSO: readInputTXT, pfd, checkpfd, checkpfdMat" 1537 1450 { 1538 1451 system("--ticks-per-sec",1000); 1539 1452 short = 0; 1540 int dotest,ignore_nonlin,output_mode,parallelize = -1,1, 0,1;1453 int dotest,ignore_nonlin,output_mode,parallelize = -1,1,1,1; 1541 1454 if(size(#)>0) {dotest = #[1];} 1542 1455 if(size(#)>1) {ignore_nonlin = #[2];} … … 1547 1460 list arguments,results; 1548 1461 1549 print(newline+"reading data "); int t0=rtimer;1462 dbprint(newline+"reading data "); int t0=rtimer; 1550 1463 1551 1464 if(typeof(infile)=="list") 1552 1465 { 1466 printlevel = printlevel+1; 1553 1467 if(output_mode>2) {list mat = readInputTXT(infile,3);} 1554 1468 else {list mat = readInputTXT(infile,2);} 1469 printlevel = printlevel-1; 1555 1470 int pos=find(infile[1],".txt"); 1556 1471 string filename = infile[1][1,pos-1]; … … 1561 1476 if(pos!=0) 1562 1477 { 1478 printlevel = printlevel+1; 1563 1479 if(output_mode>2) {list mat = readInputTXT(infile,3);} 1564 1480 else {list mat = readInputTXT(infile,2);} 1481 printlevel = printlevel-1; 1565 1482 } 1566 1483 else 1567 1484 { 1568 1485 pos=find(infile,".ssi"); 1569 if(pos!=0) 1570 { 1571 list mat = read("ssi:r "+infile); 1572 } 1486 if(pos!=0) {list mat = read("ssi:r "+infile);} 1573 1487 else {ERROR("invalid file type, expected ssi or txt");} 1574 1488 } … … 1579 1493 int n = size(mat); 1580 1494 int m = size(mat[1]); 1581 printf("done! (%s ms)",rtimer-t0); 1582 1583 if(typeof(mat[1][1][2])!="list") // denominators are not yet factorized 1584 { 1585 print("factorizing the denominators "); t0=rtimer; 1495 dbprint(sprintf("done! (%s ms)",rtimer-t0)); 1496 1497 if(typeof(mat[1][1][2])!="list") // denominators are not yet factorized 1498 { 1499 dbprint("factorizing the denominators "); t0=rtimer; 1500 printlevel = printlevel+1; 1586 1501 mat = FactDenom(mat); 1502 printlevel = printlevel-1; 1587 1503 if(output_mode>2) 1588 1504 {write("ssi:w "+filename+"_factorized_denominators.ssi",mat);} 1589 printf("done! (%s ms)",rtimer-t0);1505 dbprint(sprintf("done! (%s ms)",rtimer-t0)); 1590 1506 } 1591 1507 1592 1508 if(ignore_nonlin) 1593 1509 { 1594 print("removing nonlinear denominator factors before pfd is applied"); 1595 // +" (%s_non_linear_factors[_indexed].txt):", filename); 1596 list denom_factors; 1510 dbprint("removing nonlinear denominator factors before pfd is applied"); 1511 list nonlin_denom_factors; 1597 1512 ideal p; 1598 mat, denom_factors, p = removeNonlinearFactors(mat, filename); 1513 printlevel = printlevel+1; 1514 mat, nonlin_denom_factors, p = removeNonlinearFactors(mat, filename); 1515 printlevel = printlevel-1; 1599 1516 if(output_mode>2) 1600 1517 { 1601 print("saving nonlinear factors to "+filename+"_non_linear_factors.ssi "); t0 = rtimer; 1602 write("ssi:w "+filename+"_non_linear_factors.ssi", denom_factors); 1603 printf("done! (%s ms)",rtimer-t0); 1604 print("saving input matrix without the nonlinear factors to "+filename+"_linear_part.ssi "); t0 = rtimer; 1518 dbprint("saving nonlinear factors to "+filename+"_non_linear_factors.ssi "); 1519 t0 = rtimer; 1520 write("ssi:w "+filename+"_non_linear_factors.ssi", list(p,nonlin_denom_factors)); 1521 dbprint(sprintf("done! (%s ms)",rtimer-t0)); 1522 dbprint("saving input matrix without the nonlinear factors to " 1523 +filename+"_linear_part.ssi "); 1524 t0 = rtimer; 1605 1525 write("ssi:w "+filename+"_linear_part.ssi", mat); 1606 printf("done! (%s ms)",rtimer-t0); 1607 } 1608 //ideal p = saveNonlinearFactorsTXT(denom_factors, filename); 1609 //filename = filename + "_linear_part"; 1526 dbprint(sprintf("done! (%s ms)",rtimer-t0)); 1527 } 1610 1528 } 1611 1529 1612 1530 if(parallelize) 1613 1531 { 1614 print("creating tasks "); t0=rtimer;1532 dbprint("creating tasks "); t0=rtimer; 1615 1533 write(":w "+filename+"_pfdMat_logfile.txt","finished matrix entries with runtimes " 1616 1534 +"(calculated in parallel on "+string(getcores())+" cores):"); … … 1624 1542 } 1625 1543 } 1626 printf("done! (%s ms)",rtimer-t0);1627 1628 print("applying pfd to each matrix entry "); t0 = rtimer;1544 dbprint(sprintf("done! (%s ms)",rtimer-t0)); 1545 1546 dbprint("applying pfd to each matrix entry "); t0 = rtimer; 1629 1547 results = parallelWaitAll("pfdWrap",arguments); 1630 1548 arguments = list(); 1631 printf("done! (%s ms)",rtimer-t0);1549 dbprint(sprintf("done! (%s ms)",rtimer-t0)); 1632 1550 1633 1551 write(logfile,"decomposition: "+string(rtimer-t0)+" ms and "+string(memory(2)) 1634 1552 +" Byte Memory max. (after calling pfd on each matrix entry)"); 1635 1553 1636 print("writing results in matrix shape "); t0 = rtimer;1554 dbprint("writing results in matrix shape "); t0 = rtimer; 1637 1555 list dec_mat; 1638 1556 for(i=1; i<=n; i++) … … 1646 1564 } 1647 1565 results = list(); 1648 printf("done! (%s ms)",rtimer-t0);1566 dbprint(sprintf("done! (%s ms)",rtimer-t0)); 1649 1567 } 1650 1568 else 1651 1569 { 1652 print("applying pfd to each matrix entry "); t0=rtimer;1570 dbprint("applying pfd to each matrix entry "); t0=rtimer; 1653 1571 write(":w "+filename+"_pfdMat_logfile.txt", 1654 1572 "finished matrix entries with runtimes (no parallelization):"); … … 1663 1581 } 1664 1582 } 1665 printf("done! (%s ms)",rtimer-t0);1583 dbprint(sprintf("done! (%s ms)",rtimer-t0)); 1666 1584 write(logfile,"decomposition: "+string(rtimer-t0)+" ms and "+string(memory(2)) 1667 1585 +" Byte Memory max. (after calling pfd on each matrix entry)"); 1668 1586 } 1669 1587 1670 print("making one single list of denominator factors "); t0 = rtimer;1588 dbprint("making one single list of denominator factors "); t0 = rtimer; 1671 1589 ideal q,new_q; 1672 1590 intvec dict; … … 1691 1609 { 1692 1610 if(size(dec_mat[i][j][k][2])>0) 1693 { 1694 dec_mat[i][j][k][2] = intvec(dict[dec_mat[i][j][k][2]]); 1695 } 1696 } 1697 } 1698 printf(" row %s complete!",i); 1699 } 1700 printf("done! (%s ms)",rtimer-t0); 1611 {dec_mat[i][j][k][2] = intvec(dict[dec_mat[i][j][k][2]]);} 1612 } 1613 } 1614 dbprint(sprintf(" row %s complete!",i)); 1615 } 1616 dbprint(sprintf("done! (%s ms)",rtimer-t0)); 1701 1617 1702 1618 if(output_mode>2) 1703 1619 { 1704 print("saving result to "+filename+"_pfd_only_linear_factors.ssi "); t0 = rtimer;1705 write("ssi:w "+filename+"_pfd _only_linear_factors.ssi", list(q,dec_mat));1706 printf("done! (%s ms)",rtimer-t0);1620 dbprint("saving result to "+filename+"_pfd.ssi "); t0 = rtimer; 1621 write("ssi:w "+filename+"_pfd.ssi", list(q,dec_mat)); 1622 dbprint(sprintf("done! (%s ms)",rtimer-t0)); 1707 1623 } 1708 1624 … … 1714 1630 { 1715 1631 for(j=1; j<=m; j++) 1716 1717 denom_factors[i][j][1] =denom_factors[i][j][1]+ind; // adjust indices1718 1632 { 1633 nonlin_denom_factors[i][j][1] = nonlin_denom_factors[i][j][1]+ind; // adjust indices 1634 } 1719 1635 } 1720 1636 } 1721 1637 1722 1638 if(ignore_nonlin) 1723 { printf("creating readable .txt-files (including the nonlinear factors again)");}1639 {dbprint(sprintf("creating readable .txt-files (including the nonlinear factors again)"));} 1724 1640 else 1725 { printf("creating readable .txt-files ");}1641 {dbprint(sprintf("creating readable .txt-files "));} 1726 1642 t0 = rtimer; 1727 print(" indexed ("+filename+"_pfd_indexed.txt):"); 1643 dbprint(" indexed ("+filename+"_pfd_indexed.txt):"); 1644 printlevel = printlevel+1; 1728 1645 if(ignore_nonlin) 1729 {saveResultTXT_indexed(dec_mat, filename+"_pfd_indexed", denom_factors);}1646 {saveResultTXT_indexed(dec_mat, filename+"_pfd_indexed", nonlin_denom_factors);} 1730 1647 else {saveResultTXT_indexed(dec_mat, filename+"_pfd_indexed");} 1648 printlevel = printlevel-1; 1731 1649 for(i=1; i<=size(q); i++) {fprintf(qfile, "q%s = %s;", i, q[i]);} 1732 1650 if(output_mode>1) 1733 1651 { 1734 print(" denominators written out ("+filename+"_pfd.txt):"); 1652 dbprint(" denominators written out ("+filename+"_pfd.txt):"); 1653 printlevel = printlevel+1; 1735 1654 if(ignore_nonlin) 1736 {saveResultTXT(dec_mat, q, filename+"_pfd", denom_factors);}1655 {saveResultTXT(dec_mat, q, filename+"_pfd", nonlin_denom_factors);} 1737 1656 else {saveResultTXT(dec_mat, q, filename+"_pfd");} 1738 } 1739 printf("done! (%s ms)",rtimer-t0); 1657 printlevel = printlevel-1; 1658 } 1659 dbprint(sprintf("done! (%s ms)",rtimer-t0)); 1740 1660 1741 1661 if(dotest) 1742 1662 { 1743 1663 if(dotest<0) 1744 { print("checking for correctness (exact test) ");}1664 {dbprint("checking for correctness (exact test) ");} 1745 1665 else 1746 { printf("checking for correctness (%s random evaluations per entry) ",dotest);}1666 {dbprint(sprintf("checking for correctness (%s random evaluations per entry) ",dotest));} 1747 1667 t0 = rtimer; 1748 1668 if(parallelize) … … 1767 1687 } 1768 1688 } 1769 1770 printf("%s out of %s = %sx%s decompositions are correct! (%s ms)%n", 1771 sum(results),n*m,n,m,rtimer-t0,0); 1689 dbprint(sprintf("%s out of %s = %sx%s decompositions are correct! (%s ms)%n", 1690 sum(results),n*m,n,m,rtimer-t0,0)); 1772 1691 write(logfile,"checking for correctness: "+string(rtimer-t0)+" ms and " 1773 1692 +string(memory(2))+" Byte Memory max. (at the end of pfdMat), " … … 1802 1721 { 1803 1722 t = rtimer; 1804 for(i=1; i<=n; i++) // create argument list1723 for(i=1; i<=n; i++) // create argument list 1805 1724 { 1806 1725 for(j=1; j<=m; j++) 1807 { 1808 arguments[m*(i-1)+j] = list(denom[i][j]); 1809 } 1726 {arguments[m*(i-1)+j] = list(denom[i][j]);} 1810 1727 } 1811 1728 … … 1813 1730 arguments = list(); 1814 1731 1815 for(i=1; i<=n; i++) // update q, mat, denom1732 for(i=1; i<=n; i++) // update q, mat, denom 1816 1733 { 1817 1734 for(j=1; j<=m; j++) … … 1882 1799 } 1883 1800 } 1884 1885 1801 timeout = timeout*2; 1886 1887 printf(" %s out of %s denominators factorized completely (%s ms)",counter,n*m,rtimer-t);1802 dbprint(sprintf(" %s out of %s denominators factorized completely (%s ms)", 1803 counter,n*m,rtimer-t)); 1888 1804 } 1889 1805 return(mat); … … 1924 1840 k = size(s); 1925 1841 s[k-1] = "}"; s[k] = ","; s = s+" "; 1926 printf(" row %s done! (%s ms)",i,rtimer-t);1842 dbprint(sprintf(" row %s done! (%s ms)",i,rtimer-t)); 1927 1843 } 1928 1844 k = size(s); … … 1965 1881 k = size(s); 1966 1882 s[k-1] = "}"; s[k] = ","; s = s+" "; 1967 printf(" row %s done! (%s ms)",i,rtimer-t);1883 dbprint(sprintf(" row %s done! (%s ms)",i,rtimer-t)); 1968 1884 } 1969 1885 k = size(s); … … 1978 1894 int i,j,k,t0,ind; 1979 1895 ideal p; 1980 list denom_factors;1896 list nonlin_denom_factors; 1981 1897 intvec factors,exponents; 1982 1898 list fac; … … 1984 1900 { 1985 1901 t0 = rtimer; 1986 denom_factors[i] = list();1902 nonlin_denom_factors[i] = list(); 1987 1903 for(j=1; j<=m; j++) 1988 1904 { … … 2000 1916 for(ind=1;1;ind++) 2001 1917 { 2002 //print(" "+string(ind));2003 1918 if(p[ind]==fac[1][k]) {break;} 2004 1919 if(ind==size(p)) … … 2018 1933 } 2019 1934 fractions[i][j][2] = fac; 2020 2021 denom_factors[i][j] = list(factors,exponents); 2022 } 2023 printf(" row %s done! (%s ms)",i,rtimer-t0); 2024 } 2025 return(fractions, denom_factors, p); 1935 nonlin_denom_factors[i][j] = list(factors,exponents); 1936 } 1937 dbprint(sprintf(" row %s done! (%s ms)", i, rtimer-t0)); 1938 } 1939 return(fractions, nonlin_denom_factors, p); 2026 1940 } 2027 1941 … … 2029 1943 "USAGE: checkpfdMat(input, output, denomFactors[, N, parallelize]); 2030 1944 input,output,denomFactors string, N,parallelize int 2031 checkpfdMat(input, output, denomFactors, nonlinFactors[, N, parallelize]);2032 input,output,denomFactors,nonlinFactors string, N,parallelize int2033 1945 PURPOSE: test the output files of @code{pfdMat} for correctness. Input and 2034 output (indexed) txt-file names are given as strings. If nonlinear2035 denominator factors were ignored in the calculation (as described in2036 @ref{pfdMat}), the txt-file containing the nonlinear factors must be2037 given as fourth argument @code{nonlinFactors}. The output and2038 nonlinear Factors should be indexed and @code{denomFactors} is the2039 file containing the denominator factors @code{q1},@code{q2},....1946 output (indexed) txt-files have to be given as strings in the form 1947 \"@code{<path-to-file>/<filename>.txt}\". The output should be indexed 1948 (that is the output file ending in @code{..._pfd_indexed.txt}) and 1949 @code{denomFactors} has to be the file containing the denominator 1950 factors @code{q1}, @code{q2}, ... (the txt-file ending in 1951 @code{..._denominator_factors.txt}). 2040 1952 @* As for @code{readInputTXT} and @code{pfdMat}, the basering has to 2041 1953 match the variable names used in the input file, which has to be in … … 2043 1955 than 2 GB have to be split as described for @code{readInputTXT} and a 2044 1956 list of filenames can be given as first argument instead. 2045 @* If a n integer N is given, the test is done probabilistically by N2046 random substitutions for each entry of the matrix. If N is omitted,2047 the fractions in the decomposition will be expanded symbolically and2048 compared to the input (may be slower).2049 @* If @code{parallelize} is nonzero , the tests are run in parallel using2050 @ref{parallel_lib}. (default: @code{parallelize=1})1957 @* If a positive integer N is given, the test is done probabilistically by 1958 evaluation at N random points for each entry of the matrix. If N is 1959 nonpositive (default), the fractions in the decompositions will be 1960 expanded symbolically and compared to the input (may be slower). 1961 @* If @code{parallelize} is nonzero (default), the tests are run in 1962 parallel using @ref{parallel_lib}. 2051 1963 @* The result is printed and as in @code{pfdMat} a logfile is created 2052 1964 showing the results for each matrix entry. … … 2071 1983 if(size(#)>2) {ERROR("too many arguments");} 2072 1984 2073 print(newline+"reading input file:"); 1985 dbprint(newline+"reading input file:"); 1986 printlevel = printlevel+1; 2074 1987 list frac = readInputTXT(input,2); 1988 printlevel = printlevel-1; 2075 1989 if(typeof(input)=="string") {string filename=input;} 2076 1990 if(typeof(input)=="list") {string filename=input[1];} 2077 1991 filename = filename[1,find(filename,".txt")-1]; 2078 1992 2079 print("factorizing the denominators "); t0=rtimer; 1993 dbprint("factorizing the denominators "); t0=rtimer; 1994 printlevel = printlevel+1; 2080 1995 frac = FactDenom(frac); 2081 printf("done! (%s ms)",rtimer-t0); 2082 2083 print("reading output files:"); t0=rtimer; 2084 print(" reading list of denominator factors from "+qfile); 1996 printlevel = printlevel-1; 1997 dbprint(sprintf("done! (%s ms)",rtimer-t0)); 1998 1999 dbprint("reading output files:"); t0=rtimer; 2000 dbprint(" reading list of denominator factors from "+qfile); 2085 2001 ideal q; 2086 2002 q = readQfileTXT(qfile); 2087 //printf("denominator factors: %s",q); 2088 print(" done!"); 2089 2090 print(" reading (indexed) output decompositions "); 2003 dbprint(" done!"); 2004 2005 dbprint(" reading (indexed) output decompositions "); 2091 2006 list dec,nonlin; 2007 printlevel = printlevel+1; 2092 2008 dec,nonlin = readOutputTXT(output); 2009 printlevel = printlevel-1; 2093 2010 2094 2011 if(parallelize) 2095 { printf("done! (%s ms)%n%ncreating tasks",rtimer-t0,0); t0=rtimer;}2012 {dbprint(sprintf("done! (%s ms)%n%ncreating tasks",rtimer-t0,0)); t0=rtimer;} 2096 2013 else 2097 2014 { 2098 printf("done! (%s ms)%n%n",rtimer-t0,0);2015 dbprint(sprintf("done! (%s ms)%n",rtimer-t0,0)); 2099 2016 if(N<=0) 2100 { print("checking for correctness (exact test) ");}2017 {dbprint("checking for correctness (exact test) ");} 2101 2018 else 2102 { printf("checking for correctness (%s random evaluations per entry) ",N);}2019 {dbprint(sprintf("checking for correctness (%s random evaluations per entry) ",N));} 2103 2020 t0=rtimer; 2104 2021 } 2105 2022 2106 fprintf(":w "+filename+"_checkpfd _logfile.txt","Input file (matrix of rational functions): %s"2107 +" %nOutput file (decompositions): %s%nlist of all denominator factors: %s"2108 +" %n%nResults of checkpfdMat:",input,output,qfile,0);2109 link logfile = ":a "+filename+"_checkpfd _logfile.txt";2023 fprintf(":w "+filename+"_checkpfdMat_logfile.txt","Input file (matrix of rational functions):" 2024 +" %s%nOutput file (decompositions): %s%nlist of all denominator factors:" 2025 +" %s%n%nResults of checkpfdMat:",input,output,qfile,0); 2026 link logfile = ":a "+filename+"_checkpfdMat_logfile.txt"; 2110 2027 2111 2028 int n=size(frac); … … 2135 2052 if(parallelize) 2136 2053 { 2137 printf("done! (%s ms)%n%n",rtimer-t0,0);2054 dbprint(sprintf("done! (%s ms)%n",rtimer-t0,0)); 2138 2055 if(N<=0) 2139 { print("checking for correctness (exact test) ");}2056 {dbprint("checking for correctness (exact test) ");} 2140 2057 else 2141 { printf("checking for correctness (%s random evaluations per entry) ",N);}2058 {dbprint(sprintf("checking for correctness (%s random evaluations per entry) ",N));} 2142 2059 t0=rtimer; 2143 2060 list results = parallelWaitAll("testEntry",arguments); 2144 2061 } 2145 else { printf("done! (%s ms)",rtimer-t0);}2146 printf("%s out of %s = %sx%s decompositions are correct! (%s ms)%n",2147 sum(results),n*m,n,m,rtimer-t0,0) ;2062 else {dbprint(sprintf("done! (%s ms)",rtimer-t0));} 2063 dbprint(sprintf("%s out of %s = %sx%s decompositions are correct! (%s ms)%n", 2064 sum(results),n*m,n,m,rtimer-t0,0)); 2148 2065 fprintf(logfile,"%s out of %s = %sx%s decompositions are correct! (%s ms)%n", 2149 2066 sum(results),n*m,n,m,rtimer-t0,0); … … 2153 2070 { 2154 2071 system("--ticks-per-sec",1000); 2155 print(" reading matrix of decompositions from file "+filename__);2072 dbprint(" reading matrix of decompositions from file "+filename__); 2156 2073 int t__=rtimer; 2157 2074 int tt__; 2158 2075 string data__ = read(":r "+filename__); 2159 printf(" done! (%s ms)", rtimer-t__);2160 2161 print(" processing input "); t__ = rtimer;2076 dbprint(sprintf(" done! (%s ms)", rtimer-t__)); 2077 2078 dbprint(" processing input "); t__ = rtimer; 2162 2079 list mat__; 2163 2080 list nonlin__=list(); … … 2190 2107 mat__[i__][j__] = list(); 2191 2108 2192 2193 2109 factors__ = intvec(0:0); 2194 2110 exponents__ = intvec(0:0); … … 2196 2112 if(l__>0) 2197 2113 { 2198 ss__ = s__[1,l__-1]; //ss__ contains the nonlinear factors2114 ss__ = s__[1,l__-1]; //ss__ contains the nonlinear factors 2199 2115 s__ = s__[l__+3,size(s__)-l__-2]; 2200 2116 for(p1__=1;s__[p1__]!="(";p1__++) {} … … 2212 2128 p2__ = find(ss__,"^",p1__); 2213 2129 tmp__ = find(ss__,"*",p1__); 2214 if((p2__>tmp__ && tmp__>0) || (p2__==0)) //exponent is 12130 if((p2__>tmp__ && tmp__>0) || (p2__==0)) //exponent is 1 2215 2131 { 2216 2132 execute("factors__[l__]="+ss__[p1__,tmp__-p1__]); … … 2238 2154 if((s__[k__]=="+" && depth__==0) || k__==max__) 2239 2155 { 2240 if(tmp2__==0) // no denominator2156 if(tmp2__==0) // no denominator 2241 2157 { 2242 2158 execute("numerator__="+s__[tmp__,k__-tmp__]); … … 2249 2165 p1__ = find(ss__,"("); 2250 2166 p2__ = find(ss__,")"); 2251 ss__ = ss__[p1__+1,p2__-p1__-1]; // now ss__ is only the denominator2167 ss__ = ss__[p1__+1,p2__-p1__-1]; // now ss__ is only the denominator 2252 2168 2253 2169 ss__ = ss__+"*"; … … 2262 2178 p2__ = find(ss__,"^",p1__); 2263 2179 tmp__ = find(ss__,"*",p1__); 2264 if((p2__>tmp__ && tmp__>0) || (p2__==0)) //exponent is 12180 if((p2__>tmp__ && tmp__>0) || (p2__==0)) //exponent is 1 2265 2181 { 2266 2182 execute("factors__[l__]="+ss__[p1__,tmp__-p1__]); … … 2279 2195 } 2280 2196 } 2281 printf(" row %s done! (%s ms)",i__,rtimer-tt__);2282 } 2283 printf(" done! (%s ms)", rtimer-t__);2197 dbprint(sprintf(" row %s done! (%s ms)",i__,rtimer-tt__)); 2198 } 2199 dbprint(sprintf(" done! (%s ms)", rtimer-t__)); 2284 2200 2285 2201 return(mat__,nonlin__); … … 2291 2207 ideal q__; 2292 2208 int pos1__,pos2__=1,1; 2293 2294 2209 while(1) 2295 2210 { -
Tst/Short/pfd.tst
r6924d4 rcc07066 6 6 example pfd; 7 7 8 / * pfd(poly,poly) */8 // pfd(poly,poly) ////////////////////////////////////////////////////////////// 9 9 ring r1 = 0, x(1..5), dp; 10 10 poly f = -2*x(1)*x(3)+3*x(3)*x(4)+x(2)*x(5)+x(3)*x(5)-x(4)*x(5); … … 32 32 kill dec; 33 33 34 35 // pfd(poly,list) (denominator = ideal of factors & intvec of exponents) /////// 34 36 ring r3 = 5, (x,y,z), dp; 35 poly f = 2*x^3-x^2*y+x*y^2+y^3-2*x^2*z-2*x*y*z-y^2*z+2*x*z^2+y*z^2-2*z^3; 36 poly g = (x^2+3*y)^2*(x-2*y^2)*(x+y+1)^3*(2*x+3*y+4)*(x^2-x*y+y^2); 37 poly f = x+y+z+1; 38 list g = list(ideal((x^2+y^2+z^2),(x+y^2),(y+z^2),(z+x^2)), intvec(2,1,1,1)); 39 list dec = pfd(f,g); 40 displaypfd(dec); 41 checkpfd(list(f,g), dec); 42 checkpfd(list(f,g), dec, 10); 43 kill dec; 44 45 // different ordering, same polynomials: 46 ring r4 = 5, (x,y,z), lp; 47 poly f = fetch(r3,f); 48 list g = fetch(r3,g); 37 49 list dec = pfd(f,g); 38 50 displaypfd(dec); … … 41 53 kill dec; 42 54 43 ring r4 = 5, (x,y,z), lp;44 poly f = fetch(r3,f);45 poly g = fetch(r3,g);46 list dec = pfd(f,g);47 displaypfd(dec);48 checkpfd(list(f,g), dec);49 checkpfd(list(f,g), dec,10);50 kill dec;51 55 52 / * pfd(list) */56 // pfd(list) /////////////////////////////////////////////////////////////////// 53 57 ring s1 = 0, (x,y,z), dp; 54 58 poly f1 = x*y+y*z+z*x-x-y-z+1; … … 68 72 kill dec; 69 73 70 / * pfd(matrix) */74 // pfd(matrix) ///////////////////////////////////////////////////////////////// 71 75 ring s2 = 3, (x,y,z), dp; 72 76 poly f11 = (x+y+z+1)^3; … … 85 89 displaypfd(dec[2][1]); 86 90 displaypfd(dec[2][2]); 87 checkpfd(list(f1 1,g11), dec[1][1]);88 checkpfd(list(f 12,g12), dec[1][2]);89 checkpfd(list(f 21,g21), dec[2][1]);90 checkpfd(list(f 22,g22), dec[2][2]);91 checkpfd(list(f1,g1), dec[1][1]); 92 checkpfd(list(f2,g2), dec[1][2]); 93 checkpfd(list(f3,g3), dec[2][1]); 94 checkpfd(list(f4,g4), dec[2][2]); 91 95 kill dec; 92 96
Note: See TracChangeset
for help on using the changeset viewer.