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

[William Sit]

Hi Martin:

You raised many issues. I am separating this reply into sections with dividing
lines, to make it look better. I think I found how to fix dvdsum.

--------------------------- Use of new()$Symbol in fspace.spad ---------------

"m.rubey" wrote:
> 
> On Thu, 10 Jun 2004, William Sit wrote:
> 
> > Note: the same kind of fix should probably be done in fspace.spad.
> 
> probably. TIM: Where should we note this? Maybe it would be good to have a
> list where we could note such things on savannah.

I take this back after studying fspace.spad more. There the %%0i symbols are not
used as dummy variables in loops. gendiff (%%0) is used to generate a list of
new symbols that are needed and Bronstein generates them as %%0i where i can be
from 1 to n using symsub to concatenate gendiff with i. Whether the list of new
symbols need to be "of the same family" (in the sense of say arbitrary constants
C_1, ..., C_n) or not, I don't know. But they way it is, there should not be any
bug.

--------------------------- Display with parenthesis in OutputForm -----------
--------------------------- local patch in combfunc.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)
> 
> I noticed this to. Bug report on the way

The fix is not as easy as it seems. It is not clear whether the "bug" should be
fixed locally in combfunc.spad (a simple, temporary, but unsatisfactory fix is
to modify ddprod and ddsum with paren(...)) with limited effect, or globally in
OutputForm (outform.spad), where the bug was perhaps recognized:

    -- Todo:
    ...
    -- bug in product, paren blankSeparate []

The temperary fix is unsatisfactory because we should not SIMPLY add parenthesis
to EVERY summation or product.  This is a general problem and requires an
examination of the context to decide whether the parenthesis is necessary and
that may be subjective. (I find controlling the output in Mathematica equally
difficult, but generally, Mathematica did a much better job -- except in cases I
don't like :-). Since OutputForm is built bottom up (the responsibility of each
domain via coerce), it may be difficult to make the decision even if it can be
done objectively.

So may be when you report the bug, you can supply the temporary patch? Thanks.

---------
Temporary local patch in combfunc.spad:

  ddprod(l, x) ==
--  prod(summand(l)::O, third(l)::O = fourth(l)::O, fourth(rest l)::O)
++  paren(prod(summand(l)::O, third(l)::O = fourth(l)::O, fourth(rest l)::O))

  ddsum(l, x) ==
--  sum(summand(l)::O, third(l)::O = fourth(l)::O, fourth(rest l)::O)
++  paren(sum(summand(l)::O, third(l)::O = fourth(l)::O, fourth(rest l)::O))

--------------------------- Patch for dvdsum in combfunc.spad ---------------
> 
> > (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.
> 
> I thought about this a little further. At first I thought that
> differentiating a sum with respect to one of its bound should yield an
> error!
> 
> The reason was, that this construct is not well-defined, since sum as a
> function of one of the bounds (say the upper) is (should be, I'll comment
> on this a little later) a function from the integers into some space:
> 
> n+->sum(f(i,a,n),i=a..n)
> 
> Therefore, there is no such thing as a differential!

The meaning of differenting a sum with variable bounds can be understood this
way, by treating the RESULT as a function of the variable. For example, assuming
x takes integer values, 

  x +-> sum(f(i,x), i=x^2..x^3)

defines a function of x as f(x^2,x)+ f((x^2+1),x) + ... + f(x^3,x)
whose domain can be generalized to the real numbers or even complex numbers.
Then we may ask what is its derivative and we can actually write it down with a
summation. More generally, we want to compute

  D(sum(f(i,x),i=a(x)..b(x)),x)

which should give

  sum(D(f(a(x) + i,x),x), i=0..(b(x) - a(x)))

So if you accept this interpretation, here is my new dvdsum (which evaluates in
all cases!) Note: before calling dvdsum, the system already checks that the
lower and upper limits are functions of x, so no checking is needed here.

  dvdsum(l,x)==
    x = retract(y := third l)@SE => 0
    k := retract(d := second l)@K
    f := first l
    g := third rest l
    h := third rest rest l
    opdsum [differentiate(eval(f,k,g+d),x),d,y,0,h-g]

Test this:

  f := operator 'f
  D(sum(f(i,x), i=a..b),x)

I tried other differentiation and they seem fine.

  D(sum(x^i, i=0..n),x)
  D(sum(f(i,x), i=x^2..x^3)
  D(sum(x/i, i=x^2..x^3)

The patch leaves the summation unevaluated if it cannot be done (this is
automatic due to Bronstein's original implementation).

------------------------- Other forms of summation results ---------------------
> 
> For example, usually we say that
> 
> sum(i,i=1..n) should give n*(n+1)/2
> 
> but I don't see any good reason, why it shouldn't be
> 
> n*(n+1)/2*cos(2*n*%pi)
> 
> which has a very different derivative with respect to n...

Sure, but what you are doing is having two functions with domain the real
numbers (say) having the same restriction on integers. As long as you are using
integer inputs, they are the same. When you start differentiating, you CHOOSE a
generalization (there are many, and there may or may not be a canonical one)
that enlarges the domain. These two generalisations are different functions (of
the continuous variable n) even though they coincide on integers. You may
arbitrarily interpolate any function for non-integral inputs. Their derivatives
don't have to coincide for integer values of n (or even differentiable at
integer values, but in your example, they are and do agree, but only for
integers).

Also, equality (not just for functions) is a messy business.

---------------------- Mathematica's handling of Summation -------------------

> However: both Maple 8 and Mathematica 5 return the result unevaluated, so
> we could do so, too. (However I don't know how to do it...)

This is true only if in 

  Sum[f[i,x], {i,a[x],b[x]}]

all the functions f, a, b are undefined. In many cases, Mathematica gives
results
that Axiom (with the patch above) cannot sum. For example, 
Mathematica even gives a result when f[i,x] = x/i in terms of polygamma
functions. Unfortunately, I cannot verify the result using integration. However,
I tried other simpler cases when a[x] and b[x] are known, and results are given
(Axiom can do those too), like Sum[i^2 x, {i,x^2, x^3}].

--------------------Discussion on Expression Constructor -------------------
-------------------- related: Segment, SegmentBinding ----------------------

> The one thing which really got me going on all this is the following: I
> wrote a program that *wants* to spit out things like
> 
> sum(product(f(i),i=1..j),j=1..n),
> 
> where f is a rational function in i, with coefficients in some (specified)
> infinite field.
> 
> Unfortunately, it turned out that this is not possible in axiom at the
> moment: sum (and product too) take an "Expression" as their first
> argument, and 2*x^3::POLY PF 3 is *not* an expression, since EXPR takes an
> OrderedSet as argument, which PF 5 is not. (curiously, a::EXPR CHAR is
> valid. I was also surprised to find that Complex ORDSET has ORDSET, namely
> lexicographically.)
> 
> Now, there are several things that I do not quite understand:
> 
> * Why does EXPR ask for an OrderedSet?

Good question. I am not sure about the answer. But here is my guess.
If we are to use EXPR R for the purpose of summation, then we need EXPR R to be
ordered so we can have segments. Even though construction of a segment does NOT
need an ordered set (so Segment PF 5 is perfectly legal), the real need is the
function expand from Segment(S), which requires S to be an OrderedRing (now in
Axiom, OrderedRing simply means the domain belongs to the category of Ring and
of OrderSet -- the order has nothing to do with the ring structure otherwise; so
since R is already a Ring in most cases, EXPR only need OrderedSet).  I believe
that is the origin. 

> 
> * Why don't we have EXPR POLY INT but only EXPR INT? Maybe the reason for
> this is that there is an operation ** in EXPR?

I don't understand. Sure, I typed EXPR POLY INT and it is perfectly ok. Even
EXPR FRAC POLY INT, EXPR FRAC INT.

> * How can I adapt EXPR so that I can have sums of POLY PF 5?

The question is, why isn't PF 5 an OrderedSet? The answer is: it is cylcic, so
one would have 0 < 1 < 2 < 3 < 4 < 5 = 0. FiniteField does have some rudimentary
ordering based on SegmentBinding and StepThrough:

s:SegmentBinding PF 5:= i=1..5

SegmentBinding, like Segment, does NOT require an ordered set and logically it
may even seem that for the purpose of summation, one only needs to be able to
"StepThrough" the segment with a "nextItem" function (all these ARE available,
even for PF) But unfortunately, it is not even possible to expand a segment into
a list for PF 5 because expand requires OrderedSet. 

  t:Segment PF 5:=1..3
  [i for i in t]

because the built-in "for" only allows segments of integer values for the
iterator i (it will step through with any list though, including List PF 5).

So you may have to create something like OrderedPrimeField where the order
function is explicitly given as say, for p = 5: 0 < 1 < 2 < 3 < 4


> 
> Consider
> 
> 2*x^3::POLY PF 3
> 
> should the 3 in x^3 be an element of PF 3 (i.e., be equal to zero) or
> should it be an integer? Currently there is no exponentiation where the
> exponent could be from PF 3, except of a highly suspicious ...

If you are working with polynomials, then the 3 in x^3 should be a NNI, not
integer, and certainly not PF 3, even if the coefficients are from PF 3. It
makes perfectly good sense to compute 2^3 = 8 = 2 (mod 3). Currently, it is not
possible to compute exponents mod 3 because the PolynomialCategory only allows
OrderedAbelianMonoidSup domains as exponents and PF 3 is not ordered.

> 
> * I think that the modeline of sum *should* be
> 
> (D1->EXPR Something, SEGBIND D1)->EXPR Something  where D1 has ORDSET.
> 
> Sorry for so many questions...

That Something unfortunately needs to be an ordered set in order to expand
SEGBIND D1.

---------------------------- Documentation on box in combfunc.spad ----------
> 
> > > > ----------------
> > > > 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))  (corrected typo from previous message)
> > > >
> > > > 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.
> 
> Thanks for this explanation! This should go into the docs!

I think the confusion was because the current description (see fspace.spad)
says:

box(f) returns f with a 'box' around it that prevents f from being evaluated
when operators are applied to it.

should have been

box(f) returns f with a 'box' around it that prevents the evaluation of an
operator when the operator is applied to f.

and similarly when box is applied to a list of arguments. However, box is
designed for operators (which are unexported and unary) and so it is difficult
in the interpreter to apply box to functions like sum that takes two arguments
(instead of a list to two elements). 
You certainly have my permission to adapt the above to the pamphlet.

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