RE: [Axiom-mail] Another beginner's question: Lagrange interpolation

[Bill Page]

Alasdair,

Thanks for the question!

On Sent: April 12, 2007 9:45 PM Alasdair McAndrew wrote:

> Just to explore Axiom a little further, I was trying to do a
> Lagrange interpolation on two lists:
>
> LagrangeInterpolation([1,2,3],[4,5,-6])
>
> This doesn't work - I get an error message which I (as a beginner)
> don't understand.

This is what I get:

(1) -> LagrangeInterpolation([1,2,3],[4,5,-6])
   There are no exposed library operations named LagrangeInterpolation
      but there is one unexposed operation with that name. Use HyperDoc
      Browse or issue
                      )display op LagrangeInterpolation
      to learn more about the available operation.

   Cannot find a definition or applicable library operation named
      LagrangeInterpolation with argument type(s)
                            List PositiveInteger
                                List Integer

      Perhaps you should use "@" to indicate the required return type,
      or "$" to specify which version of the function you need.


> So then: 
>
> 1)  How do I do this?

I recommend first of all using the suggested command:

(1) -> )display op LagrangeInterpolation

There is one unexposed function called LagrangeInterpolation :
   [1] (List D3,List D3) -> D1 from PolynomialInterpolationAlgorithms(
            D3,D1)
            if D3 has FIELD and D1 has UPOLYC D3

> 2)  Without asking this mail group, where can I find this
> information in a manner which I can understand?

Hmmm.... what should I assume about what you can understand?
That is a difficult problem and the answer will change over
time.

With the above information you should be able to see that
at least that:

- the LagrangeInterpolation function does exist in Axiom
  but it is "not exposed"
- that it expects two lists as arguments
- and these lists must consist of elements of a Field
- and that it will return a polynomial over this Field
- and finally that what you passed as arguments to this
  function did not meet these requirements.

The first thing I would try is to look in the Axiom book by
using the search function in your PDF reader. No luck. But
at least you should be able to read there about what the phrase
"not exposed" means. Essentially it means that you must call
the function explicitly from the package where it is defined.

The second thing you might try is to look this function up
in HyperDoc using Browse, enter the name LagrangeInterpolation
and click Operations. Unfortunately this only tells you what
you know already plus the fact that it is "not documented
yet"...

You should also ask Axiom about 'PolynomialInterpolationAlgorithms'

(1) -> )show PolynomialInterpolationAlgorithms
 PolynomialInterpolationAlgorithms(F: Field,P: UnivariatePolynomialCategory
F)
is a package constructor
 Abbreviation for PolynomialInterpolationAlgorithms is PINTERPA
 This constructor is not exposed in this frame.
 Issue )edit c:/Program
Files/axiom/mnt/windows/../../src/algebra/PINTERPA.spad
to see algebra source code for PINTERPA

------------------------------- Operations --------------------------------
 LagrangeInterpolation : (List F,List F) -> P

So now take another look at the last point and what you
actually passed to the function:

                            List PositiveInteger
                                List Integer

So Axiom thought the first list [1,2,3] was a list of
PositiveInteger and the second list [4,5,-6] was a list of
Integer. Lets ask Axiom if these are Fields:

(1) -> PositiveInteger has Field

   (1)  false
                                          Type: Boolean


(2) -> Integer has Field

   (2)  false
                                          Type: Boolean

but at least

(3) -> Fraction Integer has Field

   (3)  true
                                          Type: Boolean

Oh! So how can I tell Axiom that these numbers are actually
from some Field? Try this:

(4) -> [1,2,3]::List Fraction Integer

   (4)  [1,2,3]
                               Type: List Fraction Integer

Ok, now lets try LagrangeInterpolation again:

(9) -> A:=[1,2,3]::List Fraction Integer

   (9)  [1,2,3]
                                 Type: List Fraction Integer

(10) -> B:=[4,5,-6]::List Fraction Integer

   (10)  [4,5,- 6]
                                 Type: List Fraction Integer

But how does one call the specific package where this function
is located? The command ')show PolynomialInterpolationAlgorithms'
gave use this information:

 PolynomialInterpolationAlgorithms(F: Field,P: UnivariatePolynomialCategory
F)

So now in my example above the Field is 'Fraction Integer' but it
may not be entirely obvious to you what

 UnivariatePolynomialCategory F

means. To make a long story short: Axiom groups it's mathematical
domains (types) into groups called "Category". This is almost (but
not quite) like a category in mathematics. Anyway let me suggest
one such domain that will work for us: The univariate polynomials
over the variable x with coefficients in the same field:

  UnivariatePolynomial(x,Fraction Integer)

so the final result looks like this:

(11) ->
LagrangeInterpolation(A,B)$PolynomialInterpolationAlgorithms(Fraction
Integer,UnivariatePolynomial(x,Fraction Integer))

             2
   (11)  - 6x  + 19x - 9
                            Type: UnivariatePolynomial(x,Fraction Integer)

There, that was easy! Right? ;) Well, perhaps not for the
beginner...

As you get a little more experienced with Axiom you might look
in Axiom's source code itself. In principle this source code is
supposed to consist of both code and the complete documentation
of the code. here:

http://wiki.axiom-developer.org/axiom--test--1/src/algebra

This is a searchable online copy of the source code for the
entire Axiom Library. Enter 'LagrangeInterpolation' in the
search box at the top right of the page, hit enter and you will
see all the places where LagrangeInterpolation is used in the
code. Usually there is also at least some minimal documentation
included in the code, but in this case we are out of luck - the
original Axiom developer did not even leave us any hints. The
best we get is a quite readable algorithm written in SPAD -
the Axiom library high level programming language.

        LagrangeInterpolation(lx, ly) ==
            #lx ^= #ly =>
                error "Different number of points and values."
            ip: P := 0
            for xi in lx for yi in ly for i in 0.. repeat
                pp: P := 1
                xp: F := 1
                for xj in lx for j in 0.. | i ^= j repeat
                    pp := pp * (monomial(1,1) - monomial(xj,0))
                    xp := xp * (xi - xj)
                ip := ip + (yi/xp) * pp
            ip

---------

Which I think is actually quite readable, but perhaps too much
for the beginner.

Another important place to look online is the Axiom Input test
files. There are a lot of useful (but sparsely documented)
examples here:

http://wiki.axiom-developer.org/axiom--test--1/src/input

But searching for 'LagrangeInterpolation' here we are again
out of luck.

---------

I am sorry that this email reply has gotten rather long, but
I wanted to illustrate the kind of help that I wish that one
day will be available on the Axiom Wiki and other Axiom
publications.

If you have any questions or comments about this please
continue.

Regards,
Bill Page.