csum:= proc(ep, np::name)
# Test convergence of series
# usage:
#  csum(ep, np)
# ep an expression involving name np
# determine convergence of sum(ep, np=a..infinity) 
# (for sufficiently large a that ep is defined for integers np >= a)
# return true (converges), false (diverges), FAIL (undecided)
# with infolevel[csum] >= 1, tell whether convergence is absolute
# with infolevel[csum] >= 2, tell results of each test
option `Maple Advisor Database 1.01 for Maple 6`,
  `Copyright (c) 2000 by Robert B. Israel.  All rights reserved`;
local e,n, ipars, k, limzero, eg1, eg2, ea, qpos, b;
global `csum/a`, `property/object`;
# assume n a large enough integer
  ipars:= map(x -> x^2, indets(ep, integer));
  `csum/a`:= max(10000, op(ipars));
  `property/object`[n]:= AndProp(RealRange(`csum/a`,infinity),integer);
  e:= eval(ep, np=n);
  userinfo(1, csum, `Checking sum of`,e,`over`,n);
# sum of 0 converges absolutely; 0 would cause trouble in some tests
  if Testzero(e) then 
    userinfo(2,csum, `Terms are identically 0`);
    userinfo(1,csum, `Converges absolutely`);
    return true 
  fi;
  if hastype(e, negative^linear(n)) then
# split up positive and negative
   `property/object`[k]:= AndProp(RealRange(`csum/a`,infinity),integer);   
    eg1:= simplify(eval(e,n=2*k));
    eg2:= simplify(eval(e,n=2*k+1));
    if Testzero(eg1) or Testzero(eg2) then
      userinfo(3,csum,`Removing terms of 0`);
      return csum(eval(eg1+eg2),k);
    fi;
    userinfo(2,csum,`Alternate terms`,eg1,`and`,eg2);
    return `csum/alternate`(eg1,eg2,k);
  fi;
# check limit of terms is 0
  limzero:= `csum/limitzero`(e,n);
  if limzero = false then 
    userinfo(1,csum,`Series diverges`);
    return false;
  fi;
# check if terms are constant sign
  if is(e <= 0) then 
# if negative, change sign
    ea:= -e;
    qpos:= true;
    userinfo(3,csum,`Changing sign to make terms >= 0`);
  elif is(e >= 0) then  
    ea:= e;
    qpos:= true;
    userinfo(3,csum,`Terms >= 0`);
  else
    ea:= `csum/abs`(e);
    qpos:= evalb(ea=e);
    if qpos then
      userinfo(3,csum,`Terms >= 0`);
    fi;
  fi;
# Test absolute convergence
  b:= `csum/abscon`(ea,n);
  if b = true then
    userinfo(1,csum,`Series`,e,`converges absolutely`);
    return true;
  elif b = false then
    if qpos then 
      userinfo(1,csum,`Series diverges`);
      return false;
    else
      userinfo(1,csum,`Series does not converge absolutely`);
    fi
  elif qpos then
    userinfo(1,csum,`Could not decide convergence`)
  else
    userinfo(1,csum,`Could not determine absolute convergence`);
  fi; 
# Test conditional convergence
  `csum/condit`(e,n);
end;    

`csum/alternate`:= proc(e1, e2, n)
# test convergence of series where terms are alternately e1 and e2
# usage: 
#   `csum/alternate`(e1, e2, n)
#   e1 and e2: expressions in n
#   n: a name, assumed to be a sufficiently large positive integer
#   series is e1(n=a) + e2(n=a) + e1(n=a+1) + e2(n=a+1) + ...
#   return true (series converges), false (diverges), FAIL (undecided)
option `Maple Advisor Database 1.01 for Maple 6`,
  `Copyright (c) 2000 by Robert B. Israel.  All rights reserved`;
local es, ea, limzero, i, qpos, qch, b;
  es:= [e1, e2];
# check limits of terms as n -> infinity
# if either is nonzero, diverges
  limzero:= [FAIL,FAIL];
  for i from 1 to 2 do
    limzero[i]:= `csum/limitzero`(es[i],n);
    if limzero[i] = false then
      userinfo(1,csum,`Series diverges`);
      return false;
    fi;
  od;
# check if terms are constant sign
  ea:= [0,0]; qpos:= [FAIL,FAIL]; qch:= [false,false];
  for i from 1 to 2 do
    if is(es[i] <= 0) then 
# if negative, change sign
      ea[i]:= -es[i];
      qpos[i]:= true;
      qch[i]:= true;
      userinfo(3,csum,`Changing sign in`,es[i],`to make terms >= 0`);
    elif is(es[i] >= 0) then  
      ea[i]:= es[i];
      qpos[i]:= true;
      userinfo(3,csum,`Terms`,es[i],`>= 0`);
    else
      ea[i]:= `csum/abs`(es[i]);
      qpos[i]:= evalb(ea[i]=es[i]);
      if qpos[i] then
        userinfo(3,csum,`Terms`,es[i],`>= 0`);
      fi;
    fi;
  od;
