Enhancement to solve_rec.mac

How about returning an explicit executable definition with a conditional,
e.g.

     c[n] := if mod(n,2)=0 then ... elseif mod(n,2)=1 then ...


On Fri, Jan 7, 2011 at 09:39, Andrej Vodopivec
<andrej.vodopivec at gmail.com>wrote:

> 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
> >
> _______________________________________________
> Maxima mailing list
> Maxima at math.utexas.edu
> http://www.math.utexas.edu/mailman/listinfo/maxima
>