polynomials with negative integer exponents
perrinmeyer
I'm new to axiom, which so far seems great.
How do I get axiom to display polynomials with negative integer
exponents, like
--> a + b * z^-1 + c * z^-2
(DSP style). See below for a transcript.
(I am using axiom on my Intel macbook 10.5 (I got sbcl and texmacs
from macports, and I compiled friCAS using sbcl. I made a soft link
between /usr/local/bin/axiom and AXIOMsys, to get texmacs to work)
Thanks,
Perrin Meyer
→
f := a*z^2 + b*z + c
→
a z2 + b z + c
(1)
Type: Polynomial Integer
→
f / z^2
→
a z2 + b z + c
z2
(2)
Type: Fraction Polynomial Integer
→
% :: UP(z,?)
→
a
z2 z 2 +
b
z2 z +
c
z2
(3)
Type: UnivariatePolynomial(z,Fraction Polynomial Integer)
→
Martin Rubey
> I'm new to axiom, which so far seems great.
Welcome!
> How do I get axiom to display polynomials with negative integer
> exponents, like
>
> --> a + b * z^-1 + c * z^-2
You want to use the LaurentPolynomial domain. Unfortunately, it is a little
braindead, it seems.
(1) -> L ==> LAUPOL(POLY INT, UP(z, POLY INT))
Type: Void
(2) -> (a + b * z^-1 + c * z^-2)::L
Cannot convert from type Fraction Polynomial Integer to
LaurentPolynomial(Polynomial Integer,UnivariatePolynomial(z,
Polynomial Integer)) for value
2
a z + b z + c
--------------
2
z
(2) -> a + b * monomial(1,-1)$L + c * monomial(1,-2)$L
- 1 - 2
(2) a + b z + c z
Type: LaurentPolynomial(Polynomial Integer,UnivariatePolynomial(z,Polynomial Integer))
I think I'd like to see the following things changed:
1) LaurentPolynomial should take parameters just as any other polynomial
domain, or, maybe better, there should be a frontend
UnivariateLaurentPolynomial(Symbol, Ring)
2) I'd like to have support for multivariate Laurent polynomials.
3) a ^ b and a ** b should work as shorthand for monomial(a, b). Actually,
that should be the case in general, shouldn't it?
I think (3) would be very worthwhile to check.
Martin
Ralf Hemmecke
>> How do I get axiom to display polynomials with negative integer
>> exponents, like
>>
>> --> a + b * z^-1 + c * z^-2
> You want to use the LaurentPolynomial domain. Unfortunately, it is a little
> braindead, it seems.
> (1) -> L ==> LAUPOL(POLY INT, UP(z, POLY INT))
Don't give your opinion too easily. Look at gpol.spad.pamphlet why there
are two parameters.
LaurentPolynomial(R, UP): Exports == Implementation where
R : IntegralDomain
UP: UnivariatePolynomialCategory R
...
Two parameters give a little more flexibility. In fact, *this*
LaurentPolynomial constructor is some kind of meta constructor. Namely,
it allows to construct sparse and dense Laurent polynomials (just give
SUP or DUP as its second parameter.
But in some sense you are right. Nobody is used to such complicated
stuff, so introducing something like
SparseUnivariateLaurentPolynomial(z: Symbol, R: IntegralDomain): ... ==
LaurentPolynomial(R, SUP(z, R)
seems like a natural "abbreviation".
> (2) -> (a + b * z^-1 + c * z^-2)::L
>
> Cannot convert from type Fraction Polynomial Integer to
> LaurentPolynomial(Polynomial Integer,UnivariatePolynomial(z,
> Polynomial Integer)) for value
> 2
> a z + b z + c
> --------------
> 2
> z
>
> (2) -> a + b * monomial(1,-1)$L + c * monomial(1,-2)$L
>
> - 1 - 2
> (2) a + b z + c z
> Type: LaurentPolynomial(Polynomial Integer,UnivariatePolynomial(z,Polynomial Integer))
Maybe "POLY INT" is getting in the way again?
(1) -> Q:=Fraction Integer
(1) Fraction Integer
Type:
Domain
(2) -> P:=DMP([a,b,c],Q)
(2) DistributedMultivariatePolynomial([a,b,c],Fraction Integer)
Type:
Domain
(3) -> L:=LAUPOL(P, UP(z,P))
(3)
LaurentPolynomial(DistributedMultivariatePolynomial([a,b,c],Fraction
Integer)
,UnivariatePolynomial(z,DistributedMultivariatePolynomial([a,b,c],Fraction
In
teger)))
Type:
Domain
(4) -> x:L := z::L
(4) z
Type:
LaurentPolynomial(DistributedMultivariatePolynomial([a,b,c],Fraction
Integer),UnivariatePolynomial(z,DistributedMultivariatePolynomial([a,b,c],Fraction
Integer)))
(5) -> l: L := a + b * x^(-1) + c * x^(-2)
- 1 - 2
(5) a + b z + c z
Type:
LaurentPolynomial(DistributedMultivariatePolynomial([a,b,c],Fraction
Integer),UnivariatePolynomial(z,DistributedMultivariatePolynomial([a,b,c],Fraction
Integer)))
> 2) I'd like to have support for multivariate Laurent polynomials.
In a recursive fashion that should be easy. Distributed ones are
probably also easy if one just takes the code of GDMP.
> 3) a ^ b and a ** b should work as shorthand for monomial(a, b). Actually,
> that should be the case in general, shouldn't it?
But that works (see above). It is only that POLY(INT) also contains z
and so there are actually two ways to embed
a + b * z^-1 + c * z^-2
into *your* L. Note that thing could just be considered to be a constant.
I hope some day people understand why POLY(INT) is just bad in some
situations (like this).
Ralf
Martin Rubey
> On 08/11/2008 09:39 AM, Martin Rubey wrote:
> >> How do I get axiom to display polynomials with negative integer
> >> exponents, like
> >>
> >> --> a + b * z^-1 + c * z^-2
>
> > You want to use the LaurentPolynomial domain. Unfortunately, it is a little
> > braindead, it seems.
>
> > (1) -> L ==> LAUPOL(POLY INT, UP(z, POLY INT))
>
> Don't give your opinion too easily.
I didn't in this particular case - you deleted something:
or, maybe better, there should be a frontend
UnivariateLaurentPolynomial(Symbol, Ring)
> Maybe "POLY INT" is getting in the way again?
no.
> (1) -> Q:=Fraction Integer
>
> (1) Fraction Integer
> Type:
> Domain
> (2) -> P:=DMP([a,b,c],Q)
>
> (2) DistributedMultivariatePolynomial([a,b,c],Fraction Integer)
> Type:
> Domain
> (3) -> L:=LAUPOL(P, UP(z,P))
>
> (3)
> LaurentPolynomial(DistributedMultivariatePolynomial([a,b,c],Fraction
> Integer)
>
> ,UnivariatePolynomial(z,DistributedMultivariatePolynomial([a,b,c],Fraction
> In
> teger)))
> Type:
> Domain
> (4) -> x:L := z::L
this is cheating. What I want is that
(1) -> L ==> LAUPOL(INT, UP(z, INT))
Type: Void
(2) -> (1 + 2 * z^-1 + 3 * z^-2)::L
Cannot convert from type Fraction Polynomial Integer to
LaurentPolynomial(Integer,UnivariatePolynomial(z,Integer)) for
value
2
z + 2z + 3
-----------
2
z
"works".
> > 3) a ^ b and a ** b should work as shorthand for monomial(a, b). Actually,
> > that should be the case in general, shouldn't it?
>
> But that works (see above). It is only that POLY(INT) also contains z
> and so there are actually two ways to embed
>
> a + b * z^-1 + c * z^-2
>
> into *your* L. Note that thing could just be considered to be a constant.
No it doesn't work:
(2) -> (z^-2)$L
There are 4 exposed and 0 unexposed library operations named ^
having 2 argument(s) but none was determined to be applicable.
Use HyperDoc Browse, or issue
)display op ^
to learn more about the available operations. Perhaps
package-calling the operation or using coercions on the arguments
will allow you to apply the operation.
Cannot find a definition or applicable library operation named ^
with argument type(s)
Symbol
Integer
Perhaps you should use "@" to indicate the required return type,
or "$" to specify which version of the function you need.
> I hope some day people understand why POLY(INT) is just bad in some
> situations (like this).
Don't worry, I understand the problematic points of POLY INT fairly well. They
are not involved here. (Although they could be)
Martin
Ralf Hemmecke
> Ralf Hemmecke <ra...@hemmecke.de> writes:
>
>> On 08/11/2008 09:39 AM, Martin Rubey wrote:
>>>> How do I get axiom to display polynomials with negative integer
>>>> exponents, like
>>>>
>>>> --> a + b * z^-1 + c * z^-2
>>> You want to use the LaurentPolynomial domain. Unfortunately, it is a little
>>> braindead, it seems.
>>> (1) -> L ==> LAUPOL(POLY INT, UP(z, POLY INT))
>> Don't give your opinion too easily.
>
> I didn't in this particular case - you deleted something:
>
> or, maybe better, there should be a frontend
>
> UnivariateLaurentPolynomial(Symbol, Ring)
Yes. Sorry, it only came at the end of your mail.
>> Maybe "POLY INT" is getting in the way again?
>
> no.
>
>> (1) -> Q:=Fraction Integer
>>
>> (1) Fraction Integer
>> Type:
>> Domain
>> (2) -> P:=DMP([a,b,c],Q)
>>
>> (2) DistributedMultivariatePolynomial([a,b,c],Fraction Integer)
>> Type:
>> Domain
>> (3) -> L:=LAUPOL(P, UP(z,P))
>>
>> (3)
>> LaurentPolynomial(DistributedMultivariatePolynomial([a,b,c],Fraction
>> Integer)
>>
>> ,UnivariatePolynomial(z,DistributedMultivariatePolynomial([a,b,c],Fraction
>> In
>> teger)))
>> Type:
>> Domain
>> (4) -> x:L := z::L
>
> this is cheating.
Really? I just give the interpreter some hint. If you simply type z then
that can be anything and it depends on the interpreter to do exactly as
you wish. But of course it already happens with people that they don't
understand each other. Why should a computer program be more clever in
general? I was just "friendly" to the interpreter. (Some people call
that cheating. *smile*)
> What I want is that
>
> (1) -> L ==> LAUPOL(INT, UP(z, INT))
> Type: Void
> (2) -> (1 + 2 * z^-1 + 3 * z^-2)::L
>
> Cannot convert from type Fraction Polynomial Integer to
> LaurentPolynomial(Integer,UnivariatePolynomial(z,Integer)) for
> value
> 2
> z + 2z + 3
> -----------
> 2
> z
>
> "works".
But how should it work? Once the interpreter has wrongly constructed an
element of type Fraction Polynomial Integer, it is, of course, a bit
hard to go back. Note that you would have to provide a function
Fraction Polynomial Integer -> L
which fails in general. For example,
q 1
--- = ----------
q-1 1 - q^(-1)
is not a LaurentPolynomial.
What I find strange, though, is
(12) -> (z:: L)^(-1)::L
1
(12) -
z
Type: Fraction LaurentPolynomial(Fraction
Integer,UnivariatePolynomial(z,Fraction Integer))
(13) -> ((z::L)^(-1))::L
- 1
(13) z
Type: LaurentPolynomial(Fraction Integer,UnivariatePolynomial(z,Fraction
Integer))
Looks like :: is binding stronger than I expected.
But the following is OK.
(14) -> h:L := (z::L)^(-1)
- 1
(14) z
Type: LaurentPolynomial(Fraction Integer,UnivariatePolynomial(z,Fraction
Integer))
>>> 3) a ^ b and a ** b should work as shorthand for monomial(a, b). Actually,
>>> that should be the case in general, shouldn't it?
>> But that works (see above). It is only that POLY(INT) also contains z
>> and so there are actually two ways to embed
>>
>> a + b * z^-1 + c * z^-2
>>
>> into *your* L. Note that thing could just be considered to be a constant.
>
> No it doesn't work:
>
> (2) -> (z^-2)$L
> There are 4 exposed and 0 unexposed library operations named ^
> having 2 argument(s) but none was determined to be applicable.
> Use HyperDoc Browse, or issue
> )display op ^
> to learn more about the available operations. Perhaps
> package-calling the operation or using coercions on the arguments
> will allow you to apply the operation.
>
> Cannot find a definition or applicable library operation named ^
> with argument type(s)
> Symbol
> Integer
>
> Perhaps you should use "@" to indicate the required return type,
> or "$" to specify which version of the function you need.
Hmm, than it seems that one should add a function
^: (Symbol, Integer) -> %
to LaurentPolynomial(R, UP(z, R))
which always fails for any first argument x except x=z. That's a stupid
function... and moreover it's a partial function.
Ralf
Martin Rubey
> > What I want is that
> >
> > (1) -> L ==> LAUPOL(INT, UP(z, INT))
> > Type: Void
> > (2) -> (1 + 2 * z^-1 + 3 * z^-2)::L
> >
> > Cannot convert from type Fraction Polynomial Integer to
> > LaurentPolynomial(Integer,UnivariatePolynomial(z,Integer)) for
> > value
> > 2
> > z + 2z + 3
> > -----------
> > 2
> > z
> >
> > "works".
>
> But how should it work?
Good question.
> Once the interpreter has wrongly constructed an element of type Fraction
> Polynomial Integer, it is, of course, a bit hard to go back. Note that you
> would have to provide a function
>
> Fraction Polynomial Integer -> L
>
> which fails in general. For example,
>
> q 1
> --- = ----------
> q-1 1 - q^(-1)
>
> is not a LaurentPolynomial.
I guess that would be:
FRAC R has RetractableTo LaurentPolynomial R
when R is some polynomial domain... Then we'd have retract and retractIfCan.
> (13) -> ((z::L)^(-1))::L
>
> - 1
> (13) z
> Type: LaurentPolynomial(Fraction Integer,UnivariatePolynomial(z,Fraction
> Integer))
I believe that the interpreter supplies some magic here. I cannot find a
coerce or convert from FRAC LAUPOL to LAUPOL:
(7) -> p := ((z::L)^(-1))
Function Selection for ^
Arguments: (LAUPOL(INT,UP(z,INT)),INT)
Default target type: Fraction LaurentPolynomial(Integer,UnivariatePolynomial(z,Integer))
[1] signature: (FRAC LAUPOL(INT,UP(z,INT)),INT) -> FRAC LAUPOL(INT,UP(z,INT))
implemented: slot $$(Integer) from FRAC LAUPOL(INT,UP(z,INT))
[2] signature: (FRAC LAUPOL(INT,UP(z,INT)),INT) -> FRAC LAUPOL(INT,UP(z,INT))
implemented: slot $$(Integer) from FRAC LAUPOL(INT,UP(z,INT))
1
(7) -
z
Type: Fraction LaurentPolynomial(Integer,UnivariatePolynomial(z,Integer))
(8) -> p::L
- 1
(8) z
Type: LaurentPolynomial(Integer,UnivariatePolynomial(z,Integer))
> > (2) -> (z^-2)$L
> > There are 4 exposed and 0 unexposed library operations named ^
> > having 2 argument(s) but none was determined to be applicable.
> > Use HyperDoc Browse, or issue
> > )display op ^
> > to learn more about the available operations. Perhaps
> > package-calling the operation or using coercions on the arguments
> > will allow you to apply the operation.
> >
> > Cannot find a definition or applicable library operation named ^
> > with argument type(s)
> > Symbol
> > Integer
> >
> > Perhaps you should use "@" to indicate the required return type,
> > or "$" to specify which version of the function you need.
>
> Hmm, than it seems that one should add a function
>
> ^: (Symbol, Integer) -> %
>
> to LaurentPolynomial(R, UP(z, R))
>
> which always fails for any first argument x except x=z. That's a stupid
> function... and moreover it's a partial function.
Yes, you are right. In UP(z, R) the command z^2 works, because it is simply
using z::UP(z, R) ^ 2.
But the signature doesn't have to be (Symbol, INT) -> %, if we supply the
variable as parameter to LaurentPolynomial as in
)abb package TEST Test
Test(x: Symbol): with
foo: Variable x -> Symbol
== add
foo z == variable()
-------------------------------------------------------------------------------
(1) -> foo(z)
(1) z
Type: Symbol
(2) -> foo(z)$Test(z)
(2) z
Type: Symbol
(3) -> foo(z)$Test(x)
There are 1 exposed and 0 unexposed library operations named foo
having 1 argument(s) but none was determined to be applicable.
Use HyperDoc Browse, or issue
)display op foo
to learn more about the available operations. Perhaps
package-calling the operation or using coercions on the arguments
will allow you to apply the operation.
Cannot find a definition or applicable library operation named foo
with argument type(s)
Variable z
Perhaps you should use "@" to indicate the required return type,
or "$" to specify which version of the function you need.
Martin
Ralf Hemmecke
>> Polynomial Integer, it is, of course, a bit hard to go back. Note that you
>> would have to provide a function
>>
>> Fraction Polynomial Integer -> L
>>
>> which fails in general. For example,
>>
>> q 1
>> --- = ----------
>> q-1 1 - q^(-1)
>>
>> is not a LaurentPolynomial.
>
> I guess that would be:
>
> FRAC R has RetractableTo LaurentPolynomial R
>
> when R is some polynomial domain... Then we'd have retract and retractIfCan.
Can you give me a (terminating) algorithm that finds for two univariate
polynomials a, b whether their quotient is a Laurent series or not, so
that one can implement 'retractIfCan'?
>> (13) -> ((z::L)^(-1))::L
>>
>> - 1
>> (13) z
>> Type: LaurentPolynomial(Fraction Integer,UnivariatePolynomial(z,Fraction
>> Integer))
> I believe that the interpreter supplies some magic here. I cannot find a
> coerce or convert from FRAC LAUPOL to LAUPOL:
Yep, that is one of the reasons why I dislike the interpreter. There is
some magic that I don't easily have access to. Maybe the interpreter
interprets ::L as @L ... no... but
(16) -> retract((z:: L)^(-1))@L
- 1
(16) z
Type: LaurentPolynomial(Fraction Integer,UnivariatePolynomial(z,Fraction
Integer))
Ralf
Martin Rubey
at the very end you'll find a nice way to generate Laurentpolynomials.
Ralf, what do you think about
(Variable x) ** Integer -> %
for LaurentPolynomial?
Ralf Hemmecke <ra...@hemmecke.de> writes:
> >> (13) -> ((z::L)^(-1))::L
> >>
> >> - 1
> >> (13) z
> >> Type: LaurentPolynomial(Fraction Integer,UnivariatePolynomial(z,Fraction
> >> Integer))
>
> > I believe that the interpreter supplies some magic here. I cannot find a
> > coerce or convert from FRAC LAUPOL to LAUPOL:
>
> Yep, that is one of the reasons why I dislike the interpreter. There is some
> magic that I don't easily have access to. Maybe the interpreter interprets
> ::L as @L ... no... but
>
> (16) -> retract((z:: L)^(-1))@L
>
> - 1
> (16) z
> Type: LaurentPolynomial(Fraction Integer,UnivariatePolynomial(z,Fraction
> Integer))
Oh, very nice! retract is implemented in FRAC, which employs exquo from
LAUPOL. I like it.
I'll have to think a little why
retract(z^(-1))@L
doesn't work. Hm,
(z^(-1))::FRAC(L)::L
does. That's at least usable:
(27) -> (1 + 2 * z^-1 + 3 * z^-2)::FRAC L::L
- 1 - 2
(27) 1 + 2z + 3z
Type: LaurentPolynomial(Integer,UnivariatePolynomial(z,Integer))
And in this setting, LAUPOL(POLY INT, UP(z, POLY INT)) is indeed problematic.
However:
(33) -> M ==> LAUPOL(DMP([a,b,c], INT), UP(z, DMP([a,b,c], INT)))
Type: Void
(34) -> (a + b * z^-1 + c * z^-2)::FRAC M::M
- 1 - 2
(34) a + b z + c z
Type: LaurentPolynomial(DistributedMultivariatePolynomial([a,b,c],Integer),UnivariatePolynomial(z,DistributedMultivariatePolynomial([a,b,c],Integer)))
Great work, Ralf!
Martin
Ralf Hemmecke
> Great work, Ralf!
I couldn't have done anything without *your* initial work. So I would
like to share your praise with you. ;-)
Ralf
perrinmeyer
I'm having trouble manipulating Laurent Polynomials (I would think
that multiplying by z^(-2), for example, would be straightforward.) Is
there a Laurent equivalent to UP? Please see transcript below. I am
using texmacs html output, is there a better way for email
discussions?
Thank you very much for your help.
Perrin Meyer
→
f(s) == (b2 * s^2 + b1 * s + b0)/(a2 * s^2 + a1 * s + a0)
Type: Void
→
f((2/t)*( (1 - z^(-1))/(1 + z^(-1))))
Compiling function f with type Fraction Polynomial Integer ->
Fraction Polynomial Integer
(b0t2 + 2b1t + 4b2)z2 + (2b0t2 - 8b2)z + b0t2 - 2b1t + 4b2
(a0t2 + 2a1t + 4a2)z2 + (2a0t2 - 8a2)z + a0t2 - 2a1t + 4a2
Type: Fraction Polynomial Integer
→
(b0t2 + 2b1t + 4b2)z2 + (2b0t2 - 8b2)z + b0t2 - 2b1t + 4b2
(3)
Type: Polynomial Integer
→
M ==> LAUPOL(DMP([a0,a1,a2,b0,b1,b2,t], INT), UP(z,
DMP([a0,a1,a2,b0,b1,b2,t], INT)))
Type: Void
→
fnl := fn :: M
(b0t2 + 2b1t + 4b2)z2 + (2b0t2 - 8b2)z + b0t2 - 2b1t + 4b2
(7)
Type:
LaurentPolynomial(DistributedMultivariatePolynomial([a0,a1,a2,b0,b1,b2,t],Integer),UnivariatePolynomial(z,DistributedMultivariatePolynomial([a0,a1,a2,b0,b1,b2,t],Integer)))
→
fnd := fnl / (a0 * t^2 + 2 * a1 * t + 4 * a2)
There are 12 exposed and 12 unexposed library operations named /
having 2 argument(s) but none was determined to be applicable.
Use HyperDoc Browse, or issue
to learn more about the available operations. Perhaps
package-calling the operation or using coercions on the arguments
will allow you to apply the operation.
with argument type(s)
LaurentPolynomial(DistributedMultivariatePolynomial([a0,a1,a2,b0,b1,b2,t],Integer),UnivariatePolynomial(z,DistributedMultivariatePolynomial([a0,a1,a2,b0,b1,b2,t],Integer)))
Polynomial Integer
Perhaps you should use "@" to indicate the required return type,
or "$" to specify which version of the function you need.
fnz := fnl * z^(-2)
There are 34 exposed and 23 unexposed library operations named *
having 2 argument(s) but none was determined to be applicable.
Use HyperDoc Browse, or issue
to learn more about the available operations. Perhaps
package-calling the operation or using coercions on the arguments
will allow you to apply the operation.
with argument type(s)
LaurentPolynomial(DistributedMultivariatePolynomial([a0,a1,a2,b0,b1,b2,t],Integer),UnivariatePolynomial(z,DistributedMultivariatePolynomial([a0,a1,a2,b0,b1,b2,t],Integer)))
Fraction Polynomial Integer
Perhaps you should use "@" to indicate the required return type,
or "$" to specify which version of the function you need.
Martin Rubey
> Thanks, but...
>
> I'm having trouble manipulating Laurent Polynomials (I would think
> that multiplying by z^(-2), for example, would be straightforward.) Is
> there a Laurent equivalent to UP? Please see transcript below. I am
> using texmacs html output, is there a better way for email
> discussions?
Yes: use eaxiom as provided in the contrib directory of the svn repository of
fricas. If you are used to emacs, use emacs + axiom.el (eaxiom is just a
shortcut for that)
Then the key sequence
ctrl-u 3 alt-k
copies the last 3 input - output combinations into the clipboard.
What concerns the trouble with z^-2: for the moment you'll have to
1) convert z^-2 into a fraction of laurent polynomials
2) retract it into a laurent polynomial:
(z^-2)::FRAC M::M
the parenthesis is important: "::" binds stronger than anything else.
Martin
(HyperDoc) Cannot connect to the X11 server!
GCL (GNU Common Lisp) 2.6.8 CLtL1 Jan 29 2008 15:05:23
Source License: LGPL(gcl,gmp), GPL(unexec,bfd,xgcl)
Binary License: GPL due to GPL'ed components: (READLINE BFD UNEXEC)
Modifications of this banner must retain notice of a compatible license
Dedicated to the memory of W. Schelter
Use (help) to get some basic information on how to use GCL.
Temporary directory for compiler files set to /tmp/
openServer result 0
FriCAS (AXIOM fork) Computer Algebra System
Version: FriCAS 2008-05-28
Timestamp: Thursday July 31, 2008 at 20:02:08
-----------------------------------------------------------------------------
Issue )copyright to view copyright notices.
Issue )summary for a summary of useful system commands.
Issue )quit to leave FriCAS and return to shell.
-----------------------------------------------------------------------------
(1) -> )se me au of
(1) -> f(s) == (b2 * s^2 + b1 * s + b0)/(a2 * s^2 + a1 * s + a0)
Type: Void
(2) -> f((2/t)*( (1 - z^(-1))/(1 + z^(-1))))
Compiling function f with type Fraction Polynomial Integer ->
Fraction Polynomial Integer
2 2 2 2
(b0 t + 2b1 t + 4b2)z + (2b0 t - 8b2)z + b0 t - 2b1 t + 4b2
(2) ---------------------------------------------------------------
2 2 2 2
(a0 t + 2a1 t + 4a2)z + (2a0 t - 8a2)z + a0 t - 2a1 t + 4a2
Type: Fraction Polynomial Integer
(3) -> fn := numer(%)
2 2 2 2
(3) (b0 t + 2b1 t + 4b2)z + (2b0 t - 8b2)z + b0 t - 2b1 t + 4b2
Type: Polynomial Integer
(4) -> M ==> LAUPOL(DMP([a0,a1,a2,b0,b1,b2,t], INT), UP(z,DMP([a0,a1,a2,b0,b1,b2,t], INT)))
Type: Void
(5) -> fnl := fn :: M
2 2 2 2
(5) (b0 t + 2b1 t + 4b2)z + (2b0 t - 8b2)z + b0 t - 2b1 t + 4b2
Type: LaurentPolynomial(DistributedMultivariatePolynomial([a0,a1,a2,b0,b1,b2,t],Integer),UnivariatePolynomial(z,DistributedMultivariatePolynomial([a0,a1,a2,b0,b1,b2,t],Integer)))
(8) -> fnl * (z^-2)::FRAC M::M
2 2 - 1 2 - 2
(8) b0 t + 2b1 t + 4b2 + (2b0 t - 8b2)z + (b0 t - 2b1 t + 4b2)z
Type: LaurentPolynomial(DistributedMultivariatePolynomial([a0,a1,a2,b0,b1,b2,t],Integer),UnivariatePolynomial(z,DistributedMultivariatePolynomial([a0,a1,a2,b0,b1,b2,t],Integer)))
(6) ->
Ralf Hemmecke
>> is there a better way for email discussions?
> Yes: use eaxiom as provided in the contrib directory of the svn
> repository of fricas. If you are used to emacs, use emacs + axiom.el
> (eaxiom is just a shortcut for that)
If you are not used to emacs, start fricas on a konsole and retype the
commands that you want to send. It is much easier for us, if we read it
in plain ASCII.
Here is my suggestion...
Ralf
(1) -> P := DMP([a0,a1,a2,b0,b1,b2,t], INT)
(1) DistributedMultivariatePolynomial([a0,a1,a2,b0,b1,b2,t],Integer)
Type:
Domain
(2) -> M := LAUPOL(P, UP(z,P))
(2)
LaurentPolynomial(DistributedMultivariatePolynomial([a0,a1,a2,b0,b1,b2,t],Int
eger),UnivariatePolynomial(z,DistributedMultivariatePolynomial([a0,a1,a2,b0,b
1,b2,t],Integer)))
Type:
Domain
(3) -> m: M := z::M
(3) z
Type:
LaurentPolynomial(DistributedMultivariatePolynomial([a0,a1,a2,b0,b1,b2,t],Integer),UnivariatePolynomial(z,DistributedMultivariatePolynomial([a0,a1,a2,b0,b1,b2,t],Integer)))
(4) -> v(i:Integer): M == m^i
Function declaration v : Integer -> LaurentPolynomial(
DistributedMultivariatePolynomial([a0,a1,a2,b0,b1,b2,t],Integer),
UnivariatePolynomial(z,DistributedMultivariatePolynomial([a0,a1,
a2,b0,b1,b2,t],Integer))) has been added to workspace.
Type: Void
(5) -> 3*v(3) + 4*v(-2)
Your expression cannot be fully compiled because it contains an
integer expression (for i ) whose sign cannot be determined (in
general) and so must be specified by you. Perhaps you can try
substituting something like
(i :: PI)
or
(i :: NNI)
into your expression for i .
FriCAS will attempt to step through and interpret the code.
Compiling function v with type Integer -> LaurentPolynomial(
DistributedMultivariatePolynomial([a0,a1,a2,b0,b1,b2,t],Integer),
UnivariatePolynomial(z,DistributedMultivariatePolynomial([a0,a1,
a2,b0,b1,b2,t],Integer)))
3 - 2
(5) 3z + 4z
Type:
LaurentPolynomial(DistributedMultivariatePolynomial([a0,a1,a2,b0,b1,b2,t],Integer),UnivariatePolynomial(z,DistributedMultivariatePolynomial([a0,a1,a2,b0,b1,b2,t],Integer)))
(6) -> 3*v(3) * 4*v(-2)
(6) 12z
Type:
LaurentPolynomial(DistributedMultivariatePolynomial([a0,a1,a2,b0,b1,b2,t],Integer),UnivariatePolynomial(z,DistributedMultivariatePolynomial([a0,a1,a2,b0,b1,b2,t],Integer)))
(7) -> (3*v(3) +1)* 4*v(-2)
- 2
(7) 12z + 4z
Ralf Hemmecke
v(i:Integer): M == monomial(1, i)
Ralf
perrinmeyer
display as
a c e
--- + --- z^-1 + --- z^-2
b d g
(i.e. and then how do I collect terms in LaurentPolynomials, like
with :: UP(z,?) that takes a Fraction Polynomial Integer and converts
it to a UnivariatePolynomial(z,Fraction Polynomial Integer). Is there
a
Again, many, many thanks for your help. I'll try to get emacs +
axiom.el working before my next questions...
Perrin
→
fnlz := fnl * (z^(-2))::FRAC M::M
b0t2 + 2b1t + 4b2 + (2b0t2 - 8b2)z( - 1) + (b0t2 - 2b1t + 4b2)z( - 2)
(6)
LaurentPolynomial(DistributedMultivariatePolynomial([a0,a1,a2,b0,b1,b2,t],Integer),UnivariatePolynomial(z,DistributedMultivariatePolynomial([a0,a1,a2,b0,b1,b2,t],Integer)))
fnlz2 := fnlz / (a0 * t^2 + 2 * a1 * t + 4 * a2)::FRAC M::M
b0t2 + 2b1t + 4b2 + (2b0t2 - 8b2)z( - 1) + (b0t2 - 2b1t + 4b2)z( - 2)
a0t2 + 2a1t + 4a2
(8)
Type: Fraction
On Aug 11, 4:51 pm, Martin Rubey <martin.ru...@univie.ac.at> wrote:
Martin Rubey
> Thanks!, I'm getting closer, but, how do I get fnlz2 (see below) to
> display as
>
> a c e
> --- + --- z^-1 + --- z^-2
> b d g
>
> (i.e. and then how do I collect terms in LaurentPolynomials, like
> with :: UP(z,?) that takes a Fraction Polynomial Integer and converts
> it to a UnivariatePolynomial(z,Fraction Polynomial Integer).
You want to use
R ==> FRAC DMP([a0,a1,a2,b0,b1,b2,t], INT)
M ==> LAUPOL(R, UP(z, R))
instead, since the coefficients of your Laurent Polynomial are rational
functions, not polynomials.
> Again, many, many thanks for your help. I'll try to get emacs +
> axiom.el working before my next questions...
Please send email if you experience any trouble. I'll be back in roughly 8
hours...
Martin
perrinmeyer
Now that I have the Laurent Polynomial that I want, I'm having trouble
evaluating it numerically. I suppose I have to use "@" or "$"??
Below is the error message and the laupolerror.input script. I've
compiled FriCAS from the SVN source using GCL, so I'm using the latest
version.
Thanks again,
Perrin
----- laupolerror.input -----------------
)clear completely
R ==> FRAC DMP([a0,a1,a2,b0,b1,b2,p],INT)
M ==> LAUPOL(R, UP(z,R))
f(s) == (b2 * s^2 + b1 * s + b0) / (a2 * s^2 + a1 * s + a0)
fz := f(p * ((z-1)/(z+1)))
fzn := numer(fz)
fzn2 := fzn / (a2 * p^2 + a1 * p + a0)
fzn3 := fzn2::UP(z,?)
fzn4 := fzn3::FRAC M::M
fzn5 := fzn4 * (z^(-2))::FRAC M::M
fzd := denom(fz)
fzd2 := fzd / (a2 * p^2 + a1 * p + a0)
fzd3 := fzd2::UP(z,?)
fzd4 := fzd3::FRAC M::M
fzd5 := fzd4 * (z^(-2))::FRAC M::M
ffinal := fzn5::FRAC M::M / fzd5::FRAC M::M
eval(ffinal,[a0=1.0,a1=1.0,a2=1.0,b0=1.0,b1=1.0,p=1.0,z=1.0])
[FriCAS output ...]
ffinal := fzn5::FRAC M::M / fzd5::FRAC M::M
Loading
/home/perrin/cas/ax-install//lib/fricas/target/x86_64-unknown-
linux/algebra/FRETRCT-.o
for domain FullyRetractableTo&
Loading
/home/perrin/cas/ax-install//lib/fricas/target/x86_64-unknown-
linux/algebra/PRS.o
for package PseudoRemainderSequence
2 2
2
b0 + b1 p + b2 p 2b0 - 2b2 p - 1 b0 - b1 p + b2
p - 2
----------------- + ----------------- z +
----------------- z
2 2
2
a0 + a1 p + a2 p a0 + a1 p + a2 p a0 + a1 p + a2 p
(15)
-------------------------------------------------------------------
2 2
2a0 - 2a2 p - 1 a0 - a1 p + a2 p - 2
1 + ----------------- z + ----------------- z
2 2
a0 + a1 p + a2 p a0 + a1 p + a2 p
Type: Fraction LaurentPolynomial(Fraction
DistributedMultivariatePolynomial([a0,a1,a2,b0,b1,b2,p],Integer),UnivariatePolynomial(z,Fraction
DistributedMultivariatePolynomial([a0,a1,a2,b0,b1,b2,p],Integer)))
eval(ffinal,[a0=1.0,a1=1.0,a2=1.0,b0=1.0,b1=1.0,p=1.0,z=1.0])
Loading
/home/perrin/cas/ax-install//lib/fricas/target/x86_64-unknown-
linux/algebra/FLOAT.o
for domain Float
Loading
/home/perrin/cas/ax-install//lib/fricas/target/x86_64-unknown-
linux/algebra/EQ.o
for domain Equation
Loading
/home/perrin/cas/ax-install//lib/fricas/target/x86_64-unknown-
linux/algebra/LIST2.o
for package ListFunctions2
Loading
/home/perrin/cas/ax-install//lib/fricas/target/x86_64-unknown-
linux/algebra/FLAGG2.o
for package FiniteLinearAggregateFunctions2
Loading
/home/perrin/cas/ax-install//lib/fricas/target/x86_64-unknown-
linux/algebra/EQ2.o
for package EquationFunctions2
There are 10 exposed and 4 unexposed library operations named eval
Use HyperDoc Browse, or issue
package-calling the operation or using coercions on the
arguments
will allow you to apply the operation.
with argument type(s)
Fraction LaurentPolynomial(Fraction
DistributedMultivariatePolynomial([a0,a1,a2,b0,b1,b2,p],Integer),UnivariatePolynomial(z,Fraction
DistributedMultivariatePolynomial([a0,a1,a2,b0,b1,b2,p],Integer)))
List Equation Polynomial Float
Perhaps you should use "@" to indicate the required return type,
or "$" to specify which version of the function you need.
Martin Rubey
are you intending to work with fractions of Laurent polynomials?
Martin
perrinmeyer
filters can be transformed into digital filters by using the "Z
transform" and the "bilinear transform", and the expression ffinal is
a typical 2nd order linear filter transformed from the analog
(Laplace) domain to the digital domain ("z transform") using the
"bilinear transform).
It's late, so I'm afraid I'm not being very precise, maybe these links
will clarify what I'm intending to do:
(see, for example eq (7.5)
http://www.dsprelated.com/dspbooks/filters/Z_Transform_Difference_Equations.html
and
http://www.dsprelated.com/dspbooks/pasp/Bilinear_Transformation.html
Does this make any sense?
Thanks,
Perrin
On Aug 12, 11:45 pm, Martin Rubey <martin.ru...@univie.ac.at> wrote:
> I'm working on it, just a quick question:
>
> are you intending to work with fractions of Laurent polynomials?
>
> Martin
>
Martin Rubey
Although the following is a hack, I hope it's useful. Note that (for certain
reasons...) an element of FRAC L can only be coerced directly to FRAC POLY INT
if it doesn't contain z. That's why all this is so difficult and ugly.
It would pay off a lot to redesign the polynomial domain hierarchy. (Not the
algorithms, but "only" the hierarchy...) However, I'm afraid nobody around has
the energy to do that.
-- First fix the domains you are going to work in:
R ==> FRAC DMP([a0,a1,a2,b0,b1,b2,p],INT)
U ==> UP(z, R)
L ==> LAUPOL(R, U)
evalLAU(pol, l) == eval(pol::FRAC U::EXPR INT, l)::FRAC POLY INT::FRAC L::L
evalFLAU(fpol, l) == evalLAU(numer fpol, l) / evalLAU(denom fpol, l)
now:
(70) -> evalFLAU(ffinal, [z=1])
4b0
-----------------
2
a0 + a1 p + a2 p
(70) -------------------
4a0
-----------------
2
a0 + a1 p + a2 p
Type: Fraction LaurentPolynomial(Fraction DistributedMultivariatePolynomial([a0,a1,a2,b0,b1,b2,p],Integer),UnivariatePolynomial(z,Fraction DistributedMultivariatePolynomial([a0,a1,a2,b0,b1,b2,p],Integer)))
(71) -> evalFLAU(ffinal, [z=1])::FRAC POLY INT
b0
(71) --
a0
Type: Fraction Polynomial Integer
(73) -> evalFLAU(ffinal, [b=1])
2 2 2
b0 + b1 p + b2 p 2b0 - 2b2 p - 1 b0 - b1 p + b2 p - 2
----------------- + ----------------- z + ----------------- z
2 2 2
a0 + a1 p + a2 p a0 + a1 p + a2 p a0 + a1 p + a2 p
(73) -------------------------------------------------------------------
2 2
2a0 - 2a2 p - 1 a0 - a1 p + a2 p - 2
1 + ----------------- z + ----------------- z
2 2
a0 + a1 p + a2 p a0 + a1 p + a2 p
Type: Fraction LaurentPolynomial(Fraction DistributedMultivariatePolynomial([a0,a1,a2,b0,b1,b2,p],Integer),UnivariatePolynomial(z,Fraction DistributedMultivariatePolynomial([a0,a1,a2,b0,b1,b2,p],Integer)))
When I get to it, I'll try to improve the LaurentPolynomial domain...
Martin
Ralf Hemmecke
Maybe a stupid question, but do you *really* need Laurent polynomials?
What I see on
http://www.dsprelated.com/dspbooks/filters/Z_Transform_Difference_Equations.html
looks like you were fine with polynomials in t and rational functions in
t. Where in the end instead of substituting a value v for z you
substitute 1/v for t.
I don't know whether that makes the problems with the evaluation process
that Martin observed any simpler, but at least, you would use ordinary
polynomials and ordinary rational functions during your calculations.
Ralf
Ralf Hemmecke
> I think so. I'm interested in digital signal processing, where analog
> filters can be transformed into digital filters by using the "Z
> transform" and the "bilinear transform", and the expression ffinal is
> a typical 2nd order linear filter transformed from the analog
> (Laplace) domain to the digital domain ("z transform") using the
> "bilinear transform).
>
> It's late, so I'm afraid I'm not being very precise, maybe these links
> will clarify what I'm intending to do:
>
> (see, for example eq (7.5)
> http://www.dsprelated.com/dspbooks/filters/Z_Transform_Difference_Equations.html
> and
> http://www.dsprelated.com/dspbooks/pasp/Bilinear_Transformation.html
>
> Does this make any sense?
Actually not really. I don't think that expression (7.5) is all you
want. Must there be some simplification on that expression (which is of
course impossible if you just use indeterminates for the a's and b's) or
what does (7.5) actually show?
Is evaluating H(b0,...bM;a0,...,aN;z) your only goal?
Ralf
perrinmeyer
No, Evaluating H(z) is not my only goal (I have Matlab for that).
What I want is be able to do the algebraic manipulations called for by
the bilinear transform (in order to transform continuous time analog
filters into linear discrete time digital filter (difference
equations).
I have uploaded two pages from the book "Understanding Digital Signal
Processing 2nd Edition" (R.G Lyons), page 266-267.
http://groups.google.com/group/fricas-devel/web/bilinear.PDF?hl=en
So I'm after obtaining symbolic equations such as (6-112). I'm
interested in higher order filters, and thus more help with the
algebra...
Perrin
On Aug 13, 2:39 am, Ralf Hemmecke <r...@hemmecke.de> wrote:
> On 08/13/2008 09:26 AM, perrinmeyer wrote:
>
>
>
> > I think so. I'm interested in digital signal processing, where analog
> > filters can be transformed into digital filters by using the "Z
> > transform" and the "bilinear transform", and the expression ffinal is
> > a typical 2nd order linear filter transformed from the analog
> > (Laplace) domain to the digital domain ("z transform") using the
> > "bilinear transform).
>
> > It's late, so I'm afraid I'm not being very precise, maybe these links
> > will clarify what I'm intending to do:
>
> > (see, for example eq (7.5)
Ralf Hemmecke
> No, Evaluating H(z) is not my only goal (I have Matlab for that).
> What I want is be able to do the algebraic manipulations called for by
> the bilinear transform (in order to transform continuous time analog
> filters into linear discrete time digital filter (difference
> equations).
> I have uploaded two pages from the book "Understanding Digital Signal
> Processing 2nd Edition" (R.G Lyons), page 266-267.
>
> http://groups.google.com/group/fricas-devel/web/bilinear.PDF?hl=en
>
> So I'm after obtaining symbolic equations such as (6-112). I'm
> interested in higher order filters, and thus more help with the
> algebra...
Well, in that case you might want to rely on the abilities of the domain
Expression Integer. ... Oh, that is me, who is suggesting that... ;-)
Anyway, it seems for what you have in mind, the types are getting in the
way. So you might be happy with the following...
Oh... "sage" only appears below, because I have used sage as an
interface to FriCAS and exported that as text.
Ralf
sage: )clear all
All user variables and function definitions have been cleared.
sage: h:Expression Integer := c/(s^2+b*s+c)
c
------------
2
s + b s + c
Type: Expression Integer
sage: k:=subst(h, s=a*(1-t)/(1+t))
2
c t + 2c t + c
---------------------------------------------
2 2 2 2
(c - a b + a )t + (2c - 2a )t + c + a b + a
Type: Expression Integer
sage: n := numerator k
2
c t + 2c t + c
Type: Expression Integer
sage: d := denominator k
2 2 2 2
(c - a b + a )t + (2c - 2a )t + c + a b + a
Type: Expression Integer
sage: dp := d::Polynomial(Integer)
2 2 2 2
(c - a b + a )t + (2c - 2a )t + c + a b + a
Type: Polynomial Integer
sage: co:=coefficient(dp,t,0)
2
c + a b + a
Type: Polynomial Integer
sage: nc:=n/co
2
c t + 2c t + c
---------------
2
c + a b + a
Type: Expression Integer
sage: dc:=d/co
2 2 2 2
(c - a b + a )t + (2c - 2a )t + c + a b + a
---------------------------------------------
2
c + a b + a
Type: Expression Integer
sage: nr:=subst(nc,[a=200, b=137.94536, c=17410.145])
2
0.2048271221 1337193847 t + 0.4096542442 2674387695 t + 0.2048271221
1337193847
Type: Expression Float
sage: dr:=subst(dc,[a=200, b=137.94536, c=17410.145])
2
0.3508393847 9103872216 t - 0.5315308963 3755096826 t + 1.0
perrinmeyer
numerically, which helps me compare to results in matlab.
Perrin
perrinmeyer
On Aug 13, 4:27 pm, Ralf Hemmecke <r...@hemmecke.de> wrote:
> On 08/14/2008 12:00 AM, perrinmeyer wrote:
>
> > No, Evaluating H(z) is not my only goal (I have Matlab for that).
> > What I want is be able to do the algebraic manipulations called for by
> > the bilinear transform (in order to transform continuous time analog
> > filters into linear discrete time digital filter (difference
> > equations).
> > I have uploaded two pages from the book "Understanding Digital Signal
> > Processing 2nd Edition" (R.G Lyons), page 266-267.
>
> >http://groups.google.com/group/fricas-devel/web/bilinear.PDF?hl=en
>
> > So I'm after obtaining symbolic equations such as (6-112). I'm
> > interested in higher order filters, and thus more help with the
> > algebra...
>
> Well, in that case you might want to rely on the abilities of the domain
> Expression Integer. ... Oh, that is me, who is suggesting that... ;-)
>
"Expression Integer" is, but it seems to make sense.
Thanks again,
Perrin
Ralf Hemmecke
>> Expression Integer. ... Oh, that is me, who is suggesting that... ;-)
> I don't understand the humor...
Sorry. For that you have to browse the panAxiom mailing lists. I am the
one who hates 'Expression Integer'.
> This is a very helpful tutorial. I guess I don't understand what an
> "Expression Integer" is, but it seems to make sense.
The "Integer" here is maybe misleading. Just think of "Expression
Integer" as something that makes FriCAS behave like a (nearly) untyped
computer algebra system like Maple or Mathematica. Roughly speaking, you
forget about all the type stuff. That's not totally true, but for a new
user it is probably the best way to do the first steps. You can only
appreciate all those different types if you do more advanced things in
FriCAS. Types especially become important if one wants to implement a
bigger library.
One big difference for example is: In 'Expression Integer' it is
undecidable whether an element is actually equal to zero. So if you say
zero?(e) for some expresssion e you might get false although e can
actually be simplified to zero.
(1) -> e: Expression Integer := sin(x)^2 -1 + cos(x)^2
2 2
(1) sin(x) + cos(x) - 1
Type: Expression Integer
(2) -> zero? e
(2) false
Type: Boolean
(3) -> simplify e
(3) 0
Type: Expression Integer
For more specialized domains like DMP([a,b], Integer), the function
zero? always returns what you expect.
Hope that helps (or maybe it confuses you). Don't be too concerned about
such issues. For what you have in mind they most certainly don't matter.
Ralf