Division by zero error in certain sums

[Kurt Pagani]

I get the following error in certain sums when starting at k=0:

(1) -> sum(4^k/binomial(2*k, k), k=0..n)

>> Error detected within library code:
catdef: division by zero


eval(4^k/binomial(2*k, k),k=0) --> 1

eval(4^k/binomial(2*k, k),k=-1) --> above error message


I suspect the error could be triggered in innersum(x, k) in
combfunc.spad (line 2785/eval k::F-1 ??) but am not sure.


-----------------------------------
--- the result is correct otherwise

(2) -> S:=1+sum(4^k/binomial(2*k, k), k=1..n)

2 n n
( ) + (2 n + 2)4
n
(2) -------------------
2 n
3 ( )
n
Type:
Expression(Integer)


(3) -> Sexpected:=4/3*4^n*(2*n+1)/binomial(2*n+2,n+1)+1/3

2 n + 2 n
( ) + (8 n + 4)4
n + 1
(67) -----------------------
2 n + 2
3 ( )
n + 1
Type:
Expression(Integer)


(4) -> test(normalize(S-Sexpected)=0)

(4) true
Type:
Boolean

[Waldek Hebisch]

On Wed, Jun 03, 2026 at 03:20:57PM +0200, Kurt Pagani wrote:
> I get the following error in certain sums when starting at k=0:
>
> (1) -> sum(4^k/binomial(2*k, k), k=0..n)
>
> >> Error detected within library code:
> catdef: division by zero
>
>
> eval(4^k/binomial(2*k, k),k=0) --> 1
>
> eval(4^k/binomial(2*k, k),k=-1) --> above error message
>
>
> I suspect the error could be triggered in innersum(x, k) in combfunc.spad
> (line 2785/eval k::F-1 ??) but am not sure.

There are similar examples like:

sum(n*factorial(n), n=0..m)


AFAICS the error is in line 2773:

eval(u@F, k, 1 + high(segment(s))) - eval(u@F, k, low(segment(s)))

Previous evals are symbolic, so should not signal error. The
above substitutes numeric value which is not always possible.
We need smarter function instead of 'eval' so that it can correctly
evaluate 0/0 things. Taking limit may help, but limits of
binomial expressions can be tricky, that is result may depend on
how exactly one approaches to limit. So it needs some analysis
before we implement something smarter.

--
Waldek Hebisch

[Kurt Pagani]

You are right, the err stems from line 2773.


I have extracted the relevant steps and it's the last one which gives
the error, obviously because of the subst k --> k-1. One solution might
be to shift the boundaries to k=1..n+1 and adjust x correspondingly?
This, however, would mean some preprocessing ...

R ==> INT
F ==> EXPR R
K ==> Kernel F
IK ==> IndexedExponents K
SMP ==> SparseMultivariatePolynomial(R, K)
GSM ==> GosperSummationMethod(IK,K,R,SMP ,F)

newk():K == kernel(new()$Symbol)

x:F:=k*factorial(k)
--x:=4^k/binomial(2*k, k)
s:=k=0..n
k := kernel(variable s)@K
x1:=eval(x,k,k::F-1)
f:=normalize(x/x1)
u0 := GospersMethod(f, k, newk)$GSM
u:=x1 * eval(u0::F, k, k::F-1 ) ---> k(k - 1)!
eval(u::F, k, 1 + high(segment(s))) - eval(u::F, k, low(segment(s)))

[Kurt Pagani]

addendum

Shifting doesn't work either because of the preceding evals (k->k-1), so
one might sum from 1..n and add eval(x,k,0). If a<0? Splitting and
adding x(0)?


a:=low segment s
b:=high segment s

if a=0 then
eval(u::F, k, 1 + b) - eval(u::F, k, 1) + eval(x,k,0)
else
if a>0 then eval(u::F, k, 1 + b) - eval(u::F, k, a)
-- else?

[Qian Yun]

Off topic a little.

I was asking LLM to review combfunc.spad, it found the
following issues:

1. hahn_p takes 5 parameters, but its documentation
lists 6.

2. For Wilson Polynomial, it is 1 for degree 0 instead of
degree 1. Function "e_wilsonW" should be fixed.

- Qian

[Waldek Hebisch]

On Wed, Jun 03, 2026 at 11:06:07PM +0200, Kurt Pagani wrote:
> addendum
>
> Shifting doesn't work either because of the preceding evals (k->k-1), so one
> might sum from 1..n and add eval(x,k,0). If a<0? Splitting and adding x(0)?
>
>
> a:=low segment s
> b:=high segment s
>
> if a=0 then
> eval(u::F, k, 1 + b) - eval(u::F, k, 1) + eval(x,k,0)
> else
> if a>0 then eval(u::F, k, 1 + b) - eval(u::F, k, a)
> -- else?

AFAICS there are two possible reasons for division by 0. One is
due to result from Gosper method. Namely in your case Gosper
method returns:

2 k + 2
-------
3

which after substitution is turned into 2*k/3. The result from
innersum is then 2*k/3*x1 where x1 is the input function shifted
by -1. Evaluating that for k = 0 gives 0/0 expression which leads
to division by 0. In principle result from Gosper method may have
zero at arbitrary point, so the correct method to discover this
kind of trouble is to look at numerator of result from Gosper
method, evaluate it a lower and upper bound and check for zeros.

But there is another possible trouble: evaluating input function
at upper or lower bound may also lead to division by 0. Input
may have pole at bound, in such case division by 0 is unavoidable.
But another case is when already input contains trouble similar to
the first case. That is rational part of input is supposed to
cancel with factorial or binomial part. AFAICS we need better
version of 'eval' to transform input (cancel zeros) first and
only evaluate after transformation.

--
Waldek Hebisch

[Waldek Hebisch]

On Thu, Jun 04, 2026 at 08:13:19AM +0800, Qian Yun wrote:
> Off topic a little.
>
> I was asking LLM to review combfunc.spad, it found the
> following issues:
>
> 1. hahn_p takes 5 parameters, but its documentation
> lists 6.

Documentation is right.

> 2. For Wilson Polynomial, it is 1 for degree 0 instead of
> degree 1. Function "e_wilsonW" should be fixed.

Commited now.

BTW: We should have more functionality, in particular reccurences,
but it takes time to work out details.

--
Waldek Hebisch

[Kurt Pagani]

On 04/06/2026 02:30, Waldek Hebisch wrote:
...

>
> But there is another possible trouble: evaluating input function
> at upper or lower bound may also lead to division by 0. Input
> may have pole at bound, in such case division by 0 is unavoidable.
> But another case is when already input contains trouble similar to
> the first case. That is rational part of input is supposed to
> cancel with factorial or binomial part. AFAICS we need better
> version of 'eval' to transform input (cancel zeros) first and
> only evaluate after transformation.
>

Indeed,as you say, it's even worse. Your example gives

>> Error detected within library code:

factorial not defined for negative integers

which clearly shows that a summand is tried to be evaluated outside of
the domain of definition.

In the book "A=B" on p. 77, where Gosper's algorithm is summarized the
following (IMO) cryptic passage is a mystery to me.

---
Once we have zn , the sum that we are looking for is sn = zn − z0 . The
lower summation bound need not be 0; for example, it may happen that the
summand in (5.1.1) is undefined for certain integer values of k, and
then we will want to start the summation at some large enough value to
skip over all the singularities. If the lower summation bound is k0,
then everything goes through as before, except that now sn = zn − zk0 .
---

I guess we may assume that in finite summation the summands are well
defined over the segment? Well, I leave it to you to rethink the whole
matter :)