[Axiom-developer] Bugs in combfunc.spad and partial patch

[William Sit]

Hi Martin:

Thanks for the detail explanations. The following discussion relates to
combfunc.spad (by Bronstein) and bug report 9218. I am changing the thread
heading to reflect the context of this discussion (previous thread: Patches)

> -- cannot use new()$Symbol because of possible re-instantiation
>       gendiff := "%%0"::SY

Could this be in the open source version, something gets re-instantiated when it
shouldn't? It seems Bronstein thinks it is ok to use one dummy within a package,
but worries when two or more are instantiated (hence the %% instead of %, which
appears in the output of Bug Report #9218).

> If I understand correctly, dummy := new()$SE::F assigns a new symbol to
> the variable "dummy" once and for all. So, if I'm not mistaken, the policy of
> axiom which you describe is violated in COMBF, and probably in other
> packages too. It is not so unlikely that this produces errors only very
> rarely, since errors will only occur if you nest procedures from one
> package ...

I agree with you on both. The first is so that the two occurrences of dummy in a
loop gets compiled into the same variable. 

The second, about the violation, is subtle. Of course, it is wrong to code:

   [[k for k in 1..5] for k in 1..6]

because it is expected that the body of the loop will involve the outer k and it
is not possible in general to figure out which k the k in the inner body refers
to. So my "local principle" cannot be automatically enforced in a nested loop:
there is no way to reassign the inner k or the outer k to a new variable to make
it work. I think, however, the Axiom compiler or interpreter should give a
warning or error flag on this one (currently, it just gives an answer).

However, it is perfectly ok to do:

  k:= 7
  foo(i)==[k for k in 1..i]
  bar(t,j)==[t for k in 1..j]

and call bar(foo(2),3). One is allowed to reuse k as long as it is not nested
even if the functions are later composed and k is also globally defined. Now
Bronstein's case is in theory similar but really not.  In combfunc.spad, if I
did not miss any count, Bronstein never used nested loops with the same dummy
identifier. There are only four functions (product, summation; each with dual
signatures) using the identifier dummy, each time in a simple "loop" (see
later). It was in the interpreter when in your example, a nested loop was used
with two functions by composing product with summation. But surely, there is
nothing wrong in the code to use the same iteration variable in different loop
constructs that are not nested. 

So where is the error? I think it lies in two facts:

(a) the identifier dummy is not a Symbol, but coerced into F (a domain such as
Expression Integer), and 

(b) Bronstein did not use the identifier dummy directly in a for-loop, because
it would not be allowed if the lower or upper limits of the summation (or
product) are not retractible to integers (say functions of x and therefore
"indefinite") in a for-construction. 

A trace of the functions (thanks for your helpful pointers) iidsum (for integer
to integer definite summation, I suppose) shows that when the limits are both
integers, a simple term by term summation is used to evaluate the sum, and the
code is a simple for-loop with a symbol identifier as the loop variable. When
either one is anything else, the function idsum (for indefinite summation?) is
called which simply holds the sum unevaluated (kernel(opdsum,l)).

Thus the compiler did not treat the identifier dummy as a local variable but as
a global variable (inspection of code.lsp confirms this).

So there are four problems to be fixed:

(0) new()$Symbol should generate a new one each time it is used (this bug does
not exist in Axiom 2.3 NAG, and seems unrelated to Bug Report 9218)

(1) the identifier dummy needs to be different for EACH INVOCATION of the four
functions in combfunc.spad (two product and two summation functions), not just
for each FUNCTION, because one can still compose product with product as in
    product(product(i*j, i=a..b),j=c..d)

The fix is not to use four separate dummy identifiers, but rather to IN-LINE
new()$Symbol with a LOCAL identifier in each of the four functions, AND WHEN
COMBINED WITH the fix (0), this problem can be resolved. See proposed patch at
end (tested in NAG version). Note: the same kind of fix should probably be done
in fspace.spad.

BTW, currently, even in NAG version, such compositions where product may be
replaced by summation, gave wrong results (as well as wrong displays for lack of
needed parenthesis -- to see true returned results, )set output tex on. So:

(1a) wrap outputs in product or summation in parenthesis for correct display:
example of error (ignore error in the value):

   )set output tex on
   product(summation(i*j, i=a..b),j=c..d)

(2) the formula for differentiating a sum (in dvdsum) needs to be revised and
the call to dvdsum (or its equivalences) should return an unevaluated expression
when the summation cannot be evaluated in close form. A similar problem probably
also exist for differentiation a product.


