Hello Andrej and users of solve_rec.mac,
A few months ago I tried to solve the following recurrence with
solve_rec, but it failed:
(%i1) display2d:false$
(%i2) load(solve_rec);
(%o2) "/usr/local/share/maxima/5.23.0/share/contrib/solve_rec/solve_rec.mac"
(%i3) solve_rec( c[n+1]=-n*c[n-1]/(n+1), c[n], c[0]=1, c[1]=2 );
(%o3) false
However, I noticed that this recurrence is actually two independent
recurrences: one recurrence for odd n, and one for even n. If I break
them apart by hand, solve_rec is able to solve each one.
I added code to solve_rec to do this job automatically:
(%i4) load(solve_rec_mhw);
(%o4) "/home/mhw/.maxima/solve_rec_mhw.mac"
(%i5) solve_rec( c[n+1]=-n*c[n-1]/(n+1), c[n], c[0]=1, c[1]=2 );
(%t5) c_even[n] = -gamma(n/2+1/2)*(-1)^(n/2-1)/(sqrt(%pi)*(n/2)!)
(%t6) c_odd[n] = -sqrt(%pi)*((n-1)/2)!*(-1)^((n-1)/2-1)/gamma((n-1)/2+3/2)
(%o6) [%t5,%t6]
It can split recurrences not only into evens and odds, but into any
finite number of independent recurrences:
(%i7) solve_rec( c[n+1]=-n*c[n-2]/(n+1), c[n], c[0]=1, c[1]=2, c[2]=3 );
(%t7) c_0[n] = -gamma(n/3+2/3)*(-1)^(n/3-1)/(gamma(2/3)*(n/3)!)
(%t8) c_1[n] = -2*gamma(1/3)*((n-1)/3)!*(-1)^((n-1)/3-1)
/(3*gamma((n-1)/3+4/3))
(%t9) c_2[n] = -6*gamma(2/3)*gamma((n-2)/3+4/3)*(-1)^((n-2)/3-1)
/(gamma(1/3)*gamma((n-2)/3+5/3))
(%o9) [%t7,%t8,%t9]
My patch follows. Comments and suggestions welcome. I'd especially be
interested in possible suggestions on how best to return these results
to the caller.
Best,
Mark
*** share/contrib/solve_rec/solve_rec-orig.mac 2009-10-05 03:21:50.000000000 -0400
--- share/contrib/solve_rec/solve_rec.mac 2010-11-02 21:14:32.000000000 -0400
***************
*** 14,19 ****
--- 14,25 ----
* *** Version: 1.08 *** *
* *** Author: Andrej Vodopivec <andrejv at users.sourceforge.net> *** *
* *** *** *
+ * *** Modified by Mark H Weaver <mhw at netris.org> in November 2010 *** *
+ * *** to solve recurrences in which the gcd of the index shifts is *** *
+ * *** greater than 1, such that the recurrence is equivalent to *** *
+ * *** multiple independent recurrences of lesser order, e.g.: *** *
+ * *** solve_rec( c[n+1]=-n*c[n-1]/(n+1), c[n], c[0]=1, c[1]=2 ); *** *
+ * *** *** *
* ************************************************************************* *
* *
* *
***************
*** 96,101 ****
--- 102,112 ----
map(lambda([%u], change_to_shift_op(%u, %f, %n)), args(expr)))
else expr$
+ gcd_shift(expr) :=
+ if atom(expr) then 0
+ else if part(expr, 0)=shift_op then part(expr, 3)
+ else first(apply('ezgcd, map('gcd_shift, args(expr))))$
+
max_shift(expr) :=
if atom(expr) then -inf
else if part(expr, 0)=shift_op then part(expr, 3)
***************
*** 156,161 ****
--- 167,217 ----
to_standard_form1(tmp, %n, m)
)$
+ /**************************************************************************
+ *
+ * reduce the order of a recurrence whose shifts have a gcd > 1
+ *
+ *************************************************************************/
+
+ srec_reduce_order(expr) :=
+ if atom(expr) then (
+ if expr=%n then (%nn * %mult) + %offset
+ else expr
+ )
+ else if part(expr, 0)='shift_op then (
+ if part(expr, 2)=%n
+ then funmake('shift_op, [%ff, %nn, part(expr, 3) / %mult])
+ else expr
+ )
+ else if member(part(expr, 0), ["+", "-", "*", "/", "^"])
+ then apply(part(expr, 0), map('srec_reduce_order, args(expr)))
+ else subst(%n = (%nn * %mult) + %offset, expr)$
+
+ srec_reduce_order_cond1(cond) :=
+ if part(cond,0)#"=" or part(cond,1,0)#%f or
+ length(part(cond,1))#1 or not integerp(part(cond,1,1))
+ then error("Invalid solve_rec condition", cond)
+ else block ([newidx : (part(cond,1,1) - %offset) / %mult],
+ if integerp(newidx)
+ then [ substpart(%ff, substpart(newidx,cond,1,1), 1, 0) ]
+ else []
+ )$
+
+ srec_reduce_order_conds(conds) :=
+ apply('append, map(srec_reduce_order_cond1, conds))$
+
+ srec_finalize_reduced_order_sol(expr) :=
+ if atom(expr) then (
+ if expr=%nn then (%n - %offset) / %mult
+ else expr
+ )
+ else if part(expr,0)=%ff then (
+ srec_finalize_reduced_order_sol(part(expr,1)) * %mult + %offset,
+ substpart(%%, fn, 1),
+ substpart(%fff, %%, 0)
+ )
+ else if mapatom(expr) then expr
+ else map('srec_finalize_reduced_order_sol, expr)$
/**************************************************************************
*
***************
*** 164,170 ****
*************************************************************************/
solve_rec(eq, fn, [cond]) := block(
! [std_form, %f, %n, deg, ord, hyper_to_product : true, sol : false],
if atom(fn) then return(false),
if length(fn)#1 then return(false),
%f : part(fn, 0),
--- 220,227 ----
*************************************************************************/
solve_rec(eq, fn, [cond]) := block(
! [std_form, %f, %n, deg, ord, %mult,
! hyper_to_product : true, sol : false],
if atom(fn) then return(false),
if length(fn)#1 then return(false),
%f : part(fn, 0),
***************
*** 175,189 ****
std_form : to_standard_form(eq, %f, %n),
std_form : expand(num(ratsimp(std_form))),
! ord : differ_order(std_form),
! if ord=1 then sol : solve_rec_order_1(std_form, %f, %n, cond)
! else sol : solve_rec_gen(std_form, %f, %n, cond),
! if sol=false then false
! else (
! if listp(sol) then sol
! else fn=sol
)
)$
--- 232,272 ----
std_form : to_standard_form(eq, %f, %n),
std_form : expand(num(ratsimp(std_form))),
! %mult : gcd_shift(std_form),
! if %mult=1 then (
! ord : differ_order(std_form),
! if ord=1 then sol : solve_rec_order_1(std_form, %f, %n, cond)
! else sol : solve_rec_gen(std_form, %f, %n, cond),
! if sol=false then false
! else (
! if listp(sol) then sol
! else fn=sol
! )
! ) else block(
! [ %nn : ?gensym(), %ff : ?gensym(), %fff,
! reduced_form, reduced_cond, solns : [] ],
!
! for %offset:0 thru %mult-1 do (
! if %mult # 2 then %fff : concat(%f, "_", %offset)
! else if %offset = 0 then %fff : concat(%f, "_even")
! else %fff : concat(%f, "_odd"),
!
! reduced_form : srec_reduce_order(std_form),
! reduced_cond : srec_reduce_order_conds(cond),
!
! ord : differ_order(reduced_form),
!
! if ord=1 then sol : solve_rec_order_1(reduced_form, %ff, %nn, reduced_cond)
! else sol : solve_rec_gen(reduced_form, %ff, %nn, reduced_cond),
!
! if sol#false and not listp(sol) then sol : arraymake(%ff,[%nn])=sol,
! sol : srec_finalize_reduced_order_sol(sol),
! solns : cons(sol, solns)
! ),
! solns : reverse(solns),
! apply('ldisp, solns)
)
)$