axiom-developer
[Top][All Lists]
Advanced

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

Re: [Axiom-developer] semantics of #1


From: Martin Rubey
Subject: Re: [Axiom-developer] semantics of #1
Date: Wed, 23 Jun 2004 16:54:52 +0000

Martin Rubey writes:
 > Does anybody know about the exact semantics of #1 in compiled code?
 > 
 > For example, the following does not work:
 > 
 > )abbrev package TEST Test
 > Test(F): Cat == Body where
 >     F: Field
 > 
 >     Cat == with
 >             tst: (F, List F, List F) -> Boolean
 > 
 >     Body == add
 >       tst (e, lst1, lst2) ==
 >         any?(e=elt(lst1, #1), lst2)
 
Sorry, this was nonsense: replacing the final List F by List Integer, the
example works. However, here is the real thing which I can't get to work:

)abbrev package RINTERPA RationalInterpolationAlgorithms
++ Description:
++ This package exports rational interpolation algorithms
RationalInterpolationAlgorithms(F, P): Cat == Body   where
    F: Field 
    P: UnivariatePolynomialCategory(F)
    Cat == with
        RationalInterpolation: (List F, List F, NonNegativeInteger, 
                                NonNegativeInteger) 
                               -> Union(Record(function: Fraction P,
                                               undefined: List Integer,
                                               unattainable: List Integer),
                                        "failed")

    Body == add
        RationalInterpolation(xlist, ylist, m, k) ==
            #xlist ^= #ylist =>
                error "Different number of points and values."
            #xlist ^= m+k+1 =>
                error "wrong number of points"
            tempvec: List F := [1 for i in 1..(m+k+1)]

            collist: List List F := cons(tempvec, 
                                         [(tempvec := [tempvec.i * xlist.i _
                                                       for i in 1..(m+k+1)]) _
                                          for j in 1..max(m, k)])

            collist := append([collist.j for j in 1..(m+1)], _
                              [[- collist.j.i * ylist.i for i in 1..(m+k+1)] _
                               for j in 1..(k+1)])
            resspace: List Vector F := nullSpace((transpose matrix collist) _
                                                 ::Matrix F)

            reslist: List List P := _
                      [[monomial((resspace.1).(i+1), i) for i in 0..m], _
                      [monomial((resspace.1).(i+m+2), i) for i in 0..k]]

            den := reduce((_+), reslist.2)
            zero?(den) => "failed"

            fun := reduce((_+), reslist.1)/den

            denres := denom(fun)

            una: List Integer := []
            und: List Integer := []
            for i in 1..#xlist repeat
              if zero?(elt(denres, xlist.i)) then 
                und := cons(i, und)
                una := cons(i, una)
              else if zero?(elt(den, xlist.i)) and 
                      elt(fun, xlist.i) ~= ylist.i then 
                     una := cons(i, una)

            [fun, una, und]

)abbrev package RINTERP RationalInterpolation
++ Description:
++ This package exports interpolation algorithms
RationalInterpolation(xx, F): Cat == Body   where
    xx: Symbol
    F:  Field
    UP  ==> UnivariatePolynomial
    SUP ==> SparseUnivariatePolynomial
 
    Cat == with
        interpolate: (Fraction UP(xx, F), List F, List F, _
                      NonNegativeInteger, NonNegativeInteger) _
                      -> Union(Record(function: Fraction UP(xx, F),
                                      undefined: List Integer,
                                      unattainable: List Integer),
                               "failed") 

    Body == add
        RIA ==> RationalInterpolationAlgorithms

        tst: (F, List F, List Integer) -> Boolean
        tst (e, lst1, lst2) ==
          any?(e=elt(lst1, #1), lst2)

        interpolate(qx, lx, ly, m, k) ==
            px := RationalInterpolation(lx, ly, m, k)$RIA(F, UP(xx, F))

            (px case "failed") => "failed"

            if ((qnum := retractIfCan(numer qx)@Union(F, "failed")) case F) and
               ((qden := retractIfCan(denom qx)@Union(F, "failed")) case F) then
              q := qnum/qden
--              if tst(q, lx, px.unattainable) then _
--                return "failed"
              if any?(q::F=elt(lx::List(F), #1), 
(px.unattainable)::List(Integer)) then 
                return "failed"                  

            [elt(px.function, qx), px.undefined, px.unattainable]

using tst, it works, using any?, it doesn't...

HELP!





reply via email to

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