share/solve_rec/simplify_sum.mac

   1 /*****************************************************************************
   2  *                                                                           *
   3  * ************************************************************************* *
   4  * ***                                                                   *** *
   5  * ***                       * simplify_sum *                            *** *
   6  * ***                                                                   *** *
   7  * ***   This file implements simplify_sum function for simplification   *** *
   8  * ***   of sums.                                                        *** *
   9  * ***   The methods used by simplify_sum:                               *** *
  10  * ***     1. internal simplification using simpsum=true                 *** *
  11  * ***     2. Gosper algorithm                                           *** *
  12  * ***     3. rational summation with conversion to psi functions        *** *
  13  * ***     4. Zeilberger algorithm                                       *** *
  14  * ***     5. conversion to hypergeometrical sums, reduced by hgfred.    *** *
  15  * ***                                                                   *** *
  16  * ***   Author:  Andrej Vodopivec <andrejv@users.sourceforge.net>       *** *
  17  * ***   Licence: GPL                                                    *** *
  18  * ***                                                                   *** *
  19  * ************************************************************************* *
  20  *                                                                           *
  21  *                                                                           *
  22  * Demo:                                                                     *
  23  *                                                                           *
  24  * (%i6) s : simplify_sum(sum(((-1)^k*x*binom(n,k))/(x+k),k,0,n))$           *
  25  * (%i7) factcomb(s);                                                        *
  26  * (%o7) (n!*x!)/(x+n)!                                                      *
  27  *                                                                           *
  28  * (%i8) simplify_sum(sum((-1)^k*binom(x-2*k,n-k)*binom(x-k+1,k),k,0,n));    *
  29  * Is  x - 2 * n   positive, negative, or zero?  pos;                        *
  30  * Is  x - n + 1   positive, negative, or zero?  pos;                        *
  31  * (%o8) ((-1)^n+1)/2                                                        *
  32  *                                                                           *
  33  * (%i9) simplify_sum(sum(binom(n,k)/(k+1),k,0,n));                          *
  34  * (%o9) (2*2^n-1)/(n+1)                                                     *
  35  *                                                                           *
  36  * More examples with load("simplify_sum_test");                             *
  37  *                                                                           *
  38  *****************************************************************************/
  39
  40 load("zeilberger")$
  41 load("solve_rec/solve_rec")$
  42 load("simplifying.lisp")$
  43 load("opsubst")$
  44
  45 eval_when(batch,
  46           ttyoff : true,
  47           nolabels : true)$
  48
  49 put('simplify_sum, 1.0, 'version)$
  50
  51 define_variable(zeilberger_check, true, boolean)$
  52 define_variable(sum_min, 0, any)$
  53 define_variable(use_simpsum, true, boolean)$
  54 define_variable(use_harmonic, true, boolean)$
  55 define_variable(use_integral, true, boolean)$
  56 define_variable(use_gosper, true, boolean)$
  57 define_variable(use_ratfun, true, boolean)$
  58 define_variable(use_zeilberger, true, boolean)$
  59 define_variable(use_hypergeometric, true, boolean)$
  60 define_variable(use_telescoping, true, boolean)$
  61
  62 define_variable(simplify_sum_depth, 0, fixnum)$
  63 define_variable(simplify_sum_max_depth, 3, fixnum)$
  64
  65 declare(simplify_sum, evfun)$
  66 declare(split_sum, evfun)$
  67
  68 /*******************
  69  *
  70  * Debugging.
  71  *
  72  *******************/
  73
  74 define_variable(verbose_level, 0, fixnum)$
  75
  76 ss_print_message(level, [mess]) :=
  77   if verbose_level>=level then (
  78     apply(print, mess))$
  79
  80 /*******************
  81  *
  82  * simplify_sum(expr) : tries to simplify all sums in expr.
  83  *
  84  *******************/
  85
  86 simplify_sum(expr) :=
  87   /* Check if we have a sum - if not then simplify all arguments.       */
  88   if mapatom(expr) then expr
  89   else if part(expr, 0)#nounify(sum) then map(simplify_sum, expr)
  90
  91
  92   /* Check for sum(expr, var, min, inf) */
  93   else if (part(expr, 3)=-inf or part(expr, 3)=minf) and part(expr, 4)=inf then block(
  94     [summand, var_, lo, hi],
  95     [summand, var_, lo, hi]: args(expr),
  96     simplify_sum(apply('sum, [summand, var_, 0, inf])) +
  97     simplify_sum(apply('sum, [subst(var_=-var_, summand), var_, 1, inf])))
  98
  99   /* Check for sum(expr, var, minf, hi) */
 100   else if part(expr, 3)=-inf or part(expr, 3)=minf then (
 101     simplify_sum(apply('sum, [subst(var_=-var_, summand), var_, -hi, inf])))
 102
 103   /* Trigonometric sums */
 104   else if catch(
 105     scanmap(
 106       lambda([u], if not atom(u) and member(part(u,0), [sin, cos, sinh, cosh]) then throw(true)),
 107       expr,
 108       bottomup))=true then block(
 109     [sm1 : simplify_sum(split_sum(expand(exponentialize(expr))))],
 110     if freeof_sum(sm1) then trigsimp(rectform(sm1))
 111     else expr)
 112
 113   /* Sums with Fibonacci numbers */
 114   else if catch(
 115     scanmap(
 116       lambda([u], if not atom(u) and part(u,0)=fib then throw(true)),
 117       expr,
 118       bottomup))=true then block(
 119     [sm1: simplify_sum(split_sum(expand(fibtophi(expr))))],
 120     if freeof_sum(sm1) then block(
 121       [algebraic:true],
 122       ratsimp(sm1)))
 123
 124   /* Check if summand is a logarighm */
 125   else if not(atom(part(expr, 1))) and part(expr, 1, 0)=log and part(expr, 3)#minf then block(
 126     [%up, %sub],
 127     if part(expr, 4)=inf then (
 128       %sub: log(simplify_product(apply(product,
 129             [exp(part(expr, 1)), part(expr, 2), part(expr, 3), %up]))),
 130       %sub: limit(%sub, %up, inf),
 131       if not(freeof(inf, %sub) and freeof(minf, %sub)) then error("Sum is divergent!")
 132       elseif freeof_limit(%sub) then %sub
 133       else expr)
 134     else (
 135       log(simplify_product(apply(product, cons(exp(part(expr, 1)), rest(args(expr))))))))
 136
 137   /* Main summation code starts here */
 138   else block(
 139     [simpsum:false,
 140      summand : ratsimp(part(expr, 1)),
 141      var : part(expr, 2),
 142      lo : part(expr, 3),
 143      hi : part(expr, 4),
 144      sm1 : expr,
 145      linsolvewarn:false,
 146      simplify_sum_depth:simplify_sum_depth+1,
 147      ss_new_context: concat('ss_context, simplify_sum_depth)],
 148
 149     /*** Prevent possible infinite recursion. ***************************/
 150     if simplify_sum_depth>simplify_sum_max_depth then return(expr),
 151
 152     /*** Simplify the summand *******************************************/
 153
 154     /* make sure we dont have contexts from previous runs */
 155     if member(ss_new_context, contexts) then
 156       killcontext(ss_new_context),
 157
 158     supcontext (ss_new_context),
 159     if lo#minf and lo#-inf then assume(var>=lo),
 160     if hi#inf then assume(var<=hi),
 161     summand : simplify_sum(summand),
 162     killcontext(ss_new_context),
 163     expr : intosum(apply(sum, [summand, var, lo, hi])),
 164
 165     /*** Default maxima simplification. *********************************/
 166     if use_simpsum=true then (
 167       ss_print_message(1, "Trying with simpsum=true ..."),
 168       sm1 : ev(apply(sum, [summand, var, lo, hi]), simpsum = true),
 169       ss_print_message(1, "sum with simpsum=true returns:", sm1)),
 170
 171     /*** Sums with harmonic_number(n) by parts. *************************/
 172     if use_harmonic=true and not freeof_sum(sm1) and hi#inf then block(
 173       [a, b, harmonic_number_args: get_harmonic_number_args(summand), harmonic_part],
 174       if length(harmonic_number_args)>0 then (
 175         harmonic_part : apply(gen_harmonic_number, first(harmonic_number_args)),
 176         [a, b]:bothcoef(summand, harmonic_part),
 177         ss_print_message(1, "Trying with sum_by_parts"),
 178         ss_print_message(2, "harmonic_number by parts"),
 179         ss_print_message(2, "Coefficients", a, b),
 180         if a#0 then
 181           sm1: simplify_sum(sum_by_parts(summand-b, a, var, lo, hi)) +
 182                simplify_sum(apply('sum, [b, var, lo, hi]))),
 183         ss_print_message(1, "sum_by_parts returns", sm1)),
 184
 185     /*** Try with integral representation *******************************/
 186     if use_integral=true and not freeof_sum(sm1) then (
 187       ss_print_message(1, "Trying with integral representation."),
 188       sm1 : sum_by_integral(summand, var, lo, hi),
 189       ss_print(1, "Integral representation returns", sm1)),
 190
 191
 192     /*** Now let's try the Gosper algorithm *****************************/
 193     if errcatch(
 194       if use_gosper and not(atom(sm1)) and not freeof_sum(sm1) then block(
 195         [hi1 : ?gensym(), lo1 : ?gensym()],
 196         ss_print_message(1, "Trying with Gosper ..."),
 197         sm1 : block([ttyoff:ttyoff], nusum(summand, var, lo1, hi1)),
 198         if freeof_sum(sm1) then (
 199           ss_print_message(1, "Gosper returns:", sm1),
 200           if lo=minf or lo=-inf then
 201             sm1 : limit(sm1, lo1, lo)
 202           else
 203             sm1 : subst(lo1=lo, sm1),
 204           if hi=inf then
 205             sm1 : limit(sm1, hi1, hi)
 206           else
 207             sm1 : subst(hi1=hi, sm1),
 208           if not freeof_limit(sm1) or not freeof('und, sm1) then sm1: expr)
 209         else sm1 : expr)) = [] then sm1 : expr,
 210
 211
 212     /*** Try the extended Gosper algorithm ******************************/
 213     if errcatch(
 214       if use_gosper and not(atom(sm1)) and not freeof_sum(sm1) then block(
 215         [hi1 : ?gensym(), lo1 : ?gensym()],
 216         ss_print_message(1, "Trying with extended Gosper ..."),
 217         sm1 : block([ttyoff:ttyoff], extended_nusum(summand, var, lo1, hi1)),
 218         if freeof_sum(sm1) then (
 219           ss_print_message(1, "Extended Gosper returns:", sm1),
 220           if lo=minf or lo=-inf then
 221             sm1 : limit(sm1, lo1, lo)
 222           else
 223             sm1 : subst(lo1=lo, sm1),
 224           if hi=inf then
 225             sm1 : limit(sm1, hi1, hi)
 226           else
 227             sm1 : subst(hi1=hi, sm1),
 228           if not freeof_limit(sm1) then sm1 : expr)
 229         else sm1 : expr)) = [] then sm1 : expr,
 230
 231     /*** Try converting the sum to finite sum of psi functions **********/
 232     if use_ratfun and not(atom(sm1)) and not( freeof_sum(sm1) ) and ?ratp(summand, var)
 233     then block(
 234       [ratfun, expr1: expr],
 235       ss_print_message(1, "Trying ratfun -> psi ..."),
 236       ratfun : part(expr1, 1),
 237       if hi=inf and hipow(num(ratfun), var) > hipow(denom(ratfun), var)-2 then
 238         error("Sum is divergent!"),
 239       polypart : first(divide(num(ratfun), denom(ratfun), var)),
 240       expr1 : intosum(substpart(ratsimp(ratfun - polypart), expr1, 1)),
 241       ss_print_message(3, "Polynomial part", polypart),
 242       ss_print_message(3, "Without polynomial part", expr1),
 243       /* polynomial part */
 244       sm1 : simplify_sum(apply(sum, [polypart, var, lo, hi])),
 245       ss_print_message(4, "Polynomial part contributes", sm1),
 246       /* fractional part */
 247       ratfun : part(expr1, 1),
 248       sm1 : sm1 + ratfun_to_psi(ratfun, var, lo, hi)),
 249
 250     /*** check for sum((-1)^n*ratfun, n, lo, inf); ****************************/
 251     if use_ratfun and not(atom(sm1)) and not( freeof_sum(sm1) ) and
 252       ?ratp(ratsimp(summand/(-1)^var), var) and hi = inf
 253     then block(
 254       [ratfun, new_var : ?gensym(), expr1],
 255       ss_print_message(1, "Trying ratfun -> psi ..."),
 256       expr1 : intosum(changevar(expr, new_var=var-lo+1, new_var, var)),
 257       ratfun : part(expr1, 1)/(-1)^new_var,
 258       if hipow(num(ratfun), new_var) <= hipow(denom(ratfun), new_var)-1 then (
 259         sm1 : ratfun_to_psi(subst(new_var=2*new_var, ratfun), new_var, 1, inf) -
 260               ratfun_to_psi(subst(new_var=2*new_var-1, ratfun), new_var, 1, inf))
 261       else sm1 : und),
 262
 263     supcontext (ss_new_context),
 264
 265     /*** Zeilberger algorithm ************************************************/
 266     if errcatch(
 267       if use_zeilberger and not(atom(sm1)) and not freeof_sum(sm1) then block(
 268         [summand:summand, var:var, lo:lo, hi:hi, nv: gensym(), supp],
 269         ss_print_message(1, "Trying with Zeilberger ..."),
 270         assume(var>lo), assume(var<hi),
 271         support : ss_support(summand, var),
 272         if atom(support[1]) and support[1]#minf then
 273           [summand, var, lo, hi]: args(intosum(
 274               changevar(apply('sum, [summand, var, lo, hi]), nv=var-support[1]+sum_min, nv, var)
 275               ))
 276         else if lo#sum_min then
 277           [summand, var, lo, hi]: args(intosum(
 278               changevar(apply('sum, [summand, var, lo, hi]), nv=var-lo+sum_min, nv, var)
 279               )),
 280         sm1 : ss_zeilb(summand, var, lo, hi),
 281         ss_print_message(1, "Zeilberger method returns:", sm1)
 282       )
 283     ) = [] then sm1 : expr,
 284
 285     killcontext(ss_new_context),
 286
 287     if sm1=false then
 288       sm1 : expr,
 289
 290     /*** Check if upper bound can be inf.  **********************************
 291     block(
 292       [expr1: factor(minfactorial(factcomb(makefact(summand)))),
 293        support],
 294       support: ss_support(expr1, var),
 295       if support[2]<=hi then hi_hyper: inf),
 296     **************************************************************************/
 297
 298     /*** Convert to hypergeometrical functions. ******************************/
 299     if not(atom(sm1)) and not( freeof_sum(sm1) ) and
 300        use_hypergeometric=true and hi=inf and lo#-inf and lo#minf then (
 301          ss_print_message(1, "Converting to hypergeometrical sum ..."),
 302          sm1 : to_hypergeometric(summand, var, lo, hi),
 303          ss_print_message(1, "hgfred method returns:", sm1)),
 304
 305
 306     if sm1=false then
 307       sm1 : expr,
 308
 309     if not(atom(sm1)) and not( freeof_sum(sm1) ) and use_telescoping=true then (
 310       ss_primt_message(1, "Using telescoping ..."),
 311       sm1: telescoping_sum(summand, var, lo, hi),
 312       ss_primt_semmace(1, "telescoping method returns:", sm1)),
 313
 314     if sm1=false then
 315       sm1 : expr,
 316
 317     sm1)$
 318
 319 /*******************
 320  *
 321  * Check if we still have some sums in expr.
 322  *
 323  *******************/
 324
 325 freeof_sum(expr) :=
 326   if atom(expr) then true
 327   else if part(expr, 0)=nounify(sum) then false
 328   else if expr=[] then true
 329   else xreduce("and", map(freeof_sum, args(expr)))$
 330
 331 freeof_integrate(expr) :=
 332   if atom(expr) then true
 333   else if part(expr, 0)=nounify(integrate) then false
 334   else if expr=[] then true
 335   else xreduce("and", map(freeof_integrate, args(expr)))$
 336
 337 freeof_limit(expr) :=
 338   if atom(expr) then true
 339   else if part(expr, 0)=nounify(limit) then false
 340   else if expr=[] then true
 341   else xreduce("and", map(freeof_limit, args(expr)))$
 342
 343 /*******************
 344  *
 345  * The extended Gosper algorithm
 346  *
 347  *******************/
 348
 349 ss_linearp(expr, var) := block(
 350   [a,b],
 351   [a,b]: bothcoeff(expand(expr), var),
 352   if freeof(var, a) and freeof(var, b) then [a,b]
 353   else [])$
 354
 355 find_coeffs(expr, var) := block([coeffs],
 356   coeffs: ss_linearp(expr, var),
 357   if coeffs=[] then (
 358     if member(part(expr, 0), ["+", "-", "*", "/", gamma, "!", binomial])
 359       then apply(append, map(lambda([e], find_coeffs(e, var)), args(expr)))
 360     else if part(expr, 0)="^" and freeof(var, part(expr, 2))
 361       then find_coeffs(part(expr, 1), var)
 362     else if part(expr, 0)="^" and freeof(var, part(expr, 1))
 363       then find_coeffs(part(expr, 2), var)
 364     else
 365       error())
 366   else
 367   [coeffs[1]])$
 368
 369 find_mfold(expr, var) := block(
 370   [coeffs: find_coeffs(expr, var)],
 371   coeffs: map(denom, coeffs),
 372   xreduce(lambda([a,b], a*b/gcd(a,b)), coeffs))$
 373
 374 extended_nusum(expr, var, lo, hi) := block(
 375   [%m%: find_mfold(expr, var), exprm, tk, sk],
 376   exprm: subst(var=%m%*var, expr),
 377   tk: nusum(exprm, var, 1, var),
 378   if not freeof_sum(tk) then error(),
 379   sk: subst(var=var/%m%, tk),
 380   apply('sum, [sk, var, hi-%m%+1, hi]) - apply('sum, [sk, var, lo-%m%, lo-1]))$
 381
 382 /*******************
 383  *
 384  * This simplifies the sum using Zeilberger
 385  *
 386  *******************/
 387
 388 ss_zeilb(expr, %k%, lo, hi, [in_zr]) := block(
 389   [vars : delete(%k%, listofvars(expr)), %n%, eq, sm, zb, deg,
 390    cond, eq_rhs, cert, %i%, expr1, ihom, support, sum_min:sum_min,
 391    upper_bound_implicit : false, lower_bound_implicit : false,
 392    solve_rec_warn : false, warnings:false, Gosper_in_Zeilberger:false,
 393    cont],
 394
 395   if lo=-inf then lo:minf,
 396
 397   /* Convert binomials to factorials. */
 398   expr1 : factor(minfactorial(factcomb(makefact(expr)))),
 399   ss_print_message(2, "Summand:", expr),
 400   ss_print_message(2, "Changed to:", expr1),
 401
 402   /* We need expr to be hypergeometric in at least two variables. */
 403   /* We prefer the second variable to appear in bounds. */
 404   if length(in_zr)=0 then (
 405     if length(vars)<1 then (
 406       ss_print(3, "Not enough variables"),
 407       return(false)),
 408     if length(listofvars(hi))>0 then %n% : first(listofvars(hi))
 409     else if length(listofvars(lo))>0 then %n% : first(listofvars(lo))
 410     else %n% : vars[1])
 411   else
 412     %n% : in_zr[1],
 413
 414   /*  Find support. */
 415   assume(%k%>lo), assume(%k%<hi),
 416   support : ss_support(expr1, %k%),
 417   support : [if numberp(support[1]) then ceiling(support[1]) else support[1],
 418              if numberp(support[2]) then floor(support[2])   else support[2]],
 419
 420   if support[2]<lo or support[1]>hi then return(0),
 421
 422   /* Check if bounds are implicit. */
 423   if ss_max(lo, support[1])=support[1] then lower_bound_implicit : true,
 424   if ss_min(hi, support[2])=support[2] then upper_bound_implicit : true,
 425
 426   if lo=minf then lo : support[1],
 427   if hi=inf then hi : support[2],
 428
 429   ss_print_message(2, "Found support:", support),
 430   if numberp(support[1]) then sum_min:max(support[1], lo),
 431
 432   /* We don't handle sums over infinite support yet! */
 433   if (support[1]=minf and lower_bound_implicit) or
 434      (support[2]=inf and upper_bound_implicit) then (
 435     ss_print_message(3, "Support not finite!"),
 436     return(false)),
 437
 438   /* Find the recurrence for the sum. */
 439   zb : Zeilberger(expr1, %k%, %n%),
 440   ss_print_message(3, "Zeilberger returns:", zb),
 441   if length(zb)=0 then error(),
 442   if not(listp(zb[1])) then error(),
 443   deg : length(part(zb, 1, 2)),
 444   cert : part(zb, 1, 1),
 445   eq : part(zb, 1, 2) . makelist(sm[%n%+%i%], %i%, 0, deg-1),
 446
 447   /* Find the initial conditions for recurrence. */
 448   cont: first(content(eq)),
 449   while subst(%n%=sum_min, cont)=0 do sum_min: sum_min+1,
 450
 451   if subst(%n%=max(support[1], lo), cont)#0 or
 452      errcatch(
 453        cond : ss_zeilb_init(expr, %n%,
 454          max(support[1], lo), min(support[2], hi), max(support[1], lo), deg))=[] or
 455      not freeof_sum(cond) then (
 456        cond : ss_zeilb_init(expr, %n%, lo, hi, sum_min, deg))
 457   else (
 458     sum_min: max(support[1], lo)),
 459
 460   /* Initial conditions should not contain sums. */
 461   if not freeof_sum(cond) then (
 462     ss_print_message(3, "Wrong initial conditions:", cond),
 463     return(false)),
 464
 465   /* Find the right hand side of the sum recurrence. */
 466   ihom : minfactorial(makefact(cert*expr)),
 467   eq_rhs : 0,
 468
 469   if not(upper_bound_implicit) then (
 470     for %i%:0 thru deg-1 do (
 471       eq_rhs : eq_rhs + part(zb, 1, 2, %i%+1)*
 472                   apply(sum, [subst(%n%=%n%+%i%, expr), %k%, hi+1, subst(%n%=%n%+%i%, hi)])),
 473     eq_rhs : eq_rhs + subst(%k%=hi+1, ihom)),
 474
 475   if not(lower_bound_implicit) then (
 476     for %i%:0 thru deg-1 do (
 477       eq_rhs : eq_rhs + part(zb, 1, 2, %i%+1)*
 478                   apply(sum, [subst(%n%=%n%+%i%, expr), %k%, subst(%n%=%n%+%i%, lo), lo-1])),
 479     eq_rhs : eq_rhs - subst(%k%=lo, ihom)),
 480
 481   eq_rhs : factor(minfactorial(factcomb(makefact(eq_rhs)))),
 482
 483   /* Right hand side must be free of sums. */
 484   if not freeof_sum(eq_rhs) then (
 485     ss_print_message(3, "Recurrence contains sums!", eq_rhs),
 486     return(false)),
 487
 488   eq : eq = eq_rhs,
 489   ss_print_message(2, "Degree of recurrence:", deg-1),
 490   ss_print_message(2, "Zeilberger recurrence:", eq),
 491   ss_print_message(2, "Initial conditions:", cond),
 492
 493   if length(in_zr)>0 then return(eq),
 494
 495   /* Solve the recurrence. */
 496   eq : apply(solve_rec, append([eq, sm[%n%]], cond)),
 497   ss_print_message(2, "Solving recurrence returns:", eq),
 498   if eq=false then return(false),
 499   if not freeof_sum(eq) then return(false),
 500
 501   eq : ratsimp(minfactorial(makefact(rhs(eq)))),
 502   ss_print_message(4, "Simplified solution:", eq),
 503
 504   /* Check the solution. */
 505   if check_sum(expr, %k%, %n%, lo, hi, deg, eq) then eq
 506   else false)$
 507
 508 ss_zeilb_init(expr, %n%, lo, hi, sum_min, deg) := block(
 509   [cond],
 510   cond : makelist(
 511     sm[%i%] =
 512     subst(%n%=%i%,
 513       apply('sum, [
 514         minfactorial(expr),
 515         %k%, lo, hi
 516       ]
 517       )
 518     ),
 519     %i%, sum_min, sum_min + deg - 2),
 520   cond : factor(minfactorial(simplify_sum(cond))),
 521   cond)$
 522
 523 /*******************
 524  *
 525  * This returns the recurrence for the sum.
 526  *
 527  *******************/
 528
 529 summand_to_rec(expr, k, n) := block(
 530   [zr, linsolvewarn:false, lo:minf, hi:inf],
 531
 532   if listp(k) then (
 533     lo : k[2],
 534     hi : k[3],
 535     k : k[1]),
 536
 537   supcontext('ss_context),
 538
 539   if errcatch(
 540     zr : ss_zeilb(expr, k, lo, hi, n)
 541   ) = [] then zr : 'failed,
 542
 543   killcontext('ss_context),
 544
 545   zr)$
 546
 547 /*******************
 548  *
 549  * Check the result - we may get something wrong!
 550  *
 551  *******************/
 552
 553 check_sum(expr, %k%, %n%, lo, hi, deg, sm) :=
 554   if not zeilberger_check then true
 555   else catch(
 556     block(
 557       [%i%, tmp_sum, real_sum, simpsum:true, sum_min:sum_min+deg],
 558       for %i%:sum_min thru deg+sum_min do(
 559
 560         real_sum : minfactorial(factcomb(makefact(
 561               simplify_sum(subst(%n%=%i%, apply(sum, [expr, %k%, lo , hi])))
 562               ))),
 563
 564         tmp_sum  : minfactorial(factcomb(makefact(subst(%n%=%i%, sm)))),
 565         dif : rectform(real_sum - tmp_sum),
 566         if not freeof_sum(real_sum) or (numberp(dif) and dif#0) then (
 567           if dif#0 then (
 568             ss_print_message(2, "Sum check failed with: ", 'i=%i%, 'tmp_sum=tmp_sum, 'dif=dif),
 569             throw(false))
 570           else
 571           print("Warning: sum check could not be completed!")))),
 572     true)$
 573
 574 /*******************
 575  *
 576  * This part checks for the support of expr in %k%
 577  *
 578  *******************/
 579
 580 ss_support(expr, %k%) := block(
 581   if freeof(%k%, expr) or expr=%k% then [minf, inf]
 582   else if member(part(expr, 0), ["+", "-"]) then
 583     lreduce(ss_union, map(lambda([u], ss_support(u,%k%)), args(expr)))
 584   else if part(expr, 0)="*" then
 585     lreduce(ss_intersection, map(lambda([u], ss_support(u,%k%)), args(expr)))
 586   else if part(expr, 0)="/" then
 587     lreduce(ss_intersection, map(lambda([u], ss_support(factcomb(u),%k%)), args(expr)))
 588   else if member(part(expr, 0), ["^"]) then ss_support(part(expr, 1), %k%)
 589   else if part(expr, 0)=binomial then ss_support_binomial(expr, %k%)
 590   else if part(expr, 0)="!" then ss_support_factorial(expr, %k%)
 591   else [minf, inf])$
 592
 593 /* we assume here that k!=0 if k<0 */
 594 ss_support_factorial(expr, %k%) :=
 595   solve_lin_ineq(part(expr,1)>0, %k%)$
 596
 597 ss_support_binomial(expr, %k%) := block(
 598   [s1, s2],
 599   s1 : solve_lin_ineq(part(expr, 2)>0, %k%),
 600   s2 : solve_lin_ineq(part(expr, 2)<part(expr, 1), %k%),
 601   ss_intersection(s1, s2))$
 602
 603 /*******************
 604  *
 605  * solves inequality which is linear in k
 606  *
 607  *******************/
 608
 609 solve_lin_ineq(eq, k) := block(
 610   [eq1, bc, a, b],
 611
 612   if op(eq)=">" then eq1 : lhs(eq) - rhs(eq)
 613   else eq1 : rhs(eq) - lhs(eq),
 614
 615   bc : bothcoef(expand(eq1), k), a:bc[1], b:bc[2],
 616
 617   if a=0 then [minf, inf]
 618   else if not freeof(k, b) then [minf, inf]
 619   else if not numberp(a) then [minf, inf]
 620   else if a>0 then [-b/a, inf]
 621   else [minf, -b/a])$
 622
 623 ss_union(l1, l2) := [ss_min(first(l1), first(l2)), ss_max(second(l1), second(l2))]$
 624 ss_intersection(l1, l2) := [ss_max(first(l1), first(l2)), ss_min(second(l1), second(l2))]$
 625
 626 ss_max(e1, e2) := block(
 627   [pnz : asksign(e1-e2)],
 628   if pnz='pos or pnz='zero then e1
 629   else e2)$
 630
 631 ss_min(e1, e2) := block(
 632   [pnz : asksign(e1-e2)],
 633   if pnz='neg or pnz='zero then e1
 634   else e2)$
 635
 636 /********************
 637  *
 638  * splits sum
 639  *
 640  ********************/
 641
 642 expand_sum(expr, k%%, lo%%, hi%%) :=
 643   if not(atom(expr)) and part(expr, 0)="+" then
 644     map(lambda([u], apply(sum, [u, k%%, lo%%, hi%%])), expr)
 645   else
 646     apply(sum, [expr, k%%, lo%%, hi%%]);
 647
 648 split_sum(expr) :=
 649   block([sm%%: opsubst(nounify(sum)=expand_sum, expr)],
 650     ev(sm%%, expand_sum))$
 651
 652 /*******************
 653  *
 654  * converts sum(ratfun, var, 1, inf) to psi functions when denom(ratfun)
 655  * can be completely factored with gfactor.
 656  *
 657  *******************/
 658
 659 ratfun_to_psi(ratfun, var, lo, hi) := block(
 660   [pf, sum: 0, denom_f : factor_with_solve(denom(ratfun), var), max_root],
 661
 662   pf : partfrac(num(ratfun)/denom_f, var),
 663
 664   max_root: lmax(sublist(map(rhs, solve(denom_f, var)), numberp)),
 665   if (numberp(max_root) and numberp(lo)) then (
 666     while (lo < max_root) do (
 667       sum: sum+subst(var=lo, ratfun),
 668       lo: lo+1)),
 669
 670   if inpart(pf, 0)="+" then pf: args(pf)
 671   else pf: [pf],
 672
 673   ss_print_message(2, "Partial fractions", pf),
 674
 675   for prt in pf do block(
 676     [term: numfactor(prt), exponent, a, b],
 677     local(a),
 678
 679     ss_print_message(3, "Working on term", prt),
 680
 681     prt : prt/term,
 682     term : term*num(prt),
 683     prt : denom(prt),
 684
 685     if not(atom(prt)) and part(prt,0)="^" then (
 686       exponent: part(prt, 2),
 687       prt: part(prt, 1))
 688     else
 689       exponent: 1,
 690
 691     ss_print_message(3, "Linear part", prt),
 692
 693     a: ratsimp(bothcoef(expand(prt), var)),
 694     b: a[2], a: a[1],
 695
 696     if not( freeof(var, a) and freeof(var, b) ) then error(),
 697     if hi=inf then (
 698       if exponent#1 then
 699         term: term*(zeta(exponent) -
 700                     gen_harmonic_number(exponent, subst(var=lo-1, prt/a)))/a^exponent
 701       else
 702         term: -term*gen_harmonic_number(exponent, subst(var=lo-1, prt/a))/a^exponent
 703     )
 704     else
 705       term: term*(gen_harmonic_number(exponent, subst(var=hi, prt/a)) -
 706                   gen_harmonic_number(exponent, subst(var=lo-1, prt/a)))/a^exponent,
 707
 708     ss_print_message(3, "Corresponding term in sum", term),
 709
 710     sum: sum+term),
 711
 712   sum)$
 713
 714 factor_with_solve(expr, n) := block(
 715   [sol, fac, expr1],
 716   sol : solve(expr, n),
 717   expr : ratexpand(expr),
 718   fac : ratcoef(expr, n, hipow(expr, n)),
 719   for i:1 thru length(sol) do (
 720     if not(freeof(n, rhs(sol[i]))) then error(),
 721     fac : fac * (n - rhs(sol[i]))^multiplicities[i]),
 722   fac)$
 723
 724 /*******************
 725  *
 726  * Reduce using hgfred
 727  *
 728  *******************/
 729
 730 to_hypergeometric(expr, var, lo, hi) := block(
 731   [quo, upper, lower, a:[], b:[], x, c, warnings:false],
 732
 733   for i:1 thru 100 while subst(var=lo, expr)=0 do lo : lo+1,
 734   expr:subst(var=var+lo, expr),
 735
 736   quo : ratsimp(shiftQuo(factor(makefact(expr)*var!), var)),
 737   if not(?ratp(quo, var)) then return(false),
 738
 739   ss_print_message(2, "Shift quotient", quo),
 740
 741   upper : -map(rhs, solve(num(quo), var)),
 742   if not(every(lambda([u], freeof(var, u)), upper)) then return(false),
 743   for i:1 thru length(upper) do (
 744     for j:1 thru multiplicities[i] do
 745       a : cons(upper[i], a)),
 746   ss_print_message(2, "a=", a),
 747
 748   lower : -map(rhs, solve(denom(quo), var)),
 749   if not(every(lambda([u], freeof(var, u)), lower)) then return(false),
 750   for i:1 thru length(lower) do (
 751     for j:1 thru multiplicities[i] do
 752       b : cons(lower[i], b)),
 753
 754   ni_coeffs: sublist(append(upper, lower), lambda([ni], is(integerp(ni) and ni<0))),
 755   if ni_coeffs#[] then block(
 756     [use_simpsum:false, use_harmonic:false, use_integral:false,
 757      use_ratfun:false, use_gosper:false, use_zeilberger:false,
 758      use_telescoping:false, min_ni: lmin(ni_coeffs)],
 759     apply('sum, [expr, var, lo, -min_ni - 1]) + simplify_sum(apply(sum, [expr, var, -min_ni, hi])))
 760   else to_hypergeometric1(expr, var, lo, hi))$
 761
 762
 763 to_hypergeometric1(expr, var, lo, hi) := block(
 764   [quo, upper, lower, a:[], b:[], x, c, warnings:false, besselexpand:true],
 765
 766   for i:1 thru 100 while subst(var=lo, expr)=0 do lo : lo+1,
 767 /*  expr:subst(var=var+lo, expr),*/
 768
 769   quo : ratsimp(shiftQuo(factor(makefact(expr)*var!), var)),
 770   if not(?ratp(quo, var)) then return(false),
 771
 772   ss_print_message(2, "Shift quotient", quo),
 773   quolim: limit(quo/(var+1), var, inf),
 774   if freeof_limit(quolim) and abs(quolim)>1 then error("Sum is divergent!"),
 775
 776   upper : -map(rhs, solve(num(quo), var)),
 777   if not(every(lambda([u], freeof(var, u)), upper)) then return(false),
 778   for i:1 thru length(upper) do (
 779     for j:1 thru multiplicities[i] do
 780       a : cons(upper[i], a)),
 781   ss_print_message(2, "a=", a),
 782
 783   lower : -map(rhs, solve(denom(quo), var)),
 784   if not(every(lambda([u], freeof(var, u)), lower)) then return(false),
 785   for i:1 thru length(lower) do (
 786     for j:1 thru multiplicities[i] do
 787       b : cons(lower[i], b)),
 788   ss_print_message(2, "b=", b),
 789
 790   x : ratsimp(quo / apply("*", map(lambda([u], var+u), a)) *
 791                     apply("*", map(lambda([u], var+u), b))),
 792   ss_print_message(2, "x=", x),
 793
 794   c : subst(var=0, expr),
 795   ss_print_message(2, "c=", c),
 796   if c=0 then return(false),
 797
 798   ratsimp(c*hgfred(a,b,x)))$
 799
 800 /*******************
 801  *
 802  * harmonic_number and gen_harmonic_number
 803  *
 804  * harmonic_number(n) = sum(1/i, i, 1, n)
 805  * gen_harmonic_number(n,k) = sum(1/i^k, i, 1, n)
 806  *
 807  *******************/
 808
 809 define_variable(harmonic_number_expand, false, boolean)$
 810
 811 simp_harmonic_number(x__):=
 812   if x__=0 then 0
 813   else if integerp(x__) and x__<1 then error("Zero to negative power computed.")
 814   else if integerp(x__) then num_harmonic_number(1, x__)
 815   else if numberp(x__) or imagpart(x__)#0 then psi[0](x__+1) + %gamma
 816   else block(
 817     [a, b, var, k%],
 818     if harmonic_number_expand then (
 819       [a,b]:split_integer_part(x__),
 820       if harmonic_number_expand and b>0 then
 821         simpfuncall('harmonic_number,a) + apply('sum, [1/k%, k%, a+1, a+b])
 822       else simpfuncall('harmonic_number, x__))
 823     else simpfuncall('harmonic_number, x__))$
 824
 825 num_harmonic_number(l, h) :=
 826     if h<l then 0
 827     else if h=l then 1/l
 828     else if h-l<50 then sum(1/i, i, l, h)
 829     else block(
 830         [mid: floor((l+h)/2)],
 831         num_harmonic_number(l, mid) +
 832         num_harmonic_number(mid+1, h))$
 833
 834 simplifying('harmonic_number,'simp_harmonic_number)$
 835
 836 simp_gen_harmonic_number(exp__, x__):=
 837   if x__=0 then 0
 838   else if integerp(x__) and x__<1 then error("Zero to negative power computed.")
 839   else if exp__=1 then harmonic_number(x__)
 840   else if x__>=inf then zeta(exp__)
 841   else if integerp(x__) and integerp(exp__) then num_gen_harmonic_number(exp__, 1, x__)
 842   else if integerp(x__) then sum(1/i^exp__, i, 1, x__)
 843   else if (numberp(x__) and numberp(exp__)) or imagpart(x__)#0 then
 844     (-1)^(exp__+1)/(exp__-1)!*(psi[exp__-1](x__+1)-psi[exp__-1](1))
 845   else block(
 846     [a, b, var, k%],
 847     if harmonic_number_expand then (
 848       [a,b]:split_integer_part(x__),
 849       if harmonic_number_expand and b>0 then
 850         simpfuncall('gen_harmonic_number,exp__,a) + apply('sum, [1/k%^exp__, k%, a+1, a+b])
 851       else simpfuncall('gen_harmonic_number, exp__, x__))
 852     else simpfuncall('gen_harmonic_number, exp__, x__))$
 853
 854 num_gen_harmonic_number(a, l, h) :=
 855     if h<l then 0
 856     else if h=l then 1/l
 857     else if h-l<50 then sum(1/i^a, i, l, h)
 858     else block(
 859         [mid: floor((l+h)/2)],
 860         num_gen_harmonic_number(a, l, mid) +
 861         num_gen_harmonic_number(a, mid+1, h))$
 862
 863 simplifying('gen_harmonic_number,'simp_gen_harmonic_number)$
 864
 865 get_harmonic_number_args(expr) :=
 866   if atom(expr) then {}
 867   else if part(expr, 0)=harmonic_number then {[1, part(expr, 1)]}
 868   else if part(expr, 0)=gen_harmonic_number then {[part(expr, 1), part(expr, 2)]}
 869   else xreduce(union, map(get_harmonic_number_args, args(expr)))$
 870
 871 split_integer_part(expr) :=
 872   if integerp(expr) then [0,expr]
 873   else if atom(expr) then [expr,0]
 874   else if part(expr,0)="+" then block(
 875     [a:0,b:0],
 876     for arg in args(expr) do (
 877       if integerp(arg) then b:b+arg
 878       else a:a+arg),
 879     [a,b])
 880   else [expr, 0]$
 881
 882 harmonic_to_psi(expr) :=
 883   opsubst([
 884     harmonic_number=lambda([x__], psi[0](x__+1)+%gamma),
 885     gen_harmonic_number=lambda([exp__, x__],
 886       (-1)^(exp__+1)/(exp__-1)!*(psi[exp__-1](x__+1)-psi[exp__-1](1)))],
 887     expr)$
 888
 889 /*******************
 890  *
 891  * sum_by_parts
 892  *
 893  *******************/
 894
 895 sum_by_parts(fkgk, gk, k__, m__, n__) := block(
 896     [fk:fkgk/gk, j__:?gensym(), gj, oth, harmonic_number_expand:true],
 897     gj: subst(k__=j__, gk),
 898     oth: (subst(k__=k__+1, fk) - fk) * apply('sum, [gj, j__, m__, k__]),
 899     subst(k__=n__+1, fk) * apply('sum, [gj, j__, m__, n__]) - apply('sum, [oth, k__, m__, n__]))$
 900
 901
 902 /*******************
 903  *
 904  * sum_by_integral:
 905  * - currently only handles harmonic_number, but could be extended to other functions.
 906  *
 907  *******************/
 908
 909 define_variable(sum_by_integral_transforms, [], list)$
 910
 911 sum_by_integral(expr, var, lo, hi) :=
 912   /* Integral representation of harmonic_number */
 913   if catch(
 914     scanmap(
 915       lambda([u], if not atom(u) and member(part(u,0), [harmonic_number]) then throw(true)),
 916       expr,
 917       bottomup))=true then
 918   block(
 919     [expr1, x_:?gensym()],
 920
 921     supcontext (concat('sum_by_integral, simplify_sum_depth)),
 922     assume(x_>0, x_<1, var>=lo, var<=hi),
 923
 924     expr1 : opsubst(harmonic_number=lambda([u], (1-x_^u)/(1-x_)), expr),
 925     if errcatch(
 926       expr1 : simplify_sum(split_sum(expand(apply('sum, [expr1, var, lo, hi]))))
 927       ) = [] then expr1 : false,
 928
 929     killcontext(concat('sum_by_integral, simplify_sum_depth)),
 930
 931     for tr in sum_by_integral_transforms do (
 932       expr1 : apply(tr, [expr1])),
 933
 934     if expr1 # false then expr1: integrate(expr1, x_, 0, 1)
 935     else expr1 : apply('sum, [expr, var, lo, hi]),
 936
 937     if freeof_integrate(expr1) and freeof_limit(expr1) then expr1
 938     else apply('sum, [expr, var, lo, hi]))
 939
 940   else apply('sum, [expr, var, lo, hi])$
 941
 942 /*****
 943  * recognize sums of the form sum(f(n) - f(n+1), n, a, b) where f(n)=f1(n)/f2(n) for
 944  * simple functions f1 and f2.
 945  ****/
 946
 947 find_f2(expr) :=
 948   if mapatom(expr) then (if subvarp(expr) then [expr] else [])
 949   else if member(part(expr, 0), ["+", "-", "/", "*"]) then
 950     apply(append, map(find_f2, args(expr)))
 951   else [expr]$
 952
 953 find_f1(expr, f2, var) := block(
 954   [f21: ratsimp(subst(var=var+1, f2)), f1, f11, algebraic:true],
 955   expr: ratsimp(expr*f2*f21),
 956   f1: coeff(ratsimp(expr), f21),
 957   f11: subst(var=var+1, f1),
 958   expr: radcan(expr - f21*f1 + f2*f11),
 959   if expr=0 then f1
 960   else false)$
 961
 962 find_quotient(expr, var) := block(
 963   [f2_list, f1, quotient: false],
 964   f2_list: find_f2(expr),
 965   for f2 in f2_list while quotient=false do (
 966     f1: find_f1(expr, f2, var),
 967     if f1#false then
 968       quotient: f1/f2),
 969   quotient)$
 970
 971 telescoping_sum(expr, var, lo, hi) := block(
 972   [quotient: find_quotient(expr, var)],
 973   if quotient#false then
 974     (if lo=minf then limit(quotient, var, lo) else subst(var=lo, quotient)) -
 975     (if hi=inf then limit(quotient, var, hi) else subst(var=hi+1, quotient))
 976   else false)$
 977
 978 eval_when(batch,
 979           ttyoff : false,
 980           nolabels : false)$