[Axiom-developer] Patches

[William Sit]

Hi Martin:

Thanks for your reply. I am sorry that I reply to your post so late (actually, I
have not caught up with these posts and I only read them once in a while). I
have started verifying your bug on the NAG Axiom 2.3 version, which presumably
have the same algebra sources. There are some places where I cannot duplicate
your error messages. (I am also handicapped as far as Linux goes and so my Axiom
set up is on a separate computer (in fact a virtual computer) that makes it
clumsy to transfer text. I am using the NAG version because it helps to identify
whether the bug is from the algebra or the open source "meddling" :-)

On my system, which is not the open source one, I have this (transcribing ;-(
nonessential things omitted)

)clear all
(1) -> %A:SYMBOL
    %A

(2) -> new()$SYMBOL
    %D

Note this was because before I clear all, I had defined %B and did a
New()$SYMBOL that gives %C. 

(3) -> G1419::SYMBOL
    G1419

(4) -> symbol(GENSYM()$Lisp)
    G84337

All output are of type Symbol.

I also got the product summation display correctly without error message, unless
the upper limit is a constant. Thus I think Bronstein's use of new()$Symbol is
correct, but somewhere in the open source version, new()$Symbol get clobbered.
Of course, you can patch Bronstein's code to force a new symbol each time, but
other packages that use new()$Symbol will still fail, as you recognized. I
believe the source of the error is more subtle and probably does not come from
new()$Symbol either.

Another reason why I think the error does not come from the source code is
because any iteration in Axiom, the iterating variable is a dummy (and
supposedly unique and local in scope). (Can Tim affirm this?)

Bug Report #9216:

Regarding dvdsum, I find it hard to follow Bronstein's code
(CombinatorialFunction constructor in combfunc.spad), because I need to trace
all the packages to figure out what opdsum does (it seems to be a hold to do
recursion, but I am not sure and opdsum is not an exported function, neither is
dvdsum). Assuming it is just a summation, mathematically speaking, then my guess
is dvdsum does the following: it takes two arguments, a list with five
parameters, and a symbol. Using Bronstein notation in the code (without regard
to type, just mathematical and forget all the coercions), the first is
[f,k,y,g,h], all belonging to some function space F, such as Expression Integer
and the second is x. Here 

f is the function to be summed with dummy variable k, 
k is the dummy (for summation), 
y is some symbol
g is the lower limit of summation, presumably a function of x
h is the upper limit of summation, presumably a function of x
x is the symbol with respect to which differentiation is to be done


What the code dvdsum does, assuming odpsum is just summation, is to return the
following: d/dx means differentiation wrt x (not the d in the code)

0 if the variable y retracts to x, otherwise the formula:

d(h)/dx f(h) - d(g)/dx f(g) + summation(d(f)/dx, k = g..h)


I have no clue what this formula means, for example, what is the role of y?
however, I hesitate to say it is wrong because Bronstein is quite careful, both
as a mathematician and a wizard in Axiom and Aldor.  dvdsum occurs only three
times in the source, once in the signature, once in the definition, and once in
an obscure setProperty statement related to opdsum. It is not exported and
therefore not usable by any other package. But it is possibly used via opdsum
which is defined as operator("%defsum"), an operator with 5 parameters defined
in CommonFunctions package. I do not understand the design of that package.
(Manuel, please help!)

I don't know what dvdsum is supposed to return (so I can't say it is correct or
not), but you interpret it as differentiating the sum with respect to x (as is
clear from your bug reports). So I did the following experiment: I input the
following code, which is a "standalone" version of dvdsum that can be read into
the interpreter (assuming, opdsum is simply summation):


EI==>Expression Integer
dvdsum(l:List EI, x:Symbol):EI==
  x = retract(y:= third l)@Symbol => 0
  k := retract(d:=second l)@(Kernel EI)
  differentiate(h:=third rest rest l,x)*eval(f:=first l, k, h)
    - differentiate(g:=third rest l, x)*eval(f,k,g)
      + sum(differentiate(f,x),d::Symbol=g..h)

Now try

dvdsum([i^2, i, y, x^2, x^3],x)

     8    5
   3x - 2x

but 

D(sum(i^2, i=x^2..x^3),x)

       8    5     3     2
   (18x + 6x + 12x  + 3x  - 2x)
  ------------------------------
                6

which is correct (this can be done because we know how to sum i^2). So it seems
the story is more complicated.

In any case, to differentiate a symbolic summation, where the
independent variable x can appear ANYWHERE (in f, in g, or in h), is hopeless.
We cannot expect answers in all cases. We know it is not always possible to have
close form summation, let alone differentiating them. Things like the above can
work, but

D(sum(1/i, i=1..x),x)

would be impossible to differentiate without a closed form formula. Even when a
term by term differentiation makes sense,

  D(sum(1/i, i=x^2..x^3))

Axiom gives the wrong answer 1/x. Unfortunately, it agrees with 

  dvdsum([1/i,i,y,x^2,x^3],x)

So it seems under some circumstances (so far the no answer cases), dvdsum is
called to differentiate a sum and that's when it gives wrong answers.

----------------
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.

---------------
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).
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.


William
-- 
William Sit
Department of Mathematics..............Email: address@hidden
City College of New York..........................Tel: 212-650-5179
Convent Ave at West 138th Street..................Fax: 212-862-0004
New York, NY 10031............Axiom, A Scientific Computation Sytem
USA..........................http://www.nongnu.org/axiom/index.html