# Test absolute convergence
# series converges absolutely iff both converge absolutely
# if all terms have same sign, series diverges if either does not conv.abs.
  b:= map(`csum/abscon`,ea,n);
  if b = [true,true] then
    userinfo(1,csum,`Series converges absolutely`);
    return true;
  fi;
  if member(false,b) then
    if qpos = [true,true] and qch[1]=qch[2] then 
      userinfo(1,csum,`Series diverges`);
      return false;
    else
      userinfo(1,csum,`Series does not converge absolutely`);
    fi
  else
    userinfo(1,csum,`Could not determine absolute convergence`);
  fi;
# Test conditional convergence
# if both converge, series converges
# if one converges and other diverges, series diverges
  for i from 1 to 2 do if (b[i] = FAIL) and not qpos[i] then
    b[i]:= `csum/condit`(es[i],n);
  fi od;
  if b = [true,true] then
    userinfo(1,csum,`Series`,e,`converges`);
    return true;
  elif b = [true,false] or b = [false,true] then
    userinfo(1,csum,`Series diverges`);
    return false;
  fi;
# test convergence of e1+e2 series if at least one limit known to be 0
  if member(true,limzero) then
     userinfo(2,csum,`Regrouping`); 
     csum(e1+e2,n);
  else
     userinfo(1,csum,`Could not determine convergence`);
     FAIL
  fi
end; 

`csum/limitzero`:= proc(e,n)
# test whether limit(e,n=infinity) = 0
# returns true, false or FAIL
# note that a nonconstant result is assumed nonzero
# even if it's 0 for some values of parameters
# e: an expression in n
# n: name, assumed positive integer
option `Maple Advisor Database 1.01 for Maple 6`,
`Copyright (c) 2000 by Robert B. Israel.  All rights reserved`;
local L, ar, br;
  userinfo(2, csum, `Checking limit of terms`);
  L := traperror(limit(e, n = infinity));
  if L <> lasterror then
    if Testzero(L) then 
      userinfo(2, csum, `OK, limit is 0`);
      return true
    fi;
    if (type(L, algebraic) and not has(L,'limit')) or
       (typematch(L, (ar::constant) .. (br::constant)) and
        is(0 < ar or br < 0)) then
          userinfo(1, csum, `Diverges, limit of terms is `, L);
          return false;
    fi;
  fi;
  L := traperror(limit(`csum/abs`(e), n = infinity));
  if L <> lasterror then
    if Testzero(L) then 
      userinfo(2, csum, `OK, limit is 0`);
      return true
    fi;
    if (type(L, algebraic) and not has(L,'limit')) or
       (typematch(L, (ar::constant) .. (br::constant)) and
        is(0 < ar or br < 0)) then
          userinfo(1, csum, 
             `Diverges, limit of absolute values of terms is `, L);
          return false;
    fi;
  fi;
  userinfo(1,csum,`limit of absolute values of terms is`, L);
  FAIL
end;

`csum/abs`:= proc(x)
local r, v, k, u;
option `Maple Advisor Database 1.01 for Maple 6`,
`Copyright (c) 2000 by Robert B. Israel.  All rights reserved`;
    if is(0 <= x) then return x elif is(x <= 0) then return -x end if;
    r := abs(x);
    if has(r, abs) then
        if typematch(r, abs(v::anything)) then
            if type(v, `*`) then return map(`csum/abs`, v)
            elif typematch(v, (u::algebraic)^(k::algebraic)) and
            is(k, real) then return `csum/abs`(u)^k
            end if
        elif type(r, `*`) then return map(`csum/abs`, r)
        end if
    end if;
    r
end proc;

`csum/abscon`:= proc(e, n)
# check absolute convergence of sum
# usage:
#  `csum/abs`(e, n)
# e an expression involving name n
# where e >= 0 for sufficiently large integer n
option `Maple Advisor Database 1.01 for Maple 6`,
`Copyright (c) 2000 by Robert B. Israel.  All rights reserved`;
local b, test;
for test in [`csum/ratiotest`, `csum/loglogtest`, `csum/indef`,
`csum/split`] do
  b:= test(e,n);
  if b <> FAIL then return b fi
od;
FAIL
end;

`csum/condit`:= proc(e,n)
# check for convergence, terms may not be nonnegative
# this is not called if terms are known to be nonnegative
# e: an expression depending on name n
option `Maple Advisor Database 1.01 for Maple 6`,
`Copyright (c) 2000 by Robert B. Israel.  All rights reserved`;
local test, b;
  for test in [`csum/osctest`,`csum/indef`, `csum/split`] do
    b := test(e, n);
    if b = true then 
      userinfo(1, csum, `Series converges`);
      return true;
    elif b = false then
      userinfo(1, csum, `Series diverges`);
      return false;
    fi;
  od;
  userinfo(1, csum, `Could not determine convergence`);
  FAIL
