7.122 bug in product, Maple V to Maple 8 (18.11.98)

7.122.1 CORNIL Jack Michel
7.122.2 Robert Israel (20.11.98)
7.122.3 Helmut Kahovec (22.11.98)
7.122.4 Willard, Daniel, Dr. (30.11.98)

7.122.1 CORNIL Jack Michel

A new bug ?

> product(1/(i-k),k=0..i-1); 
 
                                         0
 

I do not think it is possible.

7.122.2 Robert Israel (20.11.98)

Actually I think it’s quite an old bug. The explanation in a nutshell is this: Maple first evaluates the indefinite product:

   product(1/(i-k), k) = (-1)^k/GAMMA(k-i).
 

This is true in the sense that if you call this F(k), then F(k+1) = F(k)/(i-k) (as a function of i). Of course GAMMA has poles at the integers <= 0, so it doesn’t really work for the case Maple wants: it wants to calculate the definite product as F(i)/F(0). Now 1/F(0) = GAMMA(-i) which is undefined (but at this point Maple doesn’t realize that fact because it doesn’t know i is supposed to be a positive integer).

On the other hand, F(i) = (-1)^i/GAMMA(0) and Maple does know that GAMMA has a pole at 0. So it evaluates limit(1/F(k),k=0) and correctly finds 0. Therefore the answer is evaluated as (1/GAMMA(-i))*0 = 0.

On one level the problem is one of recognition of singularities (or zeros): Maple doesn’t realize that GAMMA(-i) is undefined. I would have thought that it might help to explicitly assume i is a positive integer:

> assume(i, posint);
 

Unfortunately, that doesn’t work: even if you assume i is a positive integer, GAMMA(-i) just returns unevaluated.

The problem, I think, is in the code for GAMMA which checks for (actual) integers <= 0, but doesn’t worry about values that are known to be integers <= 0: it uses type(a, integer) rather than is(a, integer).

7.122.3 Helmut Kahovec (22.11.98)

it’s a bug. When calculating

> restart; 
> product(1/(i-k),k=0..i-1); 
 
                                  0
 

Maple (Release 4) executes the following statements:

> p:=product(1/(i-k),k); 
 
                                       k 
                                   (-1) 
                         p := ------------- 
                               GAMMA(-i + k) 
 
> p1:=eval(subs(k=0,p)); 
 
                                      1 
                           p1 := --------- 
                                  GAMMA(-i) 
 
> p2:=eval(subs(k=i-1,p)); 
Error, (in GAMMA) singularity encountered 
 
> p2:=limit(p,k=i-1); 
 
                               p2 := 0 
 
> p2/p1; 
 
                                  0
 

You may verify this with the help of the Maple debugger (some extensive output shortened):

> stopat(product,73); 
 
                              [product] 
 
> product(1/(i-k),k=0..i-1); 
(-1)^k/GAMMA(-i+k) 
product: 
  73*       if r = lasterror then 
              ... 
            else 
              ... 
            fi 
> step 
(-1)^k/GAMMA(-i+k) 
product: 
  75          r := `product/subscheck`(r,i,a,b); 
> step 
`product/subscheck`: 
   1    nexpr := expr; 
> step 
(-1)^k/GAMMA(-i+k) 
`product/subscheck`: 
   2    tp1 := traperror(eval(subs(i = a,nexpr))); 
> step 
1/GAMMA(-i) 
`product/subscheck`: 
   3    if tp1 = lasterror or tp1 = 0 or tp1 = FAIL then 
          ... 
        fi; 
> step 
1/GAMMA(-i) 
`product/subscheck`: 
   5    tp2 := traperror(eval(subs(i = b,nexpr))); 
> next 
Error, singularity encountered 
`product/subscheck`: 
   6    if tp2 = lasterror or tp2 = 0 or tp2 = FAIL then 
          ... 
        fi; 
> step 
Error, singularity encountered 
`product/subscheck`: 
   7      tp2 := limit(nexpr,i = b) 
> step 
limit: 
   1    if nargs < 2 or 3 < nargs or ... then 
          ... 
        fi; 
> return 
0 
`product/subscheck`: 
   8    if tp1 = undefined or ... then 
          ... 
        else 
          ... 
        fi 
> return 
0 
product: 
  76          if r = FAIL then 
                ... 
              else 
                ... 
              fi 
> step 
0 
Warning, statement number may be incorrect 
product: 
  12        RETURN(r) 
> step
 

7.122.4 Willard, Daniel, Dr. (30.11.98)

That matter of "type(a, integer)" versus "is(a, integer)" bedevils many other programs, for example, the computation of Bessel functions, which will not provide a series solution if the index is not "type(n, integer)" though there are perfectly good series expansions for arbirary real or complex n. I have been complaining about this since the first release of Maple V, without any effective response from Maple or Waterloo.