axiom-source 20170501-3 / usr / share / axiom-20170501 / src / algebra / FR.spad

)abbrev domain FR Factored
++ Author: Robert S. Sutor, J. Grabmeier
++ Date Created: 1985
++ Change History: 21 Jan 1991 (J Grabmeier) 16 Aug 1994 (R S Sutor)
++ 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) : SIG == CODE where
  R : IntegralDomain

  fUnion ==> Union("nil", "sqfr", "irred", "prime")
  FF     ==> Record(flg: fUnion, fctr: R, xpnt: Integer)
  SRFE   ==> Set(Record(factor:R, exponent:Integer))

  SIG ==> Join(IntegralDomain, DifferentialExtension R, Algebra R,
                   FullyEvalableOver R, FullyRetractableTo 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.
    ++
    ++X f:=nilFactor(y-x,3)
    ++X expand(f)

   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.
    ++
    ++X f:=nilFactor(y-x,3)
    ++X exponent(f)

   makeFR : (R, List FF) -> %
    ++ makeFR(unit,listOfFactors) creates a factored object (for
    ++ use by factoring code).
    ++
    ++X f:=nilFactor(x-y,3)
    ++X g:=factorList f
    ++X makeFR(z,g)

   factorList : % -> List FF
    ++ factorList(u) returns the list of factors with flags (for
    ++ use by factoring code).
    ++
    ++X f:=nilFactor(x-y,3)
    ++X factorList f

   nilFactor : (R, Integer) -> %
    ++ nilFactor(base,exponent) creates a factored object with
    ++ a single factor with no information about the kind of
    ++ base (flag = "nil").
    ++
    ++X nilFactor(24,2)
    ++X nilFactor(x-y,3)

   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)}).
    ++
    ++X f:=x*y^3-3*x^2*y^2+3*x^3*y-x^4
    ++X factors f
    ++X g:=makeFR(z,factorList f)
    ++X factors g

   irreducibleFactor : (R, Integer) -> %
    ++ irreducibleFactor(base,exponent) creates a factored object with
    ++ a single factor whose base is asserted to be irreducible
    ++ (flag = "irred").
    ++
    ++X a:=irreducibleFactor(3,1)
    ++X nthFlag(a,1)

   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.
    ++
    ++X a:=factor 9720000
    ++X nthExponent(a,2)

   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.
    ++
    ++X a:=factor 9720000
    ++X nthFactor(a,2)

   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.
    ++
    ++X a:=factor 9720000
    ++X nthFlag(a,2)

   numberOfFactors : %  -> NonNegativeInteger
    ++ numberOfFactors(u) returns the number of factors in \spadvar{u}.
    ++
    ++X a:=factor 9720000
    ++X numberOfFactors a

   primeFactor : (R,Integer) -> %
    ++ primeFactor(base,exponent) creates a factored object with
    ++ a single factor whose base is asserted to be prime
    ++ (flag = "prime").
    ++
    ++X a:=primeFactor(3,4)
    ++X nthFlag(a,1)

   sqfrFactor : (R,Integer) -> %
    ++ sqfrFactor(base,exponent) creates a factored object with
    ++ a single factor whose base is asserted to be square-free
    ++ (flag = "sqfr").
    ++
    ++X a:=sqfrFactor(3,5)
    ++X nthFlag(a,1)

   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.
    ++
    ++X f:=x*y^3-3*x^2*y^2+3*x^3*y-x^4
    ++X unit f
    ++X g:=makeFR(z,factorList f)
    ++X unit g

   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.

   map : (R -> R, %) -> %
    ++ map(fn,u) maps the function \userfun{fn} across the factors of
    ++ \spadvar{u} and creates a new factored object. Note: this clears
    ++ the information flags (sets them to "nil") because the effect of
    ++ \userfun{fn} is clearly not known in general.
    ++
    ++X m(a:Factored Polynomial Integer):Factored Polynomial Integer == a^2
    ++X f:=x*y^3-3*x^2*y^2+3*x^3*y-x^4
    ++X map(m,f)
    ++X g:=makeFR(z,factorList f)
    ++X map(m,g)

    -- 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, %)

  CODE ==> 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) => l
          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 )
        empty? l => convert(unit x)@InputForm
        if unit x ^= 1 then l := concat(convert(unit x)@InputForm,l)
        empty? rest l => first l
        binary(convert(_*::Symbol)@InputForm, l)@InputForm

    orderedR? := R has OrderedSet

  -- 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
                               (u.unt = 1)

    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) ==
        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

      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
      (u = 1) => v
      (v = 1) => 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)
      if orderedR? then 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) => l
        ((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) ==
      orderedR? => y.fctr < y1.fctr
      GGREATERP(y.fctr, y1.fctr)$Lisp => false
      true

    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):% ==
              (r = 1) => 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])

    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 x1.first.xpnt < 0
          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 x.xpnt < 0
          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 ((ucar.canonical) = 1) 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) ==
        (u = 1) or (v = 1) => 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
        if orderedR? then 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
                       ((l.first.xpnt) = 1) and (l.first.flg case "prime")