[Axiom-developer] Re: [Axiom-math] Curious behavior of Taylor series

[Martin Rubey]

"Igor Khavkine" <address@hidden> writes:

> On 21 Aug 2006 22:13:13 +0200, Martin Rubey <address@hidden> wrote:

> I can try. But I hope my small patch will still be accepted even if I
> can't find the time to thoroughly document that function.

I'm afraid: no. In fact, I will very likely use it myself, but it won't go into
the distribution until it is documented.

But it is certainly enough to explain roughly what is happening in every couply
of lines in powern and to explain exactly what your patch is doing. Thus,
along the lines of


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

      smult: (RN,ST A) -> ST A
      smult(rn,x) == map(rn * #1,x)
@

The following operation [[powerrn]] computes lazily ???

<<package STTAYLOR StreamTaylorSeriesOperations>>=
      powerrn:(RN,ST A,ST A) -> ST A
      powerrn(rn,x,c) == delay
        concat(1,integ(smult(rn + 1,c * deriv x)) - rst x * c)
@

[[powern(rn, x)]] computes $x^{rn}$, $x$ being a power series. Compare with
[[cRationalPower$ISUPS]].

<<package STTAYLOR StreamTaylorSeriesOperations>>=
      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"
@

First we determine the order of [[x]], i.e., the index of the first non-zero
coefficient. (Is that true ???)

Remarks:
\begin{itemize}
\item Note that usually [[leave]] takes no arguments. I don't know what it does
      here.
\item [[empty? x]] tests whether the stream [[x]] has no elements. This is
      mathematically the same as the corresponding power series being zero.
\end{itemize}

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

[[ord]] will be the order of the new power series. If it is not an integer, we
won't get a Taylor expansion.

<<package STTAYLOR StreamTaylorSeriesOperations>>=
        co := frst x
        (invCo := recip co) case "failed" =>
           error "** rational power of coefficient undefined"
-- This error message is misleading, isn't it? see sups.spad/cRationalPower
@

The leading coefficient has to be invertible???

<<package STTAYLOR StreamTaylorSeriesOperations>>=
        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"

        zeropref: ST A := new(ord::I::NNI,0::A)$LA :: ST A
        concat(zeropref, power)
@


<<package STTAYLOR StreamTaylorSeriesOperations>>=
    if A has Field then

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

> In any case, I'm not sure I completely understand what powern does.  I've
> looked at servarl power series manipulation routines already, but I get
> stumped every time I run into this object: Y from
> ParadoxicalCombinatorsForStreams. Can someone explain what it does?

Not really. The documentation says

  A  :   Type
  ST ==> Stream

    Y: (ST A -> ST A) -> ST A
      ++ Y(f) computes a fixed point of the function f.

In our case A is the ring of coefficients and powern calls

Y(powerrn(rn,x,#1))

Here, rn is the exponent, and x is the power series. Thus, Y computes a fixed
point - whatever that is - of the function

f +-> powerrn(rn, x, f)

I guess it would help to know what powerrn does. A different possibility to
find out is to look at cRationalPower in sups.spad.pamphlet. The block
(formally) corresponding to 

        power := 
          one? co => YS(powerrn(rn, x, #1))
          one? (denom rn) =>
            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"

        zeropref: ST A := new(ord::I::NNI, 0$A)$LA :: ST A
        concat(zeropref, power)


reads (ccInv == invCo, n == ord, cc == co, r == rn, uts == x)

        ccPow :=
          one? cc => cc
          one? denom r =>
            not negative?(num := numer r) => cc ** (num::NNI)
            (ccInv::Coef) ** ((-num)::NNI)

          RATPOWERS => cc ** r
          error "** rational power of coefficient undefined"

        uts1 := (ccInv::Coef) * uts
        uts2 := uts1 * monomial(1, -order)
        monomial(ccPow, (n::I) * numer(r)) * cPower(uts2, r::Coef)



So, it seems that in ISUPS, the equivalent to YS(powerrn(rn, (invCo::A) * x,
#1)) is factored out. We have the equivalence

ccPow  * z^(n*numer r) * cPower(ccInv * uts * z^(-order), r)

power  * z^ord         * YS(powerrn(rn, (invCo::A) * x, #1)) 

Note that n*numer r = order/denom r * numer r = order * r

     but  ord       = order/denom rn          = order * rn / numer rn

so that "explains" the appearance of z^(-order) in the ISUPS version.


I have to stop here, sorry,

Martin