\documentclass{article} \usepackage{open-axiom} \begin{document} \title{\$SPAD/src/algebra fr.spad} \author{Robert S. Sutor, Johnannes Grabmeier} \maketitle \begin{abstract} \end{abstract} \eject \tableofcontents \eject -- This file contains a domain and packages for manipulating objects -- in factored form. \section{domain FR Factored} <>= )abbrev domain FR Factored ++ Author: Robert S. Sutor ++ Date Created: 1985 ++ Date Last Modified: May 19, 2013. ++ Change History: ++ 21 Jan 1991 J Grabmeier Corrected a bug in exquo. ++ 16 Aug 1994 R S Sutor Improved convert to InputForm ++ Basic Operations: ++ expand, exponent, factorList, factors, flagFactor, irreducibleFactor, ++ makeFR, map, nilFactor, nthFactor, nthFlag, numberOfFactors, ++ primeFactor, sqfrFactor, unit, unitNormalize, ++ Related Constructors: FactoredFunctionUtilities, FactoredFunctions2 ++ Also See: ++ AMS Classifications: 11A51, 11Y05 ++ Keywords: factorization, prime, square-free, irreducible, factor ++ References: ++ Description: ++ \spadtype{Factored} creates a domain whose objects are kept in ++ factored form as long as possible. Thus certain operations like ++ multiplication and gcd are relatively easy to do. Others, like ++ addition require somewhat more work, and unless the argument ++ domain provides a factor function, the result may not be ++ completely factored. Each object consists of a unit and a list of ++ factors, where a factor has a member of R (the "base"), and ++ exponent and a flag indicating what is known about the base. A ++ flag may be one of "nil", "sqfr", "irred" or "prime", which respectively mean ++ that nothing is known about the base, it is square-free, it is ++ irreducible, or it is prime. The current ++ restriction to integral domains allows simplification to be ++ performed without worrying about multiplication order. Factored(R: IntegralDomain): Exports == Implementation where fUnion ==> Union("nil", "sqfr", "irred", "prime") FF ==> Record(flg: fUnion, fctr: R, xpnt: Integer) SRFE ==> Set(Record(factor:R, exponent:Integer)) Exports == Join(IntegralDomain, DifferentialExtension R, Algebra R, FullyEvalableOver R, FullyRetractableTo R,Functorial R) with expand: % -> R ++ expand(f) multiplies the unit and factors together, yielding an ++ "unfactored" object. Note: this is purposely not called \spadfun{coerce} which would ++ cause the interpreter to do this automatically. exponent: % -> Integer ++ exponent(u) returns the exponent of the first factor of ++ \spadvar{u}, or 0 if the factored form consists solely of a unit. makeFR : (R, List FF) -> % ++ makeFR(unit,listOfFactors) creates a factored object (for ++ use by factoring code). factorList : % -> List FF ++ factorList(u) returns the list of factors with flags (for ++ use by factoring code). nilFactor: (R, Integer) -> % ++ nilFactor(base,exponent) creates a factored object with ++ a single factor with no information about the kind of ++ base (flag = "nil"). factors: % -> List Record(factor:R, exponent:Integer) ++ factors(u) returns a list of the factors in a form suitable ++ for iteration. That is, it returns a list where each element ++ is a record containing a base and exponent. The original ++ object is the product of all the factors and the unit (which ++ can be extracted by \axiom{unit(u)}). irreducibleFactor: (R, Integer) -> % ++ irreducibleFactor(base,exponent) creates a factored object with ++ a single factor whose base is asserted to be irreducible ++ (flag = "irred"). nthExponent: (%, Integer) -> Integer ++ nthExponent(u,n) returns the exponent of the nth factor of ++ \spadvar{u}. If \spadvar{n} is not a valid index for a factor ++ (for example, less than 1 or too big), 0 is returned. nthFactor: (%,Integer) -> R ++ nthFactor(u,n) returns the base of the nth factor of ++ \spadvar{u}. If \spadvar{n} is not a valid index for a factor ++ (for example, less than 1 or too big), 1 is returned. If ++ \spadvar{u} consists only of a unit, the unit is returned. nthFlag: (%,Integer) -> fUnion ++ nthFlag(u,n) returns the information flag of the nth factor of ++ \spadvar{u}. If \spadvar{n} is not a valid index for a factor ++ (for example, less than 1 or too big), "nil" is returned. numberOfFactors : % -> NonNegativeInteger ++ numberOfFactors(u) returns the number of factors in \spadvar{u}. primeFactor: (R,Integer) -> % ++ primeFactor(base,exponent) creates a factored object with ++ a single factor whose base is asserted to be prime ++ (flag = "prime"). sqfrFactor: (R,Integer) -> % ++ sqfrFactor(base,exponent) creates a factored object with ++ a single factor whose base is asserted to be square-free ++ (flag = "sqfr"). flagFactor: (R,Integer, fUnion) -> % ++ flagFactor(base,exponent,flag) creates a factored object with ++ a single factor whose base is asserted to be properly ++ described by the information flag. unit: % -> R ++ unit(u) extracts the unit part of the factorization. unitNormalize: % -> % ++ unitNormalize(u) normalizes the unit part of the factorization. ++ For example, when working with factored integers, this operation will ++ ensure that the bases are all positive integers. -- the following operations are conditional on R if R has GcdDomain then GcdDomain if R has RealConstant then RealConstant if R has UniqueFactorizationDomain then UniqueFactorizationDomain if R has ConvertibleTo InputForm then ConvertibleTo InputForm if R has IntegerNumberSystem then rational? : % -> Boolean ++ rational?(u) tests if \spadvar{u} is actually a ++ rational number (see \spadtype{Fraction Integer}). rational : % -> Fraction Integer ++ rational(u) assumes spadvar{u} is actually a rational number ++ and does the conversion to rational number ++ (see \spadtype{Fraction Integer}). rationalIfCan: % -> Union(Fraction Integer, "failed") ++ rationalIfCan(u) returns a rational number if u ++ really is one, and "failed" otherwise. if R has Eltable(%, %) then Eltable(%, %) if R has Evalable(%) then Evalable(%) if R has InnerEvalable(Symbol, %) then InnerEvalable(Symbol, %) Implementation ==> add -- Representation: -- Note: exponents are allowed to be integers so that some special cases -- may be used in simplications Rep := Record(unt:R, fct:List FF) if R has ConvertibleTo InputForm then convert(x:%):InputForm == empty?(lf := reverse factorList x) => convert(unit x)@InputForm l := empty()$List(InputForm) for rec in lf repeat ((rec.fctr) = 1) => messagePrint("WARNING (convert$Factored):_ 1 should not appear as factor.")$OutputForm iFactor : InputForm := binary( convert("::" :: Symbol)@InputForm, [convert(rec.fctr)@InputForm, (devaluate R)$Lisp :: InputForm ]$List(InputForm) ) iExpon : InputForm := convert(rec.xpnt)@InputForm iFun : List InputForm := rec.flg case "nil" => [convert("nilFactor" :: Symbol)@InputForm, iFactor, iExpon]$List(InputForm) rec.flg case "sqfr" => [convert("sqfrFactor" :: Symbol)@InputForm, iFactor, iExpon]$List(InputForm) rec.flg case "prime" => [convert("primeFactor" :: Symbol)@InputForm, iFactor, iExpon]$List(InputForm) rec.flg case "irred" => [convert("irreducibleFactor" :: Symbol)@InputForm, iFactor, iExpon]$List(InputForm) nil$List(InputForm) l := concat( iFun pretend InputForm, l ) -- one?(rec.xpnt) => -- l := concat(convert(rec.fctr)@InputForm, l) -- l := concat(convert(rec.fctr)@InputForm ** rec.xpnt, l) empty? l => convert(unit x)@InputForm if not one? unit x then l := concat(convert(unit x)@InputForm,l) empty? rest l => first l binary(convert(_*::Symbol)@InputForm, l)@InputForm -- Private function signatures: reciprocal : % -> % qexpand : % -> R negexp? : % -> Boolean SimplifyFactorization : List FF -> List FF LispLessP : (FF, FF) -> Boolean mkFF : (R, List FF) -> % SimplifyFactorization1 : (FF, List FF) -> List FF stricterFlag : (fUnion, fUnion) -> fUnion nilFactor(r, i) == flagFactor(r, i, "nil") sqfrFactor(r, i) == flagFactor(r, i, "sqfr") irreducibleFactor(r, i) == flagFactor(r, i, "irred") primeFactor(r, i) == flagFactor(r, i, "prime") unit? u == (empty? u.fct) and (not zero? u.unt) factorList u == u.fct unit u == u.unt numberOfFactors u == # u.fct 0 == [1, [["nil", 0, 1]$FF]] zero? u == # u.fct = 1 and (first u.fct).flg case "nil" and zero? (first u.fct).fctr and one? u.unt 1 == [1, empty()] one? u == empty? u.fct and u.unt = 1 mkFF(r, x) == [r, x] coerce(j:Integer):% == (j::R)::% characteristic == characteristic$R i:Integer * u:% == (i :: %) * u r:R * u:% == (r :: %) * u factors u == [[fe.fctr, fe.xpnt] for fe in factorList u] expand u == retract u negexp? x == "or"/[negative?(y.xpnt) for y in factorList x] makeFR(u, l) == -- normalizing code to be installed when contents are handled better -- current squareFree returns the content as a unit part. -- if (not unit?(u)) then -- l := cons(["nil", u, 1]$FF,l) -- u := 1 unitNormalize mkFF(u, SimplifyFactorization l) if R has IntegerNumberSystem then rational? x == true rationalIfCan x == rational x rational x == convert(unit x)@Integer * _*/[(convert(f.fctr)@Integer)::Fraction(Integer) ** f.xpnt for f in factorList x] if R has Eltable(R, R) then elt(x:%, v:%) == x(expand v) if R has Evalable(R) then eval(x:%, l:List Equation %) == eval(x,[expand lhs e = expand rhs e for e in l]$List(Equation R)) if R has InnerEvalable(Symbol, R) then eval(x:%, ls:List Symbol, lv:List %) == eval(x, ls, [expand v for v in lv]$List(R)) if R has RealConstant then --! negcount and rest commented out since RealConstant doesn't support --! positive? or negative? -- negcount: % -> Integer -- positive?(x:%):Boolean == not(zero? x) and even?(negcount x) -- negative?(x:%):Boolean == not(zero? x) and odd?(negcount x) -- negcount x == -- n := count(negative?(#1.fctr), factorList x)$List(FF) -- negative? unit x => n + 1 -- n convert(x:%):Float == convert(unit x)@Float * _*/[convert(f.fctr)@Float ** f.xpnt for f in factorList x] convert(x:%):DoubleFloat == convert(unit x)@DoubleFloat * _*/[convert(f.fctr)@DoubleFloat ** f.xpnt for f in factorList x] u:% * v:% == zero? u or zero? v => 0 one? u => v one? v => u mkFF(unit u * unit v, SimplifyFactorization concat(factorList u, copy factorList v)) u:% ** n:NonNegativeInteger == mkFF(unit(u)**n, [[x.flg, x.fctr, n * x.xpnt] for x in factorList u]) SimplifyFactorization x == empty? x => empty() x := sort!(LispLessP, x) x := SimplifyFactorization1(first x, rest x) x := sort!(LispLessP, x) x SimplifyFactorization1(f, x) == empty? x => zero?(f.xpnt) => empty() list f f1 := first x f.fctr = f1.fctr => SimplifyFactorization1([stricterFlag(f.flg, f1.flg), f.fctr, f.xpnt + f1.xpnt], rest x) l := SimplifyFactorization1(first x, rest x) zero?(f.xpnt) => l concat(f, l) coerce(x:%):OutputForm == empty?(lf := reverse factorList x) => (unit x)::OutputForm l := empty()$List(OutputForm) for rec in lf repeat ((rec.fctr) = 1) => messagePrint "WARNING (coerce$Factored): 1 should not appear as factor." ((rec.xpnt) = 1) => l := concat(rec.fctr :: OutputForm, l) l := concat(rec.fctr::OutputForm ** rec.xpnt::OutputForm, l) empty? l => (unit x) :: OutputForm e := empty? rest l => first l reduce(_*, l) 1 = unit x => e (unit x)::OutputForm * e retract(u:%):R == negexp? u => error "Negative exponent in factored object" qexpand u qexpand u == unit u * _*/[y.fctr ** (y.xpnt::NonNegativeInteger) for y in factorList u] retractIfCan(u:%):Union(R, "failed") == negexp? u => "failed" qexpand u LispLessP(y, y1) == before?(y.fctr, y1.fctr) stricterFlag(fl1, fl2) == fl1 case "prime" => fl1 fl1 case "irred" => fl2 case "prime" => fl2 fl1 fl1 case "sqfr" => fl2 case "nil" => fl1 fl2 fl2 if R has IntegerNumberSystem then coerce(r:R):% == factor(r)$IntegerFactorizationPackage(R) pretend % else if R has UniqueFactorizationDomain then coerce(r:R):% == zero? r => 0 unit? r => mkFF(r, empty()) unitNormalize(squareFree(r) pretend %) else coerce(r:R):% == one? r => 1 unitNormalize mkFF(1, [["nil", r, 1]$FF]) u = v == (unit u = unit v) and # u.fct = # v.fct and set(factors u)$SRFE =$SRFE set(factors v)$SRFE - u == zero? u => u mkFF(- unit u, factorList u) recip u == not empty? factorList u => "failed" (r := recip unit u) case "failed" => "failed" mkFF(r::R, empty()) reciprocal u == mkFF((recip unit u)::R, [[y.flg, y.fctr, - y.xpnt]$FF for y in factorList u]) exponent u == -- exponent of first factor empty?(fl := factorList u) or zero? u => 0 first(fl).xpnt nthExponent(u, i) == l := factorList u zero? u or i < 1 or i > #l => 0 (l.(minIndex(l) + i - 1)).xpnt nthFactor(u, i) == zero? u => 0 zero? i => unit u l := factorList u negative? i or i > #l => 1 (l.(minIndex(l) + i - 1)).fctr nthFlag(u, i) == l := factorList u zero? u or i < 1 or i > #l => "nil" (l.(minIndex(l) + i - 1)).flg flagFactor(r, i, fl) == zero? i => 1 zero? r => 0 unitNormalize mkFF(1, [[fl, r, i]$FF]) differentiate(u:%, deriv: R -> R) == ans := deriv(unit u) * ((u exquo unit(u)::%)::%) ans + (_+/[fact.xpnt * deriv(fact.fctr) * ((u exquo nilFactor(fact.fctr, 1))::%) for fact in factorList u]) @ This operation provides an implementation of [[differentiate]] from the category [[DifferentialExtension]]. It uses the formula $$\frac{d}{dx} f(x) = \sum_{i=1}^n \frac{f(x)}{f_i(x)}\frac{d}{dx}f_i(x),$$ where $$f(x)=\prod_{i=1}^n f_i(x).$$ Note that up to [[patch--40]] the following wrong definition was used: \begin{verbatim} differentiate(u:%, deriv: R -> R) == ans := deriv(unit u) * ((u exquo (fr := unit(u)::%))::%) ans + fr * (_+/[fact.xpnt * deriv(fact.fctr) * ((u exquo nilFactor(fact.fctr, 1))::%) for fact in factorList u]) \end{verbatim} which causes wrong results as soon as units are involved, for example in <>= D(factor (-x), x) @ (Issue~\#176) <>= map(fn, u) == fn(unit u) * _*/[irreducibleFactor(fn(f.fctr),f.xpnt) for f in factorList u] u exquo v == empty?(x1 := factorList v) => unitNormal(retract v).associate * u empty? factorList u => "failed" v1 := u * reciprocal v goodQuotient:Boolean := true while (goodQuotient and (not empty? x1)) repeat if negative? x1.first.xpnt then goodQuotient := false else x1 := rest x1 goodQuotient => v1 "failed" unitNormal u == -- does a bunch of work, but more canonical (ur := recip(un := unit u)) case "failed" => [1, u, 1] as := ur::R vl := empty()$List(FF) for x in factorList u repeat ucar := unitNormal(x.fctr) e := abs(x.xpnt)::NonNegativeInteger if negative? x.xpnt then -- associate is recip of unit un := un * (ucar.associate ** e) as := as * (ucar.unit ** e) else un := un * (ucar.unit ** e) as := as * (ucar.associate ** e) if not one?(ucar.canonical) then vl := concat([x.flg, ucar.canonical, x.xpnt], vl) [mkFF(un, empty()), mkFF(1, reverse! vl), mkFF(as, empty())] unitNormalize u == uca := unitNormal u mkFF(unit(uca.unit)*unit(uca.canonical),factorList(uca.canonical)) if R has GcdDomain then u + v == zero? u => v zero? v => u v1 := reciprocal(u1 := gcd(u, v)) (expand(u * v1) + expand(v * v1)) * u1 gcd(u, v) == one? u or one? v => 1 zero? u => v zero? v => u f1 := empty()$List(Integer) -- list of used factor indices in x f2 := f1 -- list of indices corresponding to a given factor f3 := empty()$List(List Integer) -- list of f2-like lists x := concat(factorList u, factorList v) for i in minIndex x .. maxIndex x repeat if not member?(i, f1) then f1 := concat(i, f1) f2 := [i] for j in i+1..maxIndex x repeat if x.i.fctr = x.j.fctr then f1 := concat(j, f1) f2 := concat(j, f2) f3 := concat(f2, f3) x1 := empty()$List(FF) while not empty? f3 repeat f1 := first f3 if #f1 > 1 then i := first f1 y := copy x.i f1 := rest f1 while not empty? f1 repeat i := first f1 if x.i.xpnt < y.xpnt then y.xpnt := x.i.xpnt f1 := rest f1 x1 := concat(y, x1) f3 := rest f3 x1 := sort!(LispLessP, x1) mkFF(1, x1) else -- R not a GCD domain u + v == zero? u => v zero? v => u irreducibleFactor(expand u + expand v, 1) if R has UniqueFactorizationDomain then prime? u == not(empty?(l := factorList u)) and (empty? rest l) and one?(l.first.xpnt) and (l.first.flg case "prime") @ \section{package FRUTIL FactoredFunctionUtilities} <>= )abbrev package FRUTIL FactoredFunctionUtilities ++ Author: ++ Date Created: ++ Change History: ++ Basic Operations: refine, mergeFactors ++ Related Constructors: Factored ++ Also See: ++ AMS Classifications: 11A51, 11Y05 ++ Keywords: factor ++ References: ++ Description: ++ \spadtype{FactoredFunctionUtilities} implements some utility ++ functions for manipulating factored objects. FactoredFunctionUtilities(R): Exports == Implementation where R: IntegralDomain FR ==> Factored R Exports ==> with refine: (FR, R-> FR) -> FR ++ refine(u,fn) is used to apply the function \userfun{fn} to ++ each factor of \spadvar{u} and then build a new factored ++ object from the results. For example, if \spadvar{u} were ++ created by calling \spad{nilFactor(10,2)} then ++ \spad{refine(u,factor)} would create a factored object equal ++ to that created by \spad{factor(100)} or ++ \spad{primeFactor(2,2) * primeFactor(5,2)}. mergeFactors: (FR,FR) -> FR ++ mergeFactors(u,v) is used when the factorizations of \spadvar{u} ++ and \spadvar{v} are known to be disjoint, e.g. resulting from a ++ content/primitive part split. Essentially, it creates a new ++ factored object by multiplying the units together and appending ++ the lists of factors. Implementation ==> add fg: FR func: R -> FR fUnion ==> Union("nil", "sqfr", "irred", "prime") FF ==> Record(flg: fUnion, fctr: R, xpnt: Integer) mergeFactors(f,g) == makeFR(unit(f)*unit(g),append(factorList f,factorList g)) refine(f, func) == u := unit(f) l: List FF := empty() for item in factorList f repeat fitem := func item.fctr u := u*unit(fitem) ** (item.xpnt :: NonNegativeInteger) if item.xpnt = 1 then l := concat(factorList fitem,l) else l := concat([[v.flg,v.fctr,v.xpnt*item.xpnt] for v in factorList fitem],l) makeFR(u,l) @ \section{package FR2 FactoredFunctions2} <>= )abbrev package FR2 FactoredFunctions2 ++ Author: Robert S. Sutor ++ Date Created: 1987 ++ Change History: ++ Basic Operations: map ++ Related Constructors: Factored ++ Also See: ++ AMS Classifications: 11A51, 11Y05 ++ Keywords: map, factor ++ References: ++ Description: ++ \spadtype{FactoredFunctions2} contains functions that involve ++ factored objects whose underlying domains may not be the same. ++ For example, \spadfun{map} might be used to coerce an object of ++ type \spadtype{Factored(Integer)} to ++ \spadtype{Factored(Complex(Integer))}. FactoredFunctions2(R, S): Exports == Implementation where R: IntegralDomain S: IntegralDomain Exports ==> with map: (R -> S, Factored R) -> Factored S ++ map(fn,u) is used to apply the function \userfun{fn} to every ++ factor of \spadvar{u}. The new factored object will have all its ++ information flags set to "nil". This function is used, for ++ example, to coerce every factor base to another type. Implementation ==> add map(func, f) == func(unit f) * _*/[nilFactor(func(g.factor), g.exponent) for g in factors f] @ \section{License} <>= --Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd. --All rights reserved. --Copyright (C) 2007-2013, Gabriel Dos Reis. --All rights reserved. -- --Redistribution and use in source and binary forms, with or without --modification, are permitted provided that the following conditions are --met: -- -- - Redistributions of source code must retain the above copyright -- notice, this list of conditions and the following disclaimer. -- -- - Redistributions in binary form must reproduce the above copyright -- notice, this list of conditions and the following disclaimer in -- the documentation and/or other materials provided with the -- distribution. -- -- - Neither the name of The Numerical ALgorithms Group Ltd. nor the -- names of its contributors may be used to endorse or promote products -- derived from this software without specific prior written permission. -- --THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS --IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED --TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A --PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER --OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, --EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, --PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR --PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF --LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING --NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS --SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. @ <<*>>= <> <> <> <> @ \eject \begin{thebibliography}{99} \bibitem{1} nothing \end{thebibliography} \end{document}