Examples
User documentation
The file ring.H
introduces several classes used for representing
rings and their elements. A normal user of the CoCoA library will use
principally the classes ring
and RingElem
: an object of type
ring
represents a mathematical ring with unity, and objects of
type RingElem
represent values from some ring. To make the
documentation more manageable it has been split into two: this file
describes operations on a RingElem
, whereas a separate file
describes the operations directly applicable to ring
s.
Documentation about the creation and use of homomorphisms between
rings can be found in RingHom
.
An object of type RingElem
comprises two internal parts: the ring
to which the value belongs, and the value itself. For instance, this
means that the zero elements of different rings are quite different
objects.
Constructors
Normally when creating a new RingElem
we specify both the ring to
which it belongs, and its initial value in that ring. Let R
be a
ring
. Let n
be a machine integer or a BigInt
. Let q
be a rational, i.e. a value of type BigRat
. Let r2
be a
ring element.
RingElem r; |
the zero element of RingZZ() (special case) |
RingElem r(R); |
an element of R , initially 0 |
RingElem r(R, n); |
an element of R , initially the image of n |
RingElem r(R, q); |
an element of R , initially the image of q (or error) |
RingElem r(R, str); |
an element of R , initially the value of string str |
RingElem r(R, s); |
an element of R , initially the value of symbol s |
RingElem r(R, r2); |
an element of R , maps r2 into R via CanonicalHom |
RingElem r(R, MPZ); |
an element of R , initially the value of mpz_t MPZ |
RingElem r(R, MPQ); |
an element of R , initially the value of mpq_t MPQ (or error) |
Note 1: To create a RingElem
from a value of type mpz_class
you must do RingElem(R, MPZ.get_mpz_t())
; analogously for mpq_class
.
Note 2: Construction from a rational may fail, e.g. if the denominator is a zero divisor in the ring; if it does fail then an exception is thrown (with code ERR::DivByZero
).
Note 3: There is also a string constructor for making vectors:
RingElems(R, str)
-- avector<RingElem>
of elements ofR
, initially the values in stringstr
(separated by commas)
You can create a copy of a ring element in the usual way:
RingElem r(r2); |
a copy of r2 , element of the same ring |
RingElem r = r2; |
(alternative syntax, discouraged) |
Naturally the last constructor works only if the denominator of the rational q
is 1.
These are not really constructors: you can get the zero and one of a
ring
directly using the following:
zero(R) |
the zero element of R |
one(R) |
the one element of R |
Operations on RingElems
RingElem
s are designed to be easy and safe to use; the various checks
do incur a certain run-time overhead, so a faster alternative is offered
(see below in the section Fast and Ugly Code). Arithmetic operations
between RingElem
s will fail if they do not belong to the same ring
(the exception has code ERR::MixedRings
).
Assignment & Swapping
Assigning an integer or rational to a RingElem
wil automatically map
the value into the ring to which the RingElem
belongs.
r = n; |
map n into owner(r) and assign the result |
r = q; |
map q into owner(r) and assign the result |
r = r2; |
r becomes a copy of r2 -- afterwards owner(r)==owner(r2) |
swap(r,s); |
exchange the values (and the owning rings) |
Arithmetic
Arithmetic operations between RingElem
s will fail if they do not belong
to the same ring (the exception has code ERR::MixedRings
). You may
perform arithmetic between a RingElem
and a machine integer, a BigInt
value
or a BigRat
value -- the integer/rational is automatically mapped into the
same ring as the RingElem
.
Let r
be a non-const RingElem
,
and r1
, r2
be potentially const RingElem
s.
Assume they are all associated to the same ring.
Then the operations available are: (meanings are obvious)
cout << r1
-- output value ofr1
(decimal only, see notes)r1 == r2
-- equality testr1 != r2
-- not-equal test-r1
-- negation (unary minus)r1 + r2
-- sumr1 - r2
-- differencer1 * r2
-- productr1 / r2
-- quotient, division must be exact (seeIsDivisible
)r += r1
-- equivalent tor = r + r1
r -= r1
-- equivalent tor = r - r1
r *= r1
-- equivalent tor = r * r1
r /= r1
-- equivalent tor = r / r1
division must be exact (seeIsDivisible
)power(r1, n)
--n
-th power ofr1
;n
any integer, NBpower(0,0)
gives 1r^n
-- THIS DOES NOT WORK!!! it does not even compile, you must usepower
Attempting to compute a gcd
or lcm
in a ring which not an
effective GCD domain will produce an exception with code ERR::NotTrueGCDDomain
.
If r1
or r2
is a BigRat
then an error is signalled at compile time.
gcd(r1, r2)
-- an associate of the gcdIsCoprime(r1, r2)
-- true iffIsInvertible(gcd(r1,r1))
lcm(r1, r2)
-- an associate of the lcmGcdQuot(&gcd, "1, "2, r1, r2)
-- procedure computesgcd
andquot1=r1/gcd
andquot2=r2/gcd
, herer1
andr2
must beRingElem
.
There is a class called GCDFreeBasis_RingElem
for computing a factor base
of a set of RingElem
from an effective GCD domain:
GCDFreeBasis_RingElem()
default ctor; base is initially emptyGFB.myAddInfo(x)
use alsoRingElem x
when computing he factor baseGFB.myAddInfo(v)
use also the elements ofstd::vector<RingElem> v
when computing the factor baseFactorBase(GFB)
returns the factor base obtained so far
Queries
CoCoALib offers functions for querying various properties of RingElem
s,
and about relationships between RingElems
.
Let r1
and r2
be a (possibly const) RingElem
s, and
let N
be a variable of type BigInt
, and q
a variable of type BigRat
owner(r1)
-- the ring to whichr1
is associatedIsZero(r1)
-- true iffr1
is zeroIsOne(r1)
-- true iffr1
is oneIsMinusOne(r1)
-- true iff-r1
is oneIsInvertible(r1)
-- true iffr1
has a multiplicative inverseIsZeroDivisor(r1)
-- true iffr1
is zero-divisorIsDivisible(r1, r2)
-- true iffr1
is divisible byr2
(throwsERR::DivByZero
ifr2
is a zero-divisor)IsDivisible(r, r1, r2)
--r = r1/r2
and returns true iffr1
is divisible byr2
IsDivisible_AllowFields(r1, r2)
-- true iffr1
is divisible byr2
(throwsERR::DivByZero
ifr2
is a zero-divisor)IsDivisible_AllowFields(r, r1, r2)
--r = r1/r2
and returns true iffr1
is divisible byr2
IsInteger(N, r1)
-- true iffr1
is the image of an integer (if true, a preimage is placed inN
, otherwise it is left unchanged)IsRational(q, r1)
-- true iffr1
is the image of a rational (if true, a preimage is placed inq
, otherwise it is left unchanged)IsDouble(d, r1)
-- true iffr1
is the image of a rational whose approx is put into d (false if overflow andd
unchanged)
Note that IsDivisible
tests divisibility in the ring containing the
values: so 1 is not divisible by 2 in RingZZ
, but their images in
RingQQ
would be divisible. Also note that IsDivisible
throws an
exception if the test is in a field; use IsDivisible_AllowFields
if
you want to allow divisibility testing also in a field.
Ordering
If the ring is an ordered domain then these functions may also be
used. You can discover whether CoCoALib thinks that the ring R
is
arithmetically ordered by calling IsOrderedDomain(R)
: the value is
true
iff R
is arithmetically ordered.
Note that comparison operations between RingElem
s will fail if they do not belong
to the same ring (the exception has code ERR::MixedRings
). You may
perform comparisons between a RingElem
and an integer or a rational -- the integer/rational is automatically mapped into the
same ring as the RingElem
.
Let r1
and r2
belong to an ordered ring. Trying to use any
of these functions on elements belonging to a ring which is not ordered
will produce an exception with code ERR::NotOrdDomain
.
sign(r1)
-- value is -1, 0 or +1 according asr1
is negative, zero or positiveabs(r1)
-- absolute value ofr1
floor(r1)
-- greatest integer<= r1
ceil(r1)
-- least integer>= r1
NearestInteger(r1)
-- returns nearest integer (BigInt
) tor1
(halves round as inround
, seeBigRat
).cmp(r1, r2)
-- returns a value <0, =0, >0 according asr1-r2
is <0, =0, >0CmpAbs(r1, r2)
-- equiv tocmp(abs(r1),abs(r2))
CmpDouble(r1, z)
-- compare a ring elem with adouble
, result is <0, =0, >0 according asr1-z
is <0, =0, >0r1 < r2
-- standard inequalitiesr1 > r2
-- ...r1 <= r2
-- ...r1 >= r2
-- ...
More operations on RingElems of a finite field
If owner(r)
is a finite field, or a polynomial ring whose coeffs are
ini a finite field then the following functions may be used. The P
in the names of these function refers to the characteristic of the ring.
IsPthPower(r1)
-- true iffr1
has a p-th root where p is the ring characteristic (see alsoPthRoot
)PthRoot(r1)
-- returns the p-th root ofr1
(error if no p-th root exists)
More operations on RingElems of a FractionField
If owner(r)
is a fraction field then the following functions may
be used. You can find out whether CoCoALib thinks that the ring R
is a fraction field by calling IsFractionField(R)
: the result is
true
iff R
is a fraction field.
Let K
denote a FractionField
Let r
denote an element of K
.
num(r)
-- gives a copy of the numerator ofr
as an element ofBaseRing(K)
den(r)
-- gives a copy of the denominator ofr
as an element ofBaseRing(K)
Note: the numerator and denominator are defined only upto
multiples by a unit: so it is (theoretically) possible to have two
elements of FrF
which are equal but which have different
numerators and denominators, for instance, (x-1)/(x-2) = (1-x)/(2-x)
More operations on RingElems of a QuotientRing
If owner(r)
is a quotient ring then the following function may be
called. You can find out whether CoCoALib thinks that the ring R
is a quotient ring by calling IsQuotientRing(R)
: the result is
true
iff R
is a quotient ring.
In addition to the standard RingElem
operations, elements of a
QuotientRing
may be used in other functions.
Let RmodI
denote a quotient ring.
Let r
denote a non-const element of RmodI
.
More operations on RingElems of a RingTwinFloat
You can determine if an element belongs to a twin-float ring by
calling IsRingTwinFloat(owner(r))
: this yields true
iff r
belongs to a twin-float ring.
Let x
, y
be RingElem
belonging to a RingTwinFloat
MantissaAndExponent2(x)
-- produce aMantExp2
structure representing the approximate value ofx
(forMantExp2
seeFloatApprox
)DebugPrint(out, x)
-- print out both components ofx
IsPracticallyEqual(x, y)
-- returns true ifIsZero(x-y)
otherwise false.
In contrast the test x==y
may throw a
RingTwinFloat::InsufficientPrecision
while
IsPracticallyEqual
will never throw this exception.
IsPracticallyEqual
is intended for use in a termination
criterion for an iterative approximation algorithm (e.g. see
test-RingTwinFloat4.C
).
More operations on RingElems of a PolyRing
You can determine whether an element belongs to a PolyRing by
calling IsPolyRing(owner(r))
: the result is true
iff r
belongs to a poly ring.
Let P
denote a polynomial ring.
Let f
denote a non-const element of P
.
Let f1
, f2
denote const elements of P
.
Let v
denote a const vector of elements of P
.
IsMonomial(f)
-- true ifff
is non zero and of the form coeff*ppAreMonomials(v)
-- true iffv
is non zero and of the form coeff*pp (ifv
is empty it returns true)IsConstant(f)
-- true ifff
is "constant", i.e. the image of an element of the coeff ring.IsIndet(f)
-- equivalent tof
== x[i] for some index iIsIndet(index, f)
-- equivalent tof
== x[i]; and sets index = iIsIndetPosPower(f)
-- true iff there is indexi
and exponente
such thatf == x[i]^e
IsEvenPoly(f)
-- true ifff
is even as a functionIsOddPoly(f)
-- true ifff
is odd as a functionIsIrred(f)
-- true ifff
is irreducible in Powner(f1)
-- the owner off
as aring
.NumTerms(f1)
-- the number of terms inf1
.StdDeg(f1)
-- the standard degree off1
(deg(x[i])=1); error iff1
is 0.deg(f1)
-- same asStdDeg(f1)
.deg(f1, var)
-- maximum degree ofvar
-th indet inf1
wherevar
is the index of the indet inP
(result is of type long).LC(f1)
-- the leading coeff off1
; it is an element of CoeffRing(P).ConstantCoeff(f)
-- the constant coeff off
; may be zero, belongs to CoeffRing(P).monic(f1)
-- same asf1/LC(f1)
content(f1)
-- gcd of the coeffs off1
; it is an element of CoeffRing(P). Content of zero poly is 0. If coeffs are in a fraction field then
content isc
such thatf1/c
is a primitive polynomial.
FixedDivisor(f1)
-- computes gcd off1(n)
asn
ranges over all integers (aka. intrinsic content).CommonDenom(f1)
-- the simplest common denominator for the coeffs off1
; it is an element of BaseRing(CoeffRing(P)); throws if CoeffRing is not a FractionField of a GCD domain.ClearDenom(f1)
--f1
*CommonDenom(f1)
(same restrictions as above)ClearDenom(Rx, f1)
-- likeClearDenom(f1)
but puts result in (SparsePoly)ringRx
prim(f1)
-- same asg = ClearDenom(f1); return g/content(g);
deriv(f1, var)
-- formal derivative off1
wrt. indet having indexvar
.deriv(f1, x)
-- derivative off1
w.r.t.x
,x
must be an indeterminate (also works for f1 inFractionField
of aPolyRing
)
NOTE: to compute the weighted degree of a polynomial use the function
wdeg
defined for RingElem
of a SparsePolyRing
(see below).
More operations on RingElems of a SparsePolyRing
You can determine whether an element belongs to a sparse poly ring by
calling IsSparsePolyRing(owner(r))
: the result is true
iff r
belongs to a (sparse) poly ring.
Let P
denote a SparsePolyRing
.
Let f
denote a non-const element of P
.
Let f1
, f2
denote const elements of P
.
Let expv
be a vector<long>
of size equal to the number of indeterminates.
owner(f1)
-- the owner off1
as aring
NumTerms(f1)
-- the number of terms inf1
with non-zero coefficient.UnivariateIndetIndex(f)
-- iff
is univariate in j-th indet returns j, o/w returns -1; throws error iff
is constantLPP(f1)
-- the leading PP off1
; it is an element of PPM(P). Also known as LT(f) or in(f)LF(f1)
-- the leading form off1
; sum of all summands of highest weighteddegree
; gives error iff1
is zero.CutLF(f)
-- likeLF(f)
, but also modifiesf
to becomef-LF(f)
.HomogCompt(f,d)
-- return homogeneous part off
of degd
(does not modifyf
)wdeg(f1)
-- the weighteddegree
of the leading PP off1
(see [KR] Sec.4.3); error iff1
is 0. NB result is of typeCoCoA::degree
(seedegree
). (contrast withStdDeg(f1)
anddeg(f1)
defined for generalPolyRing
)CmpWDeg(f1, f2)
-- compare the weighted degrees of the LPPs off1
andf2
; result is <0 =0 >0 according as deg(f1
) < = > deg(f2
)IsHomog(f)
-- says whetherf
is homogeneous wrt weighteddegree
.homog(f, h)
-- returnsf
homogenized with indeth
(requiresGrDim
=1 andwdeg(h)
=1)NR(f, v)
-- returns the (normal) remainder of the Groebner Division Algorithm byv
. Ifv
is a GBasis this is the Normal Form (seeNF
in the doc forideal
).indet(P,v)
-- returns thev
-th indet ofP
as aRingElem
IndetPower(P,v,e)
-- returns thee
-th power of thev
-th indet ofP
as aRingElem
monomial(P,pp)
-- returnspp
as an element ofP
monomial(P,c,pp)
-- returnsc
*pp
as an element ofP
wherec
is an integer, rational or is inCoeffRing(P)
andpp
is inPPM(P)
.monomial(P,c,expv)
-- returnsc*x[0]^expv[0]*x[1]^expv[1]*...
wherec
is an integer, rational or is inCoeffRing(P)
, and x[i] are the indets ofP
.
Let X
be an indet (i.e. a RingElem
in P
)
or a vector of indices (vector<long>
)
IndetsProd(f)
-- monomial which is product of all indets actually inf
ContentWRT(f, X)
-- the content off
wrt the indet(s)X
; result is aRingElem
inP
CoefficientsWRT(f, X)
-- returns avector<CoeffPP>
: eachCoeffPP
has fieldsmyCoeff
andmyPP
wheremyCoeff
is an element ofP
andmyPP
is inPPM(P)
being a power product of the indets inX
; the entries are in decreasing order ofmyPP
.CoeffVecWRT(f, x)
--x
must be an indet; returns avector<RingElem>
whose k-th entry contains the coeff ofx^k
as an element ofP
; NB the coeff may be zero!CoeffVecWRTSupport(f, B)
-- returns avector<RingElem>
being the representation off
in the basis given bysupport(B)
; error iff
is not in the span of the basisCoeffHeight(f)
-- returns the max of the absolute values of the coeffs off
IsPalindromic(f)
-- returnstrue
ifff
is palindromic; error iff
is not univariatereverse(f)
-- returnsx^deg(f)*f(1/x)
; error iff
is not univariatereverse(f,t)
-- returnsf
with each power productPP
replacedt/PP
graeffe(f)
-- returns graeffe transformation of univariatef
; its roots are the squares of the roots off
.graeffe3(f)
-- returns cubic graeffe transformation of univariatef
; its roots are the cubes of the roots off
.MinPoly
-- [SparsePolyOps-MinPoly].
NB For running through the summands (or terms) of a
polynomial use SparsePolyIter
s (see SparsePolyRing
).
We have still doubts on the usefulness of these two functions:
CmpWDegPartial(f1, f2, i)
-- compare the firsti
weighted degrees of the LPPs off1
andf2
; result is <0 =0 >0 according asdeg(f1)
< = >deg(f2)
IsHomogPartial(f,i)
-- says whetherf
is homogeneous wrt the firsti
components of the weighted degree
Use the following two functions with great care: they throw an error
if the PPOrdering
is not respected: (the coefficient c
may be 0)
PushFront(f, c, t)
-- add tof
the termc
*t where t is a PP belonging toPPM(owner(f))
and assuming that t >LPP(f)
orf
==0PushBack(f, c, t)
-- add tof
the termc
*t where t is a PP belonging toPPM(owner(f))
and assuming that t < t' for all t' appearing inf
.PushFront(f, c, expv)
-- add tof
the termc
*t where t is the PP with exponent vectorexpv
, and assuming that t >LPP(f)
orf
==0PushBack(f, c, expv)
-- add tof
the termc
*t where t is the PP with exponent vectorexpv
, and assuming that t < t' for all t' appearing inf
.
The corresponding member functions myPushFront/myPushBack
will not
check the validity of these assumptions: they have a CoCoA_ASSERT
to check them only in DEBUG mode.
More operations on RingElems of a DenseUPolyRing
You can determine whether an element belongs to a DenseUPolyRing by
calling IsDenseUPolyRing(owner(r))
: the result is true
iff r
belongs to a poly ring.
Let P
denote a DenseUPolyRing
.
Let f
denote an element of P
.
monomial(P,c,exp)
--c
*x^exp
as an element ofP
withc
an integer or inCoeffRing(P)
exp
aMachineInt
coeff(f,d)
-- thed
-th coefficient off
(as aConstRingElem
, read-only)
WARNING Use this functions with great care: no checks on size and degree
Let f
denote a non-const element of P
.
myAssignCoeff(f,c,d)
-- assigns thed
-th coefficient inf
toc
myAssignZeroCoeff(f,d)
myAssignNonZeroCoeff(f,c,d)
Notes on operations
Operations combining elements of different rings will cause a run-time error.
In all functions involving two RingElem
s either r1
or r2
may be replaced by a machine integer, or by a big integer (an element
of the class BigInt
). The integer value is automatically mapped into
the ring owning the RingElem
in the same expression.
The exponent n
in the power function may be zero or negative, but a
run-time error will be signalled if one attempts to compute a negative
power of a non-invertible element.
NB You cannot use ^
to compute powers -- see Bugs section.
An attempt to perform an inexact division or to compute a GCD not in a GCD domain will produce a run-time error.
The printing of ring elements is always in decimal regardless of the
ostream
settings (this is supposed to be a feature rather than a bug).
At this point, if you are new to CoCoALib, you should probably look
at some of the example programs in the examples/
directory.
Writing functions with RingElems as arguments
One would normally expect to use the type const RingElem&
for
read-only arguments which are RingElem
s, and RingElem&
for
read-write arguments. Unfortunately, doing so would lead to problems
with the CoCoA library. INSTEAD you should use the types:
ConstRefRingElem x |
for read-only arguments: morally const RingElem& x |
RingElem& x |
for read-write arguments |
RingElem x |
for read-only arguments which make a local copy |
If you are curious to know why this non-standard quirk has to be used, read on.
When accessing matrix elements or coefficients in a polynomial
CoCoALib uses proxies: these are objects which should behave much
like const RingElem
values. To allow easy use of such proxies in
functions which want a read-only RingElem
we use the type
ConstRefRingElem
(which is actually const RingElemAlias&
) for
the formal parameter.
Internally, ring element values are really smart pointers to the true
value. Now the const
keyword in C++ when applied to a pointer
makes the pointer const while the pointed-to value remains alterable
-- this is not the behaviour we want for const RingElem&
. To get
the desired behaviour we have to use another type: the type we have
called ConstRefRingElem
.
ADVANCED USE OF RingElem
The rest of this section is for more advanced use of ring
s and
RingElem
s (e.g. by CoCoA library contributors). If you are new to
CoCoA, you need not read beyond here.
Fast and Ugly Code
WE DO NOT RECOMMEND that you use what is described in this section. If you are curious to know a bit more how rings are implemented, you might find this section informative.
RingElem
s are designed to be easy and pleasant to use, but this
convenience has a price: a run-time performance penalty (and a memory
space penalty too).
Both penalities may be avoided by using raw values but at a
considerable loss of programming convenience and safety. You should
consider using raw values only if you are desperate for speed; even
so, performance gains may be only marginal except perhaps for
operations on elements of a simple ring (e.g. a small finite field).
A RingElem
object contains within itself an indication of the owning
ring, and a raw value which is a pointer to where the real
representation of the ring element value lies. These raw values may be
accessed via the raw function. They may be combined arithmetically by
calling member functions of the owning ring. For instance, if x,y,z are
all RingElem objects all BELONGING TO EXACTLY THE SAME RING then we can
achieve
x = y+z;
slightly faster by calling
owner(x)->my.Add(raw(x), raw(y), raw(z));
It should now be clear that the syntax involved is cumbersome and
somewhat obscure. For the future maintainability of the code the
simpler x = y+z;
has many advantages. Furthermore, should x,y,z
somehow happen not all to lie in the same ring then x = y+z;
will act
in a reasonable way, whereas the supposedly faster call will likely lead
to many hours of debugging grief. The member functions for arithmetic
(e.g. myAdd
) DO NOT PERFORM sanity checks on their arguments:
e.g. attempting to divide by zero could well crash the program.
If you use a debugging version of the CoCoA Library then some member functions do use assertions to check their arguments. This is useful during development, but must not be relied upon since the checks are absent from the non-debugging version of the CoCoA Library. See the file config.txt for more information.
This fast, ugly, unsafe way of programming is made available for those who desperately need the speed. If you're not desperate, don't use it!
Fast, Ugly and Unsafe operations on raw values
Read the section Fast and Ugly Code before using any of these!
Let r
be a non-const raw value (e.g. raw(x)
, with x
a
RingElem
), and r1
, r2
potentially const raw values.
Assume they are all owned by the ring R
.
Then the functions available are:
R->myNew()
-- construct a new element of R, value=0R->myNew(n)
-- construct a new element of R, value=nR->myNew(N)
-- construct a new element of R, value=NR->myNew(r1)
-- construct a new element of R, value=r1R->myDelete(r)
-- destroy r, element of R (frees resources)R->mySwap(r, s)
-- swaps the two values (s is non-const raw value)R->myAssignZero(r)
-- r = 0R->myAssign(r, r1)
-- r = r1R->myAssign(r, n)
-- r = n (n is a long)R->myAssign(r, N)
-- r = n (N is aBigInt
)R->myNegate(r, r1)
-- r = -r1R->myAdd(r, r1, r2)
-- r = r1+r2R->mySub(r, r1, r2)
-- r = r1-r2R->myMul(r, r1, r2)
-- r = r1*r2R->myDiv(r, r1, r2)
-- r = r1/r2 (division must be exact)R->myIsDivisible(r, r1, r2)
-- r = r1/r2, and returns true iff division was exactR->myIsZeroDivisor(r)
-- returns true iff r is a zero-divisorR->myIsUnit(r1)
-- IsUnit(r1)R->myGcd(r, r1, r2)
-- r = gcd(r1, r2)R->myLcm(r, r1, r2)
-- r = lcm(r1, r2)R->myPower(r, r1, n)
-- r = power(r1, n) BUT n MUST be non-negative!!R->myIsZero(r1)
-- r1 == 0R->myIsZeroAddMul(r, r1, r2)
-- ((r += r1*r2) == 0)R->myIsEqual(r1, r2)
-- r1 == r2R->myIsPrintAtom(r1)
-- true iffr1
does not need brackets when a num or denom of a fractionR->myIsPrintedWithMinus(r1)
-- true iff the printed form ofr1
begins with a minus signR->myOutput(out, r1)
-- out << r1R->mySequentialPower(r, r1, n)
-- normally it is better to use R->myPower(r, r1, n)R->myBinaryPower(r, r1, n)
-- normally it is better to use R->myPower(r, r1, n)
Maintainer documentation
(NB consider consulting also QuotientRing
, FractionField
and PolyRing
)
The design underlying rings and their elements is more complex than I would have liked, but it is not as complex as the source code may make it appear. The guiding principles are that the implementation should be flexible and easy/pleasant to use while offering a good degree of safety; extreme speed of execution was not a goal (as it is usually contrary to good flexibility) though an interface offering slightly better run-time efficiency remains.
Regarding flexibility: in CoCoALib we want to handle polynomials whose
coefficients reside in (almost) any commutative ring. Furthermore, the
actual rings to be used will be decided at run-time, and cannot
restricted to a given finite set. We have chosen to use C++ inheritance
to achieve the implementation: the abstract class RingBase
defines the
interface that every concrete ring class must offer.
Regarding ease of use: since C++ allows the common arithmetic operators
to be overloaded, it is essential that these work as expected for
elements of arbitrary rings -- with the caveat that /
means exact
division, being the only reasonable interpretation. Due to problems of
ambiguity arithmetic between elements of different rings is forbidden:
e.g. let f in Q[x,y] and g in Z[y,x], where should f+g reside?
The classes in the file ring.H are closely interrelated, and there is no obvious starting point for describing them -- you may find that you need to read the following more than once to comprehend it. Here is a list of the classes:
ring |
value represents a ring; it is a smart pointer |
RingBase |
abstract class defining what a ring is |
RingElem |
value represents an element of a ring |
RingElemAlias |
reference to a RingElem belonging to someone else |
ConstRefRingElem |
C++ const-reference to a RingElemAlias |
RingElemConstRawPtr |
raw pointer to a const ring value |
RingElemRawPtr |
raw pointer to a ring value |
For the first two see ring
.
The classes RingElem
and RingElemAlias
are
related by inheritance: they are very similar but differ in one important
way. The base class RingElemAlias
defines the data members
which are inherited by RingElem
. The essential difference is that
a RingElem
owns the value whereas a RingElemAlias
does not.
The two data members are myR
and myRawValue
: the first is the
identity of ring to which the element belongs, and the second is the
value in that ring (the value is stored in a format that only the
owning ring can comprehend). All operations on ring elements are
effected by member functions of the ring to which the value belongs.
The differing ownership inherent in RingElemAlias
and RingElem
lead to several consequences. The destructor of a RingElem
will
destroy in the internal representation of the value; in contrast, the
destructor of a RingElemAlias
does nothing. A RingElemAlias
object becomes meaningless (& dangerous) if the owner of the value it
aliases is destroyed.
Why did I create RingElemAlias
? The main reason was to allow
matrices and iterators of polynomials to be implemented cleanly and
efficiently. Clearly a matrix
should be the owner of the values
appearing as its entries, but we also want a way of reading the matrix
entries without having to copy them. Furthermore, the matrix can use
a compact representation: the ring to which its elements belong is
stored just once, and not once for each element. Analogous comments
apply to the coefficients of a polynomial.
As already stated above, the internal data layouts for objects of
types RingElem
and RingElemAlias
are identical -- this is
guaranteed by the C++ inheritance mechanism. The subfield indicating
the ring to which the value belongs is simply a ring
, which is
just a reference counting smart pointer. The subfield indicating the
value is a raw pointer of type void*
; however, when the raw
pointer value is to be handled outside a ring element object then it
is wrapped up as a RingElemRawPtr
or RingElemConstRawPtr
--
these are simply wrapped copies of the void*
.
The classes RingElemRawPtr
and RingElemConstRawPtr
are used
for two reasons. One is that if a naked void*
were used outside
the ring element objects then C++ would find the call
RingElem(R,0)
ambiguous because the constant 0
can be
interpreted either as an integer constant or as a null pointer: there
are two constructors which match the call equally well. The other
reason is that it discourages accidentally creating a ring element
object from any old pointer; it makes the programmer think -- plus I
feel uneasy when there are naked void*
pointers around. Note that
the type of the data member RingElemConstRawPtr::myPtr
is simply
void*
as opposed to void const*
which one might reasonably
expect. I implemented it this way as it is simpler to add in the
missing constness in the member function
RingElemConstRawPtr::myRawPtr
than it would be to cast it away in
the myRawPtr
function of RingElemRawPtr
.
Further comments about implementation aspects of the above classes.
The class RingBase
declares a number of pure virtual functions for
computing with ring elements. Since these functions are pure they
must all be fully defined in any instantiable ring class
(e.g. RingZZImpl
or RingFpImpl
). These member functions follow
certain conventions:
- RETURN VALUES:
- most arithmetic functions return no value, instead the result is placed in one of the arguments (normally the first argument is the one in which the result is placed), but functions which return particularly simple values (e.g. booleans or machine integers) do indeed return the values by the usual function return mechanism.
- ARG TYPES:
-
ring element values are passed as raw pointers
(i.e. a wrapped
void*
pointing to the actual value). A read-only arg is of typeRingElemConstRawPtr
, while a writable arg is of typeRingElemRawPtr
. When there are writable args they normally appear first. For brevity there are typedefsConstRawPtr
andRawPtr
in the scope ofRingBase
or any derived class. - ARG CHECKS:
- sanity checks on the arguments are NOT CONDUCTED (e.g. the division function assumes the divisor is non-zero). These member functions are supposed to be fast rather than safe.
In a few cases there are non-pure virtual member functions in
RingBase
. They exist either because there is a simple universal
definition or merely to avoid having to define inappropriate member
functions (e.g. gcd functions when the ring cannot be a gcd domain).
Here is a list of them:
myIsUnit(x)
-- default checks that 1 is divisible byx
myIsZeroDivisor(x)
-- special implementation in [QuotientRing] for setting primality flag to defining idealmyGcd(lhs, x, y)
-- gives an error: eitherNotGcdDom
orNYI
myLcm(lhs, x, y)
-- gives an error: eitherNotGcdDom
orNYI
myGcdQuot(lhs, xquot, yquot, x, y)
-- gives an error: eitherNotGcdDom
orNYI
myExgcd(lhs, xcofac, ycofac, x, y)
-- gives an error: eitherNotGcdDom
orNYI
myIsPrintAtom(x)
-- defaults to falsemyIsPrintedWithMinus(x)
-- givesShouldNeverGetHere
errormyIsMinusOne(x)
-- defaults to myIsOne(-x); calculates -xmyIsZeroAddMul(lhs, y, z)
-- computes lhs += y*z in the obvious way, and callsmyIsZero
myCmp(x, y)
-- givesNotOrdDom
errormyCmpAbs(x, y)
-- tries to computecmp(abs(x),abs(y))
so may giveNotOrdDom
errormySign(x)
-- simply callsmyCmp(x, 0)
, then returns -1,0,1 accordingly There are three non-virtual member functions for calculating powers: one uses the sequential method, the other two implement the repeated squaring method (one is an entry point, the other an implementation detail). These are non-virtual since they do not need to be redefined; they are universal for all rings. For the moment I shall assume that the intended meaning of the pure virtual functions is obvious (given the comments in the source code).
Recall that arithmetic operations on objects of type
ConstRefRingElem
(which matches RingElem
too) are converted
into member function calls of the corresponding owning ring. Here is
the source code for addition of ring elements -- it typifies the
implementation of operations on ring elements.
RingElem operator+(ConstRefRingElem x, ConstRefRingElem y) { const ring& Rx = owner(x); const ring& Ry = owner(y); if (Rx != Ry) error(CoCoAError(ERR::MixedRings, "RingElem + RingElem")); RingElem ans(Rx); Rx->myAdd(raw(ans), raw(x), raw(y)); return ans; }
The arguments are of type ConstRefRingElem
since they are
read-only, and the return type is RingElem
since it is new
self-owning value (it does not refer to a value belonging to some
other structure). Inside the function we check that the rings of the
arguments are compatible, and report an error if they are not.
Otherwise a temporary local variable is created for the answer, and
the actual computation is effected via a member function call to the
ring in which the values lie. Note the use of the raw
function
for accessing the raw pointer of a ring element. In summary, an
operation on ring elements intended for public use should fully check
its arguments for compatibility and correctness (e.g. to avoid
division by zero); if all checks pass, the result is computed by
passing raw pointers to the appropriate member functions of the ring
involved -- this member function assumes that the values handed to it
are compatible and valid; if not, undefined behaviour will result
(i.e. a crash if you are lucky).
Most of the member functions of a ring are for manipulating raw values
from that same ring, a few permit one to query properties of the ring.
The type of a raw value is RingBase::RawValue
, which helpfully
abbreviates to RawValue inside the namespace of RingBase
. Wherever
possible the concrete implementations should be exception safe, i.e. they
should offer either the strong exception guarantee or the no-throw
guarantee (according to the definitions in Exceptional C++ by Sutter).
Bugs, Shortcomings and other ideas
I have chosen not to use operator^
for computing powers because of a
significant risk of misunderstanding between programmer and compiler.
The syntax/grammar of C++ cannot be changed, and operator^
binds less
tightly than (binary) operator*
, so any expression of the form a*b^c
will be parsed as (a*b)^c
; this is almost certainly not what the
programmer intended. To avoid such problems of misunderstanding I
have preferred not to define operator^
; it seems too dangerous.
Note about comparison operators (<,<=,>,>=, and !=). The C++ STL
does have templates which will define all the relational operators
efficiently assuming the existence of operator<
and operator==
.
These are defined in the namespace std::rel_ops
in the standard
header file <utility>
. I have chosen NOT to use these because they can
define only homogeneous comparisons; so the comparisons between
ConstRefRingElem
and int
or BigInt
would still have to be written out
manually, and I prefer the symmetry of writing them all out.
See p.69ff of Josuttis for details.
The function myAssignZero
was NECESSARY because myAssign(x, 0)
was
ambiguous (ambiguated by the assignment from an mpz_t
). It is no longer
necessary, but I prefer to keep it (for the time being).
The requirement to use the type ConstRefRingElem
for function arguments
(which should normally be const RingElem&
is not ideal, but it seems hard
to find a better way. It is not nice to expect users to use a funny type
for their function arguments. How else could I implement (noncopying) access to
coefficients in a polynomial via an iterator, or access to matrix elements?
Would we want ++ and -- operators for RingElem
s???
Should (some of) the query functions return bool3
values?
What about properties which are hard to determine?
How to generate random elements from a ring?
Anna thinks that NearestInteger
could handle specially elements of
RingZZ
rather than doing the full wasteful computation. Not sure
if the extra code and complication would really make a difference in
practice.
gcd
and lcm
: there is no guarantee on sign/monic because it may
be costly to compute and generally useless.
Main changes
2013
- May (v0.9953):
- added IsZeroDivisor -