[Axiom-developer] [Guess] (new)

[billpage]

Changes http://page.axiom-developer.org/zope/mathaction/Guess/diff
--
\begin{axiom}
)abbrev package GUESS Guess
++ Description:
++ This package implements guessing of sequences
Guess(F, S, EXPRR, R, retract, coerce): Exports == Implementation where
    F: Field                                 -- zB.: FRAC POLY PF 5
-- in F wird interpoliert und gerechnet 

    S: GcdDomain

-- in guessExpRat m�chte ich die Wurzeln von Polynomen �ber F bestimmen. Wenn F
-- ein Quotientenk�rper ist, dann kann ich den Nenner loswerden.
-- F w�re dann also QFCAT S

    R: Join(OrderedSet, IntegralDomain)      -- zB.: FRAC POLY INT

-- die Ergebnisse werden aber in EXPRR dargestellt
    EXPRR: Join(ExpressionSpace, IntegralDomain, 
                RetractableTo R, RetractableTo Symbol, 
                RetractableTo Integer, CombinatorialOpsCategory) with
              _* : (%,%) -> %
              _/ : (%,%) -> %
              _*_* : (%,%) -> %
              numerator : % -> %
              denominator : % -> %

                                             -- zB.: EXPR INT
-- EXPR FRAC POLY INT ist ja verboten. Deshalb kann ich nicht einfach EXPR R
-- verwenden
-- EXPRR existiert, falls irgendwann einmal EXPR PF 5 unterst�tzt wird.


-- das folgende m�chte ich eigentlich loswerden.
  
    retract: R -> F                          -- zB.: i+->i
    coerce: F -> EXPRR                       -- zB.: i+->i
-- Achtung EXPRR ~= EXPR R

    NNI ==> NonNegativeInteger
    PI ==> PositiveInteger
    EXPRI ==> Expression Integer
    GUESSRESULT ==> List Record(function: EXPRR, order: PI)
    GUESSER ==> (Symbol, List F, BASIS) -> GUESSRESULT
    BASIS ==> EXPRI -> EXPRR                 -- zB.: i+->q^i
-- brauche hier EXPRI, weil das Ergebnis ja die allgemeine Form enth�lt

    Exports == with
      guess: (Symbol, List F, BASIS, List GUESSER, List Symbol, NNI) -> 