end proc;

`csum/ratiotest`:= proc(e,n)
# ratio test
# e an expression in variable n, which should be nonnegative for large n
local L, ep, b;
option `Maple Advisor Database 1.01 for Maple 6`,
`Copyright (c) 2000 by Robert B. Israel.  All rights reserved`;
  userinfo(2, csum, `Trying ratio test`);

  L := traperror(
        limit(subs(n = n + 1, e)/e, n = infinity));
  if L = lasterror or has(L, limit) then
     if typematch(e,abs(ep::algebraic)) then 
       L:= traperror(limit(subs(n=n+1,ep)/ep, n=infinity));
     else ep:= e
     fi;
     if L = lasterror or has(L, limit) then
       b:= simplify(subs(n=n+1,ep)/ep);
       L := traperror(limit(b, n = infinity))
     fi
  fi;
  if type(L, range) then
    b := traperror(is(-1 < op(1, L)) and is(op(2, L) < 1));
    if b = true then
      userinfo(2, csum, `Ratio test succeeds with ratio limit `, L);
      return true
    end if;
    b := traperror(is(1 < op(1, L)) or is(op(2, L) < -1));
    if b = true then
      userinfo(2, csum, `Ratio test fails with ratio limit `, L);
      return false
    end if
  else # not a range
    if L = lasterror or L = 1 or L = -1 then
      userinfo(2, csum,
        `Ratio test inconclusive with ratio limit`, L);
      return FAIL
    end if;
    b := traperror(is(L < 1) and is(-1 < L));
    if b = true then
      userinfo(2, csum, `Ratio test succeeds with ratio limit `, L);
      return true
    end if;
    b := traperror(is(1 < L) or is(L < -1));
    if b = true then
      userinfo(2, csum,
                `Ratio test shows divergence with ratio limit`, L);
      return false
    end if
  end if;
  userinfo(2, csum, `Ratio test inconclusive with ratio limit`, L);
  FAIL
end proc;

`csum/loglogtest`:= proc(e,n)
# loglog test
# e an expression in variable n, which should be nonnegative for large n
local L, b, lev, j, x, a;
option `Maple Advisor Database 1.01 for Maple 6`,
`Copyright (c) 2000 by Robert B. Israel.  All rights reserved`;
  userinfo(2, csum, `Trying loglog test`);

  for lev to 10 do
        x:= ln(e*n*mul((ln@@j)(n), j = 1 .. lev - 1))/
            (ln@@(lev + 1))(n);
        L := traperror(limit(x, n = infinity));
        if L = lasterror or has(L,limit) then
          L:= traperror(limit(asympt(x,n), n=infinity));
        fi;
        if L = lasterror or has(L,limit) then
           if is(x < -1) then 
             L:= -2;
             userinfo(4,csum, ln(a[n]*n*mul((ln@@j)(n), j = 1 .. lev - 1))/
            (ln@@(lev + 1))(n),`< -1 as `,n -> infinity);
           elif is(x > -1) then
             L:= 0;
             userinfo(4,csum, ln(a[n]*n*mul((ln@@j)(n), j = 1 .. lev - 1))/
            (ln@@(lev + 1))(n),`> -1 as `,n -> infinity);       
           fi
        else userinfo(4, csum, `limit of `,
            ln(a[n]*n*mul((ln@@j)(n), j = 1 .. lev - 1))/
            (ln@@(lev + 1))(n), is, L);
        fi;
        if L <> -1 then break end if
  end do;
  if type(L, range) then
        b := traperror(is(op(2, L) < -1));
        if b = true then
            userinfo(2, csum, `Loglog test succeeds`); return true
        end if;
        b := traperror(is(-1 < op(1, L)));
        if b = true then
            userinfo(2, csum, `Loglog test shows divergence`);
            return false
        end if
  elif L = -1 then
        userinfo(2, csum, `Loglog test fails`); 
        return FAIL
  elif L <> lasterror then
        b := traperror(is(L < -1));
        if b = true then
            userinfo(2, csum, `Loglog test succeeds`); return true
        end if;
        b := traperror(is(-1 < L));
        if b = true then
            userinfo(2, csum, `Loglog test shows divergence`);
            return false
        end if
  end if;
  userinfo(2, csum, `Loglog test inconclusive`);
  FAIL
end proc;

