The Joy of CAS

[Peter Luschny]

Simple question: What is

sum_{i=0..n} binomial(n-i-1, n-i) ?

I compiled some answers:

** Wolfram Alpha : 0

** Maple : returns unevaluated

** Sage 6.2 combined with mpmath really excels:
Depending on the order of the execution it gives
either a TypeError or the answer 1.0.

So, what is sum_{i=0..n} binomial(n-i-1, n-i) and why
can the writers of CASs not agree on the answer?

Here you can see the outputs:
http://luschny.de/temp/JoyOfCAS.png

Peter

[Axel Vogt]

On 25.07.2014 12:11, Peter Luschny wrote:
> sum_{i=0..n} binomial(n-i-1, n-i)

I think the reason might be

binomial(n-i-1,n-i);
subs(n-i=m, %); convert(%, GAMMA); subs(m=n-i, %);

GAMMA(n - i)
-------------------------
GAMMA(1 + n - i) GAMMA(0)

and ignoring 1/GAMMA(0) = 0 but computing GAMMA(0).

[Peter Luschny]

What about this approach:

Replace the last two occurrences of n by k in
sum(binomial(n-i-1,n-i), i=0..n);

sum(binomial(n-i-1,k-i), i=0..k);

Now this sum is (for 0 <= k < n)

-n*binomial(n-1, k)/(k-n)

a := (n,k) -> n*binomial(n-1, k)/(n-k);

for n from 0 to 5 do seq(a(n,k),k=0..n-1) od;

1
1, 2
1, 3, 3
1, 4, 6, 4
1, 5, 10, 10, 5

So how do you want this (quite famous) triangle to be
extended on the right hand side? OK, I agree.

The most simple trick will do: look at the limit.

for n from 0 to 5 do seq(limit(a(n,m),m=k),k=0..n) od;

0
1, 1
1, 2, 1
1, 3, 3, 1
1, 4, 6, 4, 1
1, 5, 10, 10, 5, 1

Fine. Only the case (0,0) is razor-ed here by Maple.
But in this case we have

sum(binomial(0-i-1,0-i), i=0..0) = 1

So there is no problem here and the general answer is:

sum(binomial(n-i-1,n-i), i=0..n) = 1

by extension by limit of sum(binomial(n-i-1,k-i),i=0..k).

Still I am worried that Mathematica insists on

sum(binomial(n-i-1,n-i), i=0..n) = 1.

P.S.:
Caveat: When I say 'Mathematica' I really mean 'Alpha'.
I just received a message of Hans Haverman who
alerts me that in some of these equations
"Alpha's result disagrees with Mathematica's."

[Peter Luschny]

Correction:
> Still I am worried that Mathematica [Alpha] insists on
sum(binomial(n-i-1,n-i), i=0..n) = 0.

[clicl...@freenet.de]

Peter Luschny schrieb:

>
> Correction:
> > Still I am worried that Mathematica [Alpha] insists on
> sum(binomial(n-i-1,n-i), i=0..n) = 0.
>

So does Derive 6.10 in fact:

SUM(COMB(n-i-1,n-i),i,0,n)

" COMB(m,n) -> COMB(m,m-n) "

SUM(COMB(n-i-1,-1),i,0,n)

" If n is a negative integer and m is not, COMB(m,n) -> 0 "

SUM(0,i,0,n)

Wait a moment, wait a moment. This rule is misapplied, since n-i-1 is
negative for i=n. And COMB(-1,-1) evaluates to 1. Looks like we have
caught a real bug in Derive 6.10 here!

" SUM(F(x),x,a,b) -> SUBST_DIFF(SUM(F(x),x),x,a,b+1) "

SUBST_DIFF(SUM(0,i),i,0,n+1)

" SUM(a,x) -> a*x "

SUBST_DIFF(0,i,0,n+1)

" SUBST_DIFF(F(x),x,a,b) -> F(b)-F(a) "

0

