[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Axiom-developer] Re: Bugs in combfunc.spad and partial patch
From: |
m.rubey |
Subject: |
[Axiom-developer] Re: Bugs in combfunc.spad and partial patch |
Date: |
Fri, 11 Jun 2004 12:32:10 +0200 (CEST) |
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.
> 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
> (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!
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...
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...)
Furthermore: If we decide to yield an error, sums that can be evaluated in
closed form will still give the "expected" (in the naive sense) result.
Thus, this might be somewhat inconsistent...
-----------------------------------------------------------------------------
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?
* 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?
* How can I adapt EXPR so that I can have sums of POLY PF 5?
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
[21] (D,D1) -> D from D
if D has UTSCAT D1 and D1 has RING and D1 has FIELD
(As an aside: Indeed, strange things happen when we use this ^:
(x::UPXS(PF 5,x,0))^(4::PF 5)
>> Error detected within library code:
"log of function with singularity"
)
Well, thinking about it, this doesn't seem to be a problem. The user can
always specify which operation should be used.
* 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...
> > > ----------------
> > > 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.
Thanks for this explanation! This should go into the docs!
All the best, Martin
I did not yet have time to think about the following, hence no comment...
> > > ---------------
> > > 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]]
>
- [Axiom-developer] Bugs in combfunc.spad and partial patch, William Sit, 2004/06/10
- Re: [Axiom-developer] Bugs in combfunc.spad and partial patch, William Sit, 2004/06/10
- [Axiom-developer] Re: Bugs in combfunc.spad and partial patch,
m.rubey <=
- [Axiom-developer] Re: Bugs in combfunc.spad and partial patch, William Sit, 2004/06/11
- Message not available
- Message not available
- Message not available
- Message not available
- [Axiom-developer] Re: EXPR POLY INT, Martin Rubey, 2004/06/15
- Re: [Axiom-developer] Re: EXPR POLY INT, root, 2004/06/15
- [Axiom-developer] Re: EXPR POLY INT, William Sit, 2004/06/15
- [Axiom-developer] Re: EXPR POLY INT, Martin Rubey, 2004/06/15
- [Axiom-developer] Re: EXPR POLY INT, Martin Rubey, 2004/06/16
- [Axiom-developer] Re: EXPR POLY INT, root, 2004/06/16
- [Axiom-developer] Re: EXPR POLY INT, William Sit, 2004/06/17