`csum/indef`:= proc(e,n)
# indefinite summation test: if you can find an indefinite
# sum, see if it has a limit
local F, L;
option `Maple Advisor Database 1.01 for Maple 6`,
`Copyright (c) 2000 by Robert B. Israel.  All rights reserved`;
    userinfo(2, csum, `Trying indefinite summation`);
    F := sum(e, n);
    if has(F, 'sum') then
        userinfo(2, csum, `Indefinite summation failed`); return FAIL
    end if;
    L := limit(F, n = infinity);
    if has(L, 'limit') or type(L, range) then
        userinfo(2, csum, `Could not determine limit: failed`);
        return FAIL
    end if;
    if hastype(L, infinity) or type(L, undefined) then
        userinfo(2, csum, `Limit of indefinite sum is`, L);
        return false
    end if;
    userinfo(2, csum, `Indefinite summation succeeded`);
    true
end proc;

`csum/split`:= proc(e,n)
# split up the summand into pieces and see if the pieces 
# converge
local i, b, bd, bf, bs, qpos;
option `Maple Advisor Database 1.01 for Maple 6`,
`Copyright (c) 2000 by Robert B. Israel.  All rights reserved`;
    if not type(e, `+`) then return FAIL end if;
    userinfo(2, csum, `Splitting up terms of`, e);
    b := convert(e, list);
    bs:= convert(map(signum,b),set);
    qpos:= evalb((bs = {1}) or (bs = {-1}));
# all positive or all negative: if any diverges the sum does
    b := map(csum, b, n);
    bd := 0;
    bf := 0;
    for i to nops(b) do
        if b[i] = false then 
          if qpos then 
            userinfo(2, csum, `Term has divergent sum`); 
            return false
          fi;
          bd := bd + 1
        elif b[i] = FAIL then bf := bf + 1
        end if
    end do;
    if bd = 1 and bf = 0 then
        userinfo(2, csum, `Sum of one term diverges`); return false
    end if;
    if bd = 0 and bf = 0 then
        userinfo(2, csum, `Sums of all terms converge`); true
    else userinfo(2, csum, `Splitting up failed`); FAIL
    end if
end proc;


`csum/osctest`:= proc(e, n, a)
local ep, u, v, f, S, L, q;
option `Maple Advisor Database 1.01 for Maple 6`,
`Copyright (c) 2000 by Robert B. Israel.  All rights reserved`;
# oscillating series test
# sum(u[n]*v[n],n) converges if
# sum(u[n],n) is bounded and v[n] -> 0 monotonically
    ep := normal(e);
    if type(ep, `*`^algebraic) then ep := normal(expand(e)) end if;
    if not type(ep, `*`) then return FAIL end if;
    userinfo(2, csum, `Trying oscillating series test on `, ep);
    u := 1;
    v := 1;
    for f in ep do
        q:= indets(f, (-1)^algebraic);
        if q <> {} then u:= u*f; next fi;
        q:= indets(f, specfunc(anything, exp));
        if (q <> {}) and 
            traperror(is(I*op(op(1,q)),real)) = true
          then u:= u*f; next 
        fi;
        q:= indets(f, specfunc(anything,{cos,sin}));
        if (q <> {}) and
            traperror(is(op(op(1,q)),real)) = true
          then u:= u*f; next
        fi;
        v := f*v;
    end do;
    userinfo(4, csum, Splitting, e, into, u, `and`, v);
    S := sum(u, n);
    if has(S, 'sum') then
        userinfo(4, csum, `Could not explicitly sum `, u); 
        userinfo(2, csum, `Oscillating series test fails`);
        return FAIL
    end if;
    L := traperror(limit(S, n = infinity));
    if type(L, infinity) then
        userinfo(4, csum, `Indefinite sum`, S, of, u,
            `has infinite limit`);
        userinfo(2, csum, `Oscillating series test fails`);
        return FAIL
    end if;
    if has(L, 'limit') or L = lasterror then
        userinfo(4, csum, `Could not find limit of indefinite sum`, S)
            ;
        userinfo(2, csum, `Oscillating series test fails`);
        return FAIL
    end if;
    userinfo(4, csum, `Indefinite sum`, S, `is bounded`);
    L := traperror(limit(v, n = infinity));
    if L <> 0 then
        userinfo(4, csum, `Could not show convergence to 0 of`, v);
        userinfo(2, csum, `Oscillating series test fails`);
        return FAIL
    end if;
    L := traperror(asympt(v, n));
    if L = lasterror then
        userinfo(4, csum, `Could not check monotonicity of`, v);
        userinfo(2, csum, `Oscillating series test fails`);
        return FAIL
    end if;
    L := traperror(limit(L/`csum/abs`(L), n = infinity));
    if L = 1 or L = -1 then
        userinfo(2, csum, `Oscillating series test succeeds`);
        return true
    end if;
    userinfo(4, csum, `Could not show monotonicity of`, v);
    FAIL
end proc;