This is consequently wrong, the correct result being 1. Can Alpha
display steps for this problem? Is the correspondence accidental?

Martin.

[Peter Luschny]

This time I am using Sage as the engine.

Step 1:

m, n = var('m,n')
for k in (0..4):
expand(add(m^(k-i)*binomial(n-i-1,k-i) for i in (0..k)))

[0] 1
[1] m*n - m + 1
[2] 1/2*m^2*n^2 - 3/2*m^2*n + m^2 + m*n - 2*m + 1
[3] 1/6*m^3*n^3 - m^3*n^2 + 11/6*m^3*n + 1/2*m^2*n^2
- m^3 - 5/2*m^2*n + 3*m^2 + m*n - 3*m + 1

Step 2:

m, n = var('m,n')
for k in (0..4):
print [k], expand(add(m^i*rising_factorial(n-k,i)/factorial(i) \
for i in (0..k)))

[0] 1
[1] m*n - m + 1
[2] 1/2*m^2*n^2 - 3/2*m^2*n + m^2 + m*n - 2*m + 1
[3] 1/6*m^3*n^3 - m^3*n^2 + 11/6*m^3*n + 1/2*m^2*n^2
- m^3 - 5/2*m^2*n + 3*m^2 + m*n - 3*m + 1

Step 3:

We conclude, for all m,n,k:

sum_{i=0..k} m^(k-i)*binomial(n-i-1,k-i) =
sum_{i=0..k} m^i*rising_factorial(n-k,i)/i!

We are interested in the case n=k.

sum_{i=0..n} m^(n-i)*binomial(n-i-1,n-i) =
sum_{i=0..n} m^i*rising_factorial(0,i)/i!

rising_factorial(0,i) = 1 if i = 0 otherwise 0.
Therefore the considered sum has to be 1.

====

So what does Sage return?
File "expression.pyx", line 1178,
(sage/symbolic/expression.cpp:8057)
TypeError: unable to simplify to float approximation

Also Harald Schilly told me yesterday on Google+:
<cite>
"so well, maxima agrees with mathematica
var("i n")
sum(binomial(n-i-1, n-i), i, 0, n) => 0"
</cite>

I did not check this.

CC> Can Alpha display steps for this problem?

Perhaps it can, but it does not offer this option,
at least not for the not-paying user.

Peter

[Axel Vogt]

On 25.07.2014 12:41, Axel Vogt wrote:
> On 25.07.2014 12:11, Peter Luschny wrote:
>> sum_{i=0..n} binomial(n-i-1, n-i)

...
More formally using Maple: student[changevar] is oldish, but usefull
and changevar(i=n-j,%, j) gives Sum(binomial(j-1,j),j = 0 .. n). Now

eval(binomial(j-1,j), j=0) + Sum(binomial(j-1,j),j = 1 .. n);

gives 1 + sum, it refuses for binomial(j-1,j) because of "1/GAMMA(0)".

So enforce that manually:

binomial(j-1,k):
%=convert(%, GAMMA);
Limit(rhs(%), k=j): %= value(%);

GAMMA(j)
binomial(j - 1, k) = -------------------------
GAMMA(1 + k) GAMMA(j - k)


GAMMA(j)
lim ------------------------- = 0
k -> j GAMMA(1 + k) GAMMA(j - k)

There is some "cheating" in the last result (I think)

[Peter Luschny]

This problem has a continuation for hypergeometric functions.
Consider

(*)
Table[Table[Sum[2^(k-i)*Binomial[n-i-1,k-i],{i,0,k}],{k,0,n}],{n,0,7}]

Who does not have Mathematica can enter this formula in Alpha.
Alpha provides the answer:

{{0},
{1, 0},
{1, 3, 0},
{1, 5, 7, 0},
{1, 7, 17, 15, 0},
{1, 9, 31, 49, 31, 0},
{1, 11, 49, 111, 129, 63, 0},
{1, 13, 71, 209, 351, 321, 127, 0}}

