[Top][All Lists]
[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
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Axiom-developer] Bug 312,
Waldek Hebisch <=