share/solve_rec/solve_rec.mac

   1 /*****************************************************************************
   2  *                                                                           *
   3  * ************************************************************************* *
   4  * ***                                                                   *** *
   5  * ***                           * Solve_rec *                           *** *
   6  * ***                                                                   *** *
   7  * ***           This file implements some functions 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  *****************************************************************************/
  43
  44 eval_when(batch,
  45           ttyoff : true,
  46           nolabels : true)$
  47
  48 put('solve_rec, 1.08, 'version)$
  49
  50 define_variable(solve_rec_warn,        true,   bool)$
  51
  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)$
  57
  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)$
  63
  64 define_variable(solve_rec_method, 'solve_rec_method, any)$
  65 define_variable(shift_op, 'shift_op, any)$
  66
  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)$
  77
  78 sr_print([arg]) :=
  79   if solve_rec_warn then apply(print, arg)$
  80
  81 /**************************************************************************
  82  *
  83  * shift operator handling
  84  *
  85  *************************************************************************/
  86
  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$
  98
  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)))$
 103
 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)))$
 108
 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$
 127
 128 differ_order(expr) := max_shift(expr)-min_shift(expr)$
 129
 130 /**************************************************************************
 131  *
 132  * change difference equation to standard form with shift operator
 133  *
 134  *************************************************************************/
 135
 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)$
 150
 151
 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 )$
 158
 159
 160 /**************************************************************************
 161  *
 162  * Solver functions
 163  *
 164  *************************************************************************/
 165
 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),
 173
 174   std_form : subst(cond, eq),
 175
 176   std_form : to_standard_form(eq, %f, %n),
 177   std_form : expand(num(ratsimp(std_form))),
 178   ord : differ_order(std_form),
 179
 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),
 182
 183   if sol=false then false
 184   else (
 185     if listp(sol) then sol
 186     else fn=sol
 187   )
 188 )$
 189
 190 solve_rec_order_1(std_form, %f, %n, cond) := block(
 191   [deg : differ_degree(std_form), sol],
 192
 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   )
 201
 202   else if deg=2 then (
 203     return(ricatti_type(std_form, %f, %n, cond))
 204   )
 205
 206   else return(false)
 207 )$
 208
 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%],
 212
 213   if deg#1 then return(false),
 214
 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   ),
 223
 224   ihom : expand(std_form),
 225
 226   if cc=true then sol : solve_rec_lin_cc(coeffs1, ihom/coeffs[1], %f, %n, cond)
 227
 228   else if pc=true then (
 229
 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     )
 248
 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     )
 263
 264   ),
 265   sol
 266 )$
 267
 268 /**************************************************************************
 269  *
 270  * Get a solution from general solution, which
 271  * satisfies conditions in cond.
 272  *
 273  *************************************************************************/
 274
 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 )$
 287
 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 )$
 294
 295 solve_rec_ic1(gen_sol, %f, %n, cond, var) := block(
 296   [i, %d : length(cond), eqns : [], sol, used_vars : [], s, simpsum : true],
 297
 298   if length(cond)<length(var) then return(false),
 299
 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 )$
 310
 311 /**************************************************************************
 312  *
 313  * First order linear equations
 314  *
 315  *************************************************************************/
 316
 317 solve_rec_lin_1(%a, %b, %f, %n, cond) := block(
 318   [gen_sol, simpsum:true, ap, %j, %l, lo, prodhack:true],
 319
 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),
 327
 328   %b : subst(%j, %n, %b),
 329   %a : subst(%l, %n, %a),
 330
 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   ),
 343
 344   gen_sol : gen_sol + nusum(%b*ap, %j, lo, %n-1),
 345   solve_rec_method : "First order linear reccurence",
 346
 347   if length(cond)>0 then
 348     solve_rec_ic1(gen_sol, %f, %n, cond, [%k[1]])
 349   else
 350     gen_sol
 351 )$
 352
 353 get_lower_bound(expr, %n) := block(
 354   [ex : ratsimp(expr), sol, s, i, bound : 0, de, nu],
 355
 356   de : denom(ex),
 357   nu : num(ex),
 358
 359   sol : solve(de, %n),
 360
 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   ),
 367
 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   ),
 375
 376   bound
 377 )$
 378
 379 /**************************************************************************
 380  *
 381  * Higher order linear equations with constant coefficients
 382  *
 383  *************************************************************************/
 384
 385 solve_rec_lin_cc(coeffs, ihom, %f, %n, cond) := block(
 386   [gen_sol, hi, i],
 387
 388   gen_sol : solve_rec_char_gen(coeffs, %f, %n),
 389
 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,
 398
 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),
 404
 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));
 443
 444 solve_rec_char_gen(coeffs, %f, %n) := block(
 445   [deg : length(coeffs), %x, pol, sol, i, j:1, gen_sol : 0, %k, i%],
 446
 447   pol : sum(coeffs[i%]*%x^(deg-i%), i%, 1, deg),
 448   sol : solve(pol, %x),
 449
 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   ),
 456
 457   gen_sol
 458 )$
 459
 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 )$
 486
 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   )$
 493
 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 )$
 511
 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")$
 527
 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].
 533
 534   Examples:
 535
 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]]
 542
 543   else if op(expr)="+" then
 544     apply(append,
 545           map(lambda([u], solve_rec_get_exps(u, var)), args(expr)))
 546
 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))
 550
 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)]]))
 559
 560   else if op(expr) = "*" then
 561     [solve_rec_get_exps_mult (args (expr), var)]
 562
 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)]$
 578
 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.
 583
 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])))$
 623
 624 /**************************************************************************
 625  *
 626  * Special types of difference equations
 627  *
 628  *************************************************************************/
 629
 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 )$
 657
 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  *************************************************************************/
 665
 666 /**************************************************************************
 667  *
 668  * Algorithm Poly
 669  *
 670  *************************************************************************/
 671
 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 )$
 691
 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 )$
 725
 726 /**************************************************************************
 727  *
 728  * Algorithm Hyper
 729  *
 730  *************************************************************************/
 731
 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 )$
 757
 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 )$
 817
 818 monic_factors(p, %n) :=
 819   if hyper_factor_solve then monic_factors_sol(p, %n)
 820   else monic_factors_fac(p, %n)$
 821
 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),
 828
 829   /* This is a dirty trick to get a product */
 830   if member(part(factors, 0), ["^", "+"]) then factors : %%n*factors,
 831
 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 )$
 852
 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 )$
 866
 867 /**************************************************************************
 868  *
 869  * solve_rec_rat implements Abramov's algorithm for rational solutions of
 870  * linear reccurences with polynomial coefficients.
 871  *
 872  *************************************************************************/
 873
 874 /**************************************************************************
 875  *
 876  * Algorithm RatSol
 877  *
 878  *************************************************************************/
 879
 880
 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 )$
 900
 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 )$
 925
 926 /**************************************************************************
 927  *
 928  * Select linearly independent solutions
 929  *
 930  *************************************************************************/
 931
 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 )$
 940
 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 )$
 953
 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 )$
 970
 971 /**************************************************************************
 972  *
 973  * Reduce the order of equation
 974  *
 975  *************************************************************************/
 976
 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 )$
1003
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  *************************************************************************/
1025
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),
1034
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)),
1040
1041   if lo#1 then return(simplify_product(changevar(prod, %j1=%n-lo+1, %j1, %n))),
1042
1043   /* Check for simple cases.   */
1044   if term=%n then return(hi!/(lo-1)!),
1045
1046   if atom(term) or freeof(%n, term) then return(term^(hi-lo+1)),
1047
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],
1094
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 )$
1118
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$
1150
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 )$
1156
1157 /*************************************************************************
1158  *
1159  * Zeilberger package redefines union - we need a workaround.
1160  *
1161  *************************************************************************/
1162
1163 remove_duplicates(lst) := remove_duplicates1(lst, [])$
1164
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))$
1169
1170 eval_when(batch,
1171           ttyoff : false,
1172           nolabels : false)$