1 | //////////////////////////////////////////////////////////////////// |
---|
2 | version="version pfd.lib 4.1.3.2 Aug_2020 "; |
---|
3 | category="??"; |
---|
4 | info=" |
---|
5 | LIBRARY: pfd.lib Multivariate Partial Fraction Decomposition |
---|
6 | |
---|
7 | AUTHOR: Marcel Wittmann, e-mail: mwittman@mathematik.uni-kl.de |
---|
8 | |
---|
9 | OVERVIEW: |
---|
10 | This Library implements an algorithm based on the work of E. K. Leinartas to |
---|
11 | write rational functions in mutiple variables as a sum of functions with |
---|
12 | \"smaller\" numerators and denominators. |
---|
13 | This can be used to shorten the IBP reduction coeffcients of multi-loop Feynman |
---|
14 | integrals as shown in [J. Boehm, M. Wittmann, Z. Wu, Y. Xu, Y. Zhang: 'IBP |
---|
15 | reduction coefficients made simple'] (preprint 2020). For this application, |
---|
16 | we also provide a procedure that applies the algorithm to all entries of a |
---|
17 | matrix of rational functions given as one (possibly very big) txt-file. |
---|
18 | |
---|
19 | KEYWORDS: partial fraction; decomposition; Leinartas |
---|
20 | |
---|
21 | PROCEDURES: |
---|
22 | pfd(); calculate a partial fraction decomposition |
---|
23 | of a rational function |
---|
24 | checkpfd(); test if a decomposition is equal to a rational |
---|
25 | function given by numerator/denominator polynomials |
---|
26 | evaluatepfd(); substitute values in a partial fraction |
---|
27 | decomposition gotten from @code{pfd} |
---|
28 | displaypfd(); print a decomposition gotten as output of @code{pfd} |
---|
29 | displaypfd_long(); like @code{display}, but denominators are written out |
---|
30 | getStringpfd(); turn a decomposition gotten from @code{pfd} into one |
---|
31 | string |
---|
32 | getStringpfd_indexed(); like @code{getStringpfd}, but writes the denominator |
---|
33 | factors just as @code{q1}, @code{q2}, ... |
---|
34 | readInputTXT(); read a matrix of rational functions from a txt-file |
---|
35 | pfdMat(); apply @code{pfd} to a matrix of rational functions |
---|
36 | in parallel (using @ref{parallel_lib}) and save result |
---|
37 | as easy-to-read txt-files. |
---|
38 | checkpfdMat(); test output files of @code{pfdMat} for correctness |
---|
39 | "; |
---|
40 | ///////////////////////////////////////////////////////////////////////////// |
---|
41 | |
---|
42 | LIB "random.lib"; |
---|
43 | LIB "ring.lib"; |
---|
44 | LIB "parallel.lib"; |
---|
45 | LIB "elim.lib"; |
---|
46 | LIB "poly.lib"; |
---|
47 | |
---|
48 | ///////////////////////////////////////////////////////////////////////////// |
---|
49 | |
---|
50 | static proc mod_init() |
---|
51 | { |
---|
52 | printlevel = 2; |
---|
53 | // export some static functions so they can be run in parallel using parallelWaitAll: |
---|
54 | if (!defined(Tasks)) {LIB "tasks.lib";} |
---|
55 | exportto(Tasks,pfdWrap); |
---|
56 | importfrom(Tasks,pfdWrap); |
---|
57 | exportto(Pfd,pfdWrap); |
---|
58 | exportto(Tasks,testEntry); |
---|
59 | importfrom(Tasks,testEntry); |
---|
60 | exportto(Pfd,testEntry); |
---|
61 | } |
---|
62 | |
---|
63 | proc pfd(list #) |
---|
64 | "USAGE: pfd(f,g[,debug]); f,g poly, debug int |
---|
65 | pfd(f,g[,debug]); f poly, g list, debug int |
---|
66 | pfd(arguments[, parallelize]); arguments list, parallelize int |
---|
67 | RETURN: a partial fraction decomposition of f/g as a list @code{l} where |
---|
68 | @code{l[1]} is an ideal generated by irreducible polynomials and |
---|
69 | @code{l[2]} is a list of fractions. |
---|
70 | Each fraction is represented by a list of |
---|
71 | @* 1) the numerator polynomial |
---|
72 | @* 2) an intvec of indices @code{i} for which @code{l[1][i]} occurs |
---|
73 | as a factor in the denominator |
---|
74 | @* 3) an intvec containing the exponents of those irreducible factors. |
---|
75 | @* Setting @code{debug} to a positive integer measures runtimes and |
---|
76 | creates a log file (default: @code{debug=0}). |
---|
77 | @* The denominator g can also be given in factorized form as a list of |
---|
78 | an ideal of irreducible non constant polynomials and an intvec of |
---|
79 | exponents. This can save time since the first step in the algorithm is |
---|
80 | to factorize g. (A list of the zero-ideal and an empty intvec |
---|
81 | represents a denominator of 1.) |
---|
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 @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. |
---|
88 | NOTE: The result depends on the monomial ordering. For \"small\" results |
---|
89 | use @code{dp}. |
---|
90 | EXAMPLE: example pfd; shows an example |
---|
91 | SEE ALSO: checkpfd, evaluatepfd, displaypfd, displaypfd_long, pfdMat" |
---|
92 | { |
---|
93 | short = 0; |
---|
94 | int i; |
---|
95 | if(typeof(#[1])=="list") |
---|
96 | { |
---|
97 | if(size(#)>1) |
---|
98 | { |
---|
99 | if(typeof(#[2])=="list") |
---|
100 | {list arguments = #; int parallelize = 1;} |
---|
101 | else{if(typeof(#[2])=="int") |
---|
102 | {list arguments = #[1]; int parallelize = #[2];} |
---|
103 | else {ERROR("wrong type for second argument, expected int");}} |
---|
104 | } |
---|
105 | if(parallelize) |
---|
106 | { |
---|
107 | for(i=1; i<=size(arguments); i++) |
---|
108 | { |
---|
109 | if(typeof(arguments[i][1])=="list") //input is list of lists |
---|
110 | {arguments[i] = list(arguments[i],1);} |
---|
111 | } |
---|
112 | return(parallelWaitAll("pfd",arguments)); |
---|
113 | } |
---|
114 | else |
---|
115 | { |
---|
116 | list results; |
---|
117 | for(i=1; i<=size(arguments); i++) |
---|
118 | { |
---|
119 | if(typeof(arguments[i][1])=="list") //input is list of lists |
---|
120 | {results[i] = pfd(arguments[i],0);} |
---|
121 | else |
---|
122 | { |
---|
123 | if(size(arguments[i])==2) |
---|
124 | {results[i] = pfd(arguments[i][1],arguments[i][2]);} |
---|
125 | else{if(size(arguments[i])==3) |
---|
126 | {results[i] = pfd(arguments[i][1],arguments[i][2],arguments[i][3]);} |
---|
127 | else {ERROR("wrong number of arguments, expected 2 or 3");}} |
---|
128 | } |
---|
129 | } |
---|
130 | return(results); |
---|
131 | } |
---|
132 | } |
---|
133 | poly f = #[1]; |
---|
134 | if(typeof(#[2])=="list") |
---|
135 | { |
---|
136 | list g=#[2]; |
---|
137 | } |
---|
138 | else |
---|
139 | { |
---|
140 | poly g=#[2]; |
---|
141 | } |
---|
142 | |
---|
143 | int debug=0; link l=":w "; |
---|
144 | if(size(#)>3) {ERROR("wrong number of arguments, expected 2 or 3");} |
---|
145 | if(size(#)==3) |
---|
146 | { |
---|
147 | debug=#[3]; |
---|
148 | l=":a "+string(debug)+"_log_"+datetime()+".txt"; |
---|
149 | system("--ticks-per-sec",1000); |
---|
150 | } |
---|
151 | if(debug) |
---|
152 | { |
---|
153 | fprintf(l,"debug: %s", debug); |
---|
154 | fprintf(l,"size(string(f)) = %s, size(string(g)) = %s %n", |
---|
155 | size(string(f)), size(string(g)), 0); |
---|
156 | int counter,tt,ttt; |
---|
157 | } |
---|
158 | |
---|
159 | if(f==0) |
---|
160 | { |
---|
161 | list dec = list(ideal(),list(list(poly(0),intvec(0:0),intvec(0:0)))); |
---|
162 | if(debug) |
---|
163 | {fprintf(l,"%ntotal: 0 ms (numerator was 0)",0); close(l);} |
---|
164 | if(voice<=printlevel) {displaypfd(dec);} |
---|
165 | return(dec); |
---|
166 | } |
---|
167 | if(typeof(g)=="poly") |
---|
168 | { |
---|
169 | if(deg(g)==0) |
---|
170 | { |
---|
171 | list dec = list(ideal(),list(list(f/g,intvec(0:0),intvec(0:0)))); |
---|
172 | if(debug) |
---|
173 | {fprintf(l,"%ntotal: 0 ms (denominator was constant)",0); close(l);} |
---|
174 | if(voice<=printlevel) {displaypfd(dec);} |
---|
175 | return(dec); |
---|
176 | } |
---|
177 | |
---|
178 | // (1) factorization of the denominator //////////////////////////////////// |
---|
179 | if(debug) {int t1 = rtimer; write(l,"factorizing ");} |
---|
180 | list factor = factorize(g); |
---|
181 | number lcoeff; |
---|
182 | for(i=2; i<=size(factor[1]); i++) |
---|
183 | { |
---|
184 | lcoeff = leadcoef(factor[1][i]); |
---|
185 | factor[1][i] = factor[1][i]/lcoeff; |
---|
186 | factor[1][1] = factor[1][1]*(lcoeff^factor[2][i]); // polynomial is monic (thus unique) |
---|
187 | lcoeff = content(factor[1][i]); |
---|
188 | factor[1][i] = factor[1][i]/lcoeff; |
---|
189 | factor[1][1] = factor[1][1]*(lcoeff^factor[2][i]); // polynomial has nice coefficients |
---|
190 | } |
---|
191 | ideal q = factor[1]; |
---|
192 | f = f/q[1]; |
---|
193 | q=delete(q,1); |
---|
194 | intvec e = factor[2]; e=delete(e,1); |
---|
195 | int m = size(q); |
---|
196 | if(debug) {t1 = rtimer-t1; fprintf(l,"done! (%s ms)", t1);} |
---|
197 | } |
---|
198 | else{if(typeof(g)=="list") |
---|
199 | { |
---|
200 | if(size(g[1])==0) |
---|
201 | { |
---|
202 | list dec = list(ideal(),list(list(f,intvec(0:0),intvec(0:0)))); |
---|
203 | if(debug) |
---|
204 | {fprintf(l,"%ntotal: 0 ms (denominator was constant)",0); close(l);} |
---|
205 | if(voice<=printlevel) {displaypfd(dec);} |
---|
206 | return(dec); |
---|
207 | } |
---|
208 | |
---|
209 | // denominator is already factorized |
---|
210 | for(i=1;i<=size(g[1]);i++) |
---|
211 | { |
---|
212 | if(size(factorize(g[1][i])[1])>2) |
---|
213 | {ERROR("factors should be irreducible");} |
---|
214 | } |
---|
215 | |
---|
216 | ideal q = g[1]; |
---|
217 | intvec e = g[2]; |
---|
218 | int m = size(q); |
---|
219 | if(debug) {int t1=0;} |
---|
220 | } |
---|
221 | else |
---|
222 | {ERROR("wrong type for second argument, expected poly or list(ideal,intvec)");}} |
---|
223 | |
---|
224 | // (2) Nullstellensatz decomposition ///////////////////////////////////////// |
---|
225 | if(debug) |
---|
226 | {int t2 = rtimer; write(l,"Nullstellensatz decomposition ");} |
---|
227 | list terms = list(list(poly(f),1..m,e)); |
---|
228 | list dec,newterms,result; |
---|
229 | int imax; |
---|
230 | while(size(terms)>0) |
---|
231 | { |
---|
232 | if(debug) {ttt = rtimer; counter++;} |
---|
233 | imax = size(terms); |
---|
234 | newterms = list(); |
---|
235 | for(i=1; i<=imax; i++) |
---|
236 | { |
---|
237 | result = NSSdecompStep(terms[i],q); |
---|
238 | if(size(result)>1) {newterms = mergepfd(newterms, result);} |
---|
239 | else {dec = mergepfd(dec, result);} |
---|
240 | } |
---|
241 | terms = newterms; |
---|
242 | if(debug) |
---|
243 | {fprintf(l," %s: %s ms, %s terms, %s unfinished", |
---|
244 | counter, rtimer-ttt, size(terms)+size(dec), size(terms));} |
---|
245 | } |
---|
246 | if(debug) |
---|
247 | {t2 = rtimer-t2; fprintf(l,"done! (%s ms, %s terms)", t2, size(dec));} |
---|
248 | |
---|
249 | // (3) short numerator decomposition ///////////////////////////////////////// |
---|
250 | if(debug) |
---|
251 | {int t3 = rtimer; write(l,"short numerator decompositions "); counter=0;} |
---|
252 | terms = dec; |
---|
253 | dec = list(); |
---|
254 | while(size(terms)>0) |
---|
255 | { |
---|
256 | if(debug) {tt = rtimer; counter++;} |
---|
257 | imax = size(terms); |
---|
258 | newterms = list(); |
---|
259 | for(i=1; i<=imax; i++) |
---|
260 | { |
---|
261 | result = shortNumeratorDecompStep(terms[i],q,debug,l); |
---|
262 | |
---|
263 | if(debug) {ttt = rtimer;} |
---|
264 | newterms = mergepfd(newterms, result[1]); |
---|
265 | dec = mergepfd(dec, result[2]); |
---|
266 | if(debug) |
---|
267 | {fprintf(l," merging: %s ms, %s new terms", |
---|
268 | rtimer-ttt, size(result[1]));} |
---|
269 | } |
---|
270 | terms = newterms; |
---|
271 | if(debug) |
---|
272 | {fprintf(l," %s: %s ms, %s terms, %s unfinished", |
---|
273 | counter, rtimer-tt, size(terms)+size(dec), size(terms));} |
---|
274 | } |
---|
275 | if(debug) |
---|
276 | {t3 = rtimer-t3; fprintf(l,"done! (%s ms, %s terms)", t3, size(dec));} |
---|
277 | |
---|
278 | // (4) algebraic dependence decomposition //////////////////////////////////// |
---|
279 | if(debug) |
---|
280 | {int t4 = rtimer; write(l,"algebraic dependence decomposition "); counter=0;} |
---|
281 | terms = dec; |
---|
282 | dec = list(); |
---|
283 | while(size(terms)>0) |
---|
284 | { |
---|
285 | if(debug) {tt = rtimer; counter++;} |
---|
286 | imax = size(terms); |
---|
287 | newterms = list(); |
---|
288 | for(i=1; i<=imax; i++) |
---|
289 | { |
---|
290 | result = algDependDecompStep(terms[i],q,debug,l); |
---|
291 | |
---|
292 | if(debug) {ttt = rtimer;} |
---|
293 | if(size(result)>1) {newterms = mergepfd(newterms, result);} |
---|
294 | else {dec = mergepfd(dec, result);} |
---|
295 | if(debug) |
---|
296 | {fprintf(l," merging: %s ms, %s new terms", rtimer-ttt, size(result));} |
---|
297 | } |
---|
298 | terms = newterms; |
---|
299 | if(debug) |
---|
300 | {fprintf(l," %s: %s ms, %s terms, %s unfinished", |
---|
301 | counter, rtimer-tt, size(terms)+size(dec),size(terms));} |
---|
302 | } |
---|
303 | if(debug) |
---|
304 | {t4 = rtimer-t4; fprintf(l,"done! (%s ms, %s terms)", t4, size(dec));} |
---|
305 | |
---|
306 | // (5) numerator decomposition /////////////////////////////////////////////// |
---|
307 | if(debug) |
---|
308 | {int t5 = rtimer; write(l,"numerator decompositions "); counter=0;} |
---|
309 | terms = dec; |
---|
310 | dec = list(); |
---|
311 | while(size(terms)>0) |
---|
312 | { |
---|
313 | if(debug) {tt = rtimer; counter++;} |
---|
314 | imax = size(terms); |
---|
315 | newterms = list(); |
---|
316 | for(i=1; i<=imax; i++) |
---|
317 | { |
---|
318 | result = numeratorDecompStep(terms[i],q,debug,l); |
---|
319 | |
---|
320 | if(debug) {ttt = rtimer;} |
---|
321 | newterms = mergepfd(newterms, result[1]); |
---|
322 | dec = mergepfd(dec, result[2]); |
---|
323 | if(debug) |
---|
324 | {fprintf(l," merging: %s ms, %s new terms", rtimer-ttt, size(result[1]));} |
---|
325 | } |
---|
326 | terms = newterms; |
---|
327 | if(debug) |
---|
328 | {fprintf(l," %s: %s ms, %s terms, %s unfinished", |
---|
329 | counter, rtimer-tt, size(terms)+size(dec), size(terms));} |
---|
330 | } |
---|
331 | if(debug) |
---|
332 | {t5 = rtimer-t5; fprintf(l,"done! (%s ms, %s terms)", t5, size(dec));} |
---|
333 | |
---|
334 | dec = list(q,dec); |
---|
335 | if(debug) {fprintf(l,"%ntotal: %s ms", t1+t2+t3+t4+t5); close(l);} |
---|
336 | if(voice<=printlevel) {displaypfd(dec);} |
---|
337 | return(dec); |
---|
338 | } |
---|
339 | example |
---|
340 | { |
---|
341 | "EXAMPLE:"; |
---|
342 | echo=voice; |
---|
343 | ring R = 0,(x,y),dp; |
---|
344 | poly f = x^3+3*x^2*y+2*y^2-x^2+4*x*y; |
---|
345 | poly g = x^2*y*(x-1)*(x-y)^2; |
---|
346 | |
---|
347 | list dec = pfd(f,g); |
---|
348 | |
---|
349 | displaypfd_long(dec); // display result |
---|
350 | checkpfd(list(f,g),dec); // check for equality to f/g |
---|
351 | |
---|
352 | // calculate decompositions of a 2x2 matrix of rational functions at once: |
---|
353 | list arguments = list(list(f, g), list(1, f) ), |
---|
354 | list(list(x*y, y+1), list(1, x^2-y^2)); |
---|
355 | |
---|
356 | dec = pfd(arguments); |
---|
357 | |
---|
358 | // the result has the same shape as the |
---|
359 | // input (2x2 matrix as list of lists): |
---|
360 | displaypfd_long(dec[1][1]); |
---|
361 | displaypfd_long(dec[1][2]); |
---|
362 | displaypfd_long(dec[2][1]); |
---|
363 | displaypfd_long(dec[2][2]); |
---|
364 | |
---|
365 | // a more complicated example |
---|
366 | ring S = 0,(s12,s15,s23,s34,s45),dp; |
---|
367 | poly f = 7*s12^4*s15^2 + 11*s12^3*s15^3 + 4*s12^2*s15^4 - 10*s12^4*s15*s23 |
---|
368 | - 14*s12^3*s15^2*s23 - 4*s12^2*s15^3*s23 + 3*s12^4*s23^2 + 3*s12^3*s15*s23^2 |
---|
369 | + 13*s12^4*s15*s34 + 12*s12^3*s15^2*s34 + 2*s12^2*s15^3*s34 |
---|
370 | - 5*s12^4*s23*s34 + 33*s12^3*s15*s23*s34 + 49*s12^2*s15^2*s23*s34 |
---|
371 | + 17*s12*s15^3*s23*s34 - 17*s12^3*s23^2*s34 - 19*s12^2*s15*s23^2*s34 |
---|
372 | - 5*s12*s15^2*s23^2*s34 - 24*s12^3*s15*s34^2 - 15*s12^2*s15^2*s34^2 |
---|
373 | + 2*s12*s15^3*s34^2 + 15*s12^3*s23*s34^2 - 34*s12^2*s15*s23*s34^2 |
---|
374 | - 31*s12*s15^2*s23*s34^2 + 2*s15^3*s23*s34^2 + 33*s12^2*s23^2*s34^2 |
---|
375 | + 29*s12*s15*s23^2*s34^2 + 5*s15^2*s23^2*s34^2 + 9*s12^2*s15*s34^3 |
---|
376 | - 4*s12*s15^2*s34^3 - 15*s12^2*s23*s34^3 + 9*s12*s15*s23*s34^3 |
---|
377 | - 4*s15^2*s23*s34^3 - 27*s12*s23^2*s34^3 - 13*s15*s23^2*s34^3 |
---|
378 | + 2*s12*s15*s34^4 + 5*s12*s23*s34^4 + 2*s15*s23*s34^4 + 8*s23^2*s34^4 |
---|
379 | - 6*s12^3*s15^2*s45 - 9*s12^2*s15^3*s45 - 2*s12*s15^4*s45 |
---|
380 | + 30*s12^3*s15*s23*s45 + 56*s12^2*s15^2*s23*s45 + 24*s12*s15^3*s23*s45 |
---|
381 | - 12*s12^3*s23^2*s45 - 23*s12^2*s15*s23^2*s45 - 10*s12*s15^2*s23^2*s45 |
---|
382 | - 30*s12^3*s15*s34*s45 - 32*s12^2*s15^2*s34*s45 - 6*s12*s15^3*s34*s45 |
---|
383 | + 7*s12^3*s23*s34*s45 - 86*s12^2*s15*s23*s34*s45 - 104*s12*s15^2*s23*s34*s45 |
---|
384 | - 15*s15^3*s23*s34*s45 + 41*s12^2*s23^2*s34*s45 + 51*s12*s15*s23^2*s34*s45 |
---|
385 | + 10*s15^2*s23^2*s34*s45 - 5*s12^3*s34^2*s45 + 33*s12^2*s15*s34^2*s45 |
---|
386 | + 14*s12*s15^2*s34^2*s45 - 2*s15^3*s34^2*s45 - 21*s12^2*s23*s34^2*s45 |
---|
387 | + 62*s12*s15*s23*s34^2*s45 + 28*s15^2*s23*s34^2*s45 - 46*s12*s23^2*s34^2*s45 |
---|
388 | - 28*s15*s23^2*s34^2*s45 + 10*s12^2*s34^3*s45 - s12*s15*s34^3*s45 |
---|
389 | + 4*s15^2*s34^3*s45 + 21*s12*s23*s34^3*s45 - 6*s15*s23*s34^3*s45 |
---|
390 | + 17*s23^2*s34^3*s45 - 5*s12*s34^4*s45 - 2*s15*s34^4*s45 - 7*s23*s34^4*s45 |
---|
391 | - 6*s12^2*s15^2*s45^2 - 5*s12*s15^3*s45^2 - 2*s15^4*s45^2 |
---|
392 | - 28*s12^2*s15*s23*s45^2 - 42*s12*s15^2*s23*s45^2 - 10*s15^3*s23*s45^2 |
---|
393 | + 9*s12^2*s23^2*s45^2 + 10*s12*s15*s23^2*s45^2 + 24*s12^2*s15*s34*s45^2 |
---|
394 | + 36*s12*s15^2*s34*s45^2 + 10*s15^3*s34*s45^2 - 11*s12^2*s23*s34*s45^2 |
---|
395 | + 31*s12*s15*s23*s34*s45^2 + 25*s15^2*s23*s34*s45^2 |
---|
396 | - 18*s12*s23^2*s34*s45^2 - 10*s15*s23^2*s34*s45^2 + 4*s12^2*s34^2*s45^2 |
---|
397 | - 29*s12*s15*s34^2*s45^2 - 17*s15^2*s34^2*s45^2 + 27*s12*s23*s34^2*s45^2 |
---|
398 | + 2*s15*s23*s34^2*s45^2 + 9*s23^2*s34^2*s45^2 - 3*s12*s34^3*s45^2 |
---|
399 | + 10*s15*s34^3*s45^2 - 16*s23*s34^3*s45^2 - s34^4*s45^2 + 6*s12*s15^2*s45^3 |
---|
400 | + 3*s15^3*s45^3 + 8*s12*s15*s23*s45^3 + 10*s15^2*s23*s45^3 |
---|
401 | - 8*s12*s15*s34*s45^3 - 10*s15^2*s34*s45^3 + 9*s12*s23*s34*s45^3 |
---|
402 | + s12*s34^2*s45^3 + 8*s15*s34^2*s45^3 - 9*s23*s34^2*s45^3 - s34^3*s45^3 |
---|
403 | - s15^2*s45^4 + s15*s34*s45^4; |
---|
404 | poly g = 4*s12*s15*(s12 + s15 - s34)*(s15 - s23 - s34)*(s12 + s23 - s45) |
---|
405 | *(s12 - s34 - s45)*(s12 + s15 - s34 - s45)*s45; |
---|
406 | |
---|
407 | list dec = pfd(f,g); |
---|
408 | |
---|
409 | displaypfd(dec); |
---|
410 | checkpfd(list(f,g),dec); |
---|
411 | |
---|
412 | // size comparison: |
---|
413 | size(string(f)) + size(string(g)); |
---|
414 | size(getStringpfd(dec)); |
---|
415 | } |
---|
416 | |
---|
417 | static proc NSSdecompStep(list l, ideal q) |
---|
418 | { |
---|
419 | poly f=l[1]; |
---|
420 | intvec indices=l[2]; |
---|
421 | intvec e=l[3]; |
---|
422 | int m = size(indices); |
---|
423 | |
---|
424 | if(m==0) // denominator is 1 |
---|
425 | {return(list(l));} // do nothing, return input |
---|
426 | |
---|
427 | ideal qe = q[indices]; |
---|
428 | for(int i=1; i<=m; i++) |
---|
429 | {qe[i] = qe[i]^e[i];} |
---|
430 | matrix T; |
---|
431 | ideal qe_std = liftstd(qe,T); |
---|
432 | |
---|
433 | if(deg(qe_std) == 0) |
---|
434 | { |
---|
435 | T = T/qe_std[1]; |
---|
436 | // now 1 = T[1,1]*qe[1] + ... + T[m,1]*qe[m] is a Nullstellensatz certificate |
---|
437 | list dec; |
---|
438 | poly h; |
---|
439 | for(i=1; i<=m; i++) |
---|
440 | { |
---|
441 | h = T[i,1]; |
---|
442 | if(h != 0) |
---|
443 | {dec[size(dec)+1] = list(f*h,delete(indices,i),delete(e,i));} |
---|
444 | } |
---|
445 | return(dec); |
---|
446 | } |
---|
447 | else |
---|
448 | { |
---|
449 | return(list(l)); // do nothing, return input |
---|
450 | } |
---|
451 | } |
---|
452 | |
---|
453 | static proc shortNumeratorDecompStep(list l, ideal q, list #) |
---|
454 | { |
---|
455 | int debug=0; |
---|
456 | if(size(#)>0) {debug=#[1]; link ll=#[2];} |
---|
457 | if(debug) {system("--ticks-per-sec",1000); int tt=rtimer;} |
---|
458 | |
---|
459 | poly f=l[1]; |
---|
460 | intvec indices=l[2]; |
---|
461 | intvec e=l[3]; |
---|
462 | int m = size(indices); |
---|
463 | |
---|
464 | if(m==0) // denominator is 1 |
---|
465 | { |
---|
466 | if(debug) |
---|
467 | {fprintf(ll," shortNumeratorDecompStep: %s ms (m=%s, e=%s) " |
---|
468 | +"--> 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 |
---|
473 | matrix T; |
---|
474 | ideal q_std = liftstd(q_denom,T); |
---|
475 | list divrem = division(f,q_std); |
---|
476 | poly r = divrem[2][1]/divrem[3][1,1]; |
---|
477 | |
---|
478 | if(r!=0) |
---|
479 | { |
---|
480 | if(debug) |
---|
481 | {fprintf(ll," shortNumeratorDecompStep: %s ms (m=%s, e=%s) " |
---|
482 | +"--> remainder is nonzero", rtimer-tt, m, e);} |
---|
483 | return(list(list(),list(l))); // if there is a rest, decomposing further would |
---|
484 | } // not help in the next step (alg. depend. decomposition) |
---|
485 | |
---|
486 | matrix a = divrem[1]/divrem[3][1,1]; // now f = r + a[1,1]*q_std[1] + ... +a[m,1]*q_std[m] |
---|
487 | a = T*a; // lift coefficients ==> now f = r + a[1,1]*q[1] + ... +a[m,1]*q[m] |
---|
488 | |
---|
489 | // reduce w.r.t. groebner basis of syz(q) to make the numerators "smaller": |
---|
490 | vector v; |
---|
491 | for(int i=1; i<=m; i++) {v = v + gen(i)*a[i,1];} |
---|
492 | v = reduce(v, std(syz(q_denom))); |
---|
493 | |
---|
494 | list fraction,dec; |
---|
495 | for(i=1; i<=m; i++) |
---|
496 | { |
---|
497 | if(v[i] == 0) |
---|
498 | { |
---|
499 | i++; |
---|
500 | continue; |
---|
501 | } |
---|
502 | fraction[1] = v[i]; |
---|
503 | if(e[i]==1) |
---|
504 | { |
---|
505 | fraction[2] = delete(indices,i); |
---|
506 | fraction[3] = delete(e,i); |
---|
507 | } |
---|
508 | else |
---|
509 | { |
---|
510 | fraction[2] = indices; |
---|
511 | fraction[3] = e; |
---|
512 | fraction[3][i] = fraction[3][i] - 1; |
---|
513 | } |
---|
514 | dec[size(dec)+1] = fraction; |
---|
515 | } |
---|
516 | |
---|
517 | if(debug) |
---|
518 | {fprintf(ll," shortNumeratorDecompStep: %s ms (m=%s, e=%s, deg(v)=%s, size(v)=%s)", |
---|
519 | rtimer-tt, m, e, deg(v), size(v));} |
---|
520 | return(list(dec,list())); |
---|
521 | } |
---|
522 | |
---|
523 | static proc algDependDecompStep(list l, ideal q, list #) |
---|
524 | { |
---|
525 | int debug=0; |
---|
526 | if(size(#)>0) {debug=#[1]; link ll=#[2];} |
---|
527 | if(debug) {system("--ticks-per-sec",1000); int tt=rtimer;} |
---|
528 | def br = basering; |
---|
529 | int d = nvars(br); |
---|
530 | intvec indices=l[2]; |
---|
531 | int m = size(indices); |
---|
532 | intvec e=l[3]; |
---|
533 | int i; |
---|
534 | |
---|
535 | if(m==0) // denominator is 1 |
---|
536 | { |
---|
537 | if(debug) |
---|
538 | {fprintf(ll," algDependDecompStep: %s ms (m=%s, e=%s)", rtimer-tt, m, e);} |
---|
539 | return(list(l)); // do nothing, return input |
---|
540 | } |
---|
541 | |
---|
542 | if(m<=d) |
---|
543 | { |
---|
544 | if(size(syz(module(transpose(jacob(ideal(q[indices]))))))==0) // jacobian criterion |
---|
545 | { |
---|
546 | if(debug) |
---|
547 | {fprintf(ll," algDependDecompStep: %s ms (m=%s, e=%s) " |
---|
548 | +"--> alg. indep.", rtimer-tt, m, e);} |
---|
549 | return(list(l)); // do nothing, return input |
---|
550 | } |
---|
551 | } |
---|
552 | |
---|
553 | def R = changeord(list(list("dp",m+d)),extendring(m, "y(", "dp", 1, changevar("x()",br))); |
---|
554 | setring(R); |
---|
555 | |
---|
556 | list l = fetch(br,l); |
---|
557 | ideal q = fetch(br,q); |
---|
558 | poly f=l[1]; |
---|
559 | |
---|
560 | ideal I; |
---|
561 | for(i=1; i<=m; i++) |
---|
562 | {I[i] = y(i)-q[indices[i]];} |
---|
563 | |
---|
564 | ideal annihilatingPolys = eliminate(I,intvec(1..d)); |
---|
565 | |
---|
566 | poly g = annihilatingPolys[1]; |
---|
567 | |
---|
568 | poly tail = g[size(g)]; // term of lowest dp-order (thus lowest degree) |
---|
569 | number tcoeff = leadcoef(tail); |
---|
570 | intvec texpon = leadexp(tail); |
---|
571 | texpon = texpon[(d+1)..(d+m)]; |
---|
572 | g = g-tail; |
---|
573 | poly term; |
---|
574 | number coeff; |
---|
575 | intvec expon; |
---|
576 | list fraction,dec; |
---|
577 | int pow; |
---|
578 | int jmax = size(g); |
---|
579 | for(int j=1; j<=jmax; j++) |
---|
580 | { |
---|
581 | term = g[j]; |
---|
582 | coeff = leadcoef(term); |
---|
583 | expon = leadexp(term); |
---|
584 | expon = expon[(d+1)..(d+m)]; |
---|
585 | fraction[1] = -f*coeff/tcoeff; |
---|
586 | fraction[2] = intvec(0:0); |
---|
587 | fraction[3] = intvec(0:0); |
---|
588 | for(i=1; i<=m; i++) |
---|
589 | { |
---|
590 | pow = expon[i]-texpon[i]-e[i]; |
---|
591 | if(pow>=0) |
---|
592 | { |
---|
593 | fraction[1] = fraction[1]*q[indices[i]]^pow; |
---|
594 | } |
---|
595 | else |
---|
596 | { |
---|
597 | fraction[2][size(fraction[2])+1] = indices[i]; |
---|
598 | fraction[3][size(fraction[3])+1] = -pow; |
---|
599 | } |
---|
600 | } |
---|
601 | dec[size(dec)+1] = fraction; |
---|
602 | } |
---|
603 | setring(br); |
---|
604 | list dec = fetch(R,dec); |
---|
605 | if(debug) |
---|
606 | {fprintf(ll," algDependDecompStep: %s ms (m=%s, e=%s)", rtimer-tt, m, e);} |
---|
607 | return(dec); |
---|
608 | } |
---|
609 | |
---|
610 | static proc numeratorDecompStep(list l, ideal q, list #) |
---|
611 | { |
---|
612 | int debug=0; |
---|
613 | if(size(#)>0) {debug=#[1]; link ll=#[2];} |
---|
614 | if(debug) {system("--ticks-per-sec",1000); int tt=rtimer;} |
---|
615 | |
---|
616 | poly f=l[1]; |
---|
617 | intvec indices=l[2]; |
---|
618 | intvec e=l[3]; |
---|
619 | int m = size(indices); |
---|
620 | |
---|
621 | if(m==0) // denominator is 1 |
---|
622 | { |
---|
623 | if(debug) |
---|
624 | {fprintf(ll," numeratorDecompStep: %s ms (m=%s, e=%s) " |
---|
625 | +"--> constant denominator", rtimer-tt, m, e);} |
---|
626 | return(list(list(),list(l))); // do nothing, return input |
---|
627 | } |
---|
628 | |
---|
629 | ideal q_denom = q[indices]; // factors in the denominator |
---|
630 | matrix T; |
---|
631 | ideal q_std = liftstd(q_denom,T); |
---|
632 | list divrem = division(f,q_std); |
---|
633 | matrix a = divrem[1]/divrem[3][1,1]; |
---|
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] |
---|
635 | a = T*a; // lift coefficients ==> now f = r + a[1,1]*q[1] + ... +a[m,1]*q[m] |
---|
636 | |
---|
637 | // reduce w.r.t. groebner basis of syz(q) to make the numerators "smaller" |
---|
638 | vector v; |
---|
639 | for(int i=1; i<=m; i++) {v = v + gen(i)*a[i,1];} |
---|
640 | v = reduce(v, std(syz(q_denom))); |
---|
641 | |
---|
642 | list fraction,dec,rest; |
---|
643 | if(r!=0) |
---|
644 | {rest[1] = list(r,indices,e);} |
---|
645 | for(i=1; i<=m; i++) |
---|
646 | { |
---|
647 | if(v[i] == 0) |
---|
648 | { |
---|
649 | i++; |
---|
650 | continue; |
---|
651 | } |
---|
652 | fraction[1] = v[i]; |
---|
653 | if(e[i]==1) |
---|
654 | { |
---|
655 | fraction[2] = delete(indices,i); |
---|
656 | fraction[3] = delete(e,i); |
---|
657 | } |
---|
658 | else |
---|
659 | { |
---|
660 | fraction[2] = indices; |
---|
661 | fraction[3] = e; |
---|
662 | fraction[3][i] = fraction[3][i] - 1; |
---|
663 | } |
---|
664 | dec[size(dec)+1] = fraction; |
---|
665 | } |
---|
666 | |
---|
667 | if(debug) |
---|
668 | {fprintf(ll," numeratorDecompStep: %s ms (m=%s, e=%s, deg(v)=%s, size(v)=%s)", |
---|
669 | rtimer-tt, m, e, deg(v), size(v));} |
---|
670 | return(list(dec,rest)); |
---|
671 | } |
---|
672 | |
---|
673 | static proc mergepfd(list dec1, list dec2) |
---|
674 | { |
---|
675 | // Note: assumes dec1 is already sorted w.r.t. dp in the denominator exponents |
---|
676 | int n1=size(dec1); |
---|
677 | int n2=size(dec2); |
---|
678 | if(n2==0) {return(dec1);} |
---|
679 | int i; |
---|
680 | int a,b,m; |
---|
681 | list entry; |
---|
682 | for(i=1; i<=n2; i++) |
---|
683 | { |
---|
684 | entry = dec2[i]; |
---|
685 | if(n1==0) |
---|
686 | { |
---|
687 | dec1=list(entry); n1++; |
---|
688 | i++; continue; |
---|
689 | } |
---|
690 | |
---|
691 | a=1; |
---|
692 | b=n1; |
---|
693 | m = (a+b) div 2; |
---|
694 | while(b>a) |
---|
695 | { |
---|
696 | if(is_dp_smaller(dec1[m][2], dec1[m][3], entry[2], entry[3])) |
---|
697 | { |
---|
698 | a=m+1; |
---|
699 | } |
---|
700 | else |
---|
701 | { |
---|
702 | b=m; |
---|
703 | } |
---|
704 | m = (a+b) div 2; |
---|
705 | } |
---|
706 | if(is_dp_smaller(dec1[a][2], dec1[a][3], entry[2], entry[3])) |
---|
707 | {dec1=insert(dec1,entry,a); n1++;} |
---|
708 | else{if(entry[2]==dec1[a][2] && entry[3]==dec1[a][3]) |
---|
709 | { |
---|
710 | dec1[a][1] = dec1[a][1] + entry[1]; //same denominator: add numerators |
---|
711 | if(dec1[a][1]==0) |
---|
712 | {dec1 = delete(dec1,a); n1--;} |
---|
713 | } |
---|
714 | else |
---|
715 | {dec1=insert(dec1,entry,a-1); n1++;}} |
---|
716 | } |
---|
717 | return(dec1); |
---|
718 | } |
---|
719 | |
---|
720 | static proc is_dp_smaller(intvec indices1, intvec e1, intvec indices2, intvec e2) |
---|
721 | { |
---|
722 | if(size(e2)==0) {return(0);} |
---|
723 | if(size(e1)==0) {return(1);} |
---|
724 | int s1,s2 = sum(e1),sum(e2); |
---|
725 | if(s1<s2) {return(1);} |
---|
726 | if(s1>s2) {return(0);} |
---|
727 | int n1,n2 = size(indices1),size(indices2); |
---|
728 | int imax = min(n1,n2); |
---|
729 | for(int i=0; i<imax; i++) |
---|
730 | { |
---|
731 | if(indices1[n1-i]>indices2[n2-i]) {return(1);} |
---|
732 | if(indices1[n1-i]<indices2[n2-i]) {return(0);} |
---|
733 | if(e1[n1-i]>e2[n2-i]) {return(1);} |
---|
734 | if(e1[n1-i]<e2[n2-i]) {return(0);} |
---|
735 | } |
---|
736 | return(0); |
---|
737 | } |
---|
738 | |
---|
739 | proc checkpfd(list fraction, list dec, list #) |
---|
740 | "USAGE: checkpfd(list(f,g),dec[,N,C]); f,g poly, dec list, N,C int |
---|
741 | RETURN: 0 or 1 |
---|
742 | PURPOSE: test for (mathematical) equality of f/g and a partial fraction |
---|
743 | decomposition dec. The list dec has to have the same structure as the |
---|
744 | output of @ref{pfd}. |
---|
745 | @* The denominator g can also be given in factorized form as a list of |
---|
746 | an ideal of irreducible non constant polynomials and an intvec of |
---|
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.) |
---|
750 | @* By default the test is done (exactly) by bringing all terms of the |
---|
751 | decomposition on the same denominator and comparing to f/g. |
---|
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. |
---|
755 | SEE ALSO: pfd |
---|
756 | EXAMPLE: example checkpfd; shows an example" |
---|
757 | { |
---|
758 | poly f = fraction[1]; |
---|
759 | def g = fraction[2]; |
---|
760 | if(size(#)>0) |
---|
761 | { |
---|
762 | if(#[1]>0) |
---|
763 | { |
---|
764 | int N = #[1]; // number of random tests |
---|
765 | int max_val=16; |
---|
766 | if(size(#)>1) {max_val = #[2];} |
---|
767 | ideal values; |
---|
768 | ideal vars = maxideal(1); |
---|
769 | int d=nvars(basering); |
---|
770 | number val1,val2; |
---|
771 | int div_by_0; |
---|
772 | int i,j; |
---|
773 | for(i=1; i<=N; i++) |
---|
774 | { |
---|
775 | values = ideal(random(max_val,1,d)); |
---|
776 | |
---|
777 | if(typeof(g)=="poly") |
---|
778 | {val1 = number(substitute(g,vars,values));} |
---|
779 | else{if(typeof(g)=="list") // denominator given in factorized form |
---|
780 | { |
---|
781 | val1 = number(1); |
---|
782 | for(j=1; j<=size(g[1]); j++) |
---|
783 | {val1 = val1 * number(substitute(g[1][j]^g[2][j],vars,values));} |
---|
784 | } |
---|
785 | else {ERROR("wrong type for second argument, expected poly or list");}} |
---|
786 | |
---|
787 | if(val1==0) {continue;} |
---|
788 | val1 = number(substitute(f,vars,values))/val1; |
---|
789 | val2, div_by_0 = evaluatepfd(dec,values,2); |
---|
790 | if(div_by_0) {continue;} |
---|
791 | if(val1 != val2) |
---|
792 | {return(0);} |
---|
793 | } |
---|
794 | return(1); |
---|
795 | } |
---|
796 | } |
---|
797 | ideal q = dec[1]; |
---|
798 | dec = dec[2]; |
---|
799 | |
---|
800 | int m = size(q); |
---|
801 | int jmax,j,ind,k; |
---|
802 | intvec e_max; e_max[m]=0; |
---|
803 | int imax = size(dec); |
---|
804 | for(int i=1; i<=imax; i++) |
---|
805 | { |
---|
806 | jmax = size(dec[i][2]); |
---|
807 | for(j=1; j<=jmax; j++) |
---|
808 | { |
---|
809 | ind = dec[i][2][j]; |
---|
810 | e_max[ind] = max(e_max[ind],dec[i][3][j]); |
---|
811 | } |
---|
812 | } |
---|
813 | poly num; |
---|
814 | poly sum_of_numerators = 0; |
---|
815 | intvec e; |
---|
816 | for(i=1; i<=imax; i++) |
---|
817 | { |
---|
818 | e = e_max; |
---|
819 | jmax = size(dec[i][2]); |
---|
820 | for(j=1; j<=jmax; j++) |
---|
821 | { |
---|
822 | ind = dec[i][2][j]; |
---|
823 | e[ind] = e[ind]-dec[i][3][j]; |
---|
824 | } |
---|
825 | num = dec[i][1]; |
---|
826 | for(j=1; j<=m; j++) |
---|
827 | {num = num * q[j]^(e[j]);} |
---|
828 | sum_of_numerators = sum_of_numerators + num; |
---|
829 | } |
---|
830 | // the decomposition is now equal to sum_of_numerators/product(q[i]^e_max[i]) (i from 1 to imax) |
---|
831 | // now: check if this equals f/g: |
---|
832 | if(typeof(g)=="poly") |
---|
833 | { |
---|
834 | list fact = factorize(g); |
---|
835 | ideal q_g = delete(fact[1],1); |
---|
836 | intvec e_g = delete(fact[2],1); |
---|
837 | num = f/fact[1][1]; |
---|
838 | } |
---|
839 | else{if(typeof(g)=="list") //denominator given in factorized form |
---|
840 | { |
---|
841 | ideal q_g = g[1]; |
---|
842 | intvec e_g = g[2]; |
---|
843 | num = f; |
---|
844 | } |
---|
845 | else {ERROR("wrong type for second argument, expected poly or list");}} |
---|
846 | |
---|
847 | int m_g = size(q_g); |
---|
848 | int expon; |
---|
849 | number c; |
---|
850 | for(i=1; i<=m_g; i++) |
---|
851 | { |
---|
852 | j=0; |
---|
853 | for(k=1; k<=m; k++) |
---|
854 | { |
---|
855 | c = leadcoef(q[k])/leadcoef(q_g[i]); |
---|
856 | if(c*q_g[i]==q[k]) {j=k; break;} |
---|
857 | } |
---|
858 | if(j==0) |
---|
859 | {sum_of_numerators = sum_of_numerators*q_g[i]^e_g[i];} |
---|
860 | else |
---|
861 | { |
---|
862 | num = num*(c^e_g[i]); //fix lead coefficient |
---|
863 | expon = e_g[i]-e_max[j]; |
---|
864 | if(expon>0) |
---|
865 | {sum_of_numerators = sum_of_numerators*q[j]^expon;} |
---|
866 | else{if(expon<0) |
---|
867 | {num = num*q[j]^(-expon);}} |
---|
868 | } |
---|
869 | } |
---|
870 | return(sum_of_numerators==num); |
---|
871 | } |
---|
872 | example |
---|
873 | { |
---|
874 | "EXAMPLE:"; |
---|
875 | echo=voice; |
---|
876 | ring R = 0,(x,y),dp; |
---|
877 | poly f = x^3+3*x^2*y+2*y^2-x^2+4*x*y; |
---|
878 | poly g = x^2*y*(x-1)*(x-y)^2; |
---|
879 | |
---|
880 | // partial fraction decomposition of f/g: |
---|
881 | list dec = pfd(f,g); |
---|
882 | // some other decomposition (not equal to f/g): |
---|
883 | list wrong_dec = pfd(f+1,g); |
---|
884 | |
---|
885 | displaypfd_long(dec); |
---|
886 | list fraction = f,g; |
---|
887 | |
---|
888 | // exact test: |
---|
889 | checkpfd(fraction,dec); |
---|
890 | checkpfd(fraction,wrong_dec); |
---|
891 | // probabilistic test (evaluation at 10 random points): |
---|
892 | checkpfd(fraction,dec,10); |
---|
893 | checkpfd(fraction,wrong_dec,10); |
---|
894 | } |
---|
895 | |
---|
896 | proc evaluatepfd(list dec, ideal values, list #) |
---|
897 | "USAGE: evaluatepfd(dec, values[, mode]); dec list, values ideal, mode int |
---|
898 | RETURN: the number gotten by substituting the numbers generating the ideal |
---|
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 |
---|
904 | becomes 0, and 0 otherwise. |
---|
905 | SEE ALSO: pfd |
---|
906 | EXAMPLE: example evaluatepfd; shows an example" |
---|
907 | { |
---|
908 | int mode = 1; |
---|
909 | if(size(#)>0) {mode = #[1];} |
---|
910 | |
---|
911 | ideal vars = maxideal(1); // ideal generated by ring variables |
---|
912 | number val=0; |
---|
913 | number denom; |
---|
914 | int m,i,j; |
---|
915 | for(i=1; i<=size(dec[2]); i++) |
---|
916 | { |
---|
917 | denom = 1; |
---|
918 | m = size(dec[2][i][2]); |
---|
919 | for(j=1; j<=m; j++) |
---|
920 | { |
---|
921 | denom = denom * (number(substitute(dec[1][dec[2][i][2][j]],vars,values)))^dec[2][i][3][j]; |
---|
922 | if(denom == 0) |
---|
923 | { |
---|
924 | if(mode==1) {ERROR("division by 0");} |
---|
925 | else {return(0,1);} |
---|
926 | } |
---|
927 | } |
---|
928 | val = val + number(substitute(dec[2][i][1],vars,values))/denom; |
---|
929 | } |
---|
930 | if(mode==1) {return(val);} |
---|
931 | else {return(val,0);} |
---|
932 | } |
---|
933 | example |
---|
934 | { |
---|
935 | "EXAMPLE:"; |
---|
936 | echo=voice; |
---|
937 | ring R = 0,(x,y),dp; |
---|
938 | poly f = x+2*y; |
---|
939 | poly g = x^2-y^2; |
---|
940 | |
---|
941 | // partial fraction decomposition of f/g: |
---|
942 | list dec = pfd(f,g); |
---|
943 | |
---|
944 | displaypfd_long(dec); |
---|
945 | // evaluation at x=2, y=1: |
---|
946 | ideal values = 2,1; |
---|
947 | evaluatepfd(dec,values); |
---|
948 | |
---|
949 | // compare: f(2,1)/g(2,1) = (2+2*1)/(2^2-1^1) = 4/3 |
---|
950 | } |
---|
951 | |
---|
952 | static proc find_entry(def l, def entry) |
---|
953 | { |
---|
954 | int n=size(l); |
---|
955 | for(int i=1; i<=n; i++) |
---|
956 | {if(entry==l[i]) {return(i);}} |
---|
957 | return(0); |
---|
958 | } |
---|
959 | |
---|
960 | proc displaypfd(list dec) |
---|
961 | "USAGE: displaypfd(dec); dec list |
---|
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}. |
---|
965 | SEE ALSO: pfd, displaypfd_long, getStringpfd, getStringpfd_indexed |
---|
966 | EXAMPLE: example displaypfd; shows an example" |
---|
967 | { |
---|
968 | ideal q = dec[1]; |
---|
969 | dec = dec[2]; |
---|
970 | int imax=size(dec); |
---|
971 | int jmax,j; |
---|
972 | string s1,s2; |
---|
973 | for(int i=1; i<=imax; i++) |
---|
974 | { |
---|
975 | s1 = "(" + string(dec[i][1]); |
---|
976 | if(i>1) {s1 = "+ " + s1;} |
---|
977 | else {s1 = " " + s1;} |
---|
978 | if(size(dec[i][2])>0) |
---|
979 | { |
---|
980 | s2 = ") / ("; |
---|
981 | jmax = size(dec[i][2]); |
---|
982 | for(j=1; j<=jmax; j++) |
---|
983 | { |
---|
984 | s2 = s2 + "q" + string(dec[i][2][j]); |
---|
985 | if(dec[i][3][j] != 1) {s2 = s2 + "^" + string(dec[i][3][j]);} |
---|
986 | if(j<jmax) {s2 = s2 + "*";} |
---|
987 | } |
---|
988 | s2 = s2 + ")"; |
---|
989 | } |
---|
990 | else {s2 = ")";} |
---|
991 | |
---|
992 | if(size(s1)+size(s2)>192) {s1=s1[1..(192-size(s2))]; s1 = s1 + "... ";} |
---|
993 | print(s1+s2); |
---|
994 | } |
---|
995 | print("where"); |
---|
996 | for(i=1; i<=size(q); i++) |
---|
997 | { |
---|
998 | printf("q%s = %s",i,q[i]); |
---|
999 | } |
---|
1000 | if(size(dec)==1) {printf("(%s terms)%n", 1, 0);} |
---|
1001 | else {printf("(%s terms)%n", size(dec), 0);} |
---|
1002 | } |
---|
1003 | example |
---|
1004 | { |
---|
1005 | "EXAMPLE:"; |
---|
1006 | echo=voice; |
---|
1007 | ring R = 0,(x,y),dp; |
---|
1008 | poly f = x^3+3*x^2*y+2*y^2-x^2+4*x*y; |
---|
1009 | poly g = x^2*y*(x-1)*(x-y)^2; |
---|
1010 | |
---|
1011 | list dec = pfd(f,g); |
---|
1012 | |
---|
1013 | displaypfd(dec); |
---|
1014 | } |
---|
1015 | |
---|
1016 | proc displaypfd_long(list dec) |
---|
1017 | "USAGE: displaypfd_long(dec); dec list |
---|
1018 | PURPOSE: like @ref{displaypfd}, but denominators are written out, not indexed. |
---|
1019 | SEE ALSO: pfd, displaypfd, getStringpfd, getStringpfd_indexed |
---|
1020 | EXAMPLE: example displaypfd_long; shows an example" |
---|
1021 | { |
---|
1022 | ideal q = dec[1]; |
---|
1023 | dec = dec[2]; |
---|
1024 | print(" "+getStringFraction(dec[1],q)); |
---|
1025 | for(int i=2; i<=size(dec); i++) |
---|
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);} |
---|
1029 | } |
---|
1030 | example |
---|
1031 | { |
---|
1032 | "EXAMPLE:"; |
---|
1033 | echo=voice; |
---|
1034 | ring R = 0,(x,y),dp; |
---|
1035 | poly f = x^3+3*x^2*y+2*y^2-x^2+4*x*y; |
---|
1036 | poly g = x^2*y*(x-1)*(x-y)^2; |
---|
1037 | |
---|
1038 | list dec = pfd(f,g); |
---|
1039 | |
---|
1040 | displaypfd_long(dec); |
---|
1041 | } |
---|
1042 | |
---|
1043 | proc getStringpfd(list dec) |
---|
1044 | "USAGE: getStringpfd(dec); dec list |
---|
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}. |
---|
1048 | SEE ALSO: pfd, getStringpfd_indexed, displaypfd, displaypfd_long |
---|
1049 | EXAMPLE: example getStringpfd; shows an example" |
---|
1050 | { |
---|
1051 | ideal q = dec[1]; |
---|
1052 | dec = dec[2]; |
---|
1053 | string s = getStringFraction(dec[1],q); |
---|
1054 | for(int i=2; i<=size(dec); i++) |
---|
1055 | {s = s+" + "+getStringFraction(dec[i],q);} |
---|
1056 | return(s); |
---|
1057 | } |
---|
1058 | example |
---|
1059 | { |
---|
1060 | "EXAMPLE:"; |
---|
1061 | echo=voice; |
---|
1062 | ring R = 0,(x,y),dp; |
---|
1063 | poly f = x^3+3*x^2*y+2*y^2-x^2+4*x*y; |
---|
1064 | poly g = x^2*y*(x-1)*(x-y)^2; |
---|
1065 | |
---|
1066 | list dec = pfd(f,g); |
---|
1067 | |
---|
1068 | displaypfd_long(dec); |
---|
1069 | getStringpfd(dec); |
---|
1070 | } |
---|
1071 | |
---|
1072 | proc getStringpfd_indexed(list dec) |
---|
1073 | "USAGE: getStringpfd_indexed(dec); dec list |
---|
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}. |
---|
1078 | SEE ALSO: pfd, getStringpfd, displaypfd, displaypfd_long |
---|
1079 | EXAMPLE: example getStringpfd_indexed; shows an example" |
---|
1080 | { |
---|
1081 | if(typeof(dec[1])=="ideal") {dec = dec[2];} |
---|
1082 | string s = getStringFraction_indexed(dec[1]); |
---|
1083 | for(int i=2; i<=size(dec); i++) |
---|
1084 | {s = s+" + "+getStringFraction_indexed(dec[i]);} |
---|
1085 | return(s); |
---|
1086 | } |
---|
1087 | example |
---|
1088 | { |
---|
1089 | "EXAMPLE:"; |
---|
1090 | echo=voice; |
---|
1091 | ring R = 0,(x,y),dp; |
---|
1092 | poly f = x^3+3*x^2*y+2*y^2-x^2+4*x*y; |
---|
1093 | poly g = x^2*y*(x-1)*(x-y)^2; |
---|
1094 | |
---|
1095 | list dec = pfd(f,g); |
---|
1096 | |
---|
1097 | displaypfd(dec); |
---|
1098 | getStringpfd_indexed(dec); |
---|
1099 | } |
---|
1100 | |
---|
1101 | static proc getStringFraction(list fraction, ideal q) |
---|
1102 | { |
---|
1103 | int jmax,j; |
---|
1104 | string s = "(" + string(fraction[1]); |
---|
1105 | if(size(fraction[2])>0) |
---|
1106 | { |
---|
1107 | s = s + ")/("; |
---|
1108 | jmax = size(fraction[2]); |
---|
1109 | for(j=1; j<=jmax; j++) |
---|
1110 | { |
---|
1111 | s = s + "(" + string(q[fraction[2][j]]) + ")"; |
---|
1112 | if(fraction[3][j] != 1) {s = s + "^" + string(fraction[3][j]);} |
---|
1113 | if(j<jmax) {s = s + "*";} |
---|
1114 | } |
---|
1115 | } |
---|
1116 | s = s + ")"; |
---|
1117 | return(s); |
---|
1118 | } |
---|
1119 | |
---|
1120 | static proc getStringFraction_indexed(list fraction) |
---|
1121 | { |
---|
1122 | int jmax,j; |
---|
1123 | string s = "(" + string(fraction[1]); |
---|
1124 | if(size(fraction[2])>0) |
---|
1125 | { |
---|
1126 | s = s + ")/("; |
---|
1127 | jmax = size(fraction[2]); |
---|
1128 | for(j=1; j<=jmax; j++) |
---|
1129 | { |
---|
1130 | s = s + "q" + string(fraction[2][j]); |
---|
1131 | if(fraction[3][j] != 1) {s = s + "^" + string(fraction[3][j]);} |
---|
1132 | if(j<jmax) {s = s + "*";} |
---|
1133 | } |
---|
1134 | } |
---|
1135 | s = s + ")"; |
---|
1136 | return(s); |
---|
1137 | } |
---|
1138 | |
---|
1139 | proc readInputTXT(def file__, list #) |
---|
1140 | "USAGE: readInputTXT(file[, mode]), file string, mode int |
---|
1141 | readInputTXT(filelist[, mode]), filelist list, mode int |
---|
1142 | PURPOSE: read matrix of rational functions from a txt-file and turn it into a |
---|
1143 | matrix (i.e. a list of lists) of pairs of polynomials (numerators and |
---|
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.) |
---|
1158 | @* The file should contain less than 2^31 characters (filesize < 2 GB). |
---|
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. |
---|
1163 | @* Also the basering has to match the variable names used in the |
---|
1164 | input file(s). |
---|
1165 | @* @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 |
---|
1168 | SEE ALSO: pfdMat" |
---|
1169 | { |
---|
1170 | system("--ticks-per-sec",1000); |
---|
1171 | short = 0; |
---|
1172 | if(!defined(basering)) |
---|
1173 | {ERROR("no basering defined!");} |
---|
1174 | int left__,right__,pos1__,mid__,pos2__,tmp__,i__,j__,k__,t__,tt__,depth__; |
---|
1175 | int mode__=1; |
---|
1176 | if(size(#)>0) {mode__=#[1];} |
---|
1177 | if(typeof(file__)=="list") // list of filenames given --> apply function to each |
---|
1178 | { // file and concatenate the resulting matrices |
---|
1179 | int n = size(file__); |
---|
1180 | list mat__ = list(); |
---|
1181 | for(i__=1;i__<=n;i__++) |
---|
1182 | { |
---|
1183 | dbprint(sprintf(" file %s of %s:",i__,n)); |
---|
1184 | if(find(file__[i__],".txt")==0) {ERROR("wrong file type, expected txt");} |
---|
1185 | printlevel = printlevel+1; |
---|
1186 | mat__ = mat__ + readInputTXT(file__[i__],2); |
---|
1187 | printlevel = printlevel-1; |
---|
1188 | } |
---|
1189 | |
---|
1190 | if(mode__==2) {return(mat__);} |
---|
1191 | |
---|
1192 | string filename__ = file__[1][1,find(file__[1],".txt")-1]; |
---|
1193 | dbprint(" saving to file "+filename__+".ssi "); t__ = rtimer; |
---|
1194 | write("ssi:w "+filename__+".ssi", mat__); |
---|
1195 | dbprint(sprintf(" done! (%s ms)", rtimer-t__)); |
---|
1196 | |
---|
1197 | if(mode__==3) {return(mat__);} |
---|
1198 | } |
---|
1199 | if(typeof(file__)!="string") |
---|
1200 | {ERROR("wrong type for first argument (expected string or list)");} |
---|
1201 | if(find(file__,".txt")==0) {ERROR("wrong file type, expected txt");} |
---|
1202 | |
---|
1203 | dbprint(" reading file "); t__=rtimer; |
---|
1204 | string data__ = read(":r "+file__); |
---|
1205 | dbprint(sprintf(" done! (%s ms)", rtimer-t__)); |
---|
1206 | |
---|
1207 | dbprint(" processing input "); t__ = rtimer; |
---|
1208 | left__ = find(data__,"{"); |
---|
1209 | right__ = find(data__,"}"); |
---|
1210 | tmp__ = find(data__,"{",left__+1); |
---|
1211 | if(left__<tmp__ && tmp__<right__) {left__ = tmp__;} |
---|
1212 | |
---|
1213 | int finished__,n__; |
---|
1214 | poly p__,q__; |
---|
1215 | list row__,mat__; |
---|
1216 | string s__,ss__; |
---|
1217 | i__=0; |
---|
1218 | while(1) |
---|
1219 | { |
---|
1220 | i__++; |
---|
1221 | tt__ = rtimer; |
---|
1222 | row__ = list(); |
---|
1223 | pos2__ = left__; |
---|
1224 | finished__ = 0; |
---|
1225 | j__=0; |
---|
1226 | while(not finished__) |
---|
1227 | { |
---|
1228 | j__++; |
---|
1229 | s__ = ""; |
---|
1230 | |
---|
1231 | pos1__ = pos2__+1; |
---|
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 |
---|
1244 | { |
---|
1245 | execute("p__=" + s__); |
---|
1246 | q__=1; |
---|
1247 | row__[j__] = list(p__,q__); |
---|
1248 | continue; |
---|
1249 | } |
---|
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__])); |
---|
1293 | |
---|
1294 | row__[j__] = list(p__,q__); |
---|
1295 | } |
---|
1296 | mat__[i__] = row__; // append row to matrix |
---|
1297 | dbprint(sprintf(" row %s done! (%s ms)",i__,rtimer-tt__)); |
---|
1298 | |
---|
1299 | left__ = find(data__,"{",right__); |
---|
1300 | if(left__==0) {break;} |
---|
1301 | right__ = find(data__,"}",left__); |
---|
1302 | } |
---|
1303 | dbprint(sprintf(" done! (%s ms)", rtimer-t__)); |
---|
1304 | |
---|
1305 | if(mode__==2) {return(mat__);} |
---|
1306 | |
---|
1307 | string filename__ = file__[1,find(file__,".txt")-1]; |
---|
1308 | dbprint(" saving to file "+filename__+".ssi "); t__ = rtimer; |
---|
1309 | write("ssi:w "+filename__+".ssi", mat__); |
---|
1310 | dbprint(sprintf(" done! (%s ms)", rtimer-t__)); |
---|
1311 | |
---|
1312 | if(mode__==3) {return(mat__);} |
---|
1313 | } |
---|
1314 | |
---|
1315 | static proc fixBrackets(string data) |
---|
1316 | { |
---|
1317 | int pos=0; |
---|
1318 | int left_brackets =0; |
---|
1319 | int right_brackets=0; |
---|
1320 | int n=size(data); |
---|
1321 | while(pos<n) |
---|
1322 | { |
---|
1323 | pos = find(data,"(",pos+1); |
---|
1324 | if(pos==0) {break;} |
---|
1325 | left_brackets++; |
---|
1326 | } |
---|
1327 | pos=0; |
---|
1328 | while(pos<n) |
---|
1329 | { |
---|
1330 | pos = find(data,")",pos+1); |
---|
1331 | if(pos==0) {break;} |
---|
1332 | right_brackets++; |
---|
1333 | } |
---|
1334 | int difference = left_brackets-right_brackets; |
---|
1335 | if(difference>0) |
---|
1336 | { |
---|
1337 | for(int i=1; i<=difference; i++) {data = data+")";} |
---|
1338 | } |
---|
1339 | if(difference<0) |
---|
1340 | { |
---|
1341 | for(int i=1; i<=-difference; i++) {data = "("+data;} |
---|
1342 | } |
---|
1343 | return(data); |
---|
1344 | } |
---|
1345 | |
---|
1346 | static proc pfdWrap(poly f, def g, int i, int j, link logfile, int output_mode) |
---|
1347 | { |
---|
1348 | system("--ticks-per-sec",1000); |
---|
1349 | if(output_mode>3) |
---|
1350 | {write("ssi:w pfd_results_"+string(i)+"_"+string(j) |
---|
1351 | +".ssi","task started, but not finished yet");} |
---|
1352 | int t0 = rtimer; |
---|
1353 | list result = pfd(f,g); |
---|
1354 | fprintf(logfile,"_[%s,%s]: %s ms",i,j,rtimer-t0); |
---|
1355 | if(output_mode>3) |
---|
1356 | {write("ssi:w pfd_results_"+string(i)+"_"+string(j)+".ssi",result);} |
---|
1357 | return(result); |
---|
1358 | } |
---|
1359 | |
---|
1360 | static proc testEntry(int i, int j, list fraction, list dec, link logfile, int N); |
---|
1361 | { |
---|
1362 | int t=rtimer; |
---|
1363 | system("--ticks-per-sec",1000); |
---|
1364 | int result = checkpfd(fraction,dec,N); |
---|
1365 | if(result==1) {fprintf(logfile, " _[%s,%s]: correct (%s ms)",i,j,rtimer-t);} |
---|
1366 | else {fprintf(logfile, " _[%s,%s]: WRONG! (%s ms)", i,j,rtimer-t);} |
---|
1367 | return(result); |
---|
1368 | } |
---|
1369 | |
---|
1370 | proc pfdMat(def infile, list #) |
---|
1371 | "USAGE: pfdMat(file[, dotest, ignore_nonlin, output_mode, parallelize]); |
---|
1372 | file string, dotest,ignore_nonlin,output_mode,parallelize int |
---|
1373 | PURPOSE: apply @code{pfd} to all entries of a matrix of rational functions |
---|
1374 | saved in a txt-file. The string @code{file} should be the |
---|
1375 | [directory +] name of the file. |
---|
1376 | @* The input file can either be a txt-file or an ssi-file created with |
---|
1377 | @code{readInputTXT}. In case of a txt-file, the base ring has to match |
---|
1378 | and the matrix has to be in the same format specified in |
---|
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. |
---|
1384 | @* Also a logfile is created, which protocols the memory used and the |
---|
1385 | runtimes of @code{pfd} for each matrix entry in real-time. |
---|
1386 | |
---|
1387 | @* There are also 4 optional arguments: |
---|
1388 | @* If @code{dotest} is nonzero, test the results with checkpfd: |
---|
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. |
---|
1449 | SEE ALSO: readInputTXT, pfd, checkpfd, checkpfdMat" |
---|
1450 | { |
---|
1451 | system("--ticks-per-sec",1000); |
---|
1452 | short = 0; |
---|
1453 | int dotest,ignore_nonlin,output_mode,parallelize = -1,1,1,1; |
---|
1454 | if(size(#)>0) {dotest = #[1];} |
---|
1455 | if(size(#)>1) {ignore_nonlin = #[2];} |
---|
1456 | if(size(#)>2) {output_mode = #[3];} |
---|
1457 | if(size(#)>3) {parallelize = #[4];} |
---|
1458 | |
---|
1459 | int i,j,k,l,ind; |
---|
1460 | list arguments,results; |
---|
1461 | |
---|
1462 | dbprint(newline+"reading data "); int t0=rtimer; |
---|
1463 | |
---|
1464 | if(typeof(infile)=="list") |
---|
1465 | { |
---|
1466 | printlevel = printlevel+1; |
---|
1467 | if(output_mode>2) {list mat = readInputTXT(infile,3);} |
---|
1468 | else {list mat = readInputTXT(infile,2);} |
---|
1469 | printlevel = printlevel-1; |
---|
1470 | int pos=find(infile[1],".txt"); |
---|
1471 | string filename = infile[1][1,pos-1]; |
---|
1472 | } |
---|
1473 | else |
---|
1474 | { |
---|
1475 | int pos=find(infile,".txt"); |
---|
1476 | if(pos!=0) |
---|
1477 | { |
---|
1478 | printlevel = printlevel+1; |
---|
1479 | if(output_mode>2) {list mat = readInputTXT(infile,3);} |
---|
1480 | else {list mat = readInputTXT(infile,2);} |
---|
1481 | printlevel = printlevel-1; |
---|
1482 | } |
---|
1483 | else |
---|
1484 | { |
---|
1485 | pos=find(infile,".ssi"); |
---|
1486 | if(pos!=0) {list mat = read("ssi:r "+infile);} |
---|
1487 | else {ERROR("invalid file type, expected ssi or txt");} |
---|
1488 | } |
---|
1489 | string filename = infile[1,pos-1]; |
---|
1490 | } |
---|
1491 | |
---|
1492 | link qfile = ":w "+filename+"_denominator_factors.txt"; |
---|
1493 | int n = size(mat); |
---|
1494 | int m = size(mat[1]); |
---|
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; |
---|
1501 | mat = FactDenom(mat); |
---|
1502 | printlevel = printlevel-1; |
---|
1503 | if(output_mode>2) |
---|
1504 | {write("ssi:w "+filename+"_factorized_denominators.ssi",mat);} |
---|
1505 | dbprint(sprintf("done! (%s ms)",rtimer-t0)); |
---|
1506 | } |
---|
1507 | |
---|
1508 | if(ignore_nonlin) |
---|
1509 | { |
---|
1510 | dbprint("removing nonlinear denominator factors before pfd is applied"); |
---|
1511 | list nonlin_denom_factors; |
---|
1512 | ideal p; |
---|
1513 | printlevel = printlevel+1; |
---|
1514 | mat, nonlin_denom_factors, p = removeNonlinearFactors(mat, filename); |
---|
1515 | printlevel = printlevel-1; |
---|
1516 | if(output_mode>2) |
---|
1517 | { |
---|
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; |
---|
1525 | write("ssi:w "+filename+"_linear_part.ssi", mat); |
---|
1526 | dbprint(sprintf("done! (%s ms)",rtimer-t0)); |
---|
1527 | } |
---|
1528 | } |
---|
1529 | |
---|
1530 | if(parallelize) |
---|
1531 | { |
---|
1532 | dbprint("creating tasks "); t0=rtimer; |
---|
1533 | write(":w "+filename+"_pfdMat_logfile.txt","finished matrix entries with runtimes " |
---|
1534 | +"(calculated in parallel on "+string(getcores())+" cores):"); |
---|
1535 | link logfile = ":a "+filename+"_pfdMat_logfile.txt"; |
---|
1536 | for(i=1; i<=n; i++) |
---|
1537 | { |
---|
1538 | for(j=1; j<=m; j++) |
---|
1539 | { |
---|
1540 | ind = m*(i-1)+j; |
---|
1541 | arguments[ind] = list(mat[i][j][1],mat[i][j][2],i,j,logfile,output_mode); |
---|
1542 | } |
---|
1543 | } |
---|
1544 | dbprint(sprintf("done! (%s ms)",rtimer-t0)); |
---|
1545 | |
---|
1546 | dbprint("applying pfd to each matrix entry "); t0 = rtimer; |
---|
1547 | results = parallelWaitAll("pfdWrap",arguments); |
---|
1548 | arguments = list(); |
---|
1549 | dbprint(sprintf("done! (%s ms)",rtimer-t0)); |
---|
1550 | |
---|
1551 | write(logfile,"decomposition: "+string(rtimer-t0)+" ms and "+string(memory(2)) |
---|
1552 | +" Byte Memory max. (after calling pfd on each matrix entry)"); |
---|
1553 | |
---|
1554 | dbprint("writing results in matrix shape "); t0 = rtimer; |
---|
1555 | list dec_mat; |
---|
1556 | for(i=1; i<=n; i++) |
---|
1557 | { |
---|
1558 | dec_mat[i] = list(); |
---|
1559 | for(j=1; j<=m; j++) |
---|
1560 | { |
---|
1561 | ind = m*(i-1)+j; |
---|
1562 | dec_mat[i][j] = results[ind]; |
---|
1563 | } |
---|
1564 | } |
---|
1565 | results = list(); |
---|
1566 | dbprint(sprintf("done! (%s ms)",rtimer-t0)); |
---|
1567 | } |
---|
1568 | else |
---|
1569 | { |
---|
1570 | dbprint("applying pfd to each matrix entry "); t0=rtimer; |
---|
1571 | write(":w "+filename+"_pfdMat_logfile.txt", |
---|
1572 | "finished matrix entries with runtimes (no parallelization):"); |
---|
1573 | link logfile = ":a "+filename+"_pfdMat_logfile.txt"; |
---|
1574 | list dec_mat; |
---|
1575 | for(i=1; i<=n; i++) |
---|
1576 | { |
---|
1577 | dec_mat[i] = list(); |
---|
1578 | for(j=1; j<=m; j++) |
---|
1579 | { |
---|
1580 | dec_mat[i][j] = pfdWrap(mat[i][j][1],mat[i][j][2],i,j,logfile,output_mode); |
---|
1581 | } |
---|
1582 | } |
---|
1583 | dbprint(sprintf("done! (%s ms)",rtimer-t0)); |
---|
1584 | write(logfile,"decomposition: "+string(rtimer-t0)+" ms and "+string(memory(2)) |
---|
1585 | +" Byte Memory max. (after calling pfd on each matrix entry)"); |
---|
1586 | } |
---|
1587 | |
---|
1588 | dbprint("making one single list of denominator factors "); t0 = rtimer; |
---|
1589 | ideal q,new_q; |
---|
1590 | intvec dict; |
---|
1591 | for(i=1; i<=n; i++) |
---|
1592 | { |
---|
1593 | for(j=1; j<=m; j++) |
---|
1594 | { |
---|
1595 | new_q = dec_mat[i][j][1]; |
---|
1596 | dec_mat[i][j] = dec_mat[i][j][2]; |
---|
1597 | dict = 0:0; |
---|
1598 | for(k=1; k<=size(new_q); k++) |
---|
1599 | { |
---|
1600 | ind = find_entry(q,new_q[k]); |
---|
1601 | if(ind==0) |
---|
1602 | { |
---|
1603 | ind = size(q)+1; |
---|
1604 | q[ind] = new_q[k]; |
---|
1605 | } |
---|
1606 | dict[k] = ind; |
---|
1607 | } |
---|
1608 | for(k=1; k<=size(dec_mat[i][j]); k++) |
---|
1609 | { |
---|
1610 | if(size(dec_mat[i][j][k][2])>0) |
---|
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)); |
---|
1617 | |
---|
1618 | if(output_mode>2) |
---|
1619 | { |
---|
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)); |
---|
1623 | } |
---|
1624 | |
---|
1625 | if(ignore_nonlin) |
---|
1626 | { |
---|
1627 | ind = size(q); |
---|
1628 | for(i=1; i<=size(p); i++) {q[ind+i]=p[i];} // add nonlin. polynomials to q |
---|
1629 | for(i=1; i<=n; i++) |
---|
1630 | { |
---|
1631 | for(j=1; j<=m; j++) |
---|
1632 | { |
---|
1633 | nonlin_denom_factors[i][j][1] = nonlin_denom_factors[i][j][1]+ind; // adjust indices |
---|
1634 | } |
---|
1635 | } |
---|
1636 | } |
---|
1637 | |
---|
1638 | if(ignore_nonlin) |
---|
1639 | {dbprint(sprintf("creating readable .txt-files (including the nonlinear factors again)"));} |
---|
1640 | else |
---|
1641 | {dbprint(sprintf("creating readable .txt-files "));} |
---|
1642 | t0 = rtimer; |
---|
1643 | dbprint(" indexed ("+filename+"_pfd_indexed.txt):"); |
---|
1644 | printlevel = printlevel+1; |
---|
1645 | if(ignore_nonlin) |
---|
1646 | {saveResultTXT_indexed(dec_mat, filename+"_pfd_indexed", nonlin_denom_factors);} |
---|
1647 | else {saveResultTXT_indexed(dec_mat, filename+"_pfd_indexed");} |
---|
1648 | printlevel = printlevel-1; |
---|
1649 | for(i=1; i<=size(q); i++) {fprintf(qfile, "q%s = %s;", i, q[i]);} |
---|
1650 | if(output_mode>1) |
---|
1651 | { |
---|
1652 | dbprint(" denominators written out ("+filename+"_pfd.txt):"); |
---|
1653 | printlevel = printlevel+1; |
---|
1654 | if(ignore_nonlin) |
---|
1655 | {saveResultTXT(dec_mat, q, filename+"_pfd", nonlin_denom_factors);} |
---|
1656 | else {saveResultTXT(dec_mat, q, filename+"_pfd");} |
---|
1657 | printlevel = printlevel-1; |
---|
1658 | } |
---|
1659 | dbprint(sprintf("done! (%s ms)",rtimer-t0)); |
---|
1660 | |
---|
1661 | if(dotest) |
---|
1662 | { |
---|
1663 | if(dotest<0) |
---|
1664 | {dbprint("checking for correctness (exact test) ");} |
---|
1665 | else |
---|
1666 | {dbprint(sprintf("checking for correctness (%s random evaluations per entry) ",dotest));} |
---|
1667 | t0 = rtimer; |
---|
1668 | if(parallelize) |
---|
1669 | { |
---|
1670 | for(i=1; i<=n; i++) |
---|
1671 | { |
---|
1672 | for(j=1; j<=m; j++) |
---|
1673 | { |
---|
1674 | arguments[m*(i-1)+j]=list(i,j,mat[i][j],list(q,dec_mat[i][j]),logfile,dotest); |
---|
1675 | } |
---|
1676 | } |
---|
1677 | results = parallelWaitAll("testEntry",arguments); |
---|
1678 | } |
---|
1679 | else |
---|
1680 | { |
---|
1681 | for(i=1; i<=n; i++) |
---|
1682 | { |
---|
1683 | for(j=1; j<=m; j++) |
---|
1684 | { |
---|
1685 | results[m*(i-1)+j]=testEntry(i,j,mat[i][j],list(q,dec_mat[i][j]),logfile,dotest); |
---|
1686 | } |
---|
1687 | } |
---|
1688 | } |
---|
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)); |
---|
1691 | write(logfile,"checking for correctness: "+string(rtimer-t0)+" ms and " |
---|
1692 | +string(memory(2))+" Byte Memory max. (at the end of pfdMat), " |
---|
1693 | +string(sum(results))+" correct out of "+string(n*m)); |
---|
1694 | } |
---|
1695 | } |
---|
1696 | |
---|
1697 | static proc FactDenom(list mat) |
---|
1698 | { |
---|
1699 | system("--ticks-per-sec",1000); |
---|
1700 | int i,j,k,ind,t,counter; |
---|
1701 | int n = size(mat); |
---|
1702 | int m = size(mat[1]); |
---|
1703 | list denom = list(); |
---|
1704 | for(i=1; i<=n; i++) |
---|
1705 | { |
---|
1706 | denom[i] = list(); |
---|
1707 | for(j=1; j<=m; j++) |
---|
1708 | { |
---|
1709 | denom[i][j] = mat[i][j][2]; |
---|
1710 | mat[i][j][2] = list(ideal(),intvec(0:0)); |
---|
1711 | } |
---|
1712 | } |
---|
1713 | int expon; |
---|
1714 | list fact; |
---|
1715 | number lcoeff; |
---|
1716 | int timeout = 60000; |
---|
1717 | int finished = 0; |
---|
1718 | list arguments,results; |
---|
1719 | ideal q; |
---|
1720 | while(!finished) |
---|
1721 | { |
---|
1722 | t = rtimer; |
---|
1723 | for(i=1; i<=n; i++) // create argument list |
---|
1724 | { |
---|
1725 | for(j=1; j<=m; j++) |
---|
1726 | {arguments[m*(i-1)+j] = list(denom[i][j]);} |
---|
1727 | } |
---|
1728 | |
---|
1729 | results = parallelWaitAll("factorize",arguments,timeout); |
---|
1730 | arguments = list(); |
---|
1731 | |
---|
1732 | for(i=1; i<=n; i++) // update q, mat, denom |
---|
1733 | { |
---|
1734 | for(j=1; j<=m; j++) |
---|
1735 | { |
---|
1736 | ind = m*(i-1)+j; |
---|
1737 | if(typeof(results[ind]) != "none") |
---|
1738 | { |
---|
1739 | fact = results[ind]; |
---|
1740 | for(k=2; k<=size(fact[1]); k++) |
---|
1741 | { |
---|
1742 | lcoeff = leadcoef(fact[1][k]); |
---|
1743 | fact[1][k] = fact[1][k]/lcoeff; |
---|
1744 | fact[1][1] = fact[1][1]*(lcoeff^fact[2][k]); // polynomial is monic (thus unique) |
---|
1745 | lcoeff = content(fact[1][k]); |
---|
1746 | fact[1][k] = fact[1][k]/lcoeff; |
---|
1747 | fact[1][1] = fact[1][1]*(lcoeff^fact[2][k]); // polynomial has nice coefficients |
---|
1748 | |
---|
1749 | // add a new factor to q: |
---|
1750 | ind = find_entry(q,fact[1][k]); |
---|
1751 | if(ind==0) {q[size(q)+1] = fact[1][k];} |
---|
1752 | // complete the factorization of the i,j-th denominator: |
---|
1753 | ind = find_entry(mat[i][j][2][1],fact[1][k]); |
---|
1754 | if(ind==0) |
---|
1755 | { |
---|
1756 | ind = size(mat[i][j][2][1])+1; |
---|
1757 | mat[i][j][2][1][ind] = fact[1][k]; |
---|
1758 | mat[i][j][2][2][ind] = fact[2][k]; |
---|
1759 | } |
---|
1760 | else |
---|
1761 | { |
---|
1762 | mat[i][j][2][2][ind] = mat[i][j][2][2][ind] + fact[2][k]; |
---|
1763 | } |
---|
1764 | } |
---|
1765 | mat[i][j][1] = mat[i][j][1]/fact[1][1]; |
---|
1766 | denom[i][j] = 1; |
---|
1767 | } |
---|
1768 | } |
---|
1769 | } |
---|
1770 | results = list(); |
---|
1771 | |
---|
1772 | finished = 1; |
---|
1773 | counter = n*m; |
---|
1774 | for(i=1; i<=n; i++) // factorize by any known factors (from q) |
---|
1775 | { |
---|
1776 | for(j=1; j<=m; j++) |
---|
1777 | { |
---|
1778 | for(k=1; k<=size(q); k++) |
---|
1779 | { |
---|
1780 | expon=0; |
---|
1781 | while(reduce(denom[i][j],q[k])==0) |
---|
1782 | { |
---|
1783 | denom[i][j] = denom[i][j]/q[k]; |
---|
1784 | expon++; |
---|
1785 | } |
---|
1786 | if(expon>0) |
---|
1787 | { |
---|
1788 | ind = size(mat[i][j][2][1])+1; |
---|
1789 | mat[i][j][2][1][ind] = q[k]; |
---|
1790 | mat[i][j][2][2][ind] = expon; |
---|
1791 | } |
---|
1792 | } |
---|
1793 | if(deg(denom[i][j])==0) |
---|
1794 | { |
---|
1795 | mat[i][j][1] = mat[i][j][1]/denom[i][j]; |
---|
1796 | denom[i][j] = poly(1); |
---|
1797 | } |
---|
1798 | if(denom[i][j]!=poly(1)) {finished = 0; counter--;} |
---|
1799 | } |
---|
1800 | } |
---|
1801 | timeout = timeout*2; |
---|
1802 | dbprint(sprintf(" %s out of %s denominators factorized completely (%s ms)", |
---|
1803 | counter,n*m,rtimer-t)); |
---|
1804 | } |
---|
1805 | return(mat); |
---|
1806 | } |
---|
1807 | |
---|
1808 | static proc saveResultTXT_indexed(list dec, string filename, list #) |
---|
1809 | { |
---|
1810 | // expect dec = list(list(dec11, dec12, ...), list(dec21, dec22, ...), ...), |
---|
1811 | // where dec11, dec12, ... are decompositions of form list(list(poly, intvec, intvec), ...) |
---|
1812 | system("--ticks-per-sec",1000); |
---|
1813 | list nonlinFactors = list(); |
---|
1814 | if(size(#)>0) {nonlinFactors = #;} |
---|
1815 | link file = ":w "+filename+".txt"; |
---|
1816 | |
---|
1817 | int i,j; |
---|
1818 | int t; |
---|
1819 | string s="{"; |
---|
1820 | int n=size(dec); |
---|
1821 | int m=size(dec[1]); |
---|
1822 | int k; |
---|
1823 | for(i=1; i<=n; i++) |
---|
1824 | { |
---|
1825 | t = rtimer; |
---|
1826 | s = s+"{"; |
---|
1827 | for(j=1; j<=m; j++) |
---|
1828 | { |
---|
1829 | if(size(nonlinFactors)>0) |
---|
1830 | { |
---|
1831 | if(size(nonlinFactors[i][j][1])>0) |
---|
1832 | { |
---|
1833 | s = s + getStringFraction_indexed(list(poly(1))+nonlinFactors[i][j]) |
---|
1834 | + " * (" + getStringpfd_indexed(dec[i][j]) + "), "; |
---|
1835 | j++; continue; |
---|
1836 | } |
---|
1837 | } |
---|
1838 | s = s + getStringpfd_indexed(dec[i][j]) + ", "; |
---|
1839 | } |
---|
1840 | k = size(s); |
---|
1841 | s[k-1] = "}"; s[k] = ","; s = s+" "; |
---|
1842 | dbprint(sprintf(" row %s done! (%s ms)",i,rtimer-t)); |
---|
1843 | } |
---|
1844 | k = size(s); |
---|
1845 | s[k-1] = "}"; s[k] = " "; |
---|
1846 | write(file,s); |
---|
1847 | } |
---|
1848 | |
---|
1849 | static proc saveResultTXT(list dec, ideal q, string filename, list #) |
---|
1850 | { |
---|
1851 | // expect dec = list(list(dec11, dec12, ...), list(dec21, dec22, ...), ...), |
---|
1852 | // where dec11, dec12, ... are decompositions of form list(list(poly, intvec, intvec), ...) |
---|
1853 | system("--ticks-per-sec",1000); |
---|
1854 | list nonlinFactors = list(); |
---|
1855 | if(size(#)>0) {nonlinFactors = #;} |
---|
1856 | link file = ":w "+filename+".txt"; |
---|
1857 | |
---|
1858 | int i,j; |
---|
1859 | int t; |
---|
1860 | string s="{"; |
---|
1861 | int n=size(dec); |
---|
1862 | int m=size(dec[1]); |
---|
1863 | int k; |
---|
1864 | for(i=1; i<=n; i++) |
---|
1865 | { |
---|
1866 | t = rtimer; |
---|
1867 | s = s+"{"; |
---|
1868 | for(j=1; j<=m; j++) |
---|
1869 | { |
---|
1870 | if(size(nonlinFactors)>0) |
---|
1871 | { |
---|
1872 | if(size(nonlinFactors[i][j][1])>0) |
---|
1873 | { |
---|
1874 | s = s + getStringFraction(list(poly(1))+nonlinFactors[i][j],q) |
---|
1875 | + " * (" + getStringpfd(list(q,dec[i][j])) + "), "; |
---|
1876 | j++; continue; |
---|
1877 | } |
---|
1878 | } |
---|
1879 | s = s + getStringpfd(list(q,dec[i][j])) + ", "; |
---|
1880 | } |
---|
1881 | k = size(s); |
---|
1882 | s[k-1] = "}"; s[k] = ","; s = s+" "; |
---|
1883 | dbprint(sprintf(" row %s done! (%s ms)",i,rtimer-t)); |
---|
1884 | } |
---|
1885 | k = size(s); |
---|
1886 | s[k-1] = "}"; s[k] = " "; |
---|
1887 | write(file,s); |
---|
1888 | } |
---|
1889 | |
---|
1890 | static proc removeNonlinearFactors(list fractions, string filename) |
---|
1891 | { |
---|
1892 | int n=size(fractions); |
---|
1893 | int m=size(fractions[1]); |
---|
1894 | int i,j,k,t0,ind; |
---|
1895 | ideal p; |
---|
1896 | list nonlin_denom_factors; |
---|
1897 | intvec factors,exponents; |
---|
1898 | list fac; |
---|
1899 | for(i=1; i<=n; i++) |
---|
1900 | { |
---|
1901 | t0 = rtimer; |
---|
1902 | nonlin_denom_factors[i] = list(); |
---|
1903 | for(j=1; j<=m; j++) |
---|
1904 | { |
---|
1905 | fac = fractions[i][j][2]; |
---|
1906 | factors = 0:0; |
---|
1907 | exponents = 0:0; |
---|
1908 | for(k=1; k<=size(fac[1]); k++) |
---|
1909 | { |
---|
1910 | if(deg(poly(fac[1][k]))>1) |
---|
1911 | { |
---|
1912 | // add the nonlin. factor fac[1][k] to p if necessary: |
---|
1913 | if(size(p)==0) {ind = 1; p[ind]=fac[1][k];} |
---|
1914 | else |
---|
1915 | { |
---|
1916 | for(ind=1;1;ind++) |
---|
1917 | { |
---|
1918 | if(p[ind]==fac[1][k]) {break;} |
---|
1919 | if(ind==size(p)) |
---|
1920 | { |
---|
1921 | ind++; |
---|
1922 | p[ind]=fac[1][k]; |
---|
1923 | break; |
---|
1924 | } |
---|
1925 | } |
---|
1926 | } |
---|
1927 | factors[size(factors)+1] = ind; |
---|
1928 | exponents[size(exponents)+1] = fac[2][k]; |
---|
1929 | fac[1] = delete(fac[1],k); |
---|
1930 | fac[2] = delete(fac[2],k); |
---|
1931 | continue; |
---|
1932 | } |
---|
1933 | } |
---|
1934 | fractions[i][j][2] = fac; |
---|
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); |
---|
1940 | } |
---|
1941 | |
---|
1942 | proc checkpfdMat(def input, string output, string qfile, list #) |
---|
1943 | "USAGE: checkpfdMat(input, output, denomFactors[, N, parallelize]); |
---|
1944 | input,output,denomFactors string, N,parallelize int |
---|
1945 | PURPOSE: test the output files of @code{pfdMat} for correctness. Input and |
---|
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}). |
---|
1952 | @* As for @code{readInputTXT} and @code{pfdMat}, the basering has to |
---|
1953 | match the variable names used in the input file, which has to be in |
---|
1954 | the same format specified in @ref{readInputTXT}. Also, files bigger |
---|
1955 | than 2 GB have to be split as described for @code{readInputTXT} and a |
---|
1956 | list of filenames can be given as first argument instead. |
---|
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}. |
---|
1963 | @* The result is printed and as in @code{pfdMat} a logfile is created |
---|
1964 | showing the results for each matrix entry. |
---|
1965 | SEE ALSO: readInputTXT, pfd, checkpfd, pfdMat" |
---|
1966 | { |
---|
1967 | system("--ticks-per-sec",1000); |
---|
1968 | short = 0; |
---|
1969 | int t0; |
---|
1970 | int N=0; |
---|
1971 | int parallelize=1; |
---|
1972 | if(size(#)==1) |
---|
1973 | { |
---|
1974 | if(typeof(#[1])=="int") {N=#[1];} |
---|
1975 | else {ERROR("invalid argument type: "+typeof(#[1]));} |
---|
1976 | } |
---|
1977 | if(size(#)==2) |
---|
1978 | { |
---|
1979 | if(typeof(#[1])=="int" && typeof(#[2])=="int") |
---|
1980 | {N=#[1]; parallelize=#[2];} |
---|
1981 | else {ERROR("invalid argument types: "+typeof(#[1])+", "+typeof(#[2]));} |
---|
1982 | } |
---|
1983 | if(size(#)>2) {ERROR("too many arguments");} |
---|
1984 | |
---|
1985 | dbprint(newline+"reading input file:"); |
---|
1986 | printlevel = printlevel+1; |
---|
1987 | list frac = readInputTXT(input,2); |
---|
1988 | printlevel = printlevel-1; |
---|
1989 | if(typeof(input)=="string") {string filename=input;} |
---|
1990 | if(typeof(input)=="list") {string filename=input[1];} |
---|
1991 | filename = filename[1,find(filename,".txt")-1]; |
---|
1992 | |
---|
1993 | dbprint("factorizing the denominators "); t0=rtimer; |
---|
1994 | printlevel = printlevel+1; |
---|
1995 | frac = FactDenom(frac); |
---|
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); |
---|
2001 | ideal q; |
---|
2002 | q = readQfileTXT(qfile); |
---|
2003 | dbprint(" done!"); |
---|
2004 | |
---|
2005 | dbprint(" reading (indexed) output decompositions "); |
---|
2006 | list dec,nonlin; |
---|
2007 | printlevel = printlevel+1; |
---|
2008 | dec,nonlin = readOutputTXT(output); |
---|
2009 | printlevel = printlevel-1; |
---|
2010 | |
---|
2011 | if(parallelize) |
---|
2012 | {dbprint(sprintf("done! (%s ms)%n%ncreating tasks",rtimer-t0,0)); t0=rtimer;} |
---|
2013 | else |
---|
2014 | { |
---|
2015 | dbprint(sprintf("done! (%s ms)%n",rtimer-t0,0)); |
---|
2016 | if(N<=0) |
---|
2017 | {dbprint("checking for correctness (exact test) ");} |
---|
2018 | else |
---|
2019 | {dbprint(sprintf("checking for correctness (%s random evaluations per entry) ",N));} |
---|
2020 | t0=rtimer; |
---|
2021 | } |
---|
2022 | |
---|
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"; |
---|
2027 | |
---|
2028 | int n=size(frac); |
---|
2029 | int m=size(frac[1]); |
---|
2030 | int i,j,k,ind; |
---|
2031 | list arguments; |
---|
2032 | for(i=1;i<=n;i++) |
---|
2033 | { |
---|
2034 | for(j=1;j<=m;j++) |
---|
2035 | { |
---|
2036 | for(k=1;k<=size(nonlin[i][j][1]);k++) |
---|
2037 | { |
---|
2038 | ind = find_entry(frac[i][j][2][1],q[nonlin[i][j][1][k]]); |
---|
2039 | if(ind==0) {ERROR("nonlinear factors are wrong");} |
---|
2040 | if(frac[i][j][2][2][ind]!=nonlin[i][j][2][k]) |
---|
2041 | {ERROR("nonlinear factors are wrong");} |
---|
2042 | frac[i][j][2][1] = delete(frac[i][j][2][1],ind); |
---|
2043 | frac[i][j][2][2] = delete(frac[i][j][2][2],ind); |
---|
2044 | } |
---|
2045 | if(parallelize) |
---|
2046 | {arguments[(i-1)*m+j] = list(i,j,frac[i][j],list(q,dec[i][j]),logfile,N);} |
---|
2047 | else |
---|
2048 | {results[(i-1)*m+j] = testEntry(i,j,frac[i][j],list(q,dec[i][j]),logfile,N);} |
---|
2049 | } |
---|
2050 | } |
---|
2051 | |
---|
2052 | if(parallelize) |
---|
2053 | { |
---|
2054 | dbprint(sprintf("done! (%s ms)%n",rtimer-t0,0)); |
---|
2055 | if(N<=0) |
---|
2056 | {dbprint("checking for correctness (exact test) ");} |
---|
2057 | else |
---|
2058 | {dbprint(sprintf("checking for correctness (%s random evaluations per entry) ",N));} |
---|
2059 | t0=rtimer; |
---|
2060 | list results = parallelWaitAll("testEntry",arguments); |
---|
2061 | } |
---|
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)); |
---|
2065 | fprintf(logfile,"%s out of %s = %sx%s decompositions are correct! (%s ms)%n", |
---|
2066 | sum(results),n*m,n,m,rtimer-t0,0); |
---|
2067 | } |
---|
2068 | |
---|
2069 | static proc readOutputTXT(string filename__) |
---|
2070 | { |
---|
2071 | system("--ticks-per-sec",1000); |
---|
2072 | dbprint(" reading matrix of decompositions from file "+filename__); |
---|
2073 | int t__=rtimer; |
---|
2074 | int tt__; |
---|
2075 | string data__ = read(":r "+filename__); |
---|
2076 | dbprint(sprintf(" done! (%s ms)", rtimer-t__)); |
---|
2077 | |
---|
2078 | dbprint(" processing input "); t__ = rtimer; |
---|
2079 | list mat__; |
---|
2080 | list nonlin__=list(); |
---|
2081 | int nonlin_factors_seperated__=0; |
---|
2082 | if(find(data__," * ")>0) {nonlin_factors_seperated__=1;} |
---|
2083 | int left__,right__=0,0; |
---|
2084 | left__ = find(data__,"{"); |
---|
2085 | int pos1__,pos2__=0,0; |
---|
2086 | int p1__,p2__; |
---|
2087 | int tmp__,tmp2__,depth__; |
---|
2088 | int i__,j__,k__,l__,max__; |
---|
2089 | intvec factors__,exponents__; |
---|
2090 | poly numerator__; |
---|
2091 | string s__,ss__; |
---|
2092 | for(i__=1;1;i__++) |
---|
2093 | { |
---|
2094 | tt__ = rtimer; |
---|
2095 | left__ = find(data__,"{",left__+1); |
---|
2096 | if(left__==0) {break;} |
---|
2097 | right__ = find(data__,"}",left__); |
---|
2098 | mat__[i__]=list(); |
---|
2099 | nonlin__[i__]=list(); |
---|
2100 | pos2__=left__; |
---|
2101 | for(j__=1;pos2__<right__&&pos2__>0;j__++) |
---|
2102 | { |
---|
2103 | pos1__ = pos2__+1; |
---|
2104 | pos2__ = find(data__,",",pos1__); |
---|
2105 | if(pos2__==0||pos2__>right__) {s__ = data__[pos1__,right__-pos1__];} |
---|
2106 | else {s__ = data__[pos1__,pos2__-pos1__];} |
---|
2107 | mat__[i__][j__] = list(); |
---|
2108 | |
---|
2109 | factors__ = intvec(0:0); |
---|
2110 | exponents__ = intvec(0:0); |
---|
2111 | l__ = find(s__," * "); |
---|
2112 | if(l__>0) |
---|
2113 | { |
---|
2114 | ss__ = s__[1,l__-1]; //ss__ contains the nonlinear factors |
---|
2115 | s__ = s__[l__+3,size(s__)-l__-2]; |
---|
2116 | for(p1__=1;s__[p1__]!="(";p1__++) {} |
---|
2117 | for(p2__=size(s__);s__[p2__]!=")";p2__--) {} |
---|
2118 | s__ = s__[p1__+1,p2__-p1__-1]; |
---|
2119 | l__ = find(ss__,"q"); |
---|
2120 | ss__ = ss__[l__,find(ss__,")",l__)-l__]; |
---|
2121 | ss__ = ss__+"*"; |
---|
2122 | p1__=0; |
---|
2123 | for(l__=1;1;l__++) |
---|
2124 | { |
---|
2125 | p1__ = find(ss__,"q",p1__+1); |
---|
2126 | if(p1__==0) {break;} |
---|
2127 | p1__++; |
---|
2128 | p2__ = find(ss__,"^",p1__); |
---|
2129 | tmp__ = find(ss__,"*",p1__); |
---|
2130 | if((p2__>tmp__ && tmp__>0) || (p2__==0)) //exponent is 1 |
---|
2131 | { |
---|
2132 | execute("factors__[l__]="+ss__[p1__,tmp__-p1__]); |
---|
2133 | exponents__[l__] = 1; |
---|
2134 | } |
---|
2135 | else |
---|
2136 | { |
---|
2137 | execute("factors__[l__]="+ss__[p1__,p2__-p1__]); |
---|
2138 | execute("exponents__[l__]="+ss__[p2__+1,tmp__-p2__-1]); |
---|
2139 | } |
---|
2140 | } |
---|
2141 | } |
---|
2142 | nonlin__[i__][j__] = list(factors__,exponents__); |
---|
2143 | |
---|
2144 | depth__ = 0; |
---|
2145 | s__=s__+" "; |
---|
2146 | max__ = size(s__); |
---|
2147 | tmp__ = 1; |
---|
2148 | tmp2__ = 0; |
---|
2149 | for(k__=1;k__<=max__;k__++) |
---|
2150 | { |
---|
2151 | if(s__[k__]=="(") {depth__++;k__++;continue;} |
---|
2152 | if(s__[k__]==")") {depth__--;k__++;continue;} |
---|
2153 | if(s__[k__]=="/" && depth__==0) {tmp2__ = k__;} |
---|
2154 | if((s__[k__]=="+" && depth__==0) || k__==max__) |
---|
2155 | { |
---|
2156 | if(tmp2__==0) // no denominator |
---|
2157 | { |
---|
2158 | execute("numerator__="+s__[tmp__,k__-tmp__]); |
---|
2159 | mat__[i__][j__][size(mat__[i__][j__])+1] = list(numerator__,intvec(0:0),intvec(0:0)); |
---|
2160 | tmp__ = k__+1; |
---|
2161 | k__++; continue; |
---|
2162 | } |
---|
2163 | execute("numerator__="+s__[tmp__,tmp2__-tmp__]); |
---|
2164 | ss__ = s__[tmp2__+1,k__-tmp2__-1]; |
---|
2165 | p1__ = find(ss__,"("); |
---|
2166 | p2__ = find(ss__,")"); |
---|
2167 | ss__ = ss__[p1__+1,p2__-p1__-1]; // now ss__ is only the denominator |
---|
2168 | |
---|
2169 | ss__ = ss__+"*"; |
---|
2170 | factors__ = intvec(0:0); |
---|
2171 | exponents__ = intvec(0:0); |
---|
2172 | p1__=0; |
---|
2173 | for(l__=1;1;l__++) |
---|
2174 | { |
---|
2175 | p1__ = find(ss__,"q",p1__+1); |
---|
2176 | if(p1__==0) {break;} |
---|
2177 | p1__++; |
---|
2178 | p2__ = find(ss__,"^",p1__); |
---|
2179 | tmp__ = find(ss__,"*",p1__); |
---|
2180 | if((p2__>tmp__ && tmp__>0) || (p2__==0)) //exponent is 1 |
---|
2181 | { |
---|
2182 | execute("factors__[l__]="+ss__[p1__,tmp__-p1__]); |
---|
2183 | exponents__[l__] = 1; |
---|
2184 | } |
---|
2185 | else |
---|
2186 | { |
---|
2187 | execute("factors__[l__]="+ss__[p1__,p2__-p1__]); |
---|
2188 | execute("exponents__[l__]="+ss__[p2__+1,tmp__-p2__-1]); |
---|
2189 | } |
---|
2190 | } |
---|
2191 | mat__[i__][j__][size(mat__[i__][j__])+1] = list(numerator__,factors__,exponents__); |
---|
2192 | tmp__ = k__+1; |
---|
2193 | tmp2__ = 0; |
---|
2194 | } |
---|
2195 | } |
---|
2196 | } |
---|
2197 | dbprint(sprintf(" row %s done! (%s ms)",i__,rtimer-tt__)); |
---|
2198 | } |
---|
2199 | dbprint(sprintf(" done! (%s ms)", rtimer-t__)); |
---|
2200 | |
---|
2201 | return(mat__,nonlin__); |
---|
2202 | } |
---|
2203 | |
---|
2204 | static proc readQfileTXT(string filename__) |
---|
2205 | { |
---|
2206 | string data__ = read(":r "+filename__); |
---|
2207 | ideal q__; |
---|
2208 | int pos1__,pos2__=1,1; |
---|
2209 | while(1) |
---|
2210 | { |
---|
2211 | pos1__=find(data__,"=",pos2__); |
---|
2212 | if(pos1__==0) {break;} |
---|
2213 | pos1__++; |
---|
2214 | pos2__=find(data__,";",pos1__); |
---|
2215 | execute("q__[size(q__)+1]="+data__[pos1__,pos2__-pos1__]); |
---|
2216 | } |
---|
2217 | return(q__); |
---|
2218 | } |
---|