GUESSRESULT

      guessRat: (Symbol, List F, BASIS) -> GUESSRESULT

      guessExpRat: (Symbol, List F, BASIS) -> GUESSRESULT

      guessPade: (Symbol, List F, BASIS) -> GUESSRESULT

      guessPRec: (Symbol, List F, BASIS) -> GUESSRESULT

    Implementation == add

      basis2: (BASIS, Integer) -> F
      basis2 (basis, i) == retract(retract(basis(i::EXPRI))@R)

      SUP ==> SparseUnivariatePolynomial
      SUPF ==> SUP F
      FSUPF ==> Fraction SUP F
      V ==> OrderedVariableList(['a1, 'A])
      POLYF ==> SparseMultivariatePolynomial(F, V)
      FPOLYF ==> Fraction POLYF
      FSUPFPOLYF ==> Fraction SUP FPOLYF
      POLYS ==> SparseMultivariatePolynomial(S, V)

      MPCSF ==> MPolyCatFunctions2(V, IndexedExponents V, 
                                      IndexedExponents V, S, F,
                                      POLYS, POLYF)
      MPCFS ==> MPolyCatFunctions2(V, IndexedExponents V, 
                                      IndexedExponents V, F, S,
                                      POLYF, POLYS)

      SUPF2EXPRR: (Symbol, SUP F) -> EXPRR
      SUPF2EXPRR(xx, p) ==
        zero? p => 0
        (coerce(leadingCoefficient p))::EXPRR * (xx::EXPRR)**degree p
           + SUPF2EXPRR(xx, reductum p)

      FSUPF2EXPRR: (Symbol, FSUPF) -> EXPRR
      FSUPF2EXPRR(xx, p) ==
        (SUPF2EXPRR(xx, numer p)) / (SUPF2EXPRR(xx, denom p))

      SUPS2SUPF: SUP S -> SUP F
      SUPS2SUPF p ==
        if F is S then 
          p pretend SUP(F)
        else if F is Fraction S then
          map(coerce(#1)$Fraction(S), p)
            $SparseUnivariatePolynomialFunctions2(S, F)
        else error "Type parameter F should be either equal to S or equal _
                    to Fraction S"

      POLYS2POLYF: POLYS -> POLYF
      POLYS2POLYF p ==
        if F is S then 
          p pretend POLYF
        else if F is Fraction S then
          map(coerce(#1)$Fraction(S), p)$MPCSF
        else error "Type parameter F should be either equal to S or equal _
                    to Fraction S"

      SUPFPOLYF2SUPF: (SUP FPOLYF, F, F) -> SUP F
      SUPFPOLYF2SUPF(p, a1v, Av) ==
        zero? p => 0
        lc: FPOLYF := leadingCoefficient(p)
        num: POLYF := eval(numer lc, [index(1)$V, index(2)$V]::List V, [a1v, 
Av])
        den: POLYF := eval(denom lc, [index(1)$V, index(2)$V]::List V, [a1v, 
Av])
        monomial(retract(num)/retract(den), degree p) 
          + SUPFPOLYF2SUPF(reductum p, a1v, Av)

      SUPPOLYF2SUPF: (SUP POLYF, F, F) -> SUP F
      SUPPOLYF2SUPF(p, a1v, Av) ==
        zero? p => 0
        lc: POLYF := leadingCoefficient(p)
        monomial(retract(eval(lc, [index(1)$V, index(2)$V]::List V, 
                                  [a1v, Av])),
                 degree p) 
          + SUPPOLYF2SUPF(reductum p, a1v, Av)

      SUPFPOLYF2FSUPPOLYF: (SUP FPOLYF) -> Fraction SUP POLYF
      SUPFPOLYF2FSUPPOLYF p ==
        cden := splitDenominator(p)
               $UnivariatePolynomialCommonDenominator(POLYF, FPOLYF,SUP FPOLYF)
        
        pnum: SUP POLYF 
             := map(retract(#1 * cden.den)$FPOLYF, p)
                   $SparseUnivariatePolynomialFunctions2(FPOLYF, POLYF)
        pden: SUP POLYF := (cden.den)::SUP POLYF

        pnum/pden


      POLYF2EXPRR: POLYF -> EXPRR
      POLYF2EXPRR p ==
        map(convert(#1)@Symbol::EXPRR, coerce(#1)@EXPRR, p)
           $PolynomialCategoryLifting(IndexedExponents V, V, 
                                      F, POLYF, EXPRR)

      checkResult: (EXPRR, Symbol, Integer, List F) -> PI
      checkResult (res, n, l, list) ==
        for i in l..1 by -1 repeat
          den := eval(denominator res, n::EXPRR, i::EXPRR)
          if den = 0 then return (i+1)::PI
          num := eval(numerator res, n::EXPRR, i::EXPRR)
          if list.i ~= retract(retract(num/den)@R)
          then return (i+1)::PI
        1

      PCD ==> PolynomialCommonDenominator(S, F, POLYF, IndexedExponents V, V)

      cden: POLYF -> POLYS
      cden p == 
        if F is S then 
          p pretend POLYS
        else if F is Fraction S then
          map(retract(#1)$Fraction(S), clearDenominator(p)$PCD)$MPCFS
        else error "Type parameter F should be either equal to S or equal _
                     to Fraction S"

-------------------------------------------------------------------------------
-- guessPade
-------------------------------------------------------------------------------

      UTSF ==> UnivariateTaylorSeries
      UP ==> UnivariatePolynomial

      list2poly: (List F, NNI) -> SUPF
      list2poly(list, order) ==
        if null(list) then 0
        else monomial(first(list), order)$SUPF + list2poly(rest(list), order+1)

      opcoeff := operator("coefficient"::Symbol)$CommonOperators

      guessPade(xx, list, basis) ==
-- basis is ignored currently

        x := new(x)$Symbol
        len := (#list-2)
        if len < 1 then return []
        c := 0$F
-- turn the list into a series
        s := coerce(unmakeSUP(list2poly(list, 0))::UP(x, F))$UTSF(F, x, c)

        res: List Fraction UP(x, F) := []

-- try pade for all possible degrees of numerator and denominator
        for i in 0..len repeat
          r := pade(i, (len-i)::NNI, s)$PadeApproximantPackage(F, x, c)

          if not (r case "failed") then
-- make r into a series
            rs := coerce(numer(r))$UTSF(F, x, c) 
                  * (coerce(denom(r))$UTSF(F, x, c)**((-1)::Fraction Integer))
                    $UTSF(F, x, c)

-- test the last coefficient and collect, if ok
            if coefficient(rs, (len+1)::NNI)$UTSF(F, x, c) = last(list)
               and not member?(r, res) then
              res := cons(r, res)
            
        map([kernel(opcoeff, 
                    [FSUPF2EXPRR(x, makeSUP(numer #1)/makeSUP(denom #1)), 
                     x::EXPRR, xx::EXPRR]), 1], 
            res)$ListFunctions2(Fraction UP(x, F),
                                Record(function: EXPRR, order: PI))


-------------------------------------------------------------------------------
-- guessPRec
-------------------------------------------------------------------------------

      oprecur := operator("rootof"::Symbol)$CommonOperators

      guessPRec(xx, list, basis) ==
        len := #list
        a := new('x)$Symbol 
        op := operator a
        
        for o in 1..(len-3) quo 3 repeat
-- versuche verschiedene Ordnungen
          d := ((len-2*o-2) quo (o+1))

          if (d >= 0) and (len >= (o+1)*(d+2)) then 
            m: List List F := [reduce(append, [[basis2(basis, 
(k+1))**i*list.(j+k) 
                                                for i in 0..d]
                                               for j in 1..(o+1)])
                               for k in 0..(o+1)*(d+1)]

            resspace: List Vector F := nullSpace((matrix m)::Matrix F)

-- nimm nur das erste Ergebnis
-- wenn der Raum mehrdimensional ist, kann der Nullvector nicht vorkommen...
            if (#resspace > 1) or (any?(not zero?(#1), resspace.1)) then

              reslist:List List SUPF := [[monomial((resspace.1).(1+i+j*(d+1)), 
i) 
                                          for i in 0..d]
                                         for j in 0..o]

              res: List SUPF := [reduce(_+, reslist.j) for j in 1..(o+1)]

              ex: List EXPRR := [eval(SUPF2EXPRR(xx, res.j),
                                      kernel(xx),
                                      basis(xx::EXPRI))*
                                  op(xx::EXPRR+(j-1)::EXPRR)::EXPRR 
                                 for j in 1..(o+1)]

-- das Ergebnis sollte noch getestet werden, auch auf Redundanz!
 
              return [[kernel(oprecur, [reduce(_+, ex)]), 1]]
        []


-------------------------------------------------------------------------------
-- guessRat
-------------------------------------------------------------------------------
-- this function applies rationalInterpolation to the list of data points
-- [(list.i,basis.i) for i in 1..#list-1]. The last data point is used to test
-- for plausibility.

      guessRat(xx, list, basis) ==
        len := (#list-1)::PositiveInteger
        xlist := [basis2(basis, i) for i in 1..len]
        x:F := basis2(basis, len+1)
        ylist: List F := first(list, len)
        y:F := last(list)
        res: List FSUPF := []

        for i in 0..(len-1) repeat
          output(hconcat("degree Rat "::OutputForm, 
i::OutputForm))$OutputPackage
          ri := interpolate(xlist, ylist, (len-1-i)::NNI, i)
                           $RationalInterpolation(xx, F)
          de: F := elt(denom ri, x)
          if (de ~= 0) and (elt(numer ri, x) = y*de) _
             and not any?(ri = #1, res) then
            res := cons(ri, res)

        [[res1 := eval(FSUPF2EXPRR(xx, guess), kernel(xx), basis(xx::EXPRI)), _
          checkResult(res1, xx, len, list)] _
         for guess in res]

-------------------------------------------------------------------------------
-- guessExpRat
-------------------------------------------------------------------------------

      GF ==> GeneralizedMultivariateFactorize(SingletonAsOrderedSet,
                                              IndexedExponents V, F, F,
                                              SparseUnivariatePolynomial F)

      resl: (POLYS, POLYS, Integer, Integer, V) -> List S
      resl(p1, p2, o, d, A) == 
        [(resultant(univariate(eval(p1, A, k::S)),
                    univariate(eval(p2, A, k::S)))$SUP(S) exquo _
         ((k::S)**(o::NNI)))::S for k in 1..d]

      p: (Integer, Integer, V, V, BASIS) -> POLYF 
      p(len, i, a1, A, basis) == 
--     a0 + a1*basis2(basis, i)
-- setting A=a0+basis(len+3)*a1, hence a0=A-(len+3)*a1 makes poly3 a lot
-- smaller
        A::POLYF + a1::POLYF*(basis2(basis, i)-basis2(basis, len+3))

      p2: (Integer, Symbol, F, F, BASIS) -> EXPRR 
      p2(len, i, a1v, Av, basis) == 
--     a0 + a1*basis2(basis, i)
-- setting A=a0+basis(len+3)*a1, hence a0=A-(len+3)*a1 makes poly3 a lot smaller
        coerce(Av) + coerce(a1v)*(basis(i::EXPRI)
                                 -basis((len+3)::EXPRI))

      ord1: Integer -> Integer
      ord1 n == (3*n**4-4*n**3+3*n*n-2*n) quo 12

      ord2: (Integer, Integer) -> Integer
      ord2(n, i) == 
        if i=0 
        then ord1 n + ((n*n-n) quo 2)
        else ord1 n
         
      deg1: (Integer, Integer) -> Integer
      deg1(n, i) == (3*n**4+12*i*n**3+2*n**3+12*i*i*n*n-6*i*n*n+15*n*n
                    -24*i*i*n+6*i*n-8*n+12*i*i-12*i-12) quo 12

      deg2: (Integer, Integer) -> Integer
      deg2(n,i) == deg1(n,i) + (n**2+2*i*n+n-2*i+2) quo 2

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

      guessExpRatAux: (Symbol, List F, BASIS, List Integer) -> List EXPRR
      guessExpRatAux(xx, list, basis, xValues) ==
        a1: V := index(1)$V
        A: V := index(2)$V

        len:NNI := (#list-3)::NNI
        if len < 1 then return []
        xlist := [basis2(basis, xValues.i)::POLYF::FPOLYF for i in 1..len]
        x1 := basis2(basis, xValues.(len+1))::POLYF::FPOLYF
        x2 := basis2(basis, xValues.(len+2))::POLYF::FPOLYF
        x3 := basis2(basis, xValues.(len+3))::POLYF::FPOLYF

        y: NNI -> FPOLYF := 
          (list.#1)::F::POLYF::FPOLYF * p(len, (xValues.#1)::Integer, a1, A, 
basis)::FPOLYF**(-(xValues.#1)::Integer)

        ylist: List FPOLYF := [y i for i in 1..len]

        y1 := y(len+1)
        y2 := y(len+2)
        y3 := y(len+3)

        res := []::List EXPRR
        for i in 0..len-1 repeat
          output(hconcat("degree ExpRat "::OutputForm, 
i::OutputForm))$OutputPackage
          ri: FSUPFPOLYF
             := interpolate(xlist, ylist, (len-1-i)::NNI, i)
                           $RationalInterpolation(xx, FPOLYF)

          poly1:POLYS := cden(numer(elt(ri, x1)$SUP(FPOLYF) - y1)$FPOLYF)
          poly2:POLYS := cden(numer(elt(ri, x2)$SUP(FPOLYF) - y2)$FPOLYF)
          poly3:POLYS := cden(numer(elt(ri, x3)$SUP(FPOLYF) - y3)$FPOLYF)

          n:Integer := len - i + 1
          if n > 5 then
            output("interpolating"::OutputForm)$OutputPackage
            o1:Integer := ord1(n)
            d1:Integer := deg1(n, i)
            o2:Integer := ord2(n, i)
            d2:Integer := deg2(n, i)

-- RINTERP seems to be a lot faster than PINTERP
-- another compiler bug: using i as iterator here makes the loop break

            res1: SUP S
                 := retract(interpolate([j::S for j in 1..d1-o1+1], 
                                        resl(poly1, poly3, o1, d1-o1+1, A), 
                                        (d1-o1)::NNI, 0)
                                       $RationalInterpolation(xx, S)
                            ::Fraction(SUP(S)))$Fraction(SUP(S))

            res2: SUP S
                 := retract(interpolate([j::S for j in 1..d2-o2+1], 
                                        resl(poly2, poly3, o2, d2-o2+1, A), 
                                        (d2-o2)::NNI, 0)
                                       $RationalInterpolation(xx, S)
                            ::Fraction(SUP(S)))$Fraction(SUP(S))
          else
            output("resultants"::OutputForm)$OutputPackage
            res1: SUP S := univariate(resultant(poly1, poly3, a1))
            res2: SUP S := univariate(resultant(poly2, poly3, a1))

-- we want to solve over F            
          res3: SUP F := SUPS2SUPF(primitivePart(gcd(res1, res2)))

-- res3 ist ein Polynom in A=a0+(len+3)*a1
-- jetzt muss ich das loesen

          res4 := factor(res3)$GF
          for f in factors res4 repeat
            if degree f.factor = 1 then
-- we are only interested in the linear factors
              Av: F := -coefficient(f.factor, 0)/leadingCoefficient f.factor
              if Av ~= 0 then
                res5 := factor(univariate(eval(POLYS2POLYF poly3, A, Av)))$GF
                for g in factors res5 repeat
                  if degree g.factor = 1 then
                    a1v: F := -coefficient(g.factor, 0)
                              /leadingCoefficient g.factor

                    t1:= eval(POLYS2POLYF poly1, [a1, A]::List V, 
                                                 [a1v, Av]::List F)
                    if t1 = 0 then  
                      t2:= eval(POLYS2POLYF poly2, [a1, A]::List V, 
                                                   [a1v, Av]::List F)
                      if t2 = 0 then

                        ri1: Fraction SUP POLYF 
                            := SUPFPOLYF2FSUPPOLYF(numer ri)
                             / SUPFPOLYF2FSUPPOLYF(denom ri)

                        num: SUP F := SUPPOLYF2SUPF(numer ri1, a1v, Av)
                        den: SUP F := SUPPOLYF2SUPF(denom ri1, a1v, Av)

                        if den ~= 0 then
                          res7: EXPRR := eval(FSUPF2EXPRR(xx, num/den),
                                              kernel(xx), 
                                              basis(xx::EXPRI))
                                        *p2(len, xx, a1v, Av, basis)**xx::EXPRR
                          res := cons(res7, res)
                        else if num = 0 then
                          output("numerator and denominator vanish!")
                                $OutputPackage
  
          if not null(res) then return res

        res

      guessExpRat(xx, list, basis) ==
-- guesses Functions of the Form (a1*n+a0)^n*rat(n)

-- remove zeros from list

        zeros: EXPRR := 1
        newlist: List F
        xValues: List Integer
        i: Integer

        i := 0
        for x in list repeat
          i := i+1
          if x = 0 then zeros := zeros * (basis(xx::EXPRI) - basis(i::EXPRI))

        i := 0
        for x in list repeat
          i := i+1
          if x ~= 0 then
            newlist := cons(x/retract(retract(eval(zeros, xx::EXPRR, 
                                                          i::EXPRR))@R),
                            newlist)
            xValues := cons(i, xValues)

        newlist := reverse newlist
        xValues := reverse xValues
        len := i

        res: List EXPRR 
            := [zeros * f for f in guessExpRatAux(xx, newlist, basis, xValues)]
        
        map([#1, checkResult(#1, xx, len, list)], res)
           $ListFunctions2(EXPRR,
                           Record(function: EXPRR, order: PI))


-------------------------------------------------------------------------------
-- guess
-------------------------------------------------------------------------------

      guess(xx, list, basis, guessers, ops, maxlevel) ==
        output(hconcat("guessing level "::OutputForm, maxlevel::OutputForm))
             $OutputPackage
        res: List Record(function: EXPRR, order: PI) := []
        len:= #list :: PositiveInteger
        if len > 1 then

          res := []
          for guesser in guessers repeat
            res := append(guesser(xx, list, basis), res)

            if member?('guessOne, ops) and not null(res) then return res

          if (maxlevel <= 0) then return res

          if member?('guessProduct, ops) and not member?(0$F, list) then
            prodList: List F := [(list.(i+1))/(list.i) for i in 1..(len-1)]

            if not every?(one?, prodList) then
              var: Symbol := subscript('p, [len::OutputForm])
              prodGuess := 
                  [[coerce(list.(guess.order)) 
                    * product(guess.function, _
                              equation(var, _
                                       (guess.order)::EXPRR..xx::EXPRR-1)), _
                    guess.order] _
                   for guess in guess(var, prodList, basis, guessers, 
                                      ops, (maxlevel-1)::NNI)$%]

              for guess in prodGuess 
                  | not any?(guess.function = #1.function, res) repeat
                res := cons(guess, res)

          if member?('guessSum, ops) then
            sumList: List F := [(list.(i+1))-(list.i) for i in 1..(len-1)]
  
            if not every?(#1=sumList.1, sumList) then
              var: Symbol := subscript('s, [maxlevel::OutputForm])
              sumGuess := 
                  [[coerce(list.(guess.order)) _
                    + summation(guess.function, _
                                equation(var, _
                                         (guess.order)::EXPRR..xx::EXPRR-1)), _
                    guess.order] _
                   for guess in guess(var, sumList, basis, guessers, 
                                      ops, (maxlevel-1)::NNI)$%]

              for guess in sumGuess 
                  | not any?(guess.function = #1.function, res) repeat
                res := cons(guess, res)
  
        res

)abbrev package GUESSINT GuessInteger
++ Description:
++ This package exports guessing of sequences of rational numbers
GuessInteger() == Guess(Fraction Integer, Integer, Expression Integer, 
                        Fraction Integer,
                        id$MappingPackage1(Fraction Integer), 
                        coerce$Expression(Integer))

)abbrev package GUESSP GuessPolynomial
++ Description:
++ This package exports guessing of sequences of rational functions
GuessPolynomial() == Guess(Fraction Polynomial Integer, Polynomial Integer, 
                           Expression Integer, 
                           Fraction Polynomial Integer,
                           id$MappingPackage1(Fraction Polynomial Integer), 
                           coerce$Expression(Integer))
\end{axiom}
--
forwarded from http://page.axiom-developer.org/zope/mathaction/address@hidden