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