Enhancement to solve_rec.mac

I think that this extension should be included. I did not look at the
patch carefully though. Maybe solve_rec should return the answer as

> (%t5) c[2*n] = ...
> (%t6) c[2*n+1] = ...

instead of

> (%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)

Andrej



On Thu, Jan 6, 2011 at 9:50 PM, Mark H Weaver <mhw at netris.org> wrote:
> 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)
> ? ?)
> ?)$
>
> _______________________________________________
> Maxima mailing list
> Maxima at math.utexas.edu
> http://www.math.utexas.edu/mailman/listinfo/maxima
>