1)abbrev category AMR AbelianMonoidRing
2++ Author:
3++ Basic Functions:
4++ Related Constructors:
5++ Also See:
6++ AMS Classifications:
7++ Keywords:
8++ References:
9++ Description:
10++ Abelian monoid ring elements (not necessarily of finite support)
11++ of this ring are of the form formal SUM (r_i * e_i)
12++ where the r_i are coefficents and the e_i, elements of the
13++ ordered abelian monoid, are thought of as exponents or monomials.
14++ The monomials commute with each other, but in general do not
15++ commute with the coefficients (which themselves may or may not
16++ be commutative).
17++ See \spadtype{FiniteAbelianMonoidRing} for the case of finite support.
18++ A useful common model for polynomials and power series.
19++ Conceptually at least, only the non-zero terms are ever operated on.
20
21AbelianMonoidRing(R : Join(SemiRng, AbelianMonoid), E : OrderedAbelianMonoid
22                 ) : Category == Join(SemiRng, BiModule(R, R),
23                     IndexedProductCategory(R, E)) with
24    degree : % -> E
25      ++ degree(p) returns the maximum of the exponents of the terms of p.
26    coefficient : (%, E) -> R
27         ++ coefficient(p, e) extracts the coefficient of the monomial with
28         ++ exponent e from polynomial p, or returns zero if exponent is not present.
29    if R has SemiRing then SemiRing
30    if R has Ring then Ring
31    if R has Field then "/": (%,R) -> %
32         ++ p/c divides p by the coefficient c.
33    if % has VariablesCommuteWithCoefficients and R has CommutativeRing then
34       CommutativeRing
35       Algebra R
36    if R has CharacteristicZero then CharacteristicZero
37    if R has CharacteristicNonZero then CharacteristicNonZero
38    if % has VariablesCommuteWithCoefficients and R has IntegralDomain then
39        IntegralDomain
40    if R has Algebra Fraction Integer then Algebra Fraction Integer
41  add
42    monomial? x == not(zero?(x)) and zero?(reductum(x))
43
44    map(fn : R -> R, x : %) ==
45          -- this default definition assumes that reductum is cheap
46       zero? x => 0
47       r := fn leadingCoefficient x
48       zero? r => map(fn, reductum x)
49       monomial(r, degree x) + map(fn, reductum x)
50
51    if R has Algebra Fraction Integer then
52      q : Fraction(Integer) * p : % == map(x1+->q*x1, p)
53
54)abbrev category FAMR FiniteAbelianMonoidRing
55++ Author:
56++ Basic Functions:
57++ Related Constructors:
58++ Also See:
59++ AMS Classifications:
60++ Keywords:
61++ References:
62++ Description: This category is
63++ similar to AbelianMonoidRing, except that the sum is assumed to be finite.
64++ It is a useful model for polynomials,
65++ but is somewhat more general.
66
67FiniteAbelianMonoidRing(R : Join(SemiRng, AbelianMonoid),
68                        E : OrderedAbelianMonoid) : Category ==
69   Join(AbelianMonoidRing(R, E), FreeModuleCategory(R, E),
70        FullyRetractableTo R) with
71    ground? : % -> Boolean
72      ++ ground?(p) tests if polynomial p is a member of the coefficient ring.
73                    -- can't be defined earlier, since a power series
74                    -- might not know if there were other terms or not
75    ground : % -> R
76      ++ ground(p) retracts polynomial p to the coefficient ring.
77    minimumDegree : % -> E
78      ++ minimumDegree(p) gives the least exponent of a non-zero term of polynomial p.
79      ++ Error: if applied to 0.
80    mapExponents : (E -> E, %) -> %
81      ++ mapExponents(fn, u) maps function fn onto the exponents
82      ++ of the non-zero monomials of polynomial u.
83    pomopo! : (%, R, E, %) -> %
84         ++ \spad{pomopo!(p1, r, e, p2)} returns \spad{p1 + monomial(r, e) * p2}
85         ++ and may use \spad{p1} as workspace. The constant \spad{r} is
86         ++ assumed to be nonzero.
87    if R has Ring then
88       fmecg : (%, E, R, %) -> %
89         ++ \spad{fmecg(p1, e, r, p2)} returns \spad{p1 - monomial(r, e) * p2}.
90    if % has CommutativeRing then
91       binomThmExpt : (%, %, NonNegativeInteger) -> %
92         ++ \spad{binomThmExpt(p, q, n)} returns \spad{(p+q)^n}
93         ++ by means of the binomial theorem trick.
94    if R has EntireRing then
95       "exquo": (%,R) -> Union(%,"failed")
96       ++ exquo(p,r) returns the exact quotient of polynomial p by r, or "failed"
97       ++ if none exists.
98       EntireRing
99    if R has GcdDomain then
100       content : % -> R
101         ++ content(p) gives the gcd of the coefficients of polynomial p.
102       primitivePart : % -> %
103         ++ primitivePart(p) returns the unit normalized form of polynomial p
104         ++ divided by the content of p.
105  add
106    pomopo!(p1, r, e, p2) == p1 + r * mapExponents(x1+->x1+e, p2)
107    if R has Ring then
108       fmecg(p1, e, r, p2) == p1 - r * mapExponents(x1+->x1+e, p2)
109
110    if % has CommutativeRing then
111       binomThmExpt(x, y, nn) ==
112               nn = 0 => 1$%
113               ans, xn, yn : %
114               bincoef : Integer
115               powl : List(%) := [x]
116               for i in 2..nn repeat powl := cons(x * powl.first, powl)
117               yn := y; ans := powl.first; i := 1; bincoef := nn
118               for xn in powl.rest repeat
119                  ans := bincoef * xn * yn + ans
120                  bincoef := (nn-i) * bincoef quo (i+1);  i := i+1
121                  -- last I and BINCOEF unused
122                  yn := y * yn
123               ans + yn
124    ground? x ==
125      retractIfCan(x)@Union(R,"failed") case "failed" => false
126      true
127    ground x == retract(x)@R
128    mapExponents (fn : E -> E, x : %) ==
129         -- this default definition assumes that reductum is cheap
130       zero? x => 0
131       monomial(leadingCoefficient x, fn degree x)+mapExponents(fn, reductum x)
132    coefficients x ==
133      zero? x => empty()
134      concat(leadingCoefficient x, coefficients reductum x)
135
136    if R has Field then
137       x/r == map(x1+->x1/r, x)
138    if R has EntireRing then
139       x : % exquo r : R ==
140          zero? r => error "Division by 0"
141          zero? x => 0
142          ans : % := 0
143          while not zero? x repeat
144              t := leadingCoefficient x exquo r
145              t case "failed" => return "failed"
146              ans := ans + monomial(t::R, degree x)
147              x := reductum x
148          ans
149    if R has GcdDomain then
150       content x ==       -- this assumes  reductum is cheap
151          zero? x => 0
152          r := leadingCoefficient x
153          x := reductum x
154          while not zero? x and not (r = 1) repeat
155            r := gcd(r, leadingCoefficient x)
156            x := reductum x
157          r
158       primitivePart x ==
159          zero? x => x
160          c := content x
161          unitCanonical((x exquo c)::%)
162
163
164
165)abbrev category GPOLCAT MaybeSkewPolynomialCategory
166++ Author:
167++ Basic Functions: Ring, monomial, coefficient, differentiate, eval
168++ Related Constructors: Polynomial, DistributedMultivariatePolynomial
169++ Also See: UnivariatePolynomialCategory
170++ AMS Classifications:
171++ Keywords:
172++ References:
173++ Description:
174++  The category for general multi-variate possibly skew polynomials
175++  over a ring R, in variables from VarSet, with exponents from the
176++  \spadtype{OrderedAbelianMonoidSup}.
177
178MaybeSkewPolynomialCategory(R : Join(SemiRng, AbelianMonoid),
179                            E : OrderedAbelianMonoidSup,
180                            VarSet : OrderedSet) : Category ==
181  FiniteAbelianMonoidRing(R, E) with
182    if R has Ring then
183        FullyLinearlyExplicitOver R
184    -- operations
185    degree : (%, VarSet) -> NonNegativeInteger
186      ++ degree(p, v) gives the degree of polynomial p with respect
187      ++ to the variable v.
188    degree : (%, List(VarSet)) -> List(NonNegativeInteger)
189      ++ degree(p, lv) gives the list of degrees of polynomial p
190      ++ with respect to each of the variables in the list lv.
191    coefficient : (%, VarSet, NonNegativeInteger) -> %
192      ++ coefficient(p, v, n) views the polynomial p as a univariate
193      ++ polynomial in v and returns the coefficient of the \spad{v^n} term.
194    coefficient : (%, List VarSet, List NonNegativeInteger) -> %
195      ++ coefficient(p, lv, ln) views the polynomial p as a polynomial
196      ++ in the variables of lv and returns the coefficient of the term
197      ++ \spad{lv^ln}, i.e. \spad{prod(lv_i ^ ln_i)}.
198    monomials : % -> List %
199      ++ monomials(p) returns the list of non-zero monomials of
200      ++ polynomial p, i.e.
201      ++ \spad{monomials(sum(a_(i) X^(i))) = [a_(1) X^(1), ..., a_(n) X^(n)]}.
202    mainVariable  : % -> Union(VarSet,"failed")
203      ++ mainVariable(p) returns the biggest variable which actually
204      ++ occurs in the polynomial p, or "failed" if no variables are
205      ++ present.
206      ++ fails precisely if polynomial satisfies ground?
207    monomial : (%, VarSet, NonNegativeInteger) -> %
208      ++ monomial(a, x, n) creates the monomial \spad{a*x^n} where \spad{a} is
209      ++ a polynomial, x is a variable and n is a nonnegative integer.
210    monomial : (%, List VarSet, List NonNegativeInteger) -> %
211      ++ monomial(a, [v1..vn], [e1..en]) returns \spad{a*prod(vi^ei)}.
212    totalDegree : % -> NonNegativeInteger
213      ++ totalDegree(p) returns the largest sum over all monomials
214      ++ of all exponents of a monomial.
215    totalDegree : (%, List VarSet) -> NonNegativeInteger
216      ++ totalDegree(p, lv) returns the maximum sum (over all monomials
217      ++ of polynomial p) of the variables in the list lv.
218    totalDegreeSorted : (%, List VarSet) -> NonNegativeInteger
219        ++ totalDegreeSorted(p, lv) returns the maximum sum (over all
220        ++ monomials of polynomial p) of the degree in variables in the
221        ++ list lv.  lv is assumed to be sorted in decreasing order.
222    variables : % -> List(VarSet)
223      ++ variables(p) returns the list of those variables actually
224      ++ appearing in the polynomial p.
225    if R has SemiRing then
226      primitiveMonomials : % -> List %
227        ++ primitiveMonomials(p) gives the list of monomials of the
228        ++ polynomial p with their coefficients removed.
229        ++ Note: \spad{primitiveMonomials(sum(a_(i) X^(i))) =
230        ++ [X^(1), ..., X^(n)]}.
231    if R has Comparable then Comparable
232
233    -- assertions
234    if R has canonicalUnitNormal then canonicalUnitNormal
235             ++ we can choose a unique representative for each
236             ++ associate class.
237             ++ This normalization is chosen to be normalization of
238             ++ leading coefficient (by default).
239 add
240    p : %
241    ln : List NonNegativeInteger
242    lv : List VarSet
243
244    monomials p ==
245--      sequential version for efficiency, by WMSIT, 7/30/90
246        ml := empty$List(%)
247        while p ~= 0 repeat
248            ml := concat(leadingMonomial p, ml)
249            p := reductum p
250        reverse! ml
251
252    monomial(p, lv, ln) ==
253        empty? lv =>
254            empty? ln => p
255            error "mismatched lists in monomial"
256        empty? ln  => error "mismatched lists in monomial"
257        monomial(monomial(p, first lv, first ln), rest lv, rest ln)
258
259    if R has SemiRing then
260
261        mkPrim(p : %) : % == monomial(1, degree p)
262
263        primitiveMonomials p ==
264            ml := empty$List(%)
265            while p ~= 0 repeat
266                ml := concat(mkPrim(leadingMonomial p), ml)
267                p := reductum p
268            reverse! ml
269
270
271)abbrev category POLYCAT PolynomialCategory
272++ Author:
273++ Basic Functions: Ring, monomial, coefficient, differentiate, eval
274++ Related Constructors: Polynomial, DistributedMultivariatePolynomial
275++ Also See: UnivariatePolynomialCategory
276++ AMS Classifications:
277++ Keywords:
278++ References:
279++ Description:
280++  The category for general multi-variate polynomials over a ring
281++  R, in variables from VarSet, with exponents from the
282++  \spadtype{OrderedAbelianMonoidSup}.  Here variables commute
283++  with the coefficients.
284
285PolynomialCategory(R : Join(SemiRng, AbelianMonoid),
286                   E : OrderedAbelianMonoidSup, VarSet : OrderedSet):
287        Category ==
288  Join(MaybeSkewPolynomialCategory(R, E, VarSet),
289       InnerEvalable(VarSet, R), InnerEvalable(VarSet, %),
290       VariablesCommuteWithCoefficients) with
291    if R has SemiRing then
292        RetractableTo VarSet
293        Evalable %
294    if R has Ring then
295        PartialDifferentialRing VarSet
296    -- operations
297    univariate   : (%, VarSet) -> SparseUnivariatePolynomial(%)
298      ++ univariate(p, v) converts the multivariate polynomial p
299      ++ into a univariate polynomial in v, whose coefficients are still
300      ++ multivariate polynomials (in all the other variables).
301    univariate   : % -> SparseUnivariatePolynomial(R)
302        ++ univariate(p) converts the multivariate polynomial p,
303        ++ which should actually involve only one variable,
304        ++ into a univariate polynomial
305        ++ in that variable, whose coefficients are in the ground ring.
306        ++ Error: if polynomial is genuinely multivariate
307    minimumDegree : (%, VarSet) -> NonNegativeInteger
308      ++ minimumDegree(p, v) gives the minimum degree of polynomial p
309      ++ with respect to v, i.e. viewed a univariate polynomial in v
310    minimumDegree : (%, List(VarSet)) -> List(NonNegativeInteger)
311      ++ minimumDegree(p, lv) gives the list of minimum degrees of the
312      ++ polynomial p with respect to each of the variables in the list lv
313    if R has Ring then
314      monicDivide : (%, %, VarSet) -> Record(quotient : %, remainder : %)
315       ++ monicDivide(a, b, v) divides the polynomial a by the polynomial b,
316       ++ with each viewed as a univariate polynomial in v returning
317       ++ both the quotient and remainder.
318       ++ Error: if b is not monic with respect to v.
319    multivariate : (SparseUnivariatePolynomial(R), VarSet) -> %
320        ++ multivariate(sup, v) converts an anonymous univariable
321        ++ polynomial sup to a polynomial in the variable v.
322    multivariate : (SparseUnivariatePolynomial(%), VarSet) -> %
323        ++ multivariate(sup, v) converts an anonymous univariable
324        ++ polynomial sup to a polynomial in the variable v.
325    isPlus : % -> Union(List %, "failed")
326        ++ isPlus(p) returns \spad{[m1, ..., mn]} if polynomial \spad{p = m1 + ... + mn} and
327        ++ \spad{n >= 2} and each mi is a nonzero monomial.
328    if R has SemiRing then
329      isTimes : % -> Union(List %, "failed")
330        ++ isTimes(p) returns \spad{[a1, ..., an]} if polynomial
331        ++ \spad{p = a1 ... an} and \spad{n >= 2}, and, for each i,
332        ++ ai is either a nontrivial constant in R or else of the
333        ++ form \spad{x^e}, where \spad{e > 0} is an integer and
334        ++ x is a member of VarSet.
335      isExpt : % -> Union(Record(var : VarSet, exponent : NonNegativeInteger),
336                       "failed")
337        ++ isExpt(p) returns \spad{[x, n]} if polynomial p has the form
338        ++ \spad{x^n} and \spad{n > 0}.
339    -- OrderedRing view removed to allow EXPR to define abs
340    --if R has OrderedRing then OrderedRing
341    if (R has ConvertibleTo InputForm) and
342       (VarSet has ConvertibleTo InputForm) then
343         ConvertibleTo InputForm
344    if R has Ring then
345      if (R has ConvertibleTo Pattern Integer) and
346       (VarSet has ConvertibleTo Pattern Integer) then
347         ConvertibleTo Pattern Integer
348      if (R has ConvertibleTo Pattern Float) and
349       (VarSet has ConvertibleTo Pattern Float) then
350         ConvertibleTo Pattern Float
351      if (R has PatternMatchable Integer) and
352       (VarSet has PatternMatchable Integer) then
353         PatternMatchable Integer
354      if (R has PatternMatchable Float) and
355       (VarSet has PatternMatchable Float) then
356         PatternMatchable Float
357    if R has CommutativeRing then
358      resultant : (%, %, VarSet) -> %
359         ++ resultant(p, q, v) returns the resultant of the polynomials
360         ++ p and q with respect to the variable v.
361      discriminant : (%, VarSet) -> %
362         ++ discriminant(p, v) returns the disriminant of the polynomial p
363         ++ with respect to the variable v.
364    if R has GcdDomain then
365      GcdDomain
366      content : (%, VarSet) -> %
367        ++ content(p, v) is the gcd of the coefficients of the polynomial p
368        ++ when p is viewed as a univariate polynomial with respect to the
369        ++ variable v.
370        ++ Thus, for polynomial 7*x^2*y + 14*x*y^2, the gcd of the
371        ++ coefficients with respect to x is 7*y.
372      primitivePart : % -> %
373        ++ primitivePart(p) returns the unitCanonical associate of the
374        ++ polynomial p with its content divided out.
375      primitivePart : (%, VarSet) -> %
376        ++ primitivePart(p, v) returns the unitCanonical associate of the
377        ++ polynomial p with its content with respect to the variable v
378        ++ divided out.
379      squareFree : % -> Factored %
380        ++ squareFree(p) returns the square free factorization of the
381        ++ polynomial p.
382      squareFreePart : % -> %
383        ++ squareFreePart(p) returns product of all the irreducible factors
384        ++ of polynomial p each taken with multiplicity one.
385
386    -- assertions
387    if R has PolynomialFactorizationExplicit then
388       PolynomialFactorizationExplicit
389 add
390
391    import from Factored(%)
392
393    p : %
394    v : VarSet
395    ln : List NonNegativeInteger
396    lv : List VarSet
397    n : NonNegativeInteger
398    pp, qq : SparseUnivariatePolynomial %
399    if R has SemiRing then
400        eval(p : %, l : List Equation %) ==
401            empty? l => p
402            for e in l repeat
403                retractIfCan(lhs e)@Union(VarSet,"failed") case "failed" =>
404                    error "cannot find a variable to evaluate"
405            lvar := [retract(lhs e)@VarSet for e in l]
406            eval(p, lvar, [rhs e for e in l]$List(%))
407
408    isPlus p ==
409      empty? rest(l := monomials p) => "failed"
410      l
411    if R has SemiRing then
412        isTimes p ==
413            empty?(lv := variables p) or not monomial? p => "failed"
414            l := [monomial(1, v, degree(p, v)) for v in lv]
415            ((r := leadingCoefficient p) = 1) =>
416                empty? rest lv => "failed"
417                l
418            concat(r::%, l)
419        isExpt p ==
420            (u := mainVariable p) case "failed" => "failed"
421            p = monomial(1, u::VarSet, d := degree(p, u::VarSet)) =>
422                [u::VarSet, d]
423            "failed"
424    coefficient(p, v, n) == coefficient(univariate(p, v), n)
425    coefficient(p, lv, ln) ==
426       empty? lv =>
427         empty? ln => p
428         error "mismatched lists in coefficient"
429       empty? ln  => error "mismatched lists in coefficient"
430       coefficient(coefficient(univariate(p, first lv), first ln),
431                   rest lv, rest ln)
432
433    retract(p : %) : VarSet ==
434      q := mainVariable(p)::VarSet
435      q::% = p => q
436      error "Polynomial is not a single variable"
437    retractIfCan(p:%):Union(VarSet, "failed") ==
438      ((q := mainVariable p) case VarSet) and (q::VarSet::% = p) => q
439      "failed"
440
441    totalDegree p ==
442        ground? p => 0
443        u := univariate(p, mainVariable(p)::VarSet)
444        d : NonNegativeInteger := 0
445        while u ~= 0 repeat
446          d := max(d, degree u + totalDegree leadingCoefficient u)
447          u := reductum u
448        d
449
450    totalDegreeSorted(p : %, lv : List VarSet) : NonNegativeInteger ==
451        ground? p => 0
452        empty?(lv) => 0
453        u := univariate(p, v := (mainVariable(p)::VarSet))
454        d : NonNegativeInteger := 0
455        w : NonNegativeInteger := 0
456        v0 := first(lv)
457        while v < v0 repeat
458            lv := rest(lv)
459            empty?(lv) => return 0
460            v0 := first(lv)
461        if v0 = v then
462            w := 1
463            lv := rest(lv)
464        while u ~= 0 repeat
465          d := max(d, w*(degree u) +
466                 totalDegreeSorted(leadingCoefficient u, lv))
467          u := reductum u
468        d
469
470    totalDegree(p, lv) ==
471        lv1 := sort((v1 : VarSet, v2 : VarSet) : Boolean +-> v2 < v1, lv)
472        totalDegreeSorted(p, lv1)
473
474    -- Condition on % is redundant, but compiler can not infer it
475    if % has CommutativeRing and R has CommutativeRing then
476        resultant(p1, p2, mvar) ==
477          resultant(univariate(p1, mvar), univariate(p2, mvar))
478        discriminant(p, var) ==
479          discriminant(univariate(p, var))
480
481    if R has IntegralDomain then
482      allMonoms(l : List %) : List(%) ==
483        removeDuplicates! concat [primitiveMonomials p for p in l]
484      P2R(p : %, b : List E, n : NonNegativeInteger) : Vector(R) ==
485        w := new(n, 0)$Vector(R)
486        for i in minIndex w .. maxIndex w for bj in b repeat
487          qsetelt!(w, i, coefficient(p, bj))
488        w
489      eq2R(l : List %, b : List E) : Matrix(R) ==
490        matrix [[coefficient(p, bj) for p in l] for bj in b]
491      reducedSystem(m : Matrix %) : Matrix(R) ==
492        l := listOfLists m
493        b := removeDuplicates!
494                           concat [allMonoms r for r in l]$List(List(%))
495        empty?(b) => new(0, ncols(m), 0)$Matrix(R)
496        d := [degree bj for bj in b]
497        mm := eq2R(first l, d)
498        l := rest l
499        while not empty? l repeat
500          mm := vertConcat(mm, eq2R(first l, d))
501          l := rest l
502        mm
503      reducedSystem(m : Matrix %, v : Vector %):
504       Record(mat : Matrix R, vec : Vector R) ==
505        l := listOfLists m
506        r := entries v
507        b : List % := removeDuplicates! concat(allMonoms r,
508                          concat [allMonoms s for s in l]$List(List(%)))
509        empty?(b) =>
510            [new(0, ncols(m), 0)$Matrix(R), new(0, 0)$Vector(R)]
511        d := [degree bj for bj in b]
512        n := #d
513        mm := eq2R(first l, d)
514        w := P2R(first r, d, n)
515        l := rest l
516        r := rest r
517        while not empty? l repeat
518          mm := vertConcat(mm, eq2R(first l, d))
519          w := concat(w, P2R(first r, d, n))
520          l := rest l
521          r := rest r
522        [mm, w]
523
524    if R has PolynomialFactorizationExplicit then
525       -- we might be in trouble if its actually only
526       -- a univariate polynomial category - have to remember to
527       -- over-ride these in UnivariatePolynomialCategory
528       PFBR ==>PolynomialFactorizationByRecursion(R, E, VarSet, %)
529       solveLinearPolynomialEquation(lpp, pp) ==
530           solveLinearPolynomialEquationByRecursion(lpp, pp)$PFBR
531       if R has FiniteFieldCategory then
532           MFFACT ==> MultFiniteFactorize(VarSet, E, R, %)
533           factorPolynomial(pp) == factor(pp)$MFFACT
534           factorSquareFreePolynomial(pp) == factor(pp)$MFFACT
535           factor(p) == factor(p)$MFFACT
536       else
537           if R has CharacteristicZero and R has EuclideanDomain then
538               MF ==> InnerMultFact(VarSet, E, R, %)
539               factorPolynomial(pp) ==
540                   factor(pp, factorPolynomial$R)$MF
541               factor(p) == factor(p, factorPolynomial$R)$MF
542               factorSquareFreePolynomial(pp) ==
543                   factor(pp, factorPolynomial$R)$MF
544           else
545               gcdPolynomial(pp, qq) ==
546                   gcdPolynomial(pp, qq
547                                )$GeneralPolynomialGcdPackage(E, VarSet, R, %)
548               factorPolynomial(pp) == factorByRecursion(pp)$PFBR
549               factorSquareFreePolynomial(pp) ==
550                   factorSquareFreeByRecursion(pp)$PFBR
551               factor p ==
552                   v : Union(VarSet, "failed") := mainVariable p
553                   v case "failed" =>
554                       ansR := factor leadingCoefficient p
555                       makeFR(unit(ansR)::%,
556                         [[w.flag, w.factor::%, w.exponent] for w in factorList ansR])
557                   up : SparseUnivariatePolynomial % := univariate(p, v)
558                   ansSUP := factorByRecursion(up)$PFBR
559                   makeFR(multivariate(unit(ansSUP), v),
560                     [[ww.flag, multivariate(ww.factor, v), ww.exponent]
561                       for ww in factorList ansSUP])
562
563       if R has CharacteristicNonZero then
564          mat : Matrix %
565          conditionP mat ==
566            ll := listOfLists transpose mat   -- hence each list corresponds to a
567                                            -- column, i.e. to one variable
568            llR : List List R := [ empty() for z in first ll]
569            monslist : List List % := empty()
570            ch : Integer := characteristic()$%
571            for l in ll repeat
572                mons := "setUnion"/[primitiveMonomials u for u in l]
573                redmons : List % := []
574                for m in mons repeat
575                    vars := variables m
576                    degs := degree(m, vars)
577                    deg1 : List NonNegativeInteger
578                    deg1 := [ ((nd := (d::Integer) exquo ch)
579                               case "failed" => return "failed" ;
580                                nd::Integer::NonNegativeInteger)
581                           for d in degs ]
582                    redmons := cons(monomial(1, vars, deg1), redmons)
583                    llR := [cons(ground coefficient(u, vars, degs), v)
584                              for u in l for v in llR]
585                monslist := cons(redmons, monslist)
586            ans := conditionP transpose matrix llR
587            ans case "failed" => "failed"
588            i : NonNegativeInteger := 0
589            [ +/[m*(ans.(i := i+1))::% for m in mons ]
590              for mons in monslist]
591
592    if R has CharacteristicNonZero then
593          charthRootlv : (%,List VarSet,NonNegativeInteger) -> Union(%,"failed")
594          charthRoot p ==
595            vars := variables p
596            empty? vars =>
597              ans := charthRoot ground p
598              ans case "failed" => "failed"
599              ans::R::%
600            ch := characteristic()$%
601            charthRootlv(p, vars, ch)
602          charthRootlv(p, vars, ch) ==
603            empty? vars =>
604              ans := charthRoot ground p
605              ans case "failed" => "failed"
606              ans::R::%
607            v := first vars
608            vars := rest vars
609            d := degree(p, v)
610            ans : % := 0
611            while (d>0) repeat
612               (dd := (d::Integer exquo ch::Integer)) case "failed" =>
613                      return "failed"
614               cp := coefficient(p, v, d)
615               p := p-monomial(cp, v, d)
616               ansx := charthRootlv(cp, vars, ch)
617               ansx case "failed" => return "failed"
618               d := degree(p, v)
619               ans := ans+monomial(ansx, v, dd::Integer::NonNegativeInteger)
620            ansx := charthRootlv(p, vars, ch)
621            ansx case "failed" => return "failed"
622            return ans+ansx
623
624    if R has Ring then
625        monicDivide(p1, p2, mvar) ==
626            result := monicDivide(univariate(p1, mvar), univariate(p2, mvar))
627            [multivariate(result.quotient, mvar),
628              multivariate(result.remainder, mvar)]
629
630
631    if R has GcdDomain then
632      if R has EuclideanDomain and R has CharacteristicZero then
633       squareFree p == squareFree(p)$MultivariateSquareFree(E, VarSet, R, %)
634      else
635        squareFree p == squareFree(p)$PolynomialSquareFree(VarSet, E, R, %)
636      squareFreePart p ==
637        unit(s := squareFree p) * */[f.factor for f in factorList s]
638      content(p, v) == content univariate(p, v)
639      primitivePart p ==
640        zero? p => p
641        unitNormal((p exquo content p) ::%).canonical
642      primitivePart(p, v) ==
643        zero? p => p
644        unitNormal((p exquo content(p, v)) ::%).canonical
645    if R has Comparable then
646        smaller?(p : %, q : %) ==
647            while p ~= 0 and q ~= 0 repeat
648                (dp := degree p) < (dq := degree q) =>
649                    return smaller?(0, leadingCoefficient q)
650                dq < dp => return smaller?(leadingCoefficient p, 0)
651                lp := leadingCoefficient(p)
652                lq := leadingCoefficient(q)
653                lp ~= lq => return smaller?(lp, lq)
654                p := reductum p
655                q := reductum q
656            p = 0 => smaller?(0, leadingCoefficient q)
657            smaller?(leadingCoefficient p, 0)
658
659        if R has Ring then
660          if (R has PatternMatchable Integer) and
661                (VarSet has PatternMatchable Integer) then
662              patternMatch(p : %, pat : Pattern Integer,
663                           l : PatternMatchResult(Integer, %)) ==
664                  patternMatch(p, pat, l
665                    )$PatternMatchPolynomialCategory(Integer, E, VarSet, R, %)
666          if (R has PatternMatchable Float) and
667                (VarSet has PatternMatchable Float) then
668              patternMatch(p : %, pat : Pattern Float,
669                           l : PatternMatchResult(Float, %)) ==
670                  patternMatch(p, pat, l
671                    )$PatternMatchPolynomialCategory(Float, E, VarSet, R, %)
672
673    if R has Ring then
674      RCX ==> Record(var : VarSet, exponent : NonNegativeInteger)
675      -- Must agree with 'patternMatch' in PatternMatchPolynomialCategory
676      -- so that x converted to pattern matches itself
677      conv_def(x, p_type) ==>
678          (r := retractIfCan(x)@Union(R, "failed")) case R =>
679              convert(r::R)
680          (vu := retractIfCan(x)@Union(VarSet, "failed")) case VarSet =>
681              convert(vu::VarSet)
682          (ex := isExpt x) case RCX =>
683              exr := ex::RCX
684              convert(exr.var)^(exr.exponent)
685          (lx := isTimes x) case List(%) =>
686              res : p_type := 1
687              for tx in reverse(lx::List(%)) repeat
688                  res := convert(tx)*res
689              res
690          (lx := isPlus x) case List(%) =>
691              res : p_type := 0
692              for mx in reverse(lx::List(%)) repeat
693                  res := convert(mx) + res
694              res
695          error "conv_def: impossible"
696
697      if (R has ConvertibleTo Pattern Integer) and
698            (VarSet has ConvertibleTo Pattern Integer) then
699          convert(x : %) : Pattern(Integer) ==
700              conv_def(x, Pattern(Integer))
701      if (R has ConvertibleTo Pattern Float) and
702            (VarSet has ConvertibleTo Pattern Float) then
703          convert(x : %) : Pattern(Float) ==
704              conv_def(x, Pattern(Float))
705    if (R has ConvertibleTo InputForm) and
706       (VarSet has ConvertibleTo InputForm) then
707         convert(p : %) : InputForm ==
708           map(convert, convert,
709                    p)$PolynomialCategoryLifting(E, VarSet, R, %, InputForm)
710
711
712
713)abbrev package POLYLIFT PolynomialCategoryLifting
714++ Author: Manuel Bronstein
715++ Basic Functions:
716++ Related Constructors:
717++ Also See:
718++ AMS Classifications:
719++ Keywords:
720++ References:
721++ Description:
722++ This package provides a very general map function, which
723++ given a set S and polynomials over R with maps from the
724++ variables into S and the coefficients into S, maps polynomials
725++ into S. S is assumed to support \spad{+}, \spad{*} and \spad{^}.
726
727PolynomialCategoryLifting(E, Vars, R, P, S) : Exports == Implementation where
728  E   : OrderedAbelianMonoidSup
729  Vars : OrderedSet
730  R   : Join(SemiRng, AbelianMonoid)
731  P   : PolynomialCategory(R, E, Vars)
732  S   : with
733           "+" : (%, %) -> %
734           "*" : (%, %) -> %
735           "^" : (%, NonNegativeInteger) -> %
736
737  Exports ==> with
738    map : (Vars -> S, R -> S, P) -> S
739      ++ map(varmap, coefmap, p) takes a
740      ++ varmap, a mapping from the variables of polynomial p into S,
741      ++ coefmap, a mapping from coefficients of p into S, and p, and
742      ++ produces a member of S using the corresponding arithmetic
743      ++ in S.
744
745  Implementation ==> add
746
747    if S has SemiRng then
748      map(fv, fc, p) ==
749          (x1 := mainVariable p) case "failed" => fc leadingCoefficient p
750          up := univariate(p, x1::Vars)
751          ml := reverse(monomials(up))
752          t  := fv(x1::Vars)
753          mon0 := first(ml)
754          ans : S
755          pow0 : S
756          i0 := degree(mon0)
757          lc0 := map(fv, fc, leadingCoefficient mon0)
758          if i0 > 0 then
759              pow0 := t^i0
760              ans := lc0*pow0
761          else
762              ans := lc0
763          for mon in rest(ml) repeat
764              i := degree(mon)
765              i = 0 => ans := ans + map(fv, fc, leadingCoefficient mon)
766              pow :=
767                  i0 <= 0 => t^i
768                  i > i0 => pow0*t^qcoerce(i - i0)@PositiveInteger
769                  t^i
770              i0 := i
771              pow0 := pow
772              ans := ans + map(fv, fc, leadingCoefficient mon)*pow
773          ans
774
775    else
776
777      -- preserve operations for use by InputForm
778      map(fv, fc, p) ==
779          (x1 := mainVariable p) case "failed" => fc leadingCoefficient p
780          up := univariate(p, x1::Vars)
781          t := fv(x1::Vars)
782          ans := map(fv, fc, leadingCoefficient up)*t^(degree up)
783          up := reductum up
784          while not zero? up repeat
785              ans := ans + map(fv, fc, leadingCoefficient up)*t^(degree up)
786              up  := reductum up
787          ans
788
789
790)abbrev category UPOLYC UnivariatePolynomialCategory
791++ Author:
792++ Basic Functions: Ring, monomial, coefficient, reductum, differentiate,
793++ elt, map, resultant, discriminant
794++ Related Constructors:
795++ Also See:
796++ AMS Classifications:
797++ Keywords:
798++ References:
799++ Description:
800++ The category of univariate polynomials over a ring R.
801++ No particular model is assumed - implementations can be either
802++ sparse or dense.
803
804UnivariatePolynomialCategory(R : Join(SemiRng, AbelianMonoid)) : Category ==
805 Join(PolynomialCategory(R, NonNegativeInteger, SingletonAsOrderedSet),
806      Eltable(R, R), Eltable(%, %)) with
807    if R has Ring then
808        DifferentialRing
809        DifferentialExtension R
810    vectorise : (%, NonNegativeInteger) -> Vector R
811      ++ vectorise(p, n) returns \spad{[a0, ..., a(n-1)]} where
812      ++ \spad{p = a0 + a1*x + ... + a(n-1)*x^(n-1)} + higher order terms.
813      ++ The degree of polynomial p can be different from \spad{n-1}.
814    unvectorise : Vector R -> %
815      ++ unvectorise(v) returns the polynomial which has for coefficients the
816      ++ entries of v in the increasing order.
817    makeSUP : % -> SparseUnivariatePolynomial R
818      ++ makeSUP(p) converts the polynomial p to be of type
819      ++ SparseUnivariatePolynomial over the same coefficients.
820    unmakeSUP : SparseUnivariatePolynomial R -> %
821      ++ unmakeSUP(sup) converts sup of type \spadtype{SparseUnivariatePolynomial(R)}
822      ++ to be a member of the given type.
823      ++ Note: converse of makeSUP.
824    multiplyExponents : (%, NonNegativeInteger) -> %
825      ++ multiplyExponents(p, n) returns a new polynomial resulting from
826      ++ multiplying all exponents of the polynomial p by the non negative
827      ++ integer n.
828    divideExponents : (%,NonNegativeInteger) -> Union(%,"failed")
829      ++ divideExponents(p, n) returns a new polynomial resulting from
830      ++ dividing all exponents of the polynomial p by the non negative
831      ++ integer n, or "failed" if some exponent is not exactly divisible
832      ++ by n.
833    shiftLeft : (%, NonNegativeInteger) -> %
834      ++ \spad{shiftLeft(p, n)} returns \spad{p * monomial(1, n)}
835    if R has Ring then
836      monicDivide : (%, %) -> Record(quotient : %, remainder : %)
837        ++ monicDivide(p, q) divide the polynomial p by the monic polynomial q,
838        ++ returning the pair \spad{[quotient, remainder]}.
839        ++ Error: if q isn't monic.
840      -- These two are for Karatsuba
841      karatsubaDivide : (%, NonNegativeInteger) -> Record(quotient : %, remainder : %)
842        ++ \spad{karatsubaDivide(p, n)} returns the same as
843        ++ \spad{monicDivide(p, monomial(1, n))}
844      shiftRight : (%, NonNegativeInteger) -> %
845        ++ \spad{shiftRight(p, n)} returns
846        ++ \spad{monicDivide(p, monomial(1, n)).quotient}
847      pseudoRemainder : (%, %) -> %
848        ++ pseudoRemainder(p, q) = r, for polynomials p and q, returns
849        ++ the remainder r when \spad{p' := p*lc(q)^(deg p - deg q + 1)}
850        ++ is pseudo right-divided by q, i.e. \spad{p' = s q + r}.
851      differentiate : (%, R -> R, %) -> %
852        ++ differentiate(p, d, x') extends the R-derivation d to an
853        ++ extension D in \spad{R[x]} where Dx is given by x', and
854        ++ returns \spad{Dp}.
855    if R has StepThrough then StepThrough
856    if R has CommutativeRing then
857      discriminant : % -> R
858        ++ discriminant(p) returns the discriminant of the polynomial p.
859      resultant : (%, %) -> R
860        ++ resultant(p, q) returns the resultant of the polynomials p and q.
861    if R has IntegralDomain then
862        Eltable(Fraction %, Fraction %)
863        elt  : (Fraction %, Fraction %) -> Fraction %
864             ++ elt(a, b) evaluates the fraction of univariate polynomials \spad{a}
865             ++ with the distinguished variable replaced by b.
866        order : (%, %) -> NonNegativeInteger
867             ++ order(p, q) returns the largest n such that \spad{q^n} divides polynomial p
868             ++ i.e. the order of \spad{p(x)} at \spad{q(x)=0}.
869        subResultantGcd : (%, %) -> %
870           ++ subResultantGcd(p, q) computes the gcd of the polynomials p and q
871           ++ using the SubResultant GCD algorithm.
872        composite : (%, %) -> Union(%, "failed")
873           ++ composite(p, q) returns h if \spad{p = h(q)}, and "failed" no such h exists.
874        composite : (Fraction %, %) -> Union(Fraction %, "failed")
875           ++ composite(f, q) returns h if f = h(q), and "failed" is no such h exists.
876        pseudoQuotient : (%, %) -> %
877           ++ pseudoQuotient(p, q) returns s, the quotient when
878           ++ \spad{p' := p*lc(q)^(deg p - deg q + 1)}
879           ++ is pseudo right-divided by q, i.e. \spad{p' = s q + r}.
880        pseudoDivide : (%, %) -> Record(coef : R, quotient : %, remainder : %)
881           ++ pseudoDivide(p, q) returns \spad{[c, s, r]}, when
882           ++ \spad{p' := p*lc(q)^(deg p - deg q + 1) = c * p}
883           ++ is pseudo right-divided by q, i.e. \spad{p' = s q + r}.
884    if R has GcdDomain then
885        separate : (%, %) -> Record(primePart : %, commonPart : %)
886           ++ separate(p, q) returns \spad{[a, b]} such that \spad{p = a b},
887           ++ \spad{a} is relatively prime to q and \spad{b} divides
888           ++ some power of \spad{q}.
889    if R has Field then
890        EuclideanDomain
891        additiveValuation
892          ++ euclideanSize(a*b) = euclideanSize(a) + euclideanSize(b)
893        elt : (Fraction %, R) -> R
894            ++ elt(a, r) evaluates the fraction of univariate polynomials \spad{a}
895            ++ with the distinguished variable replaced by the constant r.
896    if R has Algebra Fraction Integer then
897      integrate : % -> %
898        ++ integrate(p) integrates the univariate polynomial p with respect
899        ++ to its distinguished variable.
900  add
901
902    import from Integer
903    import from NonNegativeInteger
904
905    pp, qq : SparseUnivariatePolynomial %
906    variables(p) ==
907      zero? p or zero?(degree p) => []
908      [create()]
909    degree(p : %, v : SingletonAsOrderedSet) == degree p
910    totalDegree(p : %, lv : List SingletonAsOrderedSet) ==
911       empty? lv => 0
912       totalDegree p
913    degree(p : %, lv : List SingletonAsOrderedSet) ==
914       empty? lv => []
915       [degree p]
916    eval(p : %, lv : List SingletonAsOrderedSet, lq : List %) : % ==
917      empty? lv => p
918      not empty? rest lv => error "can only eval a univariate polynomial once"
919      eval(p, first lv, first lq)$%
920    eval(p : %, v : SingletonAsOrderedSet, q : %) : % == p(q)
921    eval(p : %, lv : List SingletonAsOrderedSet, lr : List R) : % ==
922      empty? lv => p
923      not empty? rest lv => error "can only eval a univariate polynomial once"
924      eval(p, first lv, first lr)$%
925    eval(p : %, v : SingletonAsOrderedSet, r : R) : % == p(r)::%
926    eval(p : %, le : List Equation %) : % ==
927      empty? le  => p
928      not empty? rest le => error "can only eval a univariate polynomial once"
929      mainVariable(lhs first le) case "failed" => p
930      p(rhs first le)
931    mainVariable(p : %) ==
932      zero? degree p => "failed"
933      create()$SingletonAsOrderedSet
934    minimumDegree(p : %, v : SingletonAsOrderedSet) == minimumDegree p
935    minimumDegree(p : %, lv : List SingletonAsOrderedSet) ==
936       empty? lv => []
937       [minimumDegree p]
938    monomial(p : %, v : SingletonAsOrderedSet, n : NonNegativeInteger) ==
939       mapExponents(x1+->x1+n, p)
940    if R has SemiRing then
941        coerce(v : SingletonAsOrderedSet) : % == monomial(1, 1)
942    makeSUP p ==
943      zero? p => 0
944      monomial(leadingCoefficient p, degree p) + makeSUP reductum p
945    unmakeSUP sp ==
946      zero? sp => 0
947      monomial(leadingCoefficient sp, degree sp) + unmakeSUP reductum sp
948
949    if R has PolynomialFactorizationExplicit then
950       PFBR ==> PolynomialFactorizationByRecursion(R, NonNegativeInteger,
951                                                   SingletonAsOrderedSet, %)
952       pp, qq : SparseUnivariatePolynomial %
953       lpp : List SparseUnivariatePolynomial %
954       SupR ==> SparseUnivariatePolynomial R
955       sp : SupR
956
957       solveLinearPolynomialEquation(lpp, pp) ==
958           solveLinearPolynomialEquationByRecursion(lpp, pp)$PFBR
959       factorPolynomial(pp) ==
960           factorByRecursion(pp)$PFBR
961       factorSquareFreePolynomial(pp) ==
962           factorSquareFreeByRecursion(pp)$PFBR
963       import from FactoredFunctions2(SupR, %)
964       factor p ==
965         zero? degree p  =>
966           ansR := factor leadingCoefficient p
967           makeFR(unit(ansR)::%,
968                  [[w.flag, w.factor::%, w.exponent] for w in factorList ansR])
969         map_preserving(unmakeSUP, factorPolynomial(makeSUP p)$R)
970
971    vectorise(p, n) ==
972      m := minIndex(v := new(n, 0)$Vector(R))
973      for i in minIndex v .. maxIndex v repeat
974        qsetelt!(v, i, coefficient(p, (i - m)::NonNegativeInteger))
975      v
976
977    unvectorise(v : Vector R) : % ==
978        p : % := 0
979        for i in 1..#v repeat
980            p := p + monomial(v(i), (i-1)::NonNegativeInteger)
981        p
982
983    retract(p : %) : R ==
984      zero? p => 0
985      zero? degree p => leadingCoefficient p
986      error "Polynomial is not of degree 0"
987    retractIfCan(p:%):Union(R, "failed") ==
988      zero? p => 0
989      zero? degree p => leadingCoefficient p
990      "failed"
991
992    if R has StepThrough then
993       init() == init()$R::%
994       nextItemInner : % -> Union(%,"failed")
995       nextItemInner(n) ==
996         zero? n => nextItem(0$R)::R::% -- assumed not to fail
997         zero? degree n =>
998           nn := nextItem leadingCoefficient n
999           nn case "failed" => "failed"
1000           nn::R::%
1001         n1 := reductum n
1002         n2 := nextItemInner n1 -- try stepping the reductum
1003         n2 case % => monomial(leadingCoefficient n, degree n) + n2
1004         1+degree n1 < degree n => -- there was a hole between lt n and n1
1005           monomial(leadingCoefficient n, degree n)+
1006             monomial(nextItem(init()$R)::R, 1+degree n1)
1007         n3 := nextItem leadingCoefficient n
1008         n3 case "failed" => "failed"
1009         monomial(n3, degree n)
1010       nextItem(n) ==
1011         n1 := nextItemInner n
1012         n1 case "failed" => monomial(nextItem(init()$R)::R,1+degree(n))
1013         n1
1014
1015    if R has GcdDomain then
1016
1017      content(p : %, v : SingletonAsOrderedSet) == content(p)::%
1018
1019      primeFactor : (%, %) -> %
1020
1021      primeFactor(p, q) ==
1022        (p1 := (p exquo gcd(p, q))::%) = p => p
1023        primeFactor(p1, q)
1024
1025      separate(p, q) ==
1026        a := primeFactor(p, q)
1027        [a, (p exquo a)::%]
1028
1029    if R has CommutativeRing then
1030      differentiate(x : %, deriv : R -> R, x': %) ==
1031        d : % := 0
1032        while (dg := degree x) > 0 repeat
1033          lc := leadingCoefficient x
1034          d := d + x' * monomial(dg * lc, (dg - 1)::NonNegativeInteger)
1035                 + monomial(deriv lc, dg)
1036          x := reductum x
1037        d + deriv(leadingCoefficient x)::%
1038    else if R has Ring then
1039      ncdiff : (NonNegativeInteger, %) -> %
1040      -- computes d(x^n) given dx = x', non-commutative case
1041      ncdiff(n, x') ==
1042        zero? n => 0
1043        zero?(n1 := (n - 1)::NonNegativeInteger) => x'
1044        x' * monomial(1, n1) + monomial(1, 1) * ncdiff(n1, x')
1045
1046      differentiate(x : %, deriv : R -> R, x': %) ==
1047        d : % := 0
1048        while (dg := degree x) > 0 repeat
1049          lc := leadingCoefficient x
1050          d := d + monomial(deriv lc, dg) + lc * ncdiff(dg, x')
1051          x := reductum x
1052        d + deriv(leadingCoefficient x)::%
1053    if R has Ring then
1054        differentiate(x : %, deriv : R -> R) == differentiate(x, deriv, 1$%)$%
1055        differentiate(x : %) : % ==
1056            d : % := 0
1057            while (dg := degree x) > 0 repeat
1058                d := d + monomial(dg*leadingCoefficient x,
1059                                 (dg - 1)::NonNegativeInteger)
1060                x := reductum x
1061            d
1062        differentiate(x : %, v : SingletonAsOrderedSet) : % == differentiate x
1063        pseudoRemainder(p, q) ==
1064           zero? q => error "PseudoDivision by Zero"
1065           zero? p => 0
1066           c2 := leadingCoefficient q
1067           e2 := degree q
1068           q := reductum q
1069           n := max((degree p)-e2+1, 0)::NonNegativeInteger
1070           while not zero? p repeat
1071               if (u := subtractIfCan(degree p, e2)) case "failed" then break
1072               p := fmecg(c2 * reductum p, u, leadingCoefficient p, q)
1073               n := (n - 1)::NonNegativeInteger
1074           n = 0 => p
1075           c2 ^ n * p
1076
1077    if R has IntegralDomain then
1078      elt(g : Fraction %, f : Fraction %) == ((numer g) f) / ((denom g) f)
1079
1080      pseudoQuotient(p, q) ==
1081        zero? q => error "PseudoDivision by Zero"
1082        zero? p => 0
1083        (degP, degQ) := (degree p, degree q)
1084        degP < degQ => 0
1085        lcQ := leadingCoefficient q
1086        q := reductum q
1087        i := (degP - degQ + 1)::NonNegativeInteger
1088        quot := 0$%
1089        while (delta : Integer := degree(p) - degQ) >= 0 repeat
1090          i := (i - 1)::NonNegativeInteger
1091          mon := monomial(leadingCoefficient p, delta::NonNegativeInteger)
1092          quot := quot + lcQ^i * mon
1093          p := lcQ * reductum(p) - mon * q
1094        quot
1095
1096      pseudoDivide(p, q) ==
1097        zero? q => error "PseudoDivision by Zero"
1098        zero? p => [1, 0, p]
1099        (degP, degQ) := (degree p, degree q)
1100        degP < degQ => [1, 0, p]
1101        lcQ := leadingCoefficient q
1102        q := reductum q
1103        i := (degP - degQ + 1)::NonNegativeInteger
1104        co := lcQ^i
1105        quot := 0$%
1106        while (delta : Integer := degree(p) - degQ) >= 0 repeat
1107          i := (i - 1)::NonNegativeInteger
1108          mon := monomial(leadingCoefficient p, delta::NonNegativeInteger)
1109          quot := quot + lcQ^i * mon
1110          p := lcQ * reductum(p) - mon * q
1111        p := lcQ^i * p
1112        [co, quot, p]
1113
1114      composite(f : Fraction %, q : %) ==
1115        (n := composite(numer f, q)) case "failed" => "failed"
1116        (d := composite(denom f, q)) case "failed" => "failed"
1117        n::% / d::%
1118
1119      composite(p : %, q : %) ==
1120        ground? p => p
1121        cqr := pseudoDivide(p, q)
1122        ground?(cqr.remainder) and
1123          ((v := cqr.remainder exquo cqr.coef) case %) and
1124            ((u := composite(cqr.quotient, q)) case %) and
1125              ((w := (u::%) exquo cqr.coef) case %) =>
1126                v::% + monomial(1, 1) * w::%
1127        "failed"
1128
1129      elt(p : %, f : Fraction %) ==
1130        zero? p => 0
1131        ans : Fraction(%) := (leadingCoefficient p)::%::Fraction(%)
1132        n := degree p
1133        while not zero?(p := reductum p) repeat
1134          ans := ans * f ^ (n - (n := degree p))::NonNegativeInteger +
1135                    (leadingCoefficient p)::%::Fraction(%)
1136        zero? n => ans
1137        ans * f ^ n
1138
1139      order(p, q) ==
1140        zero? p => error "order: arguments must be nonzero"
1141        degree(q) < 1 => error "order: place must be non-trivial"
1142        ans : NonNegativeInteger := 0
1143        repeat
1144          (u  := p exquo q) case "failed" => return ans
1145          p   := u::%
1146          ans := ans + 1
1147
1148    if R has GcdDomain then
1149      squareFree(p : %) ==
1150        squareFree(p)$UnivariatePolynomialSquareFree(R, %)
1151
1152      squareFreePart(p : %) ==
1153        squareFreePart(p)$UnivariatePolynomialSquareFree(R, %)
1154
1155    if R has PolynomialFactorizationExplicit then
1156
1157      gcdPolynomial(pp, qq) ==
1158            zero? pp => unitCanonical qq  -- subResultantGcd can't handle 0
1159            zero? qq => unitCanonical pp
1160            unitCanonical(gcd(content (pp), content(qq))*
1161                   primitivePart
1162                      subResultantGcd(primitivePart pp, primitivePart qq))
1163
1164      squareFreePolynomial pp ==
1165         squareFree(pp)$UnivariatePolynomialSquareFree(%,
1166                                    SparseUnivariatePolynomial %)
1167
1168    if R has Field then
1169      elt(f : Fraction %, r : R) == ((numer f) r) / ((denom f) r)
1170
1171      euclideanSize x ==
1172            zero? x =>
1173              error "euclideanSize called on 0 in Univariate Polynomial"
1174            degree x
1175      divide(x, y) ==
1176            zero? y => error "division by 0 in Univariate Polynomials"
1177            quot := 0
1178            lc := inv leadingCoefficient y
1179            while not zero?(x) and (degree x >= degree y) repeat
1180               f := lc*leadingCoefficient x
1181               n := (degree x - degree y)::NonNegativeInteger
1182               quot := quot+monomial(f, n)
1183               x := x-monomial(f, n)*y
1184            [quot, x]
1185    if R has Algebra Fraction Integer then
1186      integrate p ==
1187        ans : % := 0
1188        while p ~= 0 repeat
1189          l := leadingCoefficient p
1190          d := 1 + degree p
1191          ans := ans + inv(d::Fraction(Integer)) * monomial(l, d)
1192          p := reductum p
1193        ans
1194
1195
1196)abbrev package UPOLYC2 UnivariatePolynomialCategoryFunctions2
1197++ Author:
1198++ Basic Functions:
1199++ Related Constructors:
1200++ Also See:
1201++ AMS Classifications:
1202++ Keywords:
1203++ References:
1204++ Description:
1205++ Mapping from polynomials over R to polynomials over S
1206++ given a map from R to S assumed to send zero to zero.
1207
1208UnivariatePolynomialCategoryFunctions2(R, PR, S, PS) : Exports == Impl where
1209  R, S : Join(SemiRng, AbelianMonoid)
1210  PR  : UnivariatePolynomialCategory R
1211  PS  : UnivariatePolynomialCategory S
1212
1213  Exports ==> with
1214    map : (R -> S, PR) -> PS
1215     ++ map(f, p) takes a function f from R to S,
1216     ++ and applies it to each (non-zero) coefficient of a polynomial p
1217     ++ over R, getting a new polynomial over S.
1218     ++ Note: since the map is not applied to zero elements, it may map zero
1219     ++ to zero.
1220
1221  Impl ==> add
1222
1223    op_of_PS : Symbol := CAR(devaluate(PS)$Lisp)$Lisp
1224
1225    if op_of_PS = 'UnivariatePolynomial
1226         or PS is SparseUnivariatePolynomial(S) then
1227
1228        TermPS ==>  Record(k : Integer, c : S)
1229        RepPS :=  List TermPS
1230
1231        map(f, p) ==
1232            ans : RepPS := []
1233            while not(p = 0) repeat
1234                nc := f leadingCoefficient p
1235                if nc ~= 0 then
1236                    ans := cons([degree p, nc]$TermPS, ans)
1237                p := reductum p
1238            reverse!(ans) pretend PS
1239    else
1240        map(f, p) ==
1241            ans : PS := 0
1242            while not(p = 0) repeat
1243                ans := ans + monomial(f leadingCoefficient p, degree p)
1244                p   := reductum p
1245            ans
1246
1247)abbrev package COMMUPC CommuteUnivariatePolynomialCategory
1248++ Author: Manuel Bronstein
1249++ Basic Functions:
1250++ Related Constructors:
1251++ Also See:
1252++ AMS Classifications:
1253++ Keywords:
1254++ References:
1255++ Description:
1256++ A package for swapping the order of two variables in a tower of two
1257++ UnivariatePolynomialCategory extensions.
1258
1259CommuteUnivariatePolynomialCategory(R, UP, UPUP) : Exports == Impl where
1260  R   : Ring
1261  UP  : UnivariatePolynomialCategory R
1262  UPUP : UnivariatePolynomialCategory UP
1263
1264  N ==> NonNegativeInteger
1265
1266  Exports ==> with
1267    swap : UPUP -> UPUP
1268      ++ swap(p(x, y)) returns p(y, x).
1269
1270  Impl ==> add
1271    makePoly : (UP, N) -> UPUP
1272
1273-- converts P(x, y) to P(y, x)
1274    swap poly ==
1275      ans : UPUP := 0
1276      while poly ~= 0 repeat
1277        ans  := ans + makePoly(leadingCoefficient poly, degree poly)
1278        poly := reductum poly
1279      ans
1280
1281    makePoly(poly, d) ==
1282      ans : UPUP := 0
1283      while poly ~= 0 repeat
1284        ans  := ans +
1285             monomial(monomial(leadingCoefficient poly, d), degree poly)
1286        poly := reductum poly
1287      ans
1288
1289--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
1290--All rights reserved.
1291--
1292--Redistribution and use in source and binary forms, with or without
1293--modification, are permitted provided that the following conditions are
1294--met:
1295--
1296--    - Redistributions of source code must retain the above copyright
1297--      notice, this list of conditions and the following disclaimer.
1298--
1299--    - Redistributions in binary form must reproduce the above copyright
1300--      notice, this list of conditions and the following disclaimer in
1301--      the documentation and/or other materials provided with the
1302--      distribution.
1303--
1304--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
1305--      names of its contributors may be used to endorse or promote products
1306--      derived from this software without specific prior written permission.
1307--
1308--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
1309--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
1310--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
1311--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
1312--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
1313--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
1314--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
1315--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
1316--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
1317--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
1318--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
1319