| blob | aa873cd428fdba40b6505d67559f0399adc0643c |
1 /*****************************************************************************
2 * *
3 * ************************************************************************* *
4 * *** *** *
5 * *** * Solve_rec * *** *
6 * *** *** *
7 * *** This file implements some funtions for solving *** *
8 * *** recurrence equations with maxima. In particular it *** *
9 * *** implements the Petkovsek's hyper algorithm and *** *
10 * *** Abramov's algorithm for solving linear recurrences with *** *
11 * *** polynomial coefficients. *** *
12 * *** *** *
13 * *** *** *
14 * *** Version: 1.08 *** *
15 * *** Author: Andrej Vodopivec <andrejv@users.sourceforge.net> *** *
16 * *** *** *
17 * ************************************************************************* *
18 * *
19 * *
20 * This package solves recurrences using: *
21 * - characteristic equation method for solving linear recurrences with *
22 * constant coefficients *
23 * - algorithm hyper for solving linear recurrences with polynomial *
24 * coefficients (only homogeneous equations) *
25 * - algorithm ratsol for obtaining rational solutions for linear *
26 * reccurences with polynomial coefficients (also inhomogeneous) *
27 * - solves ricatti type recurrences *
28 * *
29 * Usage: *
30 * - solve_rec(f[n]=f[n-1]+f[n-2], f[n]); *
31 * - solve_rec(2*x*(x+1)*y[x] - (x^2+3*x-2)*y[x+1] + (x-1)*y[x+2], *
32 * y[x], y[1]=1, y[3]=3); *
33 * *
34 * Top level functions: *
35 * - solve_rec *
36 * - solve_rec_poly *
37 * - solve_rec_rat *
38 * - solve_rec_hyper *
39 * - solve_rec_ic *
40 * - reduce_order *
41 * *
42 *****************************************************************************/
44 eval_when(batch,
45 ttyoff : true,
46 nolabels : true)$
48 put('solve_rec, 1.08, 'version)$
50 define_variable(solve_rec_warn, true, bool)$
52 define_variable(simplify_products, true, bool)$
53 define_variable(normalize_products, true, bool)$
54 define_variable(distribute_products, true, bool)$
55 define_variable(simplify_products_deg, 10, fixnum)$
56 define_variable(product_use_gamma, true, bool)$
58 define_variable(use_hyper, true, bool)$
59 define_variable(use_ratsol, true, bool)$
60 define_variable(hyper_factor_solve, false, bool)$
61 define_variable(hyper_all_solutions, true, bool)$
62 define_variable(hyper_to_product, true, bool)$
64 define_variable(solve_rec_method, 'solve_rec_method, any)$
65 define_variable(shift_op, 'shift_op, any)$
67 define_variable(%n, '%n, any)$
68 define_variable(%f, '%f, any)$
69 define_variable(%m, '%m, any)$
70 define_variable(%n, '%n, any)$
71 define_variable(%l, '%l, any)$
72 define_variable(%j, '%j, any)$
73 define_variable(%k, '%k, any)$
74 define_variable(%x, '%x, any)$
75 define_variable(%u, '%u, any)$
76 define_variable(%z, '%z, any)$
78 sr_print([arg]) :=
79 if solve_rec_warn then apply(print, arg)$
81 /**************************************************************************
82 *
83 * shift operator handling
84 *
85 *************************************************************************/
87 change_to_shift_op(expr, %f, %n) :=
88 if atom(expr) then expr
89 else if part(expr, 0)=%f then block(
90 [k% : args(expr)],
91 if numberp(first(k%)) then expr
92 else funmake('shift_op, [%f, %n, k%[1]-%n])
93 )
94 else if member(part(expr, 0), ["+", "*", "/", "^", "-"]) then
95 apply(part(expr, 0),
96 map(lambda([%u], change_to_shift_op(%u, %f, %n)), args(expr)))
97 else expr$
99 max_shift(expr) :=
100 if atom(expr) then -inf
101 else if part(expr, 0)=shift_op then part(expr, 3)
102 else apply('max, map('max_shift, args(expr)))$
104 min_shift(expr) :=
105 if atom(expr) then inf
106 else if part(expr, 0)=shift_op then part(expr, 3)
107 else apply('min, map('min_shift, args(expr)))$
109 differ_degree(expr) :=
110 if atom(expr) then 0
111 else if part(expr, 0)='shift_op then 1
112 else if part(expr, 0)="+" or part(expr, 0)="-" then
113 apply('max, map('differ_degree, args(expr)))
114 else if part(expr, 0)="*" then apply("+", map('differ_degree, args(expr)))
115 else if part(expr, 0)="^" then block(
116 [dp : differ_degree(part(expr, 2)), db : differ_degree(part(expr, 1))],
117 if dp#0 then inf
118 else if db=0 then 0
119 else db*part(expr, 2)
120 )
121 else if part(expr, 0)="/" then block(
122 [dd : differ_degree(denom(expr)), nd : differ_degree(num(expr))],
123 if dd#0 then inf
124 else nd
125 )
126 else 0$
128 differ_order(expr) := max_shift(expr)-min_shift(expr)$
130 /**************************************************************************
131 *
132 * change difference equation to standard form with shift operator
133 *
134 *************************************************************************/
136 to_standard_form1(expr, %n, %m) :=
137 if atom(expr) then (
138 if expr=%n then %n-%m
139 else expr
140 )
141 else if part(expr, 0)='shift_op then (
142 if part(expr, 2)=%n then funmake('shift_op, [part(expr, 1), part(expr, 2),
143 part(expr, 3)-%m])
144 else expr
145 )
146 else if member(part(expr, 0), ["+", "-", "*", "/", "^"]) then
147 apply(part(expr, 0), map(lambda([u], to_standard_form1(u, %n, %m)),
148 args(expr)))
149 else subst(%n=%n-%m, expr)$
152 to_standard_form(expr, %f, %n) := block(
153 [tmp : expand(lhs(expr)-rhs(expr)), m],
154 tmp : change_to_shift_op(tmp, %f, %n),
155 m : min_shift(tmp),
156 to_standard_form1(tmp, %n, m)
157 )$
160 /**************************************************************************
161 *
162 * Solver functions
163 *
164 *************************************************************************/
166 solve_rec(eq, fn, [cond]) := block(
167 [std_form, %f, %n, deg, ord, hyper_to_product : true, sol : false],
168 if atom(fn) then return(false),
169 if length(fn)#1 then return(false),
170 %f : part(fn, 0),
171 %n : part(fn, 1),
172 if not(atom(%f)) or not(mapatom(%n)) then return(false),
174 std_form : subst(cond, eq),
176 std_form : to_standard_form(eq, %f, %n),
177 std_form : expand(num(ratsimp(std_form))),
178 ord : differ_order(std_form),
180 if ord=1 then sol : solve_rec_order_1(std_form, %f, %n, cond)
181 else sol : solve_rec_gen(std_form, %f, %n, cond),
183 if sol=false then false
184 else (
185 if listp(sol) then sol
186 else fn=sol
187 )
188 )$
190 solve_rec_order_1(std_form, %f, %n, cond) := block(
191 [deg : differ_degree(std_form), sol],
193 if deg=1 then block(
194 [c2, c1, c0],
195 c2 : ratcoeff(std_form, funmake('shift_op, [%f, %n, 1])),
196 c1 : ratcoeff(std_form, funmake('shift_op, [%f, %n, 0])),
197 c0 : ratsimp(std_form - c2*funmake('shift_op, [%f, %n, 1])
198 - c1*funmake('shift_op, [%f, %n, 0])),
199 solve_rec_lin_1(expand(-c1/c2), expand(-c0/c2), %f, %n, cond)
200 )
202 else if deg=2 then (
203 return(ricatti_type(std_form, %f, %n, cond))
204 )
206 else return(false)
207 )$
209 solve_rec_gen(std_form, %f, %n, cond) := block(
210 [deg : differ_degree(std_form), coeffs : [], ihom, cc : true, coeffs1,
211 pc : true, sol : false, ord : differ_order(std_form), var, num_sol, i%],
213 if deg#1 then return(false),
215 coeffs : makelist(coeff(std_form, funmake('shift_op, [%f, %n, ord-i%])),
216 i%, 0, ord),
217 coeffs1 : makelist(ratsimp(coeffs[i%]/coeffs[1]), i%, 1, ord+1),
218 for i%:0 thru ord do (
219 std_form : std_form - coeffs[i%+1]*funmake('shift_op, [%f, %n, ord-i%]),
220 if not(freeof(%n, coeffs1[i%+1])) then cc : false,
221 if not(rec_polyp(coeffs[i%+1], %n)) then pc : false
222 ),
224 ihom : expand(std_form),
226 if cc=true then sol : solve_rec_lin_cc(coeffs1, ihom/coeffs[1], %f, %n, cond)
228 else if pc=true then (
230 if ihom=0 and use_hyper then (
231 sol : hyper(reverse(coeffs), %n),
232 if sol#false then block(
233 [i%],
234 sol : sort(remove_duplicates(sol)),
235 sol : get_linearly_independent1(sol, %n, ord),
236 num_sol : length(sol),
237 sol : sum(%k[i%]*sol[i%]/abs(numfactor(sol[i%])), i%, 1, length(sol)),
238 if num_sol<ord then (
239 sr_print("WARNING: found some hypergeometrical solutions!"),
240 if length(cond)>0 then
241 sol : false
242 )
243 else if length(cond)>0 then
244 sol : solve_rec_ic1(sol, %f, %n, cond,
245 makelist(%k[i%], i%, 1, length(cond)))
246 )
247 )
249 else if rec_polyp(ihom, %n) and use_ratsol then (
250 sol : solve_rec_rat1(reverse(coeffs), ihom, %n),
251 if sol#false then (
252 var : get_variables(sol),
253 if length(var)<ord then (
254 sr_print("WARNING: found some rational solutions!"),
255 if length(cond)#0 then
256 sol: false
257 )
258 else if length(cond)#0 then
259 sol : solve_rec_ic1(sol, %f, %n, cond,
260 makelist(var[i%], i%, 1, length(cond)))
261 )
262 )
264 ),
265 sol
266 )$
268 /**************************************************************************
269 *
270 * Get a solution from general solution, which
271 * satisfies conditions in cond.
272 *
273 *************************************************************************/
275 solve_rec_ic(eq, fn, cond) := block(
276 [%f, %n, gen_sol, i, c, var],
277 if part(eq, 0)="=" then (
278 eq : rhs(eq)-lhs(eq),
279 c : coeff(eq, fn),
280 eq : eq - c*fn
281 ),
282 %f : part(fn, 0),
283 %n : part(fn, 1),
284 var : get_variables(eq),
285 fn = solve_rec_ic1(eq, %f, %n, cond, var)
286 )$
288 get_variables(expr) := block(
289 [var : [], v],
290 for v in listofvars(expr) do
291 if not(atom(v)) and part(v, 0)=%k then var : append(var, [v]),
292 var
293 )$
295 solve_rec_ic1(gen_sol, %f, %n, cond, var) := block(
296 [i, %d : length(cond), eqns : [], sol, used_vars : [], s, simpsum : true],
298 if length(cond)<length(var) then return(false),
300 for i:1 thru %d do block(
301 [f_l : part(cond, i, 1), f_r : part(cond, i, 2), m],
302 eqns : append(eqns, [subst(part(f_l, 1), %n, gen_sol)=f_r]),
303 used_vars : append(used_vars, [var[i]])
304 ),
305 eqns : ev(eqns, product),
306 sol : linsolve(eqns, used_vars),
307 gen_sol : subst(sol, gen_sol),
308 gen_sol
309 )$
311 /**************************************************************************
312 *
313 * First order linear equations
314 *
315 *************************************************************************/
317 solve_rec_lin_1(%a, %b, %f, %n, cond) := block(
318 [gen_sol, simpsum:true, ap, %j, %l, lo, prodhack:true],
320 if length(cond)=0 then (
321 if %b#0 then
322 lo : get_lower_bound(%a*%b, %n)
323 else
324 lo : get_lower_bound(%a, %n)
325 )
326 else lo : part(cond, 1, 1, 1),
328 %b : subst(%j, %n, %b),
329 %a : subst(%l, %n, %a),
331 if freeof(%l, %a) then
332 ap : %a^(%n-%j-1)
333 else (
334 ap : minfactorial(simplify_product(product(%a, %l, %j+1, %n-1))),
335 ap : ap/numfactor(ap)
336 ),
337 if freeof(%l, %a) then
338 gen_sol : %k[1]*%a^(%n-lo)
339 else (
340 gen_sol : minfactorial(simplify_product(product(%a, %l, lo, %n-1))),
341 gen_sol : %k[1]*gen_sol/abs(numfactor(gen_sol))
342 ),
344 gen_sol : gen_sol + nusum(%b*ap, %j, lo, %n-1),
345 solve_rec_method : "First order linear reccurence",
347 if length(cond)>0 then
348 solve_rec_ic1(gen_sol, %f, %n, cond, [%k[1]])
349 else
350 gen_sol
351 )$
353 get_lower_bound(expr, %n) := block(
354 [ex : ratsimp(expr), sol, s, i, bound : 0, de, nu],
356 de : denom(ex),
357 nu : num(ex),
359 sol : solve(de, %n),
361 if length(sol)>0 then (
362 for i in sol do (
363 s : rhs(i),
364 if integerp(s) and s>= bound then bound : s+1
365 )
366 ),
368 sol : solve(nu, %n),
369 if sol#all and length(sol)>0 then (
370 for i in sol do (
371 s : rhs(i),
372 if integerp(s) and s>= bound then bound : s+1
373 )
374 ),
376 bound
377 )$
379 /**************************************************************************
380 *
381 * Higher order linear equations with constant coefficients
382 *
383 *************************************************************************/
385 solve_rec_lin_cc(coeffs, ihom, %f, %n, cond) := block(
386 [gen_sol, hi, i],
388 gen_sol : solve_rec_char_gen(coeffs, %f, %n),
390 if ihom#0 then block(
391 [exps : solve_rec_get_exps(expand(ihom), %n), curr_ihom, pol],
392 /*
393 exps : remove_duplicates(exps),
394 */
395 for ex in exps do block(
396 [base, pow, mul, %emode : false],
397 [pol, ex] : ex,
399 if not(numberp(ex)) then
400 if part(ex, 0)="^" then
401 pow : part(ex, 2)
402 else
403 pow : part(ex, 2, 2),
405 if hipow(pow, %n)>1 then curr_ihom : false
406 else (
407 mul : coeff(pow, %n),
408 if numberp(ex) then base:1
409 else
410 if part(ex, 0)="^" then
411 base : part(ex, 1)
412 else
413 base : 1/part(ex, 2, 1),
414 base : base^mul,
415 if not(freeof(%n, expand(pow-mul*%n))) then curr_ihom : false
416 else
417 curr_ihom:
418 solve_rec_char_part(coeffs, base, pol, ex*pol, %n)),
419 if curr_ihom#false then (
420 gen_sol : gen_sol + curr_ihom,
421 ihom : ratexpand(ihom-pol*ex)
422 )
423 )
424 ),
425 /*
426 Try quite hard to see if ihom = 0, since the a^i/b^i => (a/b)^i
427 step in solve_rec_get_exps might have outwitted the simplifier.
428 */
429 block ([ihom_zero: is (0 = radcan (ihom))],
430 if not(ihom_zero) then
431 if not (rec_polyp(ihom, %n)) then
432 return(false)
433 else
434 ihom : solve_rec_char_part(coeffs, 1, ihom, ihom, %n),
435 if ihom=false then
436 return(false)
437 else
438 gen_sol: gen_sol + ihom,
439 solve_rec_method: "Characteristic equation method",
440 if length(cond)>0 then
441 solve_rec_ic1(gen_sol, %f, %n, cond, makelist(%k[i], i, 1, length(cond)))
442 else gen_sol));
444 solve_rec_char_gen(coeffs, %f, %n) := block(
445 [deg : length(coeffs), %x, pol, sol, i, j:1, gen_sol : 0, %k, i%],
447 pol : sum(coeffs[i%]*%x^(deg-i%), i%, 1, deg),
448 sol : solve(pol, %x),
450 for i:1 thru length(sol) do block(
451 [k%],
452 gen_sol : gen_sol + rhs(sol[i])^%n*
453 sum(%k[j+k%]*%n^k%, k%, 0, multiplicities[i]-1),
454 j : j+multiplicities[i]
455 ),
457 gen_sol
458 )$
460 solve_rec_char_part(coeffs, po, ihom_pol, ihom, %n) := block(
461 [pol : expand(ihom_pol), hi, u, p%, %u, eqns, i%, char, k, i,
462 %d : length(coeffs)],
463 hi : hipow(pol, %n),
464 if rec_polyp(pol, %n)=false and hi>0 then return(false),
465 p% : sum(%u[i%]*%n^i%, i%, 0, hi),
466 char : sum(coeffs[length(coeffs)-i%]*%x^i%, i%, 0, length(coeffs)-1),
467 while subst(po, %x, char)=0 and char#0 do (
468 p% : p%*%n,
469 char : diff(char, %x)
470 ),
471 p% : expand(p%*po^%n),
472 eqns : 0,
473 for i:1 thru %d do eqns : eqns + subst(%n+%d-i, %n, p%)*coeffs[i],
474 eqns : expand(eqns) + ihom,
475 eqns : expand(eqns/po^%n),
476 eqns : makelist(ratcoeff(eqns, %n, i)=0, i, lopow(eqns,%n), hipow(eqns, %n)),
477 if errcatch(eqns : solve(eqns, makelist(%u[i], i, 0, hi)))=[] then
478 false
479 else (
480 if (listp(eqns[1])) then
481 subst(eqns[1], p%) /*at(p, eqns[1])*/
482 else
483 subst(eqns, p%)
484 )
485 )$
487 rec_polyp(p, n) :=
488 if freeof(n, p) then true
489 else block(
490 [hi : hipow(p, n)],
491 rec_polyp1(p, n, hi)
492 )$
494 rec_polyp1(p%, %n, hi) := block(
495 [c : ratcoeff(p%, %n, hi)],
496 if not(freeof(%n, hi)) then false
497 else if freeof(%n, c) then (
498 p% : expand(p%-c*%n^hi),
499 if hi=hipow(p%, %n) then false
500 else (
501 hi : hipow(p%,%n),
502 if hi=0 then
503 true
504 else
505 rec_polyp1(p%, %n, hi)
506 )
507 )
508 else
509 false
510 )$
512 /*
513 If expr is 1 or of the form u^v then return a pair [u,v] such that
514 u^v=expr. If not, raise an error.
515 */
516 solve_rec_split_exponent (expr) :=
517 if expr = 1 then [1, 0]
518 else if atom (expr) then
519 error ("Found atom ", expr, " when expecting an exponent")
520 else if (op (expr) = "/" and first (expr) = 1) then
521 block ([split: solve_rec_split_exponent (second (expr))],
522 [1/first (split), second (split)])
523 else if (op (expr) = "^") then
524 [first (expr), second (expr)]
525 else
526 error ("Found ", expr, " when expecting an exponent")$
528 /*
529 Try to return a list of pairs [u,v] such that the sum of u*v across
530 the list is equal to expr. Moreover, v is always of the form
531 a^(b*var). If a term can't be factored helpfully like this, it is
532 returned in the form [u, 1].
534 Examples:
536 solve_rec_get_exps(u*v^i, i) => [[u, v^i]]
537 solve_rec_get_exps(u*v^i + w, i) => [[w, 1], [u, v^i]]
538 solve_rec_get_exps(u^i/v^i, i) => [[1, (u/v)^i]]
539 */
540 solve_rec_get_exps(expr, var) :=
541 if atom(expr) then [[expr, 1]]
543 else if op(expr)="+" then
544 apply(append,
545 map(lambda([u], solve_rec_get_exps(u, var)), args(expr)))
547 /* Should only ever see unary minus, but test for length anyway */
548 else if (op(expr) = "-" and length(expr) = 1) then
549 map(lambda([u], [-u[1], u[2]]), solve_rec_get_exps(-expr, var))
551 else if op(expr)="^" then
552 block([bs: first (expr), pw: second (expr)],
553 if freeof(var, pw) or not (freeof (var, bs)) then
554 [[expr, 1]]
555 else block([a0,a1],
556 [a1,a0]: bothcoef(pw, var),
557 if not freeof(var, a0) then [[expr, 1]]
558 else [[bs^a0, bs^(a1*var)]]))
560 else if op(expr) = "*" then
561 [solve_rec_get_exps_mult (args (expr), var)]
563 else if op(expr) = "/" then
564 /*
565 We want to call solve_rec_get_exps_mult with
566 [first(expr),1/second(expr)], but have to avoid an infinite
567 recursion here: if first(expr) = 1, do it by hand.
568 */
569 if first (expr) = 1 then
570 map (lambda ([pair],
571 [1/first (pair),
572 if second (pair) = 1 then 1
573 else
574 first(second(pair))^(-second(second(pair)))]),
575 solve_rec_get_exps (second (expr), var))
576 else
577 [solve_rec_get_exps_mult([first (expr), 1/second (expr)], var)]$
579 /*
580 This is the hard case for solve_rec_get_exps. It deals with the
581 multiplication case, taking the terms of the product as the "terms"
582 variable.
584 It will always return exactly one pair.
585 */
586 solve_rec_get_exps_mult (terms, var) :=
587 block([pair_lists:
588 map (lambda ([term], solve_rec_get_exps (term, var)),
589 terms)],
590 /* If any of these isn't a singleton, we've failed already */
591 if some (lambda ([list], length(list) # 1), pair_lists) then
592 [xreduce("*", terms), 1]
593 /* Split results into "easy" (of form [u,1]) and hard (not!) */
594 else
595 block ([easy: 1, hard: []],
596 for pair in map (first, pair_lists) do
597 block ([split: solve_rec_split_exponent (second (pair))],
598 if second (split) = 0 then
599 easy: first (pair) * easy
600 else
601 hard: cons ([first (pair), split], hard)),
602 /* If everything was easy, we're good! */
603 if length (hard) = 0 then [xreduce("*", terms), 1]
604 /* If there was just one "hard" example, we can use that */
605 else if length (hard) = 1 then
606 [easy * first (hard [1]), apply("^", second (hard [1]))]
607 /*
608 We can deal with the case when all the hard examples are "the
609 same" in the sense that second(second(x)) is constant across
610 hard.
611 */
612 else
613 block([e: second (second (first (hard))),
614 others:
615 map(lambda ([x], second (second (x))), rest (hard))],
616 if not every (lambda ([ee], is (ee = e)), others) then
617 [xreduce("*", terms), 1]
618 else
619 [easy * xreduce ("*", map (first, hard)),
620 xreduce ("*",
621 map (lambda ([x],
622 first (second(x))), hard)) ^ e])))$
624 /**************************************************************************
625 *
626 * Special types of difference equations
627 *
628 *************************************************************************/
630 ricatti_type(eqn, %f, %n, cond) := block(
631 [q0, p1, p2, p3, %z, eq],
632 q0 : coeff(eqn, funmake('shift_op, [%f, %n, 0])),
633 q0 : coeff(q0, funmake('shift_op, [%f, %n, 1])),
634 if q0=0 then return(false),
635 eqn : eqn - q0*funmake('shift_op, [%f, %n, 0])*funmake('shift_op, [%f, %n, 1]),
636 eqn : expand(eqn/q0),
637 p1 : coeff(eqn, funmake('shift_op, [%f, %n, 1])),
638 p2 : coeff(eqn, funmake('shift_op, [%f, %n, 0])),
639 p3 : expand(eqn - p1*funmake('shift_op, [%f, %n, 1])
640 - p2*funmake('shift_op, [%f, %n, 0])),
641 if ratsimp(p1*p2-p3)=0 then (
642 solve_rec_method: "Ricatty equation.",
643 return([%f[%n]=-p1, %f[%n]=-subst(%n-1, %n, p2)])
644 ),
645 eq : %z[%n+2] + (p2-ratsubst(%n+1, %n, p1))*%z[%n+1] +
646 (p3 - p1*p2)*%z[%n],
647 %z : solve_rec(ratsimp(eq), %z[%n]),
648 if %z=false then return(false),
649 %z : normalize_product(subst(%n+1, %n, rhs(%z))/rhs(%z)) - p1,
650 %z : subst(1, %k[2], %z),
651 solve_rec_method : "Ricatti equation.",
652 if length(cond)>0 then
653 solve_rec_ic1(%z, %f, %n, cond, [%k[1]])
654 else
655 %z
656 )$
658 /**************************************************************************
659 *
660 * solve_rec_poly and solve_rec_hyper implement Petkovsek's Hyper
661 * algorithm for hypergeometric solutions of linear recurrences with
662 * polynomial coeffs.
663 *
664 *************************************************************************/
666 /**************************************************************************
667 *
668 * Algorithm Poly
669 *
670 *************************************************************************/
672 solve_rec_poly(eqn, fn) := block(
673 [std_form, %f, %n, deg, ord, coeffs, c, sol],
674 if atom(fn) then return(false),
675 if length(fn)#1 then return(false),
676 %f : part(fn, 0),
677 %n : part(fn, 1),
678 std_form : to_standard_form(eqn, %f, %n),
679 std_form : expand(num(ratsimp(std_form))),
680 ord : differ_order(std_form),
681 deg : differ_degree(std_form),
682 if deg#1 then return(false),
683 coeffs : makelist(coeff(std_form, funmake('shift_op, [%f, %n, i])), i, 0, ord),
684 for i:1 thru length(coeffs) do
685 std_form : std_form-coeffs[i]*funmake('shift_op, [%f, %n, i-1]),
686 std_form : expand(std_form),
687 sol : hyper_poly(coeffs, std_form, %n),
688 if sol#false then fn = sol
689 else false
690 )$
692 hyper_poly(coeffs, ihom, %n) := block(
693 [m : apply('max, map(lambda([u], hipow(u,%n)), coeffs)), s : -1, i%, j%,
694 d : length(coeffs), j, deg_poly : 0, i, sol, deg, %x, gen_p, %u, eqn,
695 val, solutions, linsolvewarn:false, linsolve_params:false, %k, vars],
696 while true do (
697 s : s+1,
698 j : 0,
699 deg_poly : coeff(coeffs[1], %n, m-s) +
700 sum(
701 binomial(%x, j%)*
702 sum(i%^j%*coeff(coeffs[i%+1], %n, m-s+j%), i%, 1, d-1),
703 j%, 0, s
704 ),
705 deg_poly : ratsimp(deg_poly),
706 if deg_poly#0 then return(true)
707 ),
708 sol : solve(deg_poly, %x),
709 deg : -1,
710 for j in sol do
711 if integerp(part(j,2)) and part(j,2)>deg then deg : part(j,2),
712 if ihom#0 then deg : max(deg, hipow(ihom, %n) - m + s),
713 if deg=-1 then return(false)
714 else (
715 gen_p : sum(%k[i%+1]*%n^i%, i%, 0, deg),
716 eqn : sum(coeffs[j%]*subst(%n+j%-1, %n, gen_p), j%, 1, d),
717 eqn : expand(eqn+ihom),
718 eqn : makelist(coeff(eqn, %n, i)=0, i, 0, hipow(eqn, %n)),
719 if errcatch(sol : linsolve(eqn, makelist(%k[i+1], i, 0, deg)))=[] then
720 return(false),
721 gen_p : expand(subst(sol, gen_p)),
722 gen_p
723 )
724 )$
726 /**************************************************************************
727 *
728 * Algorithm Hyper
729 *
730 *************************************************************************/
732 solve_rec_hyper(eqn, fn) := block(
733 [std_form, %f, %n, deg, ord, coeffs, c],
734 if atom(fn) then return(false),
735 if length(fn)#1 then return(false),
736 %f : part(fn, 0),
737 %n : part(fn, 1),
738 std_form : to_standard_form(eqn, %f, %n),
739 std_form : expand(num(ratsimp(std_form))),
740 ord : differ_order(std_form),
741 deg : differ_degree(std_form),
742 if deg#1 then return(false),
743 coeffs : makelist(coeff(std_form, funmake('shift_op, [%f, %n, i])), i, 0, ord),
744 for i:1 thru length(coeffs) do
745 std_form : std_form-coeffs[i]*funmake('shift_op, [%f, %n, i-1]),
746 std_form : expand(std_form),
747 if std_form#0 then return(false),
748 coeffs : hyper(coeffs, %n),
749 if coeffs=false then false
750 else (
751 if hyper_to_product then
752 coeffs : map(lambda([u], u/abs(numfactor(u))), coeffs),
753 coeffs : remove_duplicates(coeffs),
754 coeffs
755 )
756 )$
758 hyper(polys, %n) := block(
759 [d : length(polys), coeffs : [], Al, A, Bl, B, j, i, m, eqn, %z, i%,
760 solutions : [], have_solution : false, forbidden_deg_diff : []],
761 Al : expand(polys[1]/coeff(polys[1], %n, hipow(polys[1], %n))),
762 Al : monic_factors(Al, %n),
763 Bl : expand(subst(%n-d+2, %n, polys[d])),
764 Bl : expand(Bl/coeff(Bl, %n, hipow(Bl, %n))),
765 Bl : monic_factors(Bl, %n),
766 for A in Al do (
767 if have_solution=true and hyper_all_solutions=false then return(1),
768 for B in Bl do (
769 if have_solution=true and hyper_all_solutions=false then return(1),
770 if not(member(hipow(A, %n)-hipow(B, %n), forbidden_deg_diff)) then (
771 coeffs : makelist(expand(polys[i]*
772 my_prod(subst(%n+j-1, %n, A), j, 0, i-1)*
773 my_prod(subst(%n+j-1, %n, B), j, i, d-1)), i, 1, d),
774 m : apply('max, map(lambda([u], hipow(u, %n)), coeffs)),
775 eqn : sum(%z^(i%-1)*coeff(coeffs[i%], %n, m), i%, 1, d),
776 eqn : solve(eqn=0, [%z]),
777 eqn : map('rhs, eqn),
778 if (length(eqn)=1 and eqn[1]=0) or length(eqn)=0 then
779 if not(member(hipow(A, %n)-hipow(B, %n), forbidden_deg_diff)) then
780 forbidden_deg_diff : append([hipow(A, %n)-hipow(B, %n)],
781 forbidden_deg_diff),
782 for %z in eqn do block(
783 [poly_eqn, gc],
784 if %z=0 then return(false),
785 gc : expand(A*subst(%n+d-2, %n, B)),
786 poly_eqn : makelist(%z^i*coeffs[i+1]/gc, i, 0, d-1),
787 poly_eqn : hyper_poly(expand(ratsimp(poly_eqn)), 0, %n),
788 if poly_eqn#false then block(
789 [basis : [], vars : listofvars(poly_eqn), v, sol, lo],
790 for v in vars do (
791 if not(atom(v)) and part(v, 0)=%k then
792 basis : append(basis, [coeff(poly_eqn, v)])
793 ),
794 for sol in basis do (
795 sol : %z*A/B*subst(%n+1, %n, sol)/sol,
796 if hyper_to_product then (
797 sol : subst(%j, %n, sol),
798 lo : get_lower_bound(sol, %j),
799 sol : apply(nounify('product), [sol, %j, lo, %n-1]),
800 sol : minfactorial(simplify_product(sol)),
801 sol : sol/numfactor(sol)
802 ),
803 solutions : append(solutions, [sol]),
804 have_solution : true
805 )
806 )
807 )
808 )
809 )
810 ),
811 if length(solutions)=0 then false
812 else (
813 solve_rec_method : "Hyper algorithm",
814 factor(solutions)
815 )
816 )$
818 monic_factors(p, %n) :=
819 if hyper_factor_solve then monic_factors_sol(p, %n)
820 else monic_factors_fac(p, %n)$
822 monic_factors_fac(p, %n) := block(
823 [zeros, factors, sol : [1], f, d, tmp, i, pol, fac, lcoeff, %%n, deg],
824 factors : factor(p),
825 if atom(factors) then return([1, factors]),
826 if part(factors, 0)="/" then factors : part(factors, 1),
827 if part(factors, 0)="-" then factors : factors*(-1),
829 /* This is a dirty trick to get a product */
830 if member(part(factors, 0), ["^", "+"]) then factors : %%n*factors,
832 for f in factors do (
833 if not(freeof(%n, f)) then (
834 tmp : sol,
835 if not(atom(f)) and part(f, 0)="^" then (
836 deg : part(f, 2),
837 f : part(f, 1)
838 )
839 else
840 deg : 1,
841 for i:1 thru deg do (
842 for s in sol do (
843 fac : expand(s*f^i),
844 tmp : append(tmp, [expand(s*f^i/coeff(fac, %n, hipow(fac, %n)))])
845 )
846 ),
847 sol : tmp
848 )
849 ),
850 sol
851 )$
853 monic_factors_sol(p, %n) := block(
854 [zeros : [], factors],
855 factors : [1],
856 zeros : solve(p, %n),
857 zeros : map('rhs, zeros),
858 for i:1 thru length(zeros) do block(
859 [tmp : []],
860 for j:0 thru multiplicities[i] do
861 for z in factors do tmp : append(tmp, [expand(z*(%n-zeros[i])^j)]),
862 factors : tmp
863 ),
864 factors
865 )$
867 /**************************************************************************
868 *
869 * solve_rec_rat implements Abramov's algorithm for rational solutions of
870 * linear reccurences with polynomial coefficients.
871 *
872 *************************************************************************/
874 /**************************************************************************
875 *
876 * Algorithm RatSol
877 *
878 *************************************************************************/
881 solve_rec_rat(eqn, fn) := block(
882 [std_form, %f, %n, deg, ord, coeffs, c, hyper_to_product : false],
883 if atom(fn) then return(false),
884 if length(fn)#1 then return(false),
885 %f : part(fn, 0),
886 %n : part(fn, 1),
887 std_form : to_standard_form(eqn, %f, %n),
888 std_form : expand(num(ratsimp(std_form))),
889 ord : differ_order(std_form),
890 deg : differ_degree(std_form),
891 if deg#1 then return(false),
892 coeffs : makelist(coeff(std_form, funmake('shift_op, [%f, %n, i])), i, 0, ord),
893 for i:1 thru length(coeffs) do
894 std_form : std_form-coeffs[i]*funmake('shift_op, [%f, %n, i-1]),
895 std_form : expand(std_form),
896 coeffs : solve_rec_rat1(coeffs, std_form, %n),
897 if coeffs=false then false
898 else fn = factor(coeffs)
899 )$
901 solve_rec_rat1(coeffs, ihom, %n) := block(
902 [A, B, ord : length(coeffs), res, %y, i, M : -1, d, u : 1, sol, eq, i%,
903 gc, coeffs1, prodhack:true],
904 coeffs1 : makelist(subst(%n-i+1, %n, coeffs[i]), i, 1, ord),
905 gc : apply('ezgcd, coeffs1)[1],
906 u : gc,
907 A : expand(subst(%n-ord+1, %n, coeffs[ord]))/gc,
908 B : expand(coeffs[1])/gc,
909 res : resultant(A, expand(subst(%n+%y, %n, B)), %n),
910 res : map('rhs, solve(res, %y)),
911 for i in res do
912 if integerp(i) then M : max(M, i),
913 for i:M thru 0 step -1 do (
914 d : gcd(A, subst(%n+i, %n, B)),
915 A : ratsimp(A/d),
916 B : ratsimp(B/subst(%n-i, %n, d)),
917 for j:0 thru i do u : u*subst(%n-j, %n, d)
918 ),
919 eq : sum(coeffs[i%]/subst(%n+i%-1, %n, u)*%y[%n+i%-1], i%, 1, ord) + ihom,
920 sol : solve_rec_poly(eq, %y[%n]),
921 if sol=false then return(false),
922 solve_rec_method : "RatSol algorithm",
923 rhs(sol)/u
924 )$
926 /**************************************************************************
927 *
928 * Select linearly independent solutions
929 *
930 *************************************************************************/
932 casoratian_determinant(l, %n) := block(
933 [ml, M, i, j, n : length(l)],
934 ml : makelist(
935 makelist(subst(%n+j-1, %n, l[i]), i, 1, n),
936 j, 1, n),
937 M : apply('matrix, ml),
938 determinant(M)
939 )$
941 get_linearly_independent(sol, %n) := block(
942 [lin_in : [sol[1]], s, tst, d, m, dd],
943 for s:2 thru length(sol) do (
944 d : casoratian_determinant(append(lin_in, [sol[s]]), %n),
945 d : makefact(d),
946 d : minfactorial(d),
947 d : factcomb(d),
948 d : fullratsimp(d),
949 if d#0 then lin_in : append(lin_in, [sol[s]])
950 ),
951 lin_in
952 )$
954 get_linearly_independent1(sol, %n, ord) := block(
955 [lin_in : [sol[1]], mat, %%i, %%j, mat1],
956 mat : makelist(subst(%n+%%j-1, %n, sol[1]), %%j, 1, ord),
957 mat : apply('matrix, [mat]),
958 mat : transpose(mat),
959 for s:2 thru length(sol) do (
960 mat1 : addcol(mat,
961 minfactorial(makefact(makelist(subst(%n+%%j-1, %n, sol[s]), %%j, 1, ord)))
962 ),
963 if rank(mat1)=length(lin_in)+1 then (
964 lin_in : cons(sol[s], lin_in),
965 mat : mat1
966 )
967 ),
968 lin_in
969 )$
971 /**************************************************************************
972 *
973 * Reduce the order of equation
974 *
975 *************************************************************************/
977 reduce_order(eqn, sol, fn) := block(
978 [%u, %z, eq, hi, lo, i, c, sumsplitfact : false, %j, %f, %n, e1, e2],
979 %f : part(fn, 0),
980 %n : part(fn, 1),
981 eq : to_standard_form(rhs(eqn)-lhs(eqn), %f, %n),
982 hi : max_shift(eq),
983 lo : min_shift(eq),
984 for i:lo thru hi do (
985 eq : subst(subst(%n+i, %n, sol)*funmake('shift_op, [%z, %n, i]),
986 funmake('shift_op, [%f, %n, i]), eq)
987 ),
988 for i:hi thru lo+1 step -1 do (
989 c : ratcoeff(eq, funmake('shift_op, [%z, %n, i])),
990 eq : expand(ratsubst(%u[%n+i-1]+funmake('shift_op, [%z, %n, i-1]),
991 funmake('shift_op, [%z, %n, i]), eq))
992 ),
993 eq : factcomb(eq),
994 c : ratcoeff(eq, funmake('shift_op, [%z, %n, 0])),
995 if c#0 then return(false),
996 eq : expand(eq - c*funmake('shift_op, [%z, %n, 0])),
997 e1 : %f[%n] = sol*%z[%n],
998 ldisplay(e1),
999 e2 : %z[%n] = 'sum(%u[%j], %j, 0, %n-1),
1000 ldisplay(e2),
1001 num(ratsimp(eq))
1002 )$
1004 /**************************************************************************
1005 *
1006 * Simplify products - sometimes we can get rational solutions from
1007 * hypergeometrical if we can simplify products well.
1008 *
1009 * Maxima does not handle products well - this simplifies product in some
1010 * cases.
1011 *
1012 * - simplify_product simplifies products like
1013 *
1014 * product(n/(n+3), n, 1, k) ==> 6/((k+1)*(k+2)(k+3))
1015 *
1016 * to disable automatic simplification, set simplify_products to false
1017 *
1018 * - normalize product changes the upper bound for products like
1019 *
1020 * product(n^2-1, n, 1, k-1) ==> product(n^2-1, n, 1, k)/(k^2-1)
1021 *
1022 * to disable automatic normalization, set normalize_products to false
1023 *
1024 *************************************************************************/
1026 simplify_product(prod) := block(
1027 [term, %n, lo, hi, %kk : 0, nu, de, comm : [], p% : 1,
1028 deg : simplify_products_deg, %j1],
1029 if simplify_products=false then return(prod),
1030 if atom(prod) then return(prod),
1031 if member(part(prod, 0), ["+", "-", "*", "/", "^"]) then
1032 return(apply(part(prod, 0), map('simplify_product, args(prod)))),
1033 if part(prod, 0)#nounify('product) then return(prod),
1035 /* Read product arguments. */
1036 %n : part(prod, 2),
1037 lo : part(prod, 3),
1038 hi : part(prod, 4),
1039 term : factor(part(prod, 1)),
1041 if lo#1 then return(simplify_product(changevar(prod, %j1=%n-lo+1, %j1, %n))),
1043 /* Check for simple cases. */
1044 if term=%n then return(hi!/(lo-1)!),
1046 if atom(term) or freeof(%n, term) then return(term^(hi-lo+1)),
1048 /* Distribute over products. */
1049 if part(term, 0)="*" then (
1050 p%*apply("*", map(lambda([u], simplify_product(apply(nounify('product),
1051 [u, %n, lo, hi]))),
1052 args(term)))
1053 )
1054 /* Take care of fractions. */
1055 else if part(term, 0)="/" then block(
1056 [nu : num(term), de : denom(term)],
1057 /* Check for cancellations. */
1058 for %kk:-deg thru deg do block(
1059 [g],
1060 g : gcd(expand(subst(%n+%kk, %n, nu)), de),
1061 if not(freeof(%n, g)) then (
1062 comm : append(comm, [[%kk, g, subst(%n-%kk, %n, g)]]),
1063 de : ratsimp(de/g),
1064 nu : ratsimp(nu/subst(%n-%kk, %n, g))
1065 )
1066 ),
1067 /* Cancel common terms. */
1068 for c in comm do block(
1069 [kk : c[1], d],
1070 if kk<0 then (
1071 p% : p%*my_prod(c[3], %n, hi+kk+1, hi),
1072 d : my_prod(c[2], %n, lo, lo-kk-1),
1073 p% : p%/d
1074 )
1075 else (
1076 p% : p%*1/my_prod(c[2], %n, hi-kk+1, hi),
1077 d : my_prod(c[3], %n, lo, lo+kk-1),
1078 p% : p%*d
1079 )
1080 ),
1081 /* Distribute over fractions. */
1082 nu : simplify_product(apply(nounify('product), [nu, %n, lo, hi])),
1083 de : simplify_product(apply(nounify('product), [de, %n, lo, hi])),
1084 p%*nu/de
1085 )
1086 else if part(term, 0)="^" and integerp(part(term, 2)) then (
1087 simplify_product(apply(nounify('product), [part(term, 1), %n, lo, hi]))^
1088 part(term, 2)
1089 )
1090 /* Assume we have a poly. */
1091 else block(
1092 [aa, bb, lcoeff, %m],
1093 bb : bothcoef(term, %n), aa : bb[1], bb : bb[2],
1095 /* Take care of linear products. */
1096 if freeof(%n, aa) and freeof(%n, bb) then (
1097 if aa=1 then (
1098 %m : term - %n,
1099 (hi+%m)!/(lo+%m-1)!
1100 )
1101 else if aa=-1 then (
1102 %m : term + %n,
1103 (%m-lo)!/(%m-hi-1)!
1104 )
1105 else if product_use_gamma and integerp(aa) then
1106 gamma(subst(hi+1, %n, expand(term/aa)))/
1107 gamma(subst(lo, %n, expand(term/aa)))*aa^(hi-lo+1)
1108 else
1109 prod
1110 )
1111 else if part(term, 0)="-" and length(args(term))=1 then
1112 (-1)^(hi-lo+1)*simplify_product(apply(nounify('product), [-term, %n, lo, hi]))
1113 /* Give up! */
1114 else
1115 p%*apply(nounify('product), [factor(term), %n, lo, hi])
1116 )
1117 )$
1119 normalize_product(%expr) :=
1120 if atom(%expr) or not(normalize_products) then %expr
1121 else if part(%expr, 0)=nounify('product) then block(
1122 [%term : part(%expr, 1),
1123 %var : part(%expr, 2),
1124 %lo : part(%expr, 3),
1125 %hi : part(%expr, 4),
1126 %vars],
1127 %vars : listofvars(%hi),
1128 /* We check if upper limit is const. + numb. - then we normalize to const. */
1129 if length(%vars)=1 and hipow(%hi, %vars[1])=1 then block(
1130 [%co : coeff(%hi, %vars[1]), %dd, %r],
1131 %dd : %hi - %co*%vars[1],
1132 if integerp(%dd) and %dd#0 then (
1133 if %dd>0 then
1134 apply(nounify('product), [%term, %var, %lo, %vars[1]]) *
1135 my_prod(subst(%vars[1]+%r, %var, %term), %r, 1, %dd)
1136 else if %dd<0 then
1137 apply(nounify('product), [%term, %var, %lo, %vars[1]]) /
1138 my_prod(subst(%vars[1]-%r, %var, %term), %r, 0, -%dd-1)
1139 )
1140 else
1141 %expr
1142 )
1143 else
1144 %expr
1145 )
1146 else if member(part(%expr, 0), ["+", "-", "*", "/", "^", "="]) then
1147 apply(part(%expr, 0), map(' normalize_product, args(%expr)))
1148 else
1149 %expr$
1151 my_prod(term, n, lo, hi) := block(
1152 [pr : 1],
1153 for i:lo thru hi do pr : pr*subst(i, n, term),
1154 pr
1155 )$
1157 /*************************************************************************
1158 *
1159 * Zeilberger package redefines union - we need a workaround.
1160 *
1161 *************************************************************************/
1163 remove_duplicates(lst) := remove_duplicates1(lst, [])$
1165 remove_duplicates1(lst1, lst2) :=
1166 if length(lst1)=0 then lst2
1167 else if member(first(lst1), lst2) then remove_duplicates1(rest(lst1), lst2)
1168 else remove_duplicates1(rest(lst1), cons(first(lst1), lst2))$
1170 eval_when(batch,
1171 ttyoff : false,
1172 nolabels : false)$