Thanks a lot for your quick help, yes I was not quite concentrated, in
the model GDMP one has
E : DirectProductCategory(#vl, NonNegativeInteger)
hence an argument of the domain constructor, while I had
E ==> TermDisjunctiveNormalForm(#vl)
hence only at runtime, also it is very subtle, that one can not specify such a
value in this way.
After having added the missing functions, the next problem immediately occured:
(24) -> monomial(1,t)$D3
>> System error:
The value
2
is not of type
CONS
Now, I tried to eliminate all the concrete 2's, adding superfluous parameter p:
PositiveInteger (and the problem, that one
can now instantiate with wrong numbers!) Which did not help at all: again:
(22) -> d : D3 := monomial(1,t)
>> System error:
The value
2
is not of type
CONS
My next suspection was my "strange" implementation with
OnePointCompletion for TermDisjunctiveNormalForm. But that is ok, as it
works by using the domain
PolynomialRIng was perfect well:
24) -> P := PolynomialRing(EV, T)
(24)
PolynomialRing(PRestrictedInteger(2),TermDisjunctiveNormalForm(2,3))
Type: Type
(25) -> monomial(1,t)$P
(25) [0,1,1]
Type:
PolynomialRing(PRestrictedInteger(2),TermDisjunctiveNormalForm(2,3))
This kind of problems really make it sometimes annoying in coding in
AXIOM/FriCAS.
I basically modified (for a more general use) three polynomial
categories to have more general exponents (p-restricted integers) and
more or less coded
a different output form and then one has to fight hours with these kind
of problems.
I send the actual files and hope for another help. In any case thanks a lot!
Am 19.09.17 um 16:17 schrieb oldk1331:
> I think the bug happens here:
>
> DisjunctiveNormalForm(vl: List Symbol) : public == private where
> EV ==> PRestrictedInteger(2)
> E ==> TermDisjunctiveNormalForm(#vl)
>
> I think the "#vl" is a runtime value and compiler can't handle that,
> change it to:
>
> DisjunctiveNormalForm(vl: List Symbol, dim : NonNegativeInteger) :
> public == private where
> EV ==> PRestrictedInteger(2)
> E ==> TermDisjunctiveNormalForm(dim)
>
> And in test file use "D3:=DisjunctiveNormalForm(vl, #vl)".
> Then this bug disappears and runs into a missing implementation.
>
> d : D3 := monomial(1,t)
>
> Function: ?.? : (%, Integer) -> PRestrictedInteger(2) is missing from
> domain: TermDisjunctiveNormalForm(3)
> Internal Error
> The function elt with signature (PRestrictedInteger 2)$(Integer) is
> missing from domain TermDisjunctiveNormalForm3
>
--
Mit freundlichen Grüßen
Johannes Grabmeier
Prof. Dr. Johannes Grabmeier
Köckstraße 1, D-94469 Deggendorf
Tel. +49-(0)-991-2979584, Tel. +49-(0)-151-681-70756
Tel. +49-(0)-991-3615-141 (d), Fax: +49-(0)-32224-192688
--
You received this message because you are subscribed to the Google Groups
"FriCAS - computer algebra system" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To post to this group, send email to [email protected].
Visit this group at https://groups.google.com/group/fricas-devel.
For more options, visit https://groups.google.com/d/optout.
)if false
)abbrev domain ORDCOMP OrderedCompletion
++ Completion with + and - infinity
++ Author: Manuel Bronstein
++ Description: Adjunction of two real infinites quantities to a set.
++ Date Created: 4 Oct 1989
OrderedCompletion(R : SetCategory) : Exports == Implementation where
B ==> Boolean
Exports ==> Join(SetCategory, FullyRetractableTo R) with
plusInfinity : () -> % ++ plusInfinity() returns +infinity.
minusInfinity : () -> % ++ minusInfinity() returns -infinity.
finite? : % -> B
++ finite?(x) tests if x is finite.
infinite? : % -> B
++ infinite?(x) tests if x is +infinity or -infinity,
whatInfinity : % -> SingleInteger
++ whatInfinity(x) returns 0 if x is finite,
++ 1 if x is +infinity, and -1 if x is -infinity.
if R has AbelianMonoid then
_+ : (%, %) -> %
if R has AbelianGroup then
_- : % -> %
if R has OrderedSet then OrderedSet
if R has IntegerNumberSystem then
rational? : % -> Boolean
++ rational?(x) tests if x is a finite rational number.
rational : % -> Fraction Integer
++ rational(x) returns x as a finite rational number.
++ Error: if x cannot be so converted.
rationalIfCan : % -> Union(Fraction Integer, "failed")
++ rationalIfCan(x) returns x as a finite rational number if
++ it is one and "failed" otherwise.
Implementation ==> add
Rep := Union(fin : R, inf : B) -- true = +infinity, false = -infinity
coerce(r : R) : % == [r]
retract(x:%):R == (x case fin => x.fin; error "Not finite")
finite? x == x case fin
infinite? x == x case inf
plusInfinity() == [true]
minusInfinity() == [false]
retractIfCan(x:%):Union(R, "failed") ==
x case fin => x.fin
"failed"
coerce(x : %) : OutputForm ==
x case fin => (x.fin)::OutputForm
e := 'infinity::OutputForm
x.inf => empty() + e
- e
whatInfinity x ==
x case fin => 0
x.inf => 1
-1
x = y ==
x case inf =>
y case inf => not xor(x.inf, y.inf)
false
y case inf => false
x.fin = y.fin
if R has AbelianGroup then
- x ==
x case inf => [not(x.inf)]
[- (x.fin)]
if R has AbelianMonoid then
x + y ==
x case inf =>
y case fin => x
xor(x.inf, y.inf) => error "Undefined sum"
x
y case inf => y
[x.fin + y.fin]
if R has OrderedSet then
x < y ==
x case inf =>
y case inf =>
xor(x.inf, y.inf) => y.inf
false
not(x.inf)
y case inf => y.inf
x.fin < y.fin
if R has IntegerNumberSystem then
rational? x == finite? x
rational x == rational(retract(x)@R)
rationalIfCan x ==
(r:= retractIfCan(x)@Union(R,"failed")) case "failed" =>"failed"
rational(r::R)
)abbrev package ORDCOMP2 OrderedCompletionFunctions2
++ Lifting of maps to ordered completions
++ Author: Manuel Bronstein
++ Description: Lifting of maps to ordered completions.
++ Date Created: 4 Oct 1989
OrderedCompletionFunctions2(R, S) : Exports == Implementation where
R, S : SetCategory
ORR ==> OrderedCompletion R
ORS ==> OrderedCompletion S
Exports ==> with
map : (R -> S, ORR) -> ORS
++ map(f, r) lifts f and applies it to r, assuming that
++ f(plusInfinity) = plusInfinity and that
++ f(minusInfinity) = minusInfinity.
map : (R -> S, ORR, ORS, ORS) -> ORS
++ map(f, r, p, m) lifts f and applies it to r, assuming that
++ f(plusInfinity) = p and that f(minusInfinity) = m.
Implementation ==> add
map(f, r) == map(f, r, plusInfinity(), minusInfinity())
map(f, r, p, m) ==
zero?(n := whatInfinity r) => (f retract r)::ORS
-- one? n => p
(n = 1) => p
m
)endif
)abbrev domain ONECOMP OnePointCompletion
++ Completion with infinity
++ Author: Manuel Bronstein
++ revisited and extended: Johannes Grabmeier
++ Description: Adjunction of a complex infinity to a set.
++ Date Created: 4 Oct 1989
OnePointCompletion(R : SetCategory) : Exports == Implementation where
B ==> Boolean
Exports ==> Join(SetCategory, FullyRetractableTo R) with
infinity : () -> %
++ infinity() returns infinity.
finite? : % -> B
++ finite?(x) tests if x is finite.
infinite? : % -> B
++ infinite?(x) tests if x is infinite.
if R has IntegerNumberSystem then
rational? : % -> Boolean
++ rational?(x) tests if x is a finite rational number.
rational : % -> Fraction Integer
++ rational(x) returns x as a finite rational number.
++ Error: if x is not a rational number.
rationalIfCan : % -> Union(Fraction Integer, "failed")
++ rationalIfCan(x) returns x as a finite rational number if
++ it is one, "failed" otherwise.
if R has AbelianSemiGroup then AbelianSemiGroup
if R has AbelianMonoid then AbelianMonoid
++ this is cheating, if there is a 0 in R, then 0 in % is
++ in competition with infinity, but 0+infinity = 0
if R has Monoid then Monoid
Implementation ==> add
Rep := Union(R, "infinity")
------------------------------------------------------------------------
-- from BasicType
------------------------------------------------------------------------
((x : %) = (y : %)): Boolean ==
x case "infinity" => y case "infinity"
y case "infinity" => false
x :: R = y :: R
------------------------------------------------------------------------
-- from CoercibleToOutputForm
------------------------------------------------------------------------
coerce(x : %) : OutputForm ==
--x case "infinity" => 'infinity :: OutputForm
-- uses UTF-code 16-8734 (+1) as index starts with 1
x case "infinity" => index(8735)$Character :: OutputForm
x :: R :: OutputForm
------------------------------------------------------------------------
-- from SetCategory
------------------------------------------------------------------------
-- missing: hashUpdate!(hs: HashState, x: %): HashState
------------------------------------------------------------------------
-- from FullyRetractableTo R
------------------------------------------------------------------------
coerce(r : R) : % == r
retract(x : %) : R == (x case R => x :: R; error "retract: rot finite")
retractIfCan(x : %): Union(R, "failed") == (x case R => x :: R; "failed")
------------------------------------------------------------------------
-- from with
------------------------------------------------------------------------
finite?(x : %): Boolean == x case R
infinite?(x : %): Boolean == x case "infinity"
infinity(): % == "infinity"
if R has IntegerNumberSystem then
rational?(x : %): Boolean == finite? x
rational(x : %): Fraction Integer == rational(retract(x)@R)
rationalIfCan(x : %): Union(Fraction Integer, "failed") ==
(r:= retractIfCan(x)@Union(R,"failed")) case "failed" =>"failed"
rational(r::R)
if R has AbelianSemiGroup then
++ commutative and associative extension by a+infinity = infinity+a = a
++ and infinity+infinity=infinity, used
((x : %) + (y : %)): % ==
x case "infinity" => y
y case "infinity" => x
x :: R + y :: R
if R has AbelianMonoid then
0 == 0$R
if R has Monoid then
1 == 1$R
((x : %) * (y : %)): % ==
x case "infinity" => infinity()
y case "infinity" => infinity()
x :: R * y :: R
)if false
)abbrev package ONECOMP2 OnePointCompletionFunctions2
++ Lifting of maps to one-point completions
++ Author: Manuel Bronstein
++ Description: Lifting of maps to one-point completions.
++ Date Created: 4 Oct 1989
OnePointCompletionFunctions2(R, S) : Exports == Implementation where
R, S : SetCategory
OPR ==> OnePointCompletion R
OPS ==> OnePointCompletion S
Exports ==> with
map : (R -> S, OPR) -> OPS
++ map(f, r) lifts f and applies it to r, assuming that
++ f(infinity) = infinity.
map : (R -> S, OPR, OPS) -> OPS
++ map(f, r, i) lifts f and applies it to r, assuming that
++ f(infinity) = i.
Implementation ==> add
map(f, r) == map(f, r, infinity())
map(f, r, i) ==
(u := retractIfCan r) case R => (f(u::R))::OPS
i
)abbrev package INFINITY Infinity
++ Top-level infinity
++ Author: Manuel Bronstein
++ Description: Default infinity signatures for the interpreter;
++ Date Created: 4 Oct 1989
Infinity() : with
infinity : () -> OnePointCompletion Integer
++ infinity() returns infinity.
plusInfinity : () -> OrderedCompletion Integer
++ plusInfinity() returns plusIinfinity.
minusInfinity : () -> OrderedCompletion Integer
++ minusInfinity() returns minusInfinity.
== add
infinity() == infinity()$OnePointCompletion(Integer)
plusInfinity() == plusInfinity()$OrderedCompletion(Integer)
minusInfinity() == minusInfinity()$OrderedCompletion(Integer)
--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
--All rights reserved.
--
--Redistribution and use in source and binary forms, with or without
--modification, are permitted provided that the following conditions are
--met:
--
-- - Redistributions of source code must retain the above copyright
-- notice, this list of conditions and the following disclaimer.
--
-- - Redistributions in binary form must reproduce the above copyright
-- notice, this list of conditions and the following disclaimer in
-- the documentation and/or other materials provided with the
-- distribution.
--
-- - Neither the name of The Numerical ALgorithms Group Ltd. nor the
-- names of its contributors may be used to endorse or promote products
-- derived from this software without specific prior written permission.
--
--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
)endif
p := 2
vl := [p1,p2,p3]
n := #vl
T := TermDisjunctiveNormalForm(p,n)
0$T
EV := PRSTINT(2)
E := PRDPINT(2,n)
t : T := directProduct vector [0,1,1]
t' : T := directProduct vector [1,-1,0]
t+t
t'+t'
t+t'
sup(t, t')
P := PolynomialRing(EV, T)
monomial(1,t)$P
D3:=DisjunctiveNormalForm(vl, p, #vl)
u : D3 := 0
d : D3 := monomial(1,t)
)fin
u : D := 1
)r testDPRMP
--F := PF 2
F := PRestricted 2
vl : List Symbol := ['p,'q,'r] --[Kpa, Kdu, Kde, Kbl, Ppa, Pdu, Pde, Pbl, Bpa,
Bdu, Bde, Bbl, Dpa, Ddu, Dde, Dbl]
--P := DNF([Kpa, Kdu, Kde, Kbl, Ppa, Pdu, Pde, Pbl, Bpa, Bdu, Bde, Bbl, Dpa,
Ddu, Dde, Dbl], F)
)read i-funsel.boot
)co iml
)co prdp
)co prdp
)co prdp
)r testDPRMP
)sh UnittestCount
)sh UnittestAux
)sh Unittest
t : T := [1,1,1] pretend T
t : T := [1,1,-1] pretend T
)r testDPRMP
--F := PF 2
F := PRestricted 2
vl : List Symbol := ['p,'q,'r] --[Kpa, Kdu, Kde, Kbl, Ppa, Pdu, Pde, Pbl, Bpa,
Bdu, Bde, Bbl, Dpa, Ddu, Dde, Dbl]
--P := DNF([Kpa, Kdu, Kde, Kbl, Ppa, Pdu, Pde, Pbl, Bpa, Bdu, Bde, Bbl, Dpa,
Ddu, Dde, Dbl], F)
T := TermDisjunctiveNormalForm(2,#vl)
s : T := random()
)r testDPRMP
--F := PF 2
F := PRestricted 2
vl : List Symbol := ['p,'q,'r] --[Kpa, Kdu, Kde, Kbl, Ppa, Pdu, Pde, Pbl, Bpa,
Bdu, Bde, Bbl, Dpa, Ddu, Dde, Dbl]
--P := DNF([Kpa, Kdu, Kde, Kbl, Ppa, Pdu, Pde, Pbl, Bpa, Bdu, Bde, Bbl, Dpa,
Ddu, Dde, Dbl], F)
T := TermDisjunctiveNormalForm(2,#vl)
s : T := random()
--3*s
t : T := random()
)r test
--F := PF 2
F := PRestricted 2
vl : List Symbol := ['p,'q,'r] --[Kpa, Kdu, Kde, Kbl, Ppa, Pdu, Pde, Pbl, Bpa,
Bdu, Bde, Bbl, Dpa, Ddu, Dde, Dbl]
--P := DNF([Kpa, Kdu, Kde, Kbl, Ppa, Pdu, Pde, Pbl, Bpa, Bdu, Bde, Bbl, Dpa,
Ddu, Dde, Dbl], F)
T := TermDisjunctiveNormalForm(2,#vl)
s : T := random()
--3*s
t : T := random()
--t*(-5) -- should - result in [-t1, -t2, ...]?
s,t
)r testDPRMP
--F := PF 2
F := PRestricted 2
vl : List Symbol := ['p,'q,'r] --[Kpa, Kdu, Kde, Kbl, Ppa, Pdu, Pde, Pbl, Bpa,
Bdu, Bde, Bbl, Dpa, Ddu, Dde, Dbl]
--P := DNF([Kpa, Kdu, Kde, Kbl, Ppa, Pdu, Pde, Pbl, Bpa, Bdu, Bde, Bbl, Dpa,
Ddu, Dde, Dbl], F)
T := TermDisjunctiveNormalForm(2,#vl)
s : T := random()
--3*s
t : T := random()
--t*(-5) -- should - result in [-t1, -t2, ...]?
s,t
)r testDPRMP
--F := PF 2
F := PRestricted 2
vl : List Symbol := ['p,'q,'r] --[Kpa, Kdu, Kde, Kbl, Ppa, Pdu, Pde, Pbl, Bpa,
Bdu, Bde, Bbl, Dpa, Ddu, Dde, Dbl]
--P := DNF([Kpa, Kdu, Kde, Kbl, Ppa, Pdu, Pde, Pbl, Bpa, Bdu, Bde, Bbl, Dpa,
Ddu, Dde, Dbl], F)
T := TermDisjunctiveNormalForm(2,#vl)
s : T := random()
--3*s
t : T := random()
--t*(-5) -- should - result in [-t1, -t2, ...]?
s,t
)r testDPRMP
--F := PF 2
F := PRestricted 2
vl : List Symbol := ['p,'q,'r] --[Kpa, Kdu, Kde, Kbl, Ppa, Pdu, Pde, Pbl, Bpa,
Bdu, Bde, Bbl, Dpa, Ddu, Dde, Dbl]
--P := DNF([Kpa, Kdu, Kde, Kbl, Ppa, Pdu, Pde, Pbl, Bpa, Bdu, Bde, Bbl, Dpa,
Ddu, Dde, Dbl], F)
T := TermDisjunctiveNormalForm(2,#vl)
s : T := random()
--3*s
t : T := random()
--t*(-5) -- should - result in [-t1, -t2, ...]?
s,t
)r testDPRMP
--F := PF 2
F := PRestricted 2
vl : List Symbol := ['p,'q,'r] --[Kpa, Kdu, Kde, Kbl, Ppa, Pdu, Pde, Pbl, Bpa,
Bdu, Bde, Bbl, Dpa, Ddu, Dde, Dbl]
--P := DNF([Kpa, Kdu, Kde, Kbl, Ppa, Pdu, Pde, Pbl, Bpa, Bdu, Bde, Bbl, Dpa,
Ddu, Dde, Dbl], F)
T := TermDisjunctiveNormalForm(2,#vl)
s : T := random()
--3*s
t : T := random()
--t*(-5) -- should - result in [-t1, -t2, ...]?
s,t
)r testDPRMP
--F := PF 2
F := PRestricted 2
vl : List Symbol := ['p,'q,'r] --[Kpa, Kdu, Kde, Kbl, Ppa, Pdu, Pde, Pbl, Bpa,
Bdu, Bde, Bbl, Dpa, Ddu, Dde, Dbl]
--P := DNF([Kpa, Kdu, Kde, Kbl, Ppa, Pdu, Pde, Pbl, Bpa, Bdu, Bde, Bbl, Dpa,
Ddu, Dde, Dbl], F)
T := TermDisjunctiveNormalForm(2,#vl)
s : T := random()
--3*s
t : T := random()
--t*(-5) -- should - result in [-t1, -t2, ...]?
s,t
--sup(s,t)
--s-t
s' : T := directProduct [0,1,1]
t' : T := directProduct [1,0,-1]
and'(s', t')
s'' : T := directProduct [-1,0,1]
)r testDPRMP
--F := PF 2
F := PRestricted 2
vl : List Symbol := ['p,'q,'r] --[Kpa, Kdu, Kde, Kbl, Ppa, Pdu, Pde, Pbl, Bpa,
Bdu, Bde, Bbl, Dpa, Ddu, Dde, Dbl]
--P := DNF([Kpa, Kdu, Kde, Kbl, Ppa, Pdu, Pde, Pbl, Bpa, Bdu, Bde, Bbl, Dpa,
Ddu, Dde, Dbl], F)
T := TermDisjunctiveNormalForm(2,#vl)
s : T := random()
--3*s
t : T := random()
--t*(-5) -- should - result in [-t1, -t2, ...]?
s,t
--sup(s,t)
--s-t
s' : T := directProduct [0,1,1]
t' : T := directProduct [1,0,-1]
--and'(s', t')
s'' : T := directProduct [-1,0,1]
--and'(s', s'')
P := DNF(vl, F, T)
u : P := monomial(1, t) + monomial(1,s)
tP := t :: P
u + tP
)r testDPRMP
--F := PF 2
F := PRestricted 2
vl : List Symbol := ['p,'q,'r] --[Kpa, Kdu, Kde, Kbl, Ppa, Pdu, Pde, Pbl, Bpa,
Bdu, Bde, Bbl, Dpa, Ddu, Dde, Dbl]
--P := DNF([Kpa, Kdu, Kde, Kbl, Ppa, Pdu, Pde, Pbl, Bpa, Bdu, Bde, Bbl, Dpa,
Ddu, Dde, Dbl], F)
T := TermDisjunctiveNormalForm(2,#vl)
s : T := random()
--3*s
t : T := random()
--t*(-5) -- should - result in [-t1, -t2, ...]?
s,t
--sup(s,t)
--s-t
s' : T := directProduct [0,1,1]
t' : T := directProduct [1,0,-1]
--and'(s', t')
s'' : T := directProduct [-1,0,1]
--and'(s', s'')
P := DNF(vl, F, T)
u : P := monomial(1, t) + monomial(1,s)
tP := t :: P
u + tP
)r testDPRMP
--F := PF 2
F := PRestricted 2
vl : List Symbol := ['p,'q,'r] --[Kpa, Kdu, Kde, Kbl, Ppa, Pdu, Pde, Pbl, Bpa,
Bdu, Bde, Bbl, Dpa, Ddu, Dde, Dbl]
--P := DNF([Kpa, Kdu, Kde, Kbl, Ppa, Pdu, Pde, Pbl, Bpa, Bdu, Bde, Bbl, Dpa,
Ddu, Dde, Dbl], F)
T := TermDisjunctiveNormalForm(2,#vl)
s : T := random()
--3*s
t : T := random()
--t*(-5) -- should - result in [-t1, -t2, ...]?
s,t
--sup(s,t)
--s-t
s' : T := directProduct [0,1,1]
t' : T := directProduct [1,0,-1]
--and'(s', t')
s'' : T := directProduct [-1,0,1]
--and'(s', s'')
P := DNF(vl, F, T)
u : P := monomial(1, t) + monomial(1,s)
)r testDPRMP
--F := PF 2
F := PRestricted 2
vl : List Symbol := ['p,'q,'r] --[Kpa, Kdu, Kde, Kbl, Ppa, Pdu, Pde, Pbl, Bpa,
Bdu, Bde, Bbl, Dpa, Ddu, Dde, Dbl]
--P := DNF([Kpa, Kdu, Kde, Kbl, Ppa, Pdu, Pde, Pbl, Bpa, Bdu, Bde, Bbl, Dpa,
Ddu, Dde, Dbl], F)
T := TermDisjunctiveNormalForm(2,#vl)
s : T := random()
--3*s
t : T := random()
--t*(-5) -- should - result in [-t1, -t2, ...]?
s,t
--sup(s,t)
--s-t
s' : T := directProduct [0,1,1]
t' : T := directProduct [1,0,-1]
--and'(s', t')
s'' : T := directProduct [-1,0,1]
--and'(s', s'')
P := DNF(vl, F, T)
u : P := monomial(1, t) + monomial(1,s)
w := tP := t :: P + tP
)co prdp
)r testDPRMP
)r testDPRMP
)r testDPRMP
)r testDPRMP
)r testDPRMP
)r testDPRMP
)co prdp
)r testDPRMP
)co prdp
)r testDPRMP
)r testDPRMP
)r testDPRMP
)r testDPRMP
)r testDPRMP
allVariables()$P
allVariables()$P
p : P := %.1
)r testDPRMP
--F := PF 2
F := PRestricted 2
vl : List Symbol := ['p,'q,'r] --[Kpa, Kdu, Kde, Kbl, Ppa, Pdu, Pde, Pbl, Bpa,
Bdu, Bde, Bbl, Dpa, Ddu, Dde, Dbl]
--P := DNF([Kpa, Kdu, Kde, Kbl, Ppa, Pdu, Pde, Pbl, Bpa, Bdu, Bde, Bbl, Dpa,
Ddu, Dde, Dbl], F)
T := TermDisjunctiveNormalForm(2,#vl)
s : T := random()
--3*s
t : T := random()
--t*(-5) -- should - result in [-t1, -t2, ...]?
s,t
--sup(s,t)
--s-t
s' : T := directProduct [0,1,1]
t' : T := directProduct [1,0,-1]
--and'(s', t')
s'' : T := directProduct [-1,0,1]
--and'(s', s'')
P := DNF(vl, F, T)
u : P := monomial(1, t) + monomial(1,s)
w := tP := t :: P + tP
)r testDPRMP
--F := PF 2
F := PRestricted 2
vl : List Symbol := ['p,'q,'r] --[Kpa, Kdu, Kde, Kbl, Ppa, Pdu, Pde, Pbl, Bpa,
Bdu, Bde, Bbl, Dpa, Ddu, Dde, Dbl]
--P := DNF([Kpa, Kdu, Kde, Kbl, Ppa, Pdu, Pde, Pbl, Bpa, Bdu, Bde, Bbl, Dpa,
Ddu, Dde, Dbl], F)
T := TermDisjunctiveNormalForm(2,#vl)
s : T := random()
--3*s
t : T := random()
--t*(-5) -- should - result in [-t1, -t2, ...]?
s,t
--sup(s,t)
--s-t
s' : T := directProduct [0,1,1]
t' : T := directProduct [1,0,-1]
--and'(s', t')
s'' : T := directProduct [-1,0,1]
--and'(s', s'')
P := DNF(vl, F, T)
u : P := monomial(1, t) + monomial(1,s)
w := tP := t :: P + tP
aV := allVariables()$P
p : P := aV.1
q : P := aV.2
p : P := %.1
p : P := %.1
)r testDPRMP
)r testDPRMP
t
)r testDPRMP
--F := PF 2
F := PRestricted 2
vl : List Symbol := ['p,'q,'r] --[Kpa, Kdu, Kde, Kbl, Ppa, Pdu, Pde, Pbl, Bpa,
Bdu, Bde, Bbl, Dpa, Ddu, Dde, Dbl]
--P := DNF([Kpa, Kdu, Kde, Kbl, Ppa, Pdu, Pde, Pbl, Bpa, Bdu, Bde, Bbl, Dpa,
Ddu, Dde, Dbl], F)
T := TermDisjunctiveNormalForm(2,#vl)
s : T := random()
--3*s
t : T := random()
--t*(-5) -- should - result in [-t1, -t2, ...]?
s,t
--sup(s,t)
--s-t
s' : T := directProduct [0,1,1]
t' : T := directProduct [1,0,-1]
--and'(s', t')
s'' : T := directProduct [-1,0,1]
--and'(s', s'')
P := DNF(vl, F, T)
u : P := monomial(1, t) + monomial(1,s)
w := tP := t :: P + tP
aV := allVariables()$P
p : P := aV.1
q : P := aV.2
r : P := aV.3
t+t
)co prdp
)co prdp
)co prdp
)co prdp
)co prdp
)co prdp
)co prdp
)co prdp
)co prdp
)r testDPRMP
p
p+p
p*p
p
not p
not' p
p*q
pp : T := [-1,0,0] pretend T
p*pp
u
u*u
u+u
s
t
s*t
s+t
)co prdp
)co prdp
t
t'
z : T := "empty"
z : T := "empty" pretend T
t+t
z+z
)set func sel on
)set mess sel on
z+z
)co prdp
)co prdp
setMessageLevel(1)$T
)set mess sel on
z+z
u
z
t
t+t
)co prdp
)set mess sel on
t+t
setMessageLevel(1)$T
t+t
)co prdp
)co prdp
)show Union(Integer, "failed")
)show Union(Integer, "empty")
)co prdp
)co prdp
)co prdp
)co prdp
)co prdp
)co prdp
setMessageLevel(1)$T
t+t
t+"empty"
t+"empty"::T
t+empty()$T
empty()$T+empty()$T
s+empty()$T
empty()$T+empty()$T
)co prdp
empty()$T+empty()$T
s+empty()$T
empty()$T+s
)co prdp
empty()$T+s
s+empty()$T
s+s
)co prdp
setMessageLevel(1)$T
s+s
)co prdp
setMessageLevel(1)$T
s+s
setMessageLevel(1)$T
s+empty()$T
empty()$T+s
)co prdp
empty()$T+s
empty()$T+s
s+empty()$T
)co prdp
s+empty()$T
empty()$T+s
s+empty()$T
)co prdp
s : T := random()
s : T := random()
s : T := random()
s : T := random()
s : T := random()
s : T := random()
s+s
)set func sel off
)set func off
)set mess sel off
s+s
)co prdp
s : T := random()
s+s
s+empty()
s+empty()$T
empty()$T+s
s+s
s : T := random()
s+s
u : T := random()
y : T := random()
-- Version: 0.0 2017-09-09
--Copyright (c) 2017-2017 Johannes Grabmeier
--All rights reserved.
--
--Redistribution and use in source and binary forms, with or without
--modification, are permitted provided that the following conditions are
--met:
--
-- - Redistributions of source code must retain the above copyright
-- notice, this list of conditions and the following disclaimer.
--
-- - Redistributions in binary form must reproduce the above copyright
-- notice, this list of conditions and the following disclaimer in
-- the documentation and/or other materials provided with the
-- distribution.
--
--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
--)abbrev domain DNF DisjunctiveNormalForm
)abbrev domain TDNF TermDisjunctiveNormalForm
++ Authors: J. Grabmeier
++ Date Created: 19.09.2017
++ Change History:
++ Version: 0.0 19.09.2017
++ Basic Operations:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ Reference:
++ Description:
++ Domain implements Exponents using p-restricted integers, i.e.
++ -1,0,1 for the disjunctive normal form
++ of Boolean functions, where internally -1 means negation.
++ Need an element infinity to represent 0 in the free module of the terms.
TermDisjunctiveNormalForm(p: PositiveInteger, n: NonNegativeInteger):
public == private where
public ==> Join(InternalMessageLevel, PRestrictedDirectProductCategory(p, n,
PRestrictedInteger(p)),
CoercibleTo Vector Integer) with
-------------------------------------------------------------------------
-- from OnePointCompletion:
-------------------------------------------------------------------------
infinity : () -> %
++ infinity() returns infinity.
finite? : % -> Boolean
++ finite?(x) tests if x is finite.
infinite? : % -> Boolean
++ infinite?(x) tests if x is infinite.
negativeUnitVector: PositiveInteger -> %
++ unitVector(p) constructs 0..010..0 with 1 at position p.
positiveUnitVector: PositiveInteger -> %
and': (%, %) -> %
private ==> OnePointCompletion PRestrictedDirectProductInteger(p, n) add
NNI ==> NonNegativeInteger
PI ==> PositiveInteger
OF ==> OutputForm
mLP ==> messageLevelPointer
msg ==> message
EV ==> PRestrictedInteger(p)
E ==> PRestrictedDirectProductInteger(p, n)
Rep := OnePointCompletion E
rep(x:%):Rep == x :: Rep
per(r:Rep):% == r :: %
import OF
-------------------------------------------------------------------------
-- constants
-------------------------------------------------------------------------
-------------------------------------------------------------------------
-- local functions
-------------------------------------------------------------------------
randomSign(): Integer ==
random()$PrimeField(p) = 0 => 1
-1
-------------------------------------------------------------------------
-- global functions
-------------------------------------------------------------------------
-------------------------------------------------------------------------
-- from AbelianGroup
-------------------------------------------------------------------------
((u: %) + (v: %)): % ==
infinite? u => v
infinite? v => u
uE : E := retract(u)
vE : E := retract(v)
for ui in entries uE for vi in entries vE repeat
((ui = 1) and (vi = -1)) or ((ui = -1) and (vi = 1)) => return
infinity()
(uE + vE) :: Rep
-------------------------------------------------------------------------
-- from OrderedAbelianMonoidSup
-------------------------------------------------------------------------
sup(u: %, v: %): % ==
infinite? u => v
infinite? v => u
uE : E := retract(u)
vE : E := retract(v)
sup(uE, vE) :: Rep
-------------------------------------------------------------------------
-- from RetractableTo DirectProduct(n,PRestrictedInteger(p))
-------------------------------------------------------------------------
retractIfCan(u: %): Union(DirectProduct(n,EV),"failed") ==
infinite? u => "failed"
retract u
-------------------------------------------------------------------------
-- from PRestrictedDirectProductInteger(p, n)
-------------------------------------------------------------------------
directProduct(v: Vector EV): % ==
empty? v => infinity()
not(#v=n) => error "directProduct: not the right vector length"
directProduct(v)$E :: Rep
elt(p: %, i: Integer): EV ==
infinite? p => error "elt: can't take components from infinity"
retract(p).i
)if false
negativeUnitVector(i : PI): % ==
v : List Integer := new(n, 0)$List(Integer)
v.i := -1
directProduct vector v
positiveUnitVector(i : PI): % ==
v : List Integer := new(n, 0)$List(Integer)
v.i := 1
directProduct vector v
and'(s : %, t : %): % ==
for si in entries s for ti in entries t repeat
si = 1 and ti = -1 => return 0
si = -1 and ti = 1 => return 0
s+t
-------------------------------------------------------------------------
-- from SetCategory
-------------------------------------------------------------------------
set(p: %): OutputForm == outputForm "u"
-------------------------------------------------------------------------
-- from Finite
-------------------------------------------------------------------------
random(): % ==
lr : Vector Integer := vector [randomSign()*random(p::Integer)$Integer
for i in 1..n]
(directProduct(lr)) :: %
--directProduct [randomSign()*random(p::Integer)$Integer for i in 1..n]
--"empty"
--lookup(x: %): PI ==
-- xI : List Integer := [xx :: Integer for xx in entries(x::Rep)]
-- if mLP 1 then print blankSeparate [msg "xRep = ", xI :: OF]
-- (1 + wholePart(wholeRadix(xI)$RadixExpansion(p))) :: PI
--index(q: PI): % ==
-- q > p^n => error "index: argument too large!"
-- lq : List NNI := [d::NNI for d in wholeRagits ((q-1) ::
RadixExpansion(p)) ]
-- directProduct vector concat( new((n-#lq)::NNI,0), lq)
-------------------------------------------------------------------------
-- from RetractableFrom DirectProduct
-------------------------------------------------------------------------
coerce(u: %): DirectProduct(n, Integer) ==
u case "empty" => error "coerce: cannot coere to DirectProduct(n,
Integer)"
u
--luplusv : Vector Integer := new(n, 0)
-------------------------------------------------------------------------
-- from CancellationAbelianMonoid
-------------------------------------------------------------------------
-------------------------------------------------------------------------
-- from Monoid
-------------------------------------------------------------------------
((u: %) * (r: Integer)): % ==
directProduct vector [reducePm1(ui*r) for ui in entries u]
((r: Integer) * (u: %)): % == u*r
-------------------------------------------------------------------------
-- from Ring
-------------------------------------------------------------------------
((u: %) * (v: %)): % ==
-- der folgende Code gibt Probleme:
-- reducePm1DP (rep(u) * rep(v))
-- das passt:
directProduct vector [reducePm1(ui*vi) for ui in entries u for vi in
entries v]
recip(z: %): Union(%, "failed") ==
pm1 : Integer := p-1
w : List Integer := new(n,0)$List(Integer)
for i in minIndex w .. maxIndex w repeat
zi : Integer := qelt(z, i) :: Integer
not gcd(zi, pm1) = 1 => return "failed"
izi : Integer := extendedEuclidean(zi, pm1).coef1
if izi < 0 then izi := izi + pm1
qsetelt!(w, i, izi)
directProduct vector w
-------------------------------------------------------------------------
-- from OrderedSet
-------------------------------------------------------------------------
((x: %) < (y: %)): Boolean ==
for i in 1..n repeat
qelt(x,i) < qelt(y,i) => return true
qelt(x,i) > qelt(y,i) => return false
false
)endif
-- SemiRng --+
-- +--- AbelianMonoidRing ---- FiniteAbelianMonoidRing --
MaybeSkewPolynomialCategory --
-- BiModule --+
)abbrev domain DNF DisjunctiveNormalForm
DisjunctiveNormalForm(vl: List Symbol, p: PositiveInteger, n:
NonNegativeInteger) : public == private where
EV ==> PRestrictedInteger(p)
E ==> TermDisjunctiveNormalForm(p, n)
OV ==> OrderedVariableList(vl)
msg ==> message
mLP ==> messageLevelPointer
OF ==> OutputForm
public == Join(PRestrictedPolynomialCategory(EV, EV, E, OV), CoercibleFrom E)
with
reorder : (%, List Integer) -> %
++ reorder(p, perm) applies the permutation perm to the variables
++ in a polynomial and returns the new correctly ordered polynomial
private == DistributedPRestrictedMultivariatePolynomial(p, vl, EV, EV, E) add
import OF
-------------------------------------------------------------------------
-- Representation
-------------------------------------------------------------------------
Term := Record(k : E, c : EV)
Rep := List Term
rep(x:%):Rep == x :: Rep
per(r:Rep):% == r :: %
-------------------------------------------------------------------------
-- local constants
-------------------------------------------------------------------------
-------------------------------------------------------------------------
-- global functions
-------------------------------------------------------------------------
-------------------------------------------------------------------------
-- from SetCategory
-------------------------------------------------------------------------
coerce(x : %) : OutputForm ==
zero?(x) => outputForm 0
l : List OutputForm := nil
lt : List OutputForm := nil
vlOF := [v::OutputForm for v in vl]
for t in reverse x repeat
l := nil
for i in 1..#vlOF repeat
t.k.i = 1 => l := cons(vlOF.i, l)
t.k.i = -1 => l := cons(overbar vlOF.i, l)
--t.k.i = 0 => "next"
--l := cons(vlOF.i ^ t.k.i ::OutputForm, l)
l := reverse l
if (t.c ~= 1) or (null l) then l := cons(t.c :: OutputForm, l)
1 = #l => lt := cons(first l, lt)
lt := cons(reduce("*",l),lt)
1 = #lt => first lt
reduce("+",lt)
)fin
-------------------------------------------------------------------------
-- modified from MaybeSkewPolynomialCategory
monomials : % -> List %
++ monomials(p) returns the list of non-zero monomials of
++ polynomial p, i.e.
++ \spad{monomials(sum(a_(i) X^(i))) = [a_(1) X^(1), ..., a_(n) X^(n)]}.
monomial : (%, OV, Integer) -> %
++ monomial(a, x, n) creates the monomial \spad{a*x^n} where \spad{a} is
++ a polynomial, x is a variable and n is an integer.
monomial : (%, List OV, List Integer) -> %
++ monomial(a, [v1..vn], [e1..en]) returns \spad{a*prod(vi^ei)}
totalDegree : % -> Integer
++ totalDegree(p) returns the largest sum over all monomials
++ of all exponents of a monomial.
variables : % -> List OV
++ variables(p) returns the list of those variables actually
++ appearing in the polynomial p.
-------------------------------------------------------------------------
reorder: (%, List Integer) -> %
++ reorder(p, perm) applies the permutation perm to the variables
++ in a polynomial and returns the new correctly ordered polynomial
allVariables: () -> List %
++ allVariables() returns all boolean variables.
-------------------------------------------------------------------------
not': E -> %
not': % -> %
++ not'(p) negation of the boolean expression p
or': (E, E) -> %
or': (%, %) -> %
private == PolynomialRing(R, E) add
import OF
-------------------------------------------------------------------------
-- Representation
-------------------------------------------------------------------------
Term := Record(k : E, c : R)
Rep := List Term
rep(x:%):Rep == x :: Rep
per(r:Rep):% == r :: %
-------------------------------------------------------------------------
-- local constants
-------------------------------------------------------------------------
n := #vl
zero?(p : %) : Boolean == null(p :: Rep)
totalDegree(p : %): Integer ==
zero? p => 0
"max"/[reduce("+",(t.k)::(Vector Integer), 0) for t in p]
-------------------------------------------------------------------------
-- modified from MaybeSkewPolynomialCategory
-------------------------------------------------------------------------
monomials(p : %): List % ==
-- sequential version for efficiency, by WMSIT, 7/30/90
ml := empty$List(%)
while p ~= 0 repeat
ml := concat(leadingMonomial p, ml)
p := reductum p
reverse! ml
monomial(p: %, lv: List OV, ln: List Integer): % ==
empty? lv =>
empty? ln => p
error "mismatched lists in monomial"
empty? ln => error "mismatched lists in monomial"
monomial(monomial(p, first lv, first ln), rest lv, rest ln)
-------------------------------------------------------------------------
-- modified from GeneralDistributedMultivariatePolynomial
-------------------------------------------------------------------------
monomial(p : %, v : OV, e : Integer) : % ==
locv := lookup v
p*monomial(1,
directProduct [if z = locv then e else 0 for z in
1..n]$Vector(Integer))
allVariables(): List % == [monomial(1, index(i::PositiveInteger)$OV,1) for
i in 1..#vl]
-------------------------------------------------------------------------
-- Boolean Functions
-------------------------------------------------------------------------
not'(t: E): % ==
lm : List % := []
uV : % := 0
for i in 1..n for ti in entries t repeat
if not zero? ti then
if ti > 0
then 1+1
--then uV := monomial(1, negativeUnitVector(i :: PositiveInteger))
--else uV := monomial(1, positiveUnitVector(i :: PositiveInteger))
lm := cons(uV, lm)
reduce("+", lm, 0)
not'(u: %): % ==
reduce("*", [not'(degree m) for m in monomials u])
or'(u : E, v: E): % ==
lu : List Integer := entries u
lv : List Integer := entries v
if mLP 1 then print blankSeparate [msg "u = ", lu :: OF]
if mLP 1 then print blankSeparate [msg "v = ", lv :: OF]
for ui in lu for vi in entries v for i in 1..n repeat
if mLP 1 then print blankSeparate [ui :: OF, vi :: OF, i :: OF]
if not (ui = vi) then break
if mLP 1 then print blankSeparate [msg "i = ", i :: OF]
i > n => u :: % -- as u = v
for ui in rest(lu, i) for vi in rest(entries v, i) for j in i+1..n repeat
if mLP 1 then print blankSeparate [ui :: OF, vi :: OF, j :: OF]
if not (ui = vi) then break
if mLP 1 then print blankSeparate [msg "j = ", j :: OF]
j <= n => (u :: %) + (v :: %)
-- know: j = n+1
1 <= i and i <= n =>
lu.i := 0
directProduct(vector(lu))$E :: %
or'(u : %, v: %): % ==
w : % := 0
for mu in monomials u repeat
for mv in monomials v repeat
w := w+or'(degree mu, degree mv)
w
-------------------------------------------------------------------------
-- CoercibleFrom E
-------------------------------------------------------------------------
coerce(t : E) : % == monomial(1,t)
-- coerce(v : OV) : % == monomial(1, v, 1)
-- listCoef(p : %) : List R ==
-- rec : Term
-- [rec.c for rec in (p :: Rep)]
-- te : List Integer := entries t
-- lE : List % := [monomial(1, -c*unitVector(i :: PositiveInteger)$E) for i
in 1..#vl for c in te]
-- reduce("+", lE)
--)endif
-- Copyright (c) 2017 J. Grabmeier
-- All rights reserved.
-- prestricted numbers:
-- )abbrev category PRNS PRestrictedNumberSystem
-- )abbrev domain PRSTINT PRestrictedInteger pINT
-- )abbrev domain PRSTNNI PRestrictedNonNegativeInteger pNNI
-- prestricted exponents:
-- )abbrev category PRDP PRestrictedDirectProductCategory
pINT^n
-- )abbrev domain PRDPINT PRestrictedDirectProductInteger
pINT^n
-- )abbrev domain PRDPNNI PRestrictedDirectProductNonNegativeInteger
pNNI^n
-- generalized polynomial categories: use pINT and pNNI instaed of INT and NNI
-- and pINT^n instead of INT^n resp. pNNI^n and NNI^n
-- many functions are commented out, reason often: univariate, so need
univariate
-- p-restricted polynomials, not done yet!
-- )abbrev category GPRPCAT MaybeSkewPRestrictedPolynomialCategory
-- )abbrev category PRPCAT PRestrictedPolynomialCategory
-- )abbrev domain GDPRMP
GeneralDistributedPRestrictedMultivariatePolynomial
-- )abbrev domain DPRMP DistributedPRestrictedMultivariatePolynomial
-- )abbrev domain HDMP HomogeneousDistributedMultivariatePolynomial
-- )abbrev domain DPRMP DistributedPRestrictedMultivariatePolynomial
)abbrev category PRNS PRestrictedNumberSystem
++ Author: J. Grabmeier
++ Date Created: 19.07.2017
++ Change History:
++ Basic Functions: pRestriction
++ Related Constructors: PRestrictedInteger, PRestrictedNonNegativeInteger
++ Also See:
++ AMS Classifications:
++ Keywords: p-restriction
++ References:
++ Description: the category provides basic setting for p-restricted integers,
in particular
++ SemiRing (not Ring as NNI has not Ring)
PRestrictedNumberSystem(p: PositiveInteger): Category ==
Join(SemiRing, CommutativeStar, OrderedAbelianMonoidSup,
CoercibleFrom NonNegativeInteger, ConvertibleTo InputForm ) with
pRestriction: () -> PositiveInteger
++ pRestriction() returns p
add
pRestriction(): PositiveInteger == p
)abbrev domain PRSTINT PRestrictedInteger
++ Authors: J. Grabmeier
++ Date Created: 19.07.2016
++ Change History:
++ Version: 0.1 19.07.2017 new setting, abbreviation is PRSTINT, was PR
++ Version: 0.0 in prdp.spad
++ Basic Operations:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ Reference:
++ Description:
++ Domain implements p-restricted integers, i.e. -p+1,...,-1,0,1,2,...,p-1
++ all larger and smaller results of multiplication and addition of the
++ corresponding integers are reduced modulo p-1. This is the additive
++ version of multiplicative terms x^i implementing the relation x^p=x,
x^0=1.
++ It is used to have sympolic elements of fields of characteristic p and
++ in particular for the distributive normal forms of Boolean expressions
++ where p = 2 and x^(-1) represents the negation of x.
PRestrictedInteger(p: PositiveInteger): Exports == Implementation where
PI ==> PositiveInteger
NNI ==> NonNegativeInteger
I ==> Integer
OF ==> OutputForm
mLP ==> messageLevelPointer
msg ==> message
Exports ==> Join(PRestrictedNumberSystem p, IntegralDomain, CoercibleTo
Integer,
CoercibleFrom Integer)
Implementation ==> Integer add
--------------------------------------------------------------------------------
-- import of domains and packages
--------------------------------------------------------------------------------
import OF
--------------------------------------------------------------------------------
-- global variables
--------------------------------------------------------------------------------
Rep := Integer
rep(s : %) : Rep == s @ Rep
per(rp : Rep) : % == rp @ %
--------------------------------------------------------------------------------
-- implementation of local functions
--------------------------------------------------------------------------------
normalize(n: Integer): Integer ==
++ if n >= p x^p = x is equivalent to x^n = x^r if n rem (p-1) = r and
++ r > 0, in case r=0, we have to keep x^(p-1).
n = 0 => 0
n < 0 => - normalize(-n)
p = 2 => 1
n < p => n
r : Integer := n rem (p-1)
r = 0 => (p-1) --:: Integer
r
-------------------------------------------------------------------------
-- implementation of global functions
-------------------------------------------------------------------------
-------------------------------------------------------------------------
-- from CoercibleTo OutputForm
-------------------------------------------------------------------------
-- the following code should be superfluous, if we can guarantee
-- that only the numbers -p+1 ... p-1 are used as internal representatives!
coerce(u: %): OutputForm == coerce(normalize(rep u))$Integer
-------------------------------------------------------------------------
-- from SemiGroup
-------------------------------------------------------------------------
((u: %) + (v: %)): % == per normalize (rep(u) + rep(v))
-------------------------------------------------------------------------
-- from Magma
-------------------------------------------------------------------------
((u: %) * (v: %)): % == per normalize (rep(u) * rep(v))
-------------------------------------------------------------------------
-- from NonAssociativeRing
-------------------------------------------------------------------------
coerce(i: Integer): % == normalize i
coerce(i: NonNegativeInteger): % == normalize i::Integer
-------------------------------------------------------------------------
-- from CoercibleTo NNI
-------------------------------------------------------------------------
coerce(u: %): Integer == rep u
-------------------------------------------------------------------------
-- from OrderedAbelianMonoid
-------------------------------------------------------------------------
sup(u: %, v: %): % == max(rep u, rep v)
-------------------------------------------------------------------------
-- from with
-------------------------------------------------------------------------
)abbrev domain PRSTNNI PRestrictedNonNegativeInteger
++ Authors: J. Grabmeier
++ Date Created: 19.07.2016
++ Change History:
++ Version: 0.1 19.07.2017 new setting
++ Version: 0.0 in prdp.spad
++ Basic Operations:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ Reference:
++ Description:
++ Domain implements p-restricted non negative integers, i.e. 0,1,2,...,p-1
++ all larger and smaller results of multiplication and addition of the
++ corresponding integers are reduced modulo p-1. This is the additive
++ version of multiplicative terms x^i implementing the relation x^p=x,
x^0=1.
++ It is used to have sympolic elements of fields of characteristic p and
++ in particular for the distributive normal forms of Boolean expressions
++ where p = 2 and x^(-1) represents the negation of x.
PRestrictedNonNegativeInteger(p: PositiveInteger): Exports == Implementation
where
PI ==> PositiveInteger
NNI ==> NonNegativeInteger
I ==> Integer
OF ==> OutputForm
mLP ==> messageLevelPointer
msg ==> message
Exports ==> Join(PRestrictedNumberSystem p, CoercibleTo
NonNegativeInteger)
Implementation ==> NonNegativeInteger add
--------------------------------------------------------------------------------
-- import of domains and packages
--------------------------------------------------------------------------------
import OF
--------------------------------------------------------------------------------
-- global variables
--------------------------------------------------------------------------------
Rep := NonNegativeInteger
rep(s : %) : Rep == s @ Rep
per(rp : Rep) : % == rp @ %
--------------------------------------------------------------------------------
-- implementation of local functions
--------------------------------------------------------------------------------
normalize(n: NonNegativeInteger): NonNegativeInteger ==
++ if n >= p x^p = x is equivalent to x^n = x^r if n rem (p-1) = r and
++ r > 0, in case r=0, we have to keep x^(p-1).
n = 0 => 0
p = 2 => 1
n < p => n
r : NonNegativeInteger := (n rem (p-1)) :: NonNegativeInteger
r = 0 => (p-1) :: NonNegativeInteger
r
--------------------------------------------------------------------------------
-- implementation of global functions
--------------------------------------------------------------------------------
-------------------------------------------------------------------------
-- from CoercibleTo OutputForm
-------------------------------------------------------------------------
-- the following code should be superfluous, if we can guarantee
-- that only the numbers -p+1 ... p-1 are used as internal representatives!
coerce(u: %): OutputForm == coerce(normalize(rep u))$NonNegativeInteger
-------------------------------------------------------------------------
-- from SemiGroup
-------------------------------------------------------------------------
((u: %) + (v: %)): % == per normalize (rep(u) + rep(v))
-------------------------------------------------------------------------
-- from Magma
-------------------------------------------------------------------------
((u: %) * (v: %)): % == per normalize (rep(u) * rep(v))
-------------------------------------------------------------------------
-- from CoercibleFrom NonNegativeInteger
-------------------------------------------------------------------------
coerce(i: NonNegativeInteger): % == normalize i
-------------------------------------------------------------------------
-- from CoercibleTo NNI
-------------------------------------------------------------------------
coerce(u: %): NonNegativeInteger == rep u
-------------------------------------------------------------------------
-- from with
-------------------------------------------------------------------------
)abbrev category PRDP PRestrictedDirectProductCategory
++ Author: J. Grabmeier
++ Date Created: 19.07.2017
++ Change History:
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
PRestrictedDirectProductCategory(p: PositiveInteger, n: NonNegativeInteger,
PR: PRestrictedNumberSystem(p)): Category == Join(InternalMessageLevel,
Finite,
OrderedAbelianMonoidSup, DirectProductCategory(n, PR),
RetractableTo DirectProduct(n, PR),
CoercibleTo DirectProduct(n, PR) )
)abbrev domain PRDPINT PRestrictedDirectProductInteger
PRestrictedDirectProductInteger(p: PositiveInteger, n: NonNegativeInteger):
public == private where
public ==> Join(PRestrictedDirectProductCategory(p, n,
PRestrictedInteger(p)),
CoercibleTo Vector Integer)
private ==> DirectProduct(n, PRestrictedInteger(p)) add
NNI ==> NonNegativeInteger
PRNNI ==> PRestrictedNonNegativeInteger(p)
PI ==> PositiveInteger
OF ==> OutputForm
mLP ==> messageLevelPointer
msg ==> message
-------------------------------------------------------------------------
-- representation
-------------------------------------------------------------------------
Rep := DirectProduct(n, PRestrictedInteger(p))
rep(x:%):Rep == x :: Rep
per(r:Rep):% == r :: %
-------------------------------------------------------------------------
-- import of domain and packages
-------------------------------------------------------------------------
import OF
-------------------------------------------------------------------------
-- global functions
-------------------------------------------------------------------------
-------------------------------------------------------------------------
-- from CoercibleTo Vector NonNegativeInteger
-------------------------------------------------------------------------
coerce(z: %): Vector Integer == vector [zi :: Integer for zi in entries rep
z]
)abbrev domain PRDPNNI PRestrictedDirectProductNonNegativeInteger
PRestrictedDirectProductNonNegativeInteger(p: PositiveInteger, n:
NonNegativeInteger):
public == private where
public ==> Join(PRestrictedDirectProductCategory(p, n,
PRestrictedNonNegativeInteger(p)),
CoercibleTo Vector NonNegativeInteger)
private ==> DirectProduct(n, PRestrictedNonNegativeInteger(p)) add
NNI ==> NonNegativeInteger
PRNNI ==> PRestrictedNonNegativeInteger(p)
PI ==> PositiveInteger
OF ==> OutputForm
mLP ==> messageLevelPointer
msg ==> message
-------------------------------------------------------------------------
-- representation
-------------------------------------------------------------------------
Rep := DirectProduct(n, PRestrictedNonNegativeInteger(p))
rep(x:%):Rep == x :: Rep
per(r:Rep):% == r :: %
-------------------------------------------------------------------------
-- import of domain and packages
-------------------------------------------------------------------------
import OF
-------------------------------------------------------------------------
-- global functions
-------------------------------------------------------------------------
-------------------------------------------------------------------------
-- from CoercibleTo Vector NonNegativeInteger
-------------------------------------------------------------------------
coerce(z: %): Vector NNI == vector [zi :: NNI for zi in entries rep z]
)abbrev category GPRPCAT MaybeSkewPRestrictedPolynomialCategory
++ Author:
++ Basic Functions: Ring, monomial, coefficient, differentiate, eval
++ Related Constructors: Polynomial, DistributedMultivariatePolynomial
++ Also See: UnivariatePolynomialCategory
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++ The category for general multi-variate possibly skew polynomials
++ over a ring R, in variables from VarSet, with exponents from the
++ \spadtype{OrderedAbelianMonoidSup}.
MaybeSkewPRestrictedPolynomialCategory(R : Join(SemiRng, AbelianMonoid),
EV: OrderedAbelianMonoidSup,
E : OrderedAbelianMonoidSup,
VarSet : OrderedSet) : Category ==
Join(InternalMessageLevel, FiniteAbelianMonoidRing(R, E)) with
if R has Ring then
FullyLinearlyExplicitOver R
-- operations
degree : (%, VarSet) -> EV
++ degree(p, v) gives the degree of polynomial p with respect
++ to the variable v.
degree : (%, List(VarSet)) -> List(EV)
++ degree(p, lv) gives the list of degrees of polynomial p
++ with respect to each of the variables in the list lv.
coefficient : (%, VarSet, EV) -> %
++ coefficient(p, v, n) views the polynomial p as a univariate
++ polynomial in v and returns the coefficient of the \spad{v^n} term.
coefficient : (%, List VarSet, List EV) -> %
++ coefficient(p, lv, ln) views the polynomial p as a polynomial
++ in the variables of lv and returns the coefficient of the term
++ \spad{lv^ln}, i.e. \spad{prod(lv_i ^ ln_i)}.
monomials : % -> List %
++ monomials(p) returns the list of non-zero monomials of
++ polynomial p, i.e.
++ \spad{monomials(sum(a_(i) X^(i))) = [a_(1) X^(1), ..., a_(n) X^(n)]}.
mainVariable : % -> Union(VarSet,"failed")
++ mainVariable(p) returns the biggest variable which actually
++ occurs in the polynomial p, or "failed" if no variables are
++ present.
++ fails precisely if polynomial satisfies ground?
monomial : (%, VarSet, EV) -> %
++ monomial(a, x, n) creates the monomial \spad{a*x^n} where \spad{a} is
++ a polynomial, x is a variable and n is a nonnegative integer.
monomial : (%, List VarSet, List EV) -> %
++ monomial(a, [v1..vn], [e1..en]) returns \spad{a*prod(vi^ei)}.
totalDegree : % -> EV
++ totalDegree(p) returns the largest sum over all monomials
++ of all exponents of a monomial.
totalDegree : (%, List VarSet) -> EV
++ totalDegree(p, lv) returns the maximum sum (over all monomials
++ of polynomial p) of the variables in the list lv.
totalDegreeSorted : (%, List VarSet) -> EV
++ totalDegreeSorted(p, lv) returns the maximum sum (over all
++ monomials of polynomial p) of the degree in variables in the
++ list lv. lv is assumed to be sorted in decreasing order.
variables : % -> List(VarSet)
++ variables(p) returns the list of those variables actually
++ appearing in the polynomial p.
if R has SemiRing then
primitiveMonomials : % -> List %
++ primitiveMonomials(p) gives the list of monomials of the
++ polynomial p with their coefficients removed.
++ Note: \spad{primitiveMonomials(sum(a_(i) X^(i))) =
++ [X^(1), ..., X^(n)]}.
if R has Comparable then Comparable
-- assertions
if R has canonicalUnitNormal then canonicalUnitNormal
++ we can choose a unique representative for each
++ associate class.
++ This normalization is chosen to be normalization of
++ leading coefficient (by default).
add
p : %
ln : List EV
lv : List VarSet
monomials p ==
-- sequential version for efficiency, by WMSIT, 7/30/90
ml := empty$List(%)
while p ~= 0 repeat
ml := concat(leadingMonomial p, ml)
p := reductum p
reverse! ml
monomial(p, lv, ln) ==
empty? lv =>
empty? ln => p
error "mismatched lists in monomial"
empty? ln => error "mismatched lists in monomial"
monomial(monomial(p, first lv, first ln), rest lv, rest ln)
if R has SemiRing then
mkPrim(p : %) : % == monomial(1, degree p)
-- primitiveMonomials p ==
-- ml := empty$List(%)
-- while p ~= 0 repeat
-- ml := concat(mkPrim(leadingMonomial p), ml)
-- p := reductum p
-- reverse! ml
)abbrev category PRPCAT PRestrictedPolynomialCategory
++ Author:
++ Basic Functions: Ring, monomial, coefficient, differentiate, eval
++ Related Constructors: Polynomial, DistributedMultivariatePolynomial
++ Also See: UnivariatePolynomialCategory
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++ The category for general multi-variate polynomials over a ring
++ R, in variables from VarSet, with exponents from the
++ \spadtype{OrderedAbelianMonoidSup}. Here variables commute
++ with the coefficients.
PRestrictedPolynomialCategory(R : Join(SemiRng, AbelianMonoid),
EV : OrderedAbelianMonoidSup,
E : OrderedAbelianMonoidSup, VarSet : OrderedSet):
Category ==
Join(MaybeSkewPRestrictedPolynomialCategory(R, EV, E, VarSet),
InnerEvalable(VarSet, R), InnerEvalable(VarSet, %),
VariablesCommuteWithCoefficients) with
if R has SemiRing then
RetractableTo VarSet
Evalable %
if R has Ring then
PartialDifferentialRing VarSet
-- operations
univariate : (%, VarSet) -> SparseUnivariatePolynomial(%)
++ univariate(p, v) converts the multivariate polynomial p
++ into a univariate polynomial in v, whose coefficients are still
++ multivariate polynomials (in all the other variables).
univariate : % -> SparseUnivariatePolynomial(R)
++ univariate(p) converts the multivariate polynomial p,
++ which should actually involve only one variable,
++ into a univariate polynomial
++ in that variable, whose coefficients are in the ground ring.
++ Error: if polynomial is genuinely multivariate
minimumDegree : (%, VarSet) -> EV
++ minimumDegree(p, v) gives the minimum degree of polynomial p
++ with respect to v, i.e. viewed a univariate polynomial in v
minimumDegree : (%, List(VarSet)) -> List(EV)
++ minimumDegree(p, lv) gives the list of minimum degrees of the
++ polynomial p with respect to each of the variables in the list lv
if R has Ring then
monicDivide : (%, %, VarSet) -> Record(quotient : %, remainder : %)
++ monicDivide(a, b, v) divides the polynomial a by the polynomial b,
++ with each viewed as a univariate polynomial in v returning
++ both the quotient and remainder.
++ Error: if b is not monic with respect to v.
multivariate : (SparseUnivariatePolynomial(R), VarSet) -> %
++ multivariate(sup, v) converts an anonymous univariable
++ polynomial sup to a polynomial in the variable v.
multivariate : (SparseUnivariatePolynomial(%), VarSet) -> %
++ multivariate(sup, v) converts an anonymous univariable
++ polynomial sup to a polynomial in the variable v.
isPlus : % -> Union(List %, "failed")
++ isPlus(p) returns \spad{[m1, ..., mn]} if polynomial \spad{p = m1 +
... + mn} and
++ \spad{n >= 2} and each mi is a nonzero monomial.
if R has SemiRing then
isTimes : % -> Union(List %, "failed")
++ isTimes(p) returns \spad{[a1, ..., an]} if polynomial
++ \spad{p = a1 ... an} and \spad{n >= 2}, and, for each i,
++ ai is either a nontrivial constant in R or else of the
++ form \spad{x^e}, where \spad{e > 0} is an integer and
++ x is a member of VarSet.
isExpt : % -> Union(Record(var : VarSet, exponent : EV),
"failed")
++ isExpt(p) returns \spad{[x, n]} if polynomial p has the form
++ \spad{x^n} and \spad{n > 0}.
-- OrderedRing view removed to allow EXPR to define abs
--if R has OrderedRing then OrderedRing
if (R has ConvertibleTo InputForm) and
(VarSet has ConvertibleTo InputForm) then
ConvertibleTo InputForm
if R has Ring then
if (R has ConvertibleTo Pattern Integer) and
(VarSet has ConvertibleTo Pattern Integer) then
ConvertibleTo Pattern Integer
if (R has ConvertibleTo Pattern Float) and
(VarSet has ConvertibleTo Pattern Float) then
ConvertibleTo Pattern Float
if (R has PatternMatchable Integer) and
(VarSet has PatternMatchable Integer) then
PatternMatchable Integer
if (R has PatternMatchable Float) and
(VarSet has PatternMatchable Float) then
PatternMatchable Float
if R has CommutativeRing then
resultant : (%, %, VarSet) -> %
++ resultant(p, q, v) returns the resultant of the polynomials
++ p and q with respect to the variable v.
discriminant : (%, VarSet) -> %
++ discriminant(p, v) returns the disriminant of the polynomial p
++ with respect to the variable v.
if R has GcdDomain then
GcdDomain
content : (%, VarSet) -> %
++ content(p, v) is the gcd of the coefficients of the polynomial p
++ when p is viewed as a univariate polynomial with respect to the
++ variable v.
++ Thus, for polynomial 7*x^2*y + 14*x*y^2, the gcd of the
++ coefficients with respect to x is 7*y.
primitivePart : % -> %
++ primitivePart(p) returns the unitCanonical associate of the
++ polynomial p with its content divided out.
primitivePart : (%, VarSet) -> %
++ primitivePart(p, v) returns the unitCanonical associate of the
++ polynomial p with its content with respect to the variable v
++ divided out.
squareFree : % -> Factored %
++ squareFree(p) returns the square free factorization of the
++ polynomial p.
squareFreePart : % -> %
++ squareFreePart(p) returns product of all the irreducible factors
++ of polynomial p each taken with multiplicity one.
-- assertions
if R has PolynomialFactorizationExplicit then
PolynomialFactorizationExplicit
add
import from Factored(%)
p : %
v : VarSet
ln : List EV
lv : List VarSet
n : EV
pp, qq : SparseUnivariatePolynomial %
if R has SemiRing then
eval(p : %, l : List Equation %) ==
empty? l => p
for e in l repeat
retractIfCan(lhs e)@Union(VarSet,"failed") case "failed" =>
error "cannot find a variable to evaluate"
lvar := [retract(lhs e)@VarSet for e in l]
eval(p, lvar, [rhs e for e in l]$List(%))
isPlus p ==
empty? rest(l := monomials p) => "failed"
l
if R has SemiRing then
isTimes p ==
empty?(lv := variables p) or not monomial? p => "failed"
l := [monomial(1, v, degree(p, v)) for v in lv]
((r := leadingCoefficient p) = 1) =>
empty? rest lv => "failed"
l
concat(r::%, l)
isExpt p ==
(u := mainVariable p) case "failed" => "failed"
p = monomial(1, u::VarSet, d := degree(p, u::VarSet)) =>
[u::VarSet, d]
"failed"
-- to make the following work we would need SUP over pRestricted Exponents
--coefficient(p, v, n) == coefficient(univariate(p, v), n)
--coefficient(p, lv, ln) ==
-- empty? lv =>
-- empty? ln => p
-- error "mismatched lists in coefficient"
-- empty? ln => error "mismatched lists in coefficient"
-- coefficient(coefficient(univariate(p, first lv), first ln),
-- rest lv, rest ln)
retract(p : %) : VarSet ==
q := mainVariable(p)::VarSet
q::% = p => q
error "Polynomial is not a single variable"
retractIfCan(p:%):Union(VarSet, "failed") ==
((q := mainVariable p) case VarSet) and (q::VarSet::% = p) => q
"failed"
-- totalDegree p ==
-- ground? p => 0
-- u := univariate(p, mainVariable(p)::VarSet)
-- d : EV := 0
-- while u ~= 0 repeat
-- d := max(d, degree u + totalDegree leadingCoefficient u)
-- u := reductum u
-- d
-- totalDegreeSorted(p : %, lv : List VarSet) : EV ==
-- ground? p => 0
-- empty?(lv) => 0
-- u := univariate(p, v := (mainVariable(p)::VarSet))
-- d : EV := 0
-- w : EV := 0
-- v0 := first(lv)
-- while v < v0 repeat
-- lv := rest(lv)
-- empty?(lv) => return 0
-- v0 := first(lv)
-- if v0 = v then
-- w := 1
-- lv := rest(lv)
-- while u ~= 0 repeat
-- d := max(d, w*(degree u) +
-- totalDegreeSorted(leadingCoefficient u, lv))
-- u := reductum u
-- d
totalDegree(p, lv) ==
lv1 := sort((v1 : VarSet, v2 : VarSet) : Boolean +-> v2 < v1, lv)
totalDegreeSorted(p, lv1)
-- Condition on % is redundant, but compiler can not infer it
if % has CommutativeRing and R has CommutativeRing then
resultant(p1, p2, mvar) ==
resultant(univariate(p1, mvar), univariate(p2, mvar))
discriminant(p, var) ==
discriminant(univariate(p, var))
if R has IntegralDomain then
-- allMonoms(l : List %) : List(%) ==
-- removeDuplicates! concat [primitiveMonomials p for p in l]
-- P2R(p : %, b : List E, n : NonNegativeInteger) : Vector(R) ==
-- w := new(n, 0)$Vector(R)
-- for i in minIndex w .. maxIndex w for bj in b repeat
-- qsetelt!(w, i, coefficient(p, bj))
-- w
-- eq2R(l : List %, b : List E) : Matrix(R) ==
-- matrix [[coefficient(p, bj) for p in l] for bj in b]
-- reducedSystem(m : Matrix %) : Matrix(R) ==
-- l := listOfLists m
-- b := removeDuplicates!
-- concat [allMonoms r for r in l]$List(List(%))
-- empty?(b) => new(0, ncols(m), 0)$Matrix(R)
-- d := [degree bj for bj in b]
-- mm := eq2R(first l, d)
-- l := rest l
-- while not empty? l repeat
-- mm := vertConcat(mm, eq2R(first l, d))
-- l := rest l
-- mm
-- reducedSystem(m : Matrix %, v : Vector %):
-- Record(mat : Matrix R, vec : Vector R) ==
-- l := listOfLists m
-- r := entries v
-- b : List % := removeDuplicates! concat(allMonoms r,
-- concat [allMonoms s for s in l]$List(List(%)))
-- empty?(b) =>
-- [new(0, ncols(m), 0)$Matrix(R), new(0, 0)$Vector(R)]
-- d := [degree bj for bj in b]
-- n := #d
-- mm := eq2R(first l, d)
-- w := P2R(first r, d, n)
-- l := rest l
-- r := rest r
-- while not empty? l repeat
-- mm := vertConcat(mm, eq2R(first l, d))
-- w := concat(w, P2R(first r, d, n))
-- l := rest l
-- r := rest r
-- [mm, w]
-- if R has PolynomialFactorizationExplicit then
-- -- we might be in trouble if its actually only
-- -- a univariate polynomial category - have to remember to
-- -- over-ride these in UnivariatePolynomialCategory
-- PFBR ==>PolynomialFactorizationByRecursion(R, E, VarSet, %)
-- solveLinearPolynomialEquation(lpp, pp) ==
-- solveLinearPolynomialEquationByRecursion(lpp, pp)$PFBR
-- if R has FiniteFieldCategory then
-- MFFACT ==> MultFiniteFactorize(VarSet, E, R, %)
-- factorPolynomial(pp) == factor(pp)$MFFACT
-- factorSquareFreePolynomial(pp) == factor(pp)$MFFACT
-- factor(p) == factor(p)$MFFACT
-- else
-- if R has CharacteristicZero and R has EuclideanDomain then
-- MF ==> InnerMultFact(VarSet, E, R, %)
-- factorPolynomial(pp) ==
-- factor(pp, factorPolynomial$R)$MF
-- factor(p) == factor(p, factorPolynomial$R)$MF
-- factorSquareFreePolynomial(pp) ==
-- factor(pp, factorPolynomial$R)$MF
-- else
-- gcdPolynomial(pp, qq) ==
-- gcdPolynomial(pp, qq
-- )$GeneralPolynomialGcdPackage(E, VarSet, R, %)
-- factorPolynomial(pp) == factorByRecursion(pp)$PFBR
-- factorSquareFreePolynomial(pp) ==
-- factorSquareFreeByRecursion(pp)$PFBR
-- factor p ==
-- v : Union(VarSet, "failed") := mainVariable p
-- v case "failed" =>
-- ansR := factor leadingCoefficient p
-- makeFR(unit(ansR)::%,
-- [[w.flg, w.fctr::%, w.xpnt] for w in factorList
ansR])
-- up : SparseUnivariatePolynomial % := univariate(p, v)
-- ansSUP := factorByRecursion(up)$PFBR
-- makeFR(multivariate(unit(ansSUP), v),
-- [[ww.flg, multivariate(ww.fctr, v), ww.xpnt]
-- for ww in factorList ansSUP])
if R has CharacteristicNonZero then
mat : Matrix %
-- conditionP mat ==
-- ll := listOfLists transpose mat -- hence each list corresponds
to a
-- -- column, i.e. to one variable
-- llR : List List R := [ empty() for z in first ll]
-- monslist : List List % := empty()
-- ch : Integer := characteristic()$%
-- for l in ll repeat
-- mons := "setUnion"/[primitiveMonomials u for u in l]
-- redmons : List % := []
-- for m in mons repeat
-- vars := variables m
-- degs := degree(m, vars)
-- deg1 : List EV
-- deg1 := [ ((nd := (d::Integer) exquo ch)
-- case "failed" => return "failed" ;
-- nd::Integer::EV)
-- for d in degs ]
-- redmons := cons(monomial(1, vars, deg1), redmons)
-- llR := [cons(ground coefficient(u, vars, degs), v)
-- for u in l for v in llR]
-- monslist := cons(redmons, monslist)
-- ans := conditionP transpose matrix llR
-- ans case "failed" => "failed"
-- i : NonNegativeInteger := 0
-- [ +/[m*(ans.(i := i+1))::% for m in mons ]
-- for mons in monslist]
--)fin
-- if R has CharacteristicNonZero then
-- charthRootlv : (%,List VarSet, EV) -> Union(%,"failed")
-- charthRoot p ==
-- vars := variables p
-- empty? vars =>
-- ans := charthRoot ground p
-- ans case "failed" => "failed"
-- ans::R::%
-- ch := characteristic()$%
-- charthRootlv(p, vars, ch)
-- charthRootlv(p, vars, ch) ==
-- empty? vars =>
-- ans := charthRoot ground p
-- ans case "failed" => "failed"
-- ans::R::%
-- v := first vars
-- vars := rest vars
-- d := degree(p, v)
-- ans : % := 0
-- while (d>0) repeat
-- (dd := (d::Integer exquo ch::Integer)) case "failed" =>
-- return "failed"
-- cp := coefficient(p, v, d)
-- p := p-monomial(cp, v, d)
-- ansx := charthRootlv(cp, vars, ch)
-- ansx case "failed" => return "failed"
-- d := degree(p, v)
-- ans := ans+monomial(ansx, v, dd::Integer::NonNegativeInteger)
-- ansx := charthRootlv(p, vars, ch)
-- ansx case "failed" => return "failed"
-- return ans+ansx
if R has Ring then
monicDivide(p1, p2, mvar) ==
result := monicDivide(univariate(p1, mvar), univariate(p2, mvar))
[multivariate(result.quotient, mvar),
multivariate(result.remainder, mvar)]
if R has GcdDomain then
-- if R has EuclideanDomain and R has CharacteristicZero then
-- squareFree p == squareFree(p)$MultivariateSquareFree(E, VarSet, R, %)
-- else
-- squareFree p == squareFree(p)$PolynomialSquareFree(VarSet, E, R, %)
-- squareFreePart p ==
-- unit(s := squareFree p) * */[f.factor for f in factors s]
content(p, v) == content univariate(p, v)
primitivePart p ==
zero? p => p
unitNormal((p exquo content p) ::%).canonical
primitivePart(p, v) ==
zero? p => p
unitNormal((p exquo content(p, v)) ::%).canonical
-- if R has Comparable then
-- smaller?(p : %, q : %) ==
-- while p ~= 0 and q ~= 0 repeat
-- (dp := degree p) < (dq := degree q) =>
-- return smaller?(0, leadingCoefficient q)
-- dq < dp => return smaller?(leadingCoefficient p, 0)
-- lp := leadingCoefficient(p)
-- lq := leadingCoefficient(q)
-- lp ~= lq => return smaller?(lp, lq)
-- p := reductum p
-- q := reductum q
-- p = 0 => smaller?(0, leadingCoefficient q)
-- smaller?(leadingCoefficient p, 0)
--
-- if R has Ring then
-- if (R has PatternMatchable Integer) and
-- (VarSet has PatternMatchable Integer) then
-- patternMatch(p : %, pat : Pattern Integer,
-- l : PatternMatchResult(Integer, %)) ==
-- patternMatch(p, pat, l
-- )$PatternMatchPolynomialCategory(Integer, E, VarSet, R, %)
-- if (R has PatternMatchable Float) and
-- (VarSet has PatternMatchable Float) then
-- patternMatch(p : %, pat : Pattern Float,
-- l : PatternMatchResult(Float, %)) ==
-- patternMatch(p, pat, l
-- )$PatternMatchPolynomialCategory(Float, E, VarSet, R, %)
-- if R has Ring then
-- if (R has ConvertibleTo Pattern Integer) and
-- (VarSet has ConvertibleTo Pattern Integer) then
-- convert(x : %) : Pattern(Integer) ==
-- map(convert, convert, x
-- )$PolynomialCategoryLifting(E, VarSet, R, %, Pattern Integer)
-- if (R has ConvertibleTo Pattern Float) and
-- (VarSet has ConvertibleTo Pattern Float) then
-- convert(x : %) : Pattern(Float) ==
-- map(convert, convert, x
-- )$PolynomialCategoryLifting(E, VarSet, R, %, Pattern Float)
-- if (R has ConvertibleTo InputForm) and
-- (VarSet has ConvertibleTo InputForm) then
-- convert(p : %) : InputForm ==
-- map(convert, convert,
-- p)$PolynomialCategoryLifting(E, VarSet, R, %, InputForm)
--
)abbrev domain GDPRMP GeneralDistributedPRestrictedMultivariatePolynomial
++ Author: Johannes Grabmeier using Barry Trager's code for GDMP
++ Basic Functions: Ring, degree, eval, coefficient, monomial, differentiate,
++ resultant, gcd, leadingCoefficient
++ Related Constructors: DistributedMultivariatePolynomial,
++ HomogeneousDistributedMultivariatePolynomial
++ Also See: Polynomial
++ AMS Classifications:
++ Keywords: polynomial, multivariate, distributed
++ References:
++ Description:
++ This type supports distributed multivariate polynomials
++ whose variables are from a user specified list of symbols.
++ The coefficient ring may be non commutative,
++ but the variables are assumed to commute.
++ The term ordering is specified by its third parameter.
++ Suggested types which define term orderings include:
\spadtype{DirectProduct},
++ \spadtype{HomogeneousDirectProduct}, \spadtype{SplitHomogeneousDirectProduct}
++ and finally \spadtype{OrderedDirectProduct} which accepts an arbitrary user
++ function to define a term ordering.
GeneralDistributedPRestrictedMultivariatePolynomial(q, vl, R, EV, E) : public
== private where
q : PositiveInteger
vl : List Symbol
R : Ring
EV : PRestrictedNumberSystem q
E : PRestrictedDirectProductCategory(q, #vl, EV)
OV ==> OrderedVariableList(vl)
SUP ==> SparseUnivariatePolynomial
NNI ==> NonNegativeInteger
public == PRestrictedPolynomialCategory(R, EV, E, OV) with
reorder : (%, List Integer) -> %
++ reorder(p, perm) applies the permutation perm to the variables
++ in a polynomial and returns the new correctly ordered polynomial
private == PolynomialRing(R, E) add
--representations
Term := Record(k : E, c : R)
Rep := List Term
n := #vl
VE ==> Vector E
zero?(p : %) : Boolean == null(p :: Rep)
-- totalDegree(p: %): EV ==
-- zero? p => 0
-- "max"/[reduce("+", entries(t.k), 0) for t in p]
monomial(p : %, v : OV, e : EV) : % ==
locv := lookup v
vEV : Vector EV := [if z = locv then e else 0 for z in 1..n]
p*monomial(1, directProduct vEV)
coerce(v : OV) : % == monomial(1, v, 1)
listCoef(p : %) : List R ==
rec : Term
[rec.c for rec in (p::Rep)]
-- mainVariable(p : %) ==
-- zero?(p) => "failed"
-- for v in vl repeat
-- vv := variable(v)::OV
-- if degree(p, vv)>0 then return vv
-- "failed"
--
-- ground?(p) == mainVariable(p) case "failed"
--
-- retract(p : %) : R ==
-- not ground? p => error "not a constant"
-- leadingCoefficient p
--
-- retractIfCan(p : %) : Union(R,"failed") ==
-- ground?(p) => leadingCoefficient p
-- "failed"
degree(p : %, v : OV): EV ==
-- degree(univariate(p, v))
zero?(p) => 0
res : EV := 0
locv := lookup v
while not empty? p repeat
t := first p
j := t.k.locv
if j > res then res := j
p := rest p
res
-- minimumDegree(p : %, v : OV) == minimumDegree(univariate(p, v))
-- differentiate(p : %, v : OV) ==
-- multivariate(differentiate(univariate(p, v)), v)
degree(p : %, lv : List OV) == [degree(p, v) for v in lv]
-- minimumDegree(p : %, lv : List OV) == [minimumDegree(p, v) for v in lv]
numberOfMonomials(p : %) ==
l : Rep := p :: Rep
null(l) => 1
#l
monomial?(p : %) : Boolean ==
l : Rep := p :: Rep
null(l) or null rest(l)
if R has OrderedRing then
maxNorm(p : %) : R ==
l : List R := nil
r, m : R
m := 0
for r in listCoef(p) repeat
if r > m then m := r
else if (-r) > m then m := -r
m
if R has Field then
(p : %) / (r : R) == inv(r) * p
variables(p : %) ==
maxdeg : Vector EV := new(n, 0)
while not zero?(p) repeat
tdeg := degree p
p := reductum p
for i in 1..n repeat
maxdeg.i := max(maxdeg.i, tdeg.i)
[index(i::PositiveInteger) for i in 1..n | maxdeg.i ~= 0]
-- reorder(p : %, perm : List Integer) : % ==
-- #perm ~= n => error "must be a complete permutation of all vars"
-- q := [[directProduct [term.k.j for j in perm]$VE, term.c]$Term
-- for term in p]
-- sort((z1, z2) +-> z1.k > z2.k, q)
--
-- univariate(p : %, v : OV) : SUP(%) ==
-- zero?(p) => 0
-- exp := degree p
-- locv := lookup v
-- deg : EV := 0
-- nexp := directProduct [if i = locv then (deg := exp.i;0) else exp.i
-- for i in 1..n]$VE
-- monomial(monomial(leadingCoefficient p, nexp), deg)+
-- univariate(reductum p, v)
eval(p : %, v : OV, val : %) : % == univariate(p, v)(val)
eval(p : %, v : OV, val : R) : % == eval(p, v, val::%)$%
eval(p : %, lv : List OV, lval : List R) : % ==
lv = [] => p
eval(eval(p, first lv, (first lval)::%)$%, rest lv, rest lval)$%
-- assume Lvar are sorted correctly
evalSortedVarlist(p : %, Lvar : List OV, Lpval : List %) : % ==
v := mainVariable p
v case "failed" => p
pv := v:: OV
Lvar=[] or Lpval=[] => p
mvar := Lvar.first
mvar > pv => evalSortedVarlist(p, Lvar.rest, Lpval.rest)
pval := Lpval.first
pts : SUP(%) := map(x+->evalSortedVarlist(x, Lvar, Lpval),
univariate(p, pv))
mvar = pv => pts(pval)
multivariate(pts, pv)
eval(p : %, Lvar : List OV, Lpval : List %) ==
nlvar : List OV := sort((x, y) +-> x > y, Lvar)
nlpval :=
Lvar = nlvar => Lpval
nlpval := [Lpval.position(mvar, Lvar) for mvar in nlvar]
evalSortedVarlist(p, nlvar, nlpval)
-- multivariate(p1 : SUP(%), v : OV) : % ==
-- 0 = p1 => 0
-- degree p1 = 0 => leadingCoefficient p1
-- leadingCoefficient(p1)*(v::%)^degree(p1) +
-- multivariate(reductum p1, v)
--
-- univariate(p : %) : SUP(R) ==
-- (v := mainVariable p) case "failed" =>
-- monomial(leadingCoefficient p, 0)
-- q := univariate(p, v:: OV)
-- ans : SUP(R) := 0
-- while q ~= 0 repeat
-- ans := ans + monomial(ground leadingCoefficient q, degree q)
-- q := reductum q
-- ans
--
-- multivariate(p : SUP(R), v : OV) : % ==
-- 0 = p => 0
-- (leadingCoefficient p)*monomial(1, v, degree p) +
-- multivariate(reductum p, v)
if R has GcdDomain then
content(p : %) : R ==
zero?(p) => 0
"gcd"/[t.c for t in p]
if R has EuclideanDomain and not(R has FloatingPointSystem) then
gcd(p : %, q : %) : % ==
gcd(p, q)$PolynomialGcdPackage(E, OV, R, %)
else gcd(p : %, q : %) : % ==
r : R
(pv := mainVariable(p)) case "failed" =>
(r := leadingCoefficient p) = 0$R => q
gcd(r, content q)::%
(qv := mainVariable(q)) case "failed" =>
(r := leadingCoefficient q) = 0$R => p
gcd(r, content p)::%
pv<qv => gcd(p, content univariate(q, qv))
qv<pv => gcd(q, content univariate(p, pv))
multivariate(gcd(univariate(p, pv), univariate(q, qv)), pv)
coerce(p : %) : OutputForm ==
zero?(p) => (0$R) :: OutputForm
l, lt : List OutputForm
lt := nil
vl1 := [v::OutputForm for v in vl]
for t in reverse p repeat
l := nil
for i in 1..#vl1 repeat
t.k.i = 0 => "next"
t.k.i = 1 => l := cons(vl1.i, l)
l := cons(vl1.i ^ t.k.i ::OutputForm, l)
l := reverse l
if (t.c ~= 1) or (null l) then l := cons(t.c :: OutputForm, l)
1 = #l => lt := cons(first l, lt)
lt := cons(reduce("*",l),lt)
1 = #lt => first lt
reduce("+",lt)
)abbrev domain DPRMP DistributedPRestrictedMultivariatePolynomial
++ Author: Barry Trager
++ Basic Functions: Ring, degree, eval, coefficient, monomial, differentiate,
++ resultant, gcd, leadingCoefficient
++ Related Constructors: GeneralDistributedMultivariatePolynomial,
++ HomogeneousDistributedMultivariatePolynomial
++ Also See: Polynomial
++ AMS Classifications:
++ Keywords: polynomial, multivariate, distributed
++ References:
++ Description:
++ This type supports distributed multivariate polynomials
++ whose variables are from a user specified list of symbols.
++ The coefficient ring may be non commutative,
++ but the variables are assumed to commute.
++ The term ordering is lexicographic specified by the variable
++ list parameter with the most significant variable first in the list.
DistributedPRestrictedMultivariatePolynomial(p, vl, R, EV, E) : public ==
private where
p : PositiveInteger
vl : List Symbol
R : Ring
EV : PRestrictedNumberSystem p
E : PRestrictedDirectProductCategory(p, #vl, EV)
OV ==> OrderedVariableList(vl)
public == PRestrictedPolynomialCategory(R, EV, E, OV) with
reorder : (%, List Integer) -> %
++ reorder(p, perm) applies the permutation perm to the variables
++ in a polynomial and returns the new correctly ordered polynomial
private ==
GeneralDistributedPRestrictedMultivariatePolynomial(p, vl, R, EV, E)
)fin
)abbrev domain HDMP HomogeneousDistributedMultivariatePolynomial
++ Author: Barry Trager
++ Basic Functions: Ring, degree, eval, coefficient, monomial, differentiate,
++ resultant, gcd, leadingCoefficient
++ Related Constructors: DistributedMultivariatePolynomial,
++ GeneralDistributedMultivariatePolynomial
++ Also See: Polynomial
++ AMS Classifications:
++ Keywords: polynomial, multivariate, distributed
++ References:
++ Description:
++ This type supports distributed multivariate polynomials
++ whose variables are from a user specified list of symbols.
++ The coefficient ring may be non commutative,
++ but the variables are assumed to commute.
++ The term ordering is total degree ordering refined by reverse
++ lexicographic ordering with respect to the position that the variables
++ appear in the list of variables parameter.
HomogeneousDistributedMultivariatePolynomial(vl, R) : public == private where
vl : List Symbol
R : Ring
E ==> HomogeneousDirectProduct(#vl, NonNegativeInteger)
OV ==> OrderedVariableList(vl)
public == PolynomialCategory(R, E, OV) with
reorder : (%, List Integer) -> %
++ reorder(p, perm) applies the permutation perm to the variables
++ in a polynomial and returns the new correctly ordered polynomial
private ==
GeneralDistributedMultivariatePolynomial(vl, R, E)
--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
--All rights reserved.
--
--Redistribution and use in source and binary forms, with or without
--modification, are permitted provided that the following conditions are
--met:
--
-- - Redistributions of source code must retain the above copyright
-- notice, this list of conditions and the following disclaimer.
--
-- - Redistributions in binary form must reproduce the above copyright
-- notice, this list of conditions and the following disclaimer in
-- the documentation and/or other materials provided with the
-- distribution.
--
-- - Neither the name of The Numerical ALgorithms Group Ltd. nor the
-- names of its contributors may be used to endorse or promote products
-- derived from this software without specific prior written permission.
--
--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
)fin
)abbrev domain DPRMP DistributedPRestrictedMultivariatePolynomial
++ Author: Johannes Grabmeier
++ Date Created:
++ Date Last Updated:
++ Basic Functions: Ring, degree, eval, coefficient, monomial, differentiate,
++ resultant, gcd, leadingCoefficient
++ Related Constructors: GeneralDistributedMultivariatePolynomial,
++ HomogeneousDistributedMultivariatePolynomial
++ Also See: Polynomial
++ AMS Classifications:
++ Keywords: polynomial, multivariate, distributed
++ References:
++ Description:
++ This type supports distributed multivariate polynomials
++ whose variables are from a user specified list of symbols x
++ reduced modulo the ideal generated by x^p-x.
++
++ The coefficient ring may be non commutative,
++ but the variables are assumed to commute.
++ The term ordering is lexicographic specified by the variable
++ list parameter with the most significant variable first in the list.
DistributedPRestrictedMultivariatePolynomial(p, vl, S): public == private where
p : PositiveInteger
vl : List Symbol
S : Ring
E ==> PRestrictedDirectProductNonNegativeInteger(p, #vl)
OV ==> OrderedVariableList(vl)
public == PolynomialCategory(S, E, OV) with
reorder: (%, List Integer) -> %
++ reorder(p, perm) applies the permutation perm to the variables
++ in a polynomial and returns the new correctly ordered polynomial
allVariables: () -> List %
private == GeneralDistributedMultivariatePolynomial(vl, S, E) --add
--allVariables(): List % == [index(i::PositiveInteger)$OV::% for i in
1..#vl]
-- überflüssiger Code:
)if false
-------------------------------------------------------------------------
-- from RetractableTo DirectProduct(n, PRNNI)
-------------------------------------------------------------------------
-- coerce(z: DirectProduct(n, NNI)): % ==
-- vzpr : Vector PRNNI := vector [zi :: PRNNI for zi in entries z]
-- directProduct(vzpr)$Rep
-- retractIfCan(z: %): Union(DirectProduct(n, NNI), "failed") ==
retractIfCan(rep z)
--
-- coerce(z: %): DirectProduct(n, NNI) ==
-- vzpr : Vector NNI := vector [zi :: NNI for zi in entries rep z]
-- directProduct vzpr
-------------------------------------------------------------------------
-- from Aggregate
-------------------------------------------------------------------------
--elt(e: %, i: Integer): NonNegativeInteger ==
-- eDP := e :: DirectProduct(n, NNI)
-- eDP.i
-- parts(e: %): List NNI == [e.i :: NNI for i in 1..n]
directProduct(z: Vector NNI): % ==
if mLP 1 then print blankSeparate [msg "directProduct entered"]
not size?(z, n) => error "directProduct: argument is not of the correct
length!"
lz : List NNI := entries z
if mLP 1 then print blankSeparate [msg "lz", lz::OF]
lzpr : List PRNNI := [zi :: PRNNI for zi in lz]
if mLP 1 then print blankSeparate [msg "lzpr", lzpr::OF]
vzpr : Vector PRNNI := vector lzpr
if mLP 1 then print blankSeparate [msg "vzpr", vzpr::OF]
directProduct(vzpr)$Rep
--per directProduct(vector [zi :: PRestrictedNonNegativeInteger(p) for zi
in entries z])$Rep
coerce(z: %): Vector NNI == coerce(rep(z))$Rep
-------------------------------------------------------------------------
-- from Finite
-------------------------------------------------------------------------
random(): % == directProduct [random(p::Integer)$Integer :: NNI for i in
1..n]
lookup(x: %): PI ==
xI : List Integer := [xx :: Integer for xx in entries(x::Rep)]
if mLP 1 then print blankSeparate [msg "xRep = ", xI :: OF]
(1 + wholePart(wholeRadix(xI)$RadixExpansion(p))) :: PI
index(q: PI): % ==
q > p^n => error "index: argument too large!"
lq : List NNI := [d::NNI for d in wholeRagits ((q-1) ::
RadixExpansion(p)) ]
directProduct vector concat( new((n#lq)::NNI,0), lq)
-------------------------------------------------------------------------
-- from RetractableFrom DirectProduct
-------------------------------------------------------------------------
coerce(u: %): DirectProduct(n, NNI) == rep u
-------------------------------------------------------------------------
-- from AbelianSemiGroup
-------------------------------------------------------------------------
--(u: % + v: %): % == reducePm1DP per( rep(u)+rep(v) )
(u: % + v: %): % ==
-- der folgende Code gibt Probleme:
-- reducePm1DP( rep(u)+rep(v) ) Control stack exhausted (no more space
for function call frames).
-- das passt:
directProduct vector [reducePm1(ui+vi) for ui in entries u for vi in
entries v]
-------------------------------------------------------------------------
-- from CancellationAbelianMonoid
-------------------------------------------------------------------------
-------------------------------------------------------------------------
-- from Monoid
-------------------------------------------------------------------------
((u: %) * (r: NNI)): % ==
directProduct vector [reducePm1(ui*r) for ui in entries u]
((r: NNI) * (u: %)): % == u*r
-------------------------------------------------------------------------
-- from Ring
-------------------------------------------------------------------------
((u: %) * (v: %)): % ==
-- der folgende Code gibt Probleme:
-- reducePm1DP (rep(u) * rep(v))
-- das passt:
directProduct vector [reducePm1(ui*vi) for ui in entries u for vi in
entries v]
recip(z: %): Union(%, "failed") ==
pm1 : Integer := p-1
w : List NNI := new(n,0)$List(NNI)
for i in minIndex w .. maxIndex w repeat
zi : Integer := qelt(z, i) :: Integer
not gcd(zi, pm1) = 1 => return "failed"
izi : Integer := extendedEuclidean(zi, pm1).coef1
if izi < 0 then izi := izi + pm1
qsetelt!(w, i, izi :: NNI)
directProduct vector w
-------------------------------------------------------------------------
-- from OrderedSet
-------------------------------------------------------------------------
((x: %) < (y: %)): Boolean ==
for i in 1..n repeat
qelt(x,i) < qelt(y,i) => return true
qelt(x,i) > qelt(y,i) => return false
false
-------------------------------------------------------------------------
-- from OrderedAbelianMonoidSup
-------------------------------------------------------------------------
sup(u: %, v: %): % ==
-- der folgende Code gibt Probleme:
-- reducePm1DP sup( rep(u), rep(v))
-- das passt:
directProduct vector [reducePm1(sup(ui,vi)) for ui in entries u for vi in
entries v]
)endif
