axiom-developer
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Axiom-developer] Bug 312


From: Waldek Hebisch
Subject: [Axiom-developer] Bug 312
Date: Mon, 12 Feb 2007 21:43:39 +0100 (CET)

I would like here to make some remarks concering bug 312.  First,
this bug is really in 'powern' function defined in 
'sttaylor.spad.pamphlet'.  Below is (corrected) version contained
in wh-sandbox:

      powern(rn,x) ==
        order : I := 0
        for n in 0.. repeat
          empty? x => return zro()
          not zero? frst x => (order := n; leave x)
          x := rst x
          n = 1000 =>
            error "**: series with many leading zero coefficients"
        (ord := (order exquo denom(rn))) case "failed" =>
          error "**: rational power does not exist"
        co := frst x
        if ord > 0 and rn < 0 then
           error "**: negative power does not exist"
        (invCo := recip co) case "failed" =>
           error "** rational power of coefficient undefined"
-- This error message is misleading, isn't it? see sups.spad/cRationalPower
        power :=
--          one? co => YS(powerrn(rn,x,#1))
          (co = 1) => YS(powerrn(rn,x,#1))
          (denom rn) = 1 =>
            not negative?(num := numer rn) =>
-- It seems that this cannot happen, but I don't know why
              (co**num::NNI) * YS(powerrn(rn,(invCo :: A) * x,#1))
            (invCo :: A)**((-num)::NNI) * YS(powerrn(rn,(invCo :: A) * x,#1))
          RATPOWERS => co**rn * YS(powerrn(rn,(invCo :: A) * x,#1))
          error "** rational power of coefficient undefined"
        monom(1,(ord :: I) * numer(rn)) * power

-------------------------

At first glance this function looks cryptic.  But look at documentation:

      powern : (RN,ST A) -> ST A
        ++ powern(r,f) raises power series f to the power r.

So 'powern' is supposed to raises Taylor series to a rational power.
How one can to such thing?  If Taylor series has form f=1+c_1*x+...
(the lowest order term is zero) one can recursively compute coefficients
using the fact that:

  (f+a*x^n)^r = f^r+r*a*x^n mod x^(n+1)

Axiom code uses YS operator to solve this reccurence: this is done
in assignment to 'power' variable.

If the series has non-zero, non one lowest order term, that is

f = c_0*(1+c_1*x+...),

we just have to divide f by c_0 and multiply the result by
c_0^r (this is also handled in the assignment to power).

So the remaining case is when the lowest order term is 0.  In this
case 

f = x^a*f_1 

where f_1 has non-zero lowest order term.  We handle f_1 by earlier
method while (x^a)^r = x^(a*r).  However, we want a Taylor series
as a result, so we must require that a*r is a natural number.

Now, the first part of 'powern' handles splitting this power of x:

        order : I := 0
        for n in 0.. repeat
          empty? x => return zro()
          not zero? frst x => (order := n; leave x)
          x := rst x
          n = 1000 =>
            error "**: series with many leading zero coefficients"

(note that 'x' in the code is what I call 'f' above, while my 'a' is
called 'order' in the code).  Next, we check that order*r is an 
integer:

        (ord := (order exquo denom(rn))) case "failed" =>
          error "**: rational power does not exist"

Remember first nonzero coefficient:

        co := frst x

Check that we get nonnegative order*r:

        if ord > 0 and rn < 0 then
           error "**: negative power does not exist"

I skip assignment to power (explained eariler).  Finally we have
to multiply the result by correct power of 'x':

        monom(1,(ord :: I) * numer(rn)) * power

(note that ord = order/denom(rn), so ord * numer(rn) = order*rn).

Checking in SVN archive shows that the last line was omitted in
revision 19.  The change in revision 19 added support for negative
integer powers, it is clear that omission of last line was
unintentional.

-- 
                              Waldek Hebisch
address@hidden 




reply via email to

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