next version...
I don't particularly like the design, improvements welcome. Documentation
non-existing currently, too.
Usage example:
(33) -> cat2 := coerce(univariate(x::UP(x, FRAC INT)*y^2-y+1,y)::SUP UP(x, FRAC
INT), 0)$UALG FRAC INT
2
(33) [x ? - ? + 1,0]
Type: UnivariateAlgebraicSeries(Fraction(Integer))
Time: 0.14 (IN) + 0.004 (EV) + 0.05 (OT) = 0.19 sec
(34) -> cat2 * cat2
3 3 2 2
(34) [x ? + (x - x)? + (- x + 1)? - 1,0]
Type: UnivariateAlgebraicSeries(Fraction(Integer))
Time: 0.01 (EV) + 0.02 (OT) = 0.03 sec
check:
(71) -> cs := series([1/(2*n+1)*binomial(2*n+1, n) for n in 0..])$UFPS FRAC INT
(71)
2 3 4 5 6 7 8 9
1 + x + 2x + 5x + 14x + 42x + 132x + 429x + 1430x + 4862x
+
10 11
16796x + O(x )
Type: UnivariateFormalPowerSeries(Fraction(Integer))
Time: 0.02 (IN) + 0.03 (OT) = 0.05 sec
(72) -> guessAlg(entries complete first(coefficients(cs2*cs2), 42),
maxPower==6)
(72)
[
[
function =
n 2 2
[[x ]f(x): x f(x) + (2x - 1)f(x) + 1= 0,
2 3 4 5
f(x)= 1 + 2x + 5x + 14x + 42x + O(x )]
,
order= 0]
]
Type: List(Record(function: Expression(Integer),order: NonNegativeInteger))
Time: 0.08 (EV) + 0.01 (OT) = 0.08 sec
Martin
\documentclass{article}
\usepackage{axiom,amsthm,amsmath,url}
\newtheorem{ToDo}{ToDo}[section]
\newcommand{\Axiom}{{\tt Axiom}}
\newcommand{\Rate}{{\tt Rate}}
\newcommand{\GFUN}{{\tt GFUN}}
\begin{document}
\title{holonomic.spad}
\author{Martin Rubey}
\maketitle
\begin{abstract}
\end{abstract}
\tableofcontents
\section{domain UALG UnivariateAlgebraicSeries}
<<domain UALG UnivariateAlgebraicSeries>>=
)abbrev domain UALG UnivariateAlgebraicSeries
UnivariateAlgebraicSeries(Coef: Algebra Fraction Integer):
Exports == Implementation where
UFPS ==> UnivariateFormalPowerSeries Coef
EQCOEF ==> UnivariatePolynomial(x, Coef)
ALGEQ ==> SparseUnivariatePolynomial EQCOEF
Exports ==> Join(AbelianMonoidRing(Coef, NonNegativeInteger),
DifferentialRing, Field) with
coerce: (ALGEQ, UFPS) -> %
Implementation ==> add
Rep := Record(eq: ALGEQ, init: UFPS)
coerce(f: %): OutputForm ==
bracket([f.eq::OutputForm,
f.init::OutputForm])$OutputForm
coerce(feq: ALGEQ, finit: UFPS): % == [feq, finit]$Rep
0 == [monomial(1,1)$ALGEQ, 0]$Rep
1 == [monomial(1,1)$ALGEQ-1, 1]$Rep
monomial(c, n) ==
[monomial(1,1)$ALGEQ-monomial(c,n)$EQCOEF::ALGEQ,
monomial(c,n)]$Rep
FEQCOEF ==> Fraction EQCOEF
OVAR ==> OrderedVariableList ['u, 'v, 'y]
SMP ==> SparseMultivariatePolynomial(FEQCOEF, OVAR)
GROE ==> GroebnerPackage(FEQCOEF, IndexedExponents OVAR, OVAR, SMP)
CDEN ==> CommonDenominator(EQCOEF,FEQCOEF,List FEQCOEF)
SUP2EF ==> SparseUnivariatePolynomialFunctions2(EQCOEF, FEQCOEF)
SUP2FE ==> SparseUnivariatePolynomialFunctions2(FEQCOEF, EQCOEF)
f + g ==
ff: SparseUnivariatePolynomial FEQCOEF
:= map(#1::FEQCOEF, f.eq)$SUP2EF
gf: SparseUnivariatePolynomial FEQCOEF
:= map(#1::FEQCOEF, g.eq)$SUP2EF
leq: List SMP :=
[index(3)$OVAR::SMP-index(1)$OVAR::SMP-index(2)$OVAR::SMP, _
multivariate(ff, index(1)$OVAR), _
multivariate(gf, index(2)$OVAR)]
neweq: SMP := third groebner(leq)$GROE
neweq2: SparseUnivariatePolynomial FEQCOEF := univariate(neweq)$SMP
cden: EQCOEF := commonDenominator(coefficients neweq2)$CDEN
algeq: ALGEQ := map(retract(cden*#1)@EQCOEF, neweq2)$SUP2FE
[algeq, f.init + g.init]$Rep
(f: %) * (g: %) ==
ff: SparseUnivariatePolynomial FEQCOEF
:= map(#1::FEQCOEF, f.eq)$SUP2EF
gf: SparseUnivariatePolynomial FEQCOEF
:= map(#1::FEQCOEF, g.eq)$SUP2EF
leq: List SMP :=
[index(3)$OVAR::SMP-index(1)$OVAR::SMP*index(2)$OVAR::SMP, _
multivariate(ff, index(1)$OVAR), _
multivariate(gf, index(2)$OVAR)]
gb := groebner(leq)$GROE
output(gb::OutputForm)$OutputPackage
neweq: SMP := last gb
output(variables(neweq)$SMP::OutputForm)$OutputPackage
neweq2: SparseUnivariatePolynomial FEQCOEF := univariate(neweq)$SMP
cden: EQCOEF := commonDenominator(coefficients neweq2)$CDEN
algeq: ALGEQ := map(retract(cden*#1)@EQCOEF, neweq2)$SUP2FE
[algeq, f.init * g.init]$Rep
@
\section{domain UHOLO1 UnivariateHolonomicSeries1}
<<domain UHOLO1 UnivariateHolonomicSeries1>>=
)abbrev domain UHOLO1 UnivariateHolonomicSeries1
UnivariateHolonomicSeries1(Coef: Algebra Fraction Integer):
Exports == Implementation where
UFPS ==> UnivariateFormalPowerSeries Coef
DEQCOEF ==> UnivariatePolynomial(x, Coef)
FDEQCOEF ==> Fraction DEQCOEF
LODO1 ==> LinearOrdinaryDifferentialOperator1 FDEQCOEF
Exports ==> Join(AbelianMonoidRing(Coef, NonNegativeInteger),
DifferentialRing) with
sin: () -> %
cos: () -> %
atan: () -> %
exp: () -> %
integrate: % -> %
series: % -> UFPS
differentialEquation: % -> LODO1
Implementation ==> add
Rep := Record(eq : LODO1, init: UFPS)
series f == f.init
differentialEquation f == f.eq
coerce(f: %): OutputForm ==
bracket([f.eq::OutputForm,
f.init::OutputForm])$OutputForm
Dx := D()$LODO1
0 == [1, 0]$Rep
1 == [Dx, 1]$Rep
-- x*Dx x^n = x*n*x^(n-1)=n*x^n
monomial(c, n) ==
c1: FDEQCOEF := monomial(1,1)$DEQCOEF :: FDEQCOEF
c0: FDEQCOEF := (-n)::DEQCOEF :: FDEQCOEF
[c1*Dx+c0*1$LODO1, monomial(c,n)]$Rep
f + g == [directSum(f.eq, g.eq), f.init + g.init]
(n: Integer) * (f: %) == [f.eq, n*f.init]
(f: %) * (g: %) == [symmetricProduct(f.eq, g.eq), f.init * g.init]$Rep
integrate f == [f.eq*Dx, integrate f.init]
D f ==
p0 := coefficient(f.eq, 0)
Deq: LODO1
if zero? p0
then Deq := f.eq
else
-- without the parenthesis below, we would use ordinary multiplication!
Deq := Dx * (inv p0 * f.eq)
lmo: List LODO1
:= [monomial(coefficient(Deq, i+1), i)$LODO1 _
for i in 0..degree Deq - 1]
[reduce(_+$LODO1, lmo, 0), D f.init]
sin() == [Dx**2+1, sin(monomial(1,1)$UFPS)$UFPS]$Rep
cos() == [Dx**2+1, cos(monomial(1,1)$UFPS)$UFPS]$Rep
exp() == [Dx-1, exp(monomial(1,1)$UFPS)$UFPS]$Rep
atan() ==
c2: FDEQCOEF := (monomial(1,2)$DEQCOEF+1)::FDEQCOEF
c1: FDEQCOEF := monomial(2*1$Coef,1)$DEQCOEF::FDEQCOEF
[c2*Dx**2+c1*Dx, atan(monomial(1,1)$UFPS)$UFPS]$Rep
@
%$
<<*>>=
<<domain UALG UnivariateAlgebraicSeries>>
<<domain UHOLO1 UnivariateHolonomicSeries1>>
@
\end{document}
--~--~---------~--~----~------------~-------~--~----~
You received this message because you are subscribed to the Google Groups
"FriCAS - computer algebra system" group.
To post to this group, send email to [email protected]
To unsubscribe from this group, send email to
[email protected]
For more options, visit this group at
http://groups.google.com/group/fricas-devel?hl=en
-~----------~----~----~----~------~----~------~--~---