> > ----------------
> > Bug Report #9215:
> >
> >  sum(box(i),i=a..b)
> >
> > returns the answer correctly, but the outputForm missed a pair of 
> > parenthesis
> > aa-1 should be a(a-1);
> >
> >  box(sum(i),i=1..n)
> >
> > also returns correctly, because after sum is evaluated to (n^2+n)/2, an
> > (invisible) box is placed around it. Note the box is useful ONLY when
> > there is an operator operating on the box content. The only way this is
> > easy to use is when the operator takes one argument, such as sin, cos,
> > log. The argument can be a list of course.
> 
> I thought that box would simply prevent axiom from doing evaluation. I
> still don't understand.

box is a function, so all its arguments are evaluated before being passed. What
box prevents evaluation is the function of which it is the argument. So in the
first example, sum(box(i),i=a..b), there is NO function of which box(i) is an
argument. It is not easy to wrap two arguments in the box without destroying the
syntax, since box takes only one argument. Now if you were to try
  factor((x-1)*(x-2))
and 
  factor(box((x-1)*(x-2)))
you will see the difference. In the first case, the answer is factored, in the
second case it is not. In each case, the two given factors were multiplied
before passing to box, but in the second case, factor does not apply.

> > ---------------
> > Bug Report 9218:
> >
> > I believe the CombinatorialFunctions package is not meant to do what may
> > be called indefinite summation, only definite summations (hence the name
> > %defsum).
> 
> I don't understand the relation of this to my bug report. Sorry.
> 
> > Note that the result from
> >
> > eval(D(f(x),x),f,y+->sum(g(i,y),i=a..y))
> >
> >    %defsum (g(%A,%%01),%A,i,a,x)
> >
> > comes from dvdsum again and is consistent with the formula I gave above,
> > except that I would not expect %%01 but rather just x.
> 

What I meant was that from the trace of the code, Bronstein wants to leave a
summation unevaluated if either of the summation limit is "indefinite" (in his
code, this means not an integer). I am not sure under what condition he intended
dvdsum to be invoked, but as you pointed out, the formula encoded in dvdsum can
only be true under certain conditions. Now it seems that every time a summation
cannot be evaluated, dvdsum is called, giving wrong results in the cases we
tried. The eval example above is mathematically equivalent (at least I think
that was your intention because f is just an unknown operator) to

  D(sum(g(i,x),i=a..x),x) which gives

summation (displayed) from i=a to x of the partial with respect to the second
argument of g at (i,x) plus g(x,x), which is following the formula in dvdsum
because the lower limit is independent of x. Normally, eval is not executed
until all its arguments are evaluated first. Now Bronstein implements diffEval
in fspace.spad, so eval probably should not perform a substitution
differentially (replacing f and differentiate). It should not have succeeded in
any substitution since f no longer "appear" in f'(x) --- at least that is the
case in Mathematica or Maple. In Axiom, it seems the substitution is still made.
So the answer should be the same as above, with two terms: a summation and
g(x,x).

> The %%01 is the dummy variable used for differentiation (see line 625 of
> fspace.spad).

Are you suggesting that 

  %defsum (g(%A,%%01),%A,i,a,x)

means sum from i=a to x the partial of g with respect the second variable
evaluated at (i,x)? What happened to g(x,x)? This is not consistent with:

eval(D(f(x),x),f,y+->sum(g(i+y),i=a..y))

which gives

  %defsum (g(%A+%%01),%A,i,a,x)

All this is garbage-in garbage-out because such operations were probably not
intended. But this is only a conjecture (that is what archeologists do). We can
follow the code, but it does not make sense to me.

William
------
Proposed fix for (1) (tested for NAG version) with the following inputs:

  product(summation(i*j), i=a..b), j=c..d)

and similar variations. However, I have NOT tested whether this will affect
other computations! (It should not, but who knows?)

Patch for combfunc.spad:

Step 1. Delete line
  dummy := new()$SE::F

Step 2. Replace the definition of product, summation by:

  product(x:F, i:SE) ==
--  opprod [eval(x, k := kernel(i)$K, dummy), dummy, k::F]
++  opprod [eval(x, k := kernel(i)$K, h:=new()$SE::F), h, k::F]

  summation(x:F, i:SE) ==
--  opsum [eval(x, k := kernel(i)$K, dummy), dummy, k::F]
++  opsum [eval(x, k := kernel(i)$K, h:=new()$SE::F), h, k::F]

  product(x:F, s:SegmentBinding F) ==
    k := kernel(variable s)$K
--  opdprod [eval(x, k, dummy), dummy, k::F, lo segment s, hi segment s]
++  opprod [eval(x, k, h:=new()$SE::F), h, k::F, lo segment s, hi segment s]]

  summation(x:F, s:SegmentBinding F) ==
    k := kernel(variable s)$K
--  opdsum [eval(x, k, dummy), dummy, k::F, lo segment s, hi segment s]
++  opdsum [eval(x, k, h:=new()$SE::F), h, k::F, lo segment s, hi segment s]]