As Hans Havermann told me Mathematica 10 yields final terms 1,
not 0. So Wolfram's left hand does not know what Wolfram's
right hand does.

But now my question is: what is the right answer for
this expression:

(**)
Table[Table[2^k*Binomial[n-1,k]*Hypergeometric2F1[1,-k,-n+1,1/2],{k,0,n}],{n,0,7}]

Consider only the final terms in the rows.

(a) 0 (as in the Alpha table above)
(b) 1 (as in the Mathematica table above)
(c) Indeterminate

If you think (c) is the right answer then please explain
why you think the answer has to be different from (*).

Peter

[Peter Luschny]

Note that I considered above the case k=n.

But as it turns out this is by no means the only case
which reveals differences.

Consider now the case k = n-1: Alpha gives for

Table[2^(n-1)*Binomial[n-1,n-1]*Hypergeometric2F1[1,-n+1,-n+1,1/2],{n,0,9}]

| 1 | 1 | 3 | 7 | 15 | 31 | 63 | 127 | 255 | 511
(this time Mathematica gives the same results), but Maple

seq(evalf(2^(n-1)*binomial(n-1,n-1)*hypergeom([1,-n+1],[-n+1],1/2)), n=0..7);
1, 2, 4, 8, 16, 32, ...

Peter

[Nasser M. Abbasi]

I think There is different definitions working here. Looking
at the function as

hypergeom([r1,r2],[d],z)

This is the same as Mathematica 2F1 since we have 2 r's and one d.
The "2" in 2F1 is the number of upper parameters and the "1" is the
number of lower parameters.

Maple says in the help, when d value is non-positive, then
hypergeom() is undefined if z is not zero (the case here), _unless_
there is a _negative_ r present of smaller absolute value than d.

Notice, it says smaller absolute value than d.

http://www.maplesoft.com/support/help/Maple/view.aspx?path=hypergeom
http://reference.wolfram.com/language/ref/Hypergeometric2F1.html

In your example, when n=0, then d=1, ok, still not negative,
when n=1 then d=0, the rule should apply, at least definitely it
applied from n=2 all the way to n=7, since d is negative now.

Lets apply it:

Looking at the r1 and r2. r1 is always positive (it is fixed
at +1), so looking at r2. It is negative for n>=1, but its
absolute values is _not_ smaller than the current d value.

Actually, |r2|=|d| for all n, but the rule applies only when
d is non-positive (should this have been negative?)

So the rule kicks in when n>=1.

So according to Maple, 2F1 is not defined for n>=1 ?

Maple returns 2.0 for all these values. This
is also a bug, it should return undefined. This is
what Maple own help says. All these values below return 2.0

hypergeom([1,0],[0],.5); ---> 2
hypergeom([1,-1],[-1],.5); ---> 2
hypergeom([1,-2],[-2],.5); ---> 2
hypergeom([1,-3],[-3],.5); ---> 2
etc....

When Mathematica

Hypergeometric2F1[1, 0, 0, .5] ---> 1
Hypergeometric2F1[1, -1, -1, .5] ---> 1.5
Hypergeometric2F1[1, -2, -2, .5] ---> 1.75
Hypergeometric2F1[1, -3, -3, .5] ---> 1.9375
etc..

So, in Maple, you were basically doing this:

seq(evalf(2^(n-1)*binomial(n-1,n-1) * 2 ), n=0..7);
1., 2., 4., 8., 16., 32., 64., 128.

Here is the one-off I think Maple is doing:

hypergeom( [1,-3],[-2], 0.5);
Float(infinity)+Float(infinity)*I

But it should have returned undefined even for this:

hypergeom( [1,-2],[-2], .5); --> 2.0

Based on the help, which says "smaller absolute value"
Either way, it is not defined the same as Mathematica.

May be there is more than one definition to 2F1, I do
not know. May be a math person would know. May be Wikipedia
has something on this.

--Nasser

[Peter Luschny]

Nasser, you are absolutely on the right track!

How many users of Maple/Mathematica are aware of the fact
that Maple and Mathematica use different definitions
for the hypergeometric functions? I was not aware of this
fact. And I wonder how many bugs are out there just because
it was assumed naively that both definitions are the same.

Moreover, do the docs give any hints to this fact? Shouldn't
there be a big warning sign? I think that the docs are
not mathematical precise enough to enable this insight.
In any case this circumstance has certainly unhealthy
consequences on the entire scope of applications of CASs.

I will come back to a more precise description of the
differences below. But for now we have seen in the above
discussion what the user is really confronted with when
he dares to use hypergeometric functions. He has to struggle
with the question: "Is this ... Maple's definition or
a Maple bug or inconsistency or Mathematica's definition
or a Mathematica bug or inconsistency." OK, you never know
and this is what I called the "Joy of CAS". We have seen all
four cases in the above discussion.

Recently I found a place where the differences seems to
be explained. Look at formulas 15.2.5 and 15.2.6 at
http://dlmf.nist.gov/15.2#ii

Peter

[clicl...@freenet.de]

"Nasser M. Abbasi" schrieb:
>
> [had to snip Nasser's text because it is classified as spam.]
>

Let m,n be nonnegative integer numbers. The condition m<n for
2F1(a,-m;-n;z) to make sense follows from the series definition: the
hypergeometric series breaks off after the term involving z^(-n), and
terms must not become infinite up to that point. The condition naturally
extends to m=n with 2F1(a,-n;-n;z) = 1F0(a;z), which remains an
infinite series, however. There is no reason for Maple not to implement
this case too, for any pFq.

Martin.

[Peter Luschny]

> (c) Indeterminate
> If you think (c) is the right answer then please explain
> why you think the answer has to be different from (*).

Well, let us assume Maple and Mathematica do not agree
on a common definition and prefer to hazard the confusion
of their users.

Thankfully there is now a third library, which will be
integrated into Sage soon, mpmath, the new shooting star
in the hypergeometric industry.

So let's see what we get with mpmath.

def T(n, k):
return 2^k*binomial(n-1, k)*hyp2f1(1, -k, -n+1, 1/2)
for n in range(8):
print [T(n, k) for k in (0..n)]

[1.0]
[2.0, nan]
[1.0, 4.0, nan]
[1.0, 5.0, 8.0, nan]
[1.0, 7.0, 17.0, 16.0, nan]
[1.0, 9.0, 31.0, 49.0, 32.0, nan]
[1.0, 11.0, 49.0, 111.0, 129.0, 64.0, nan]
[1.0, 13.0, 71.0, 209.0, 351.0, 321.0, 128.0, nan]

The case k=n-1 is here like in Maple and the case
k=n gives not-a-number.

But there is an option, which I do not know where it is
documented, although this page mentions it existence:
http://mpmath.org/doc/0.19/functions/hypergeometric.html

Fredrik Johansson explained it on G+:
"In general, cancelling duplicate parameters is useful
because the hypergeometric functions of lower order are
much easier to compute. [..] Negative integers are just
not handled specially by mpmath, so you have to use
eliminate=False."

OK, so let's try this:

def T(n, k):
return 2^k*binomial(n-1,k)*hyp2f1(1,-k,-n+1,1/2,eliminate=False)
for n in range(8):
print [T(n, k) for k in (0..n)]

[1.0]
[1.0, nan]
[1.0, 3.0, nan]
[1.0, 5.0, 7.0, nan]
[1.0, 7.0, 17.0, 15.0, nan]
[1.0, 9.0, 31.0, 49.0, 31.0, nan]
[1.0, 11.0, 49.0, 111.0, 129.0, 63.0, nan]
[1.0, 13.0, 71.0, 209.0, 351.0, 321.0, 127.0, nan]

The case k=n-1 is here like in Mathematica but the case
k=n gives not-a-number (recall that Mathematica/Alpha yields
1 resp. 0 in this case).

Although it is a good idea to let the user choose different
definitions through an option (as soon this is sufficiently
documented) we see that in effect we now have a third variant.

Peter

[Christopher J. Henrich]

In article <ff403a68-2060-4333...@googlegroups.com>,

Unless this question is meant to be a stress test for computer algebra
systems, it has a very simple answer.

First, assume that n >= 0.
Change the variable of summation: j = n-i. Then we are looking at
sum_{j = 0 .. n} binomial(j-1, j). The first term of this sum is
binomial(-1,0) = 1; the rest are all 0. So the sum is 1.

If n < 0, then I think the range 0..n must be regarded as empty, and
the sum of an empty sequence of summands is 0.

--
Chris Henrich <http://www.mathinteract.com>
Yes, one can rant about the program designs, but generally things keep getting
more and more confused as time goes on. --Sea Wasp

[clicl...@freenet.de]

clicl...@freenet.de schrieb:

>
> Let m,n be nonnegative integer numbers. The condition m<n for
> 2F1(a,-m;-n;z) to make sense follows from the series definition: the
> hypergeometric series breaks off after the term involving z^(-n), and
> terms must not become infinite up to that point. The condition
> naturally extends to m=n with 2F1(a,-n;-n;z) = 1F0(a;z), which remains
> an infinite series, however. There is no reason for Maple not to
> implement this case too, for any pFq.
>

I should have looked more carefully: Maple seems to be doing just that,
so that its hypergeom([1,-n],[-n],.5) ---> 2 would be a perfectly valid
choice, whereas Mathematica seems to be doing something else, presumably
equally valid. My guess is that analytic continuation here depends on
the direction along which the point (-m,-n) is approached in (b,c)
space.

Martin.

[Axel Vogt]

I miss the time for the discussion, but Lebedev in Special Function says:

(in § 9.4): It follows that for fixed z in the plane cut
along [1, oo], the hypergeometric function F(a,b;c; z) is an entire function
of a and b , and a meromorphic function of c, with simple poles at the points
c = 0, - 1, - 2, . . .

So it seems it is how to approach that point (b,c) = (-n,n)

[clicl...@freenet.de]

Axel Vogt schrieb:

>
> On 28.07.2014 21:10, clicl...@freenet.de wrote:
> >
> > clicl...@freenet.de schrieb:
> >>
> >> Let m,n be nonnegative integer numbers. The condition m<n for
> >> 2F1(a,-m;-n;z) to make sense follows from the series definition: the
> >> hypergeometric series breaks off after the term involving z^(-n), and
> >> terms must not become infinite up to that point. The condition
> >> naturally extends to m=n with 2F1(a,-n;-n;z) = 1F0(a;z), which remains
> >> an infinite series, however. There is no reason for Maple not to
> >> implement this case too, for any pFq.
> >>
> >
> > I should have looked more carefully: Maple seems to be doing just that,
> > so that its hypergeom([1,-n],[-n],.5) ---> 2 would be a perfectly valid
> > choice, whereas Mathematica seems to be doing something else, presumably
> > equally valid. My guess is that analytic continuation here depends on
> > the direction along which the point (-m,-n) is approached in (b,c)
> > space.
> >
>

> I miss the time for the discussion, but Lebedev in Special Function says:
>
> (in � 9.4): It follows that for fixed z in the plane cut
> along [1, oo], the hypergeometric function F(a,b;c; z) is an entire function
> of a and b , and a meromorphic function of c, with simple poles at the points
> c = 0, - 1, - 2, . . .
>
> So it seems it is how to approach that point (b,c) = (-n,n)
>

Yes, (-n,-n)? "Meromorphic" in c would really be quite simple; could
there be problems with that statement if b = c = -n? Something
logarithmic? And a correction: My text quoted above should have read:
the series for F1(a,-m;-n;z) "breaks off after the term involving
z^(-m)".

Of course I agree with the view that a CAS must carefully document its
behavior in borderline situations.

Martin.

[garguni...@gmail.com]

https://www.youtube.com/watch?v=IAfptUkGhok