axiom-developer
[Top][All Lists]
Advanced

[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]]
> 





reply via email to

[Prev in Thread] Current Thread [Next in Thread]