path: root/src/algebra/padic.spad.pamphlet

\documentclass{article}
\usepackage{axiom}
\begin{document}
\title{\$SPAD/src/algebra padic.spad}
\author{Clifton J. Williamson}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\section{category PADICCT PAdicIntegerCategory}
<<category PADICCT PAdicIntegerCategory>>=
)abbrev category PADICCT PAdicIntegerCategory
++ Author: Clifton J. Williamson
++ Date Created: 15 May 1990
++ Date Last Updated: 15 May 1990
++ Basic Operations:
++ Related Domains:
++ Also See:
++ AMS Classifications:
++ Keywords: p-adic, completion
++ Examples:
++ References:
++ Description: This is the catefory of stream-based representations of
++   the p-adic integers.
PAdicIntegerCategory(p): Category == Definition where
  p   :   Integer
  I   ==> Integer
  NNI ==> NonNegativeInteger
  ST  ==> Stream
  SUP ==> SparseUnivariatePolynomial

  Definition ==> Join(EuclideanDomain,CharacteristicZero) with
    digits: % -> ST I
      ++ \spad{digits(x)} returns a stream of p-adic digits of x.
    order: % -> NNI
      ++ \spad{order(x)} returns the exponent of the highest power of p
      ++ dividing x.
    extend: (%,I) -> %
      ++ \spad{extend(x,n)} forces the computation of digits up to order n.
    complete: % -> %
      ++ \spad{complete(x)} forces the computation of all digits.
    modulus: () -> I
      ++ \spad{modulus()} returns the value of p.
    moduloP: % -> I
      ++ \spad{modulo(x)} returns a, where \spad{x = a + b p}.
    quotientByP: % -> %
      ++ \spad{quotientByP(x)} returns b, where \spad{x = a + b p}.
    approximate: (%,I) -> I
      ++ \spad{approximate(x,n)} returns an integer y such that
      ++ \spad{y = x (mod p^n)}
      ++ when n is positive, and 0 otherwise.
    sqrt: (%,I) -> %
      ++ \spad{sqrt(b,a)} returns a square root of b.
      ++ Argument \spad{a} is a square root of b \spad{(mod p)}.
    root: (SUP I,I) -> %
      ++ \spad{root(f,a)} returns a root of the polynomial \spad{f}.
      ++ Argument \spad{a} must be a root of \spad{f} \spad{(mod p)}.

@
\section{domain IPADIC InnerPAdicInteger}
<<domain IPADIC InnerPAdicInteger>>=
)abbrev domain IPADIC InnerPAdicInteger
++ Author: Clifton J. Williamson
++ Date Created: 20 August 1989
++ Date Last Updated: 15 May 1990
++ Basic Operations:
++ Related Domains:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ Keywords: p-adic, completion
++ Examples:
++ References:
++ Description:
++   This domain implements Zp, the p-adic completion of the integers.
++   This is an internal domain.
InnerPAdicInteger(p,unBalanced?): Exports == Implementation where
  p           : Integer
  unBalanced? : Boolean
  I   ==> Integer
  NNI ==> NonNegativeInteger
  OUT ==> OutputForm
  L   ==> List
  ST  ==> Stream
  SUP ==> SparseUnivariatePolynomial

  Exports ==> PAdicIntegerCategory p

  Implementation ==> add

    PEXPR := p :: OUT

    Rep := ST I

    characteristic == 0
    euclideanSize(x) == order(x)

    stream(x:%):ST I == x pretend ST(I)
    padic(x:ST I):% == x pretend %
    digits x == stream x

    extend(x,n) == extend(x,n + 1)$Rep
    complete x == complete(x)$Rep

--     notBalanced?:() -> Boolean
--     notBalanced?() == unBalanced?

    modP:I -> I
    modP n ==
      unBalanced? or (p = 2) => positiveRemainder(n,p)
      symmetricRemainder(n,p)

    modPInfo:I -> Record(digit:I,carry:I)
    modPInfo n ==
      dv := divide(n,p)
      r0 := dv.remainder; q := dv.quotient
      if (r := modP r0) ~= r0 then q := q + ((r0 - r) quo p)
      [r,q]

    invModP: I -> I
    invModP n == invmod(n,p)

    modulus()     == p
    moduloP x     == (empty? x => 0; frst x)
    quotientByP x == (empty? x => x; rst x)

    approximate(x,n) ==
      n <= 0 or empty? x => 0
      frst x + p * approximate(rst x,n - 1)

    x = y ==
      st : ST I := stream(x - y)
      n : I := _$streamCount$Lisp
      for i in 0..n repeat
        empty? st => return true
        frst st ~= 0 => return false
        st := rst st
      empty? st

    order x ==
      st := stream x
      for i in 0..1000 repeat
        empty? st => return 0
        frst st ~= 0 => return i
        st := rst st
      error "order: series has more than 1000 leading zero coefs"

    0 == padic concat(0$I,empty())
    1 == padic concat(1$I,empty())

    intToPAdic: I -> ST I
    intToPAdic n == delay
      n = 0 => empty()
      modp := modPInfo n
      concat(modp.digit,intToPAdic modp.carry)

    intPlusPAdic: (I,ST I) -> ST I
    intPlusPAdic(n,x) == delay
      empty? x => intToPAdic n
      modp := modPInfo(n + frst x)
      concat(modp.digit,intPlusPAdic(modp.carry,rst x))

    intMinusPAdic: (I,ST I) -> ST I
    intMinusPAdic(n,x) == delay
      empty? x => intToPAdic n
      modp := modPInfo(n - frst x)
      concat(modp.digit,intMinusPAdic(modp.carry,rst x))

    plusAux: (I,ST I,ST I) -> ST I
    plusAux(n,x,y) == delay
      empty? x and empty? y => intToPAdic n
      empty? x => intPlusPAdic(n,y)
      empty? y => intPlusPAdic(n,x)
      modp := modPInfo(n + frst x + frst y)
      concat(modp.digit,plusAux(modp.carry,rst x,rst y))

    minusAux: (I,ST I,ST I) -> ST I
    minusAux(n,x,y) == delay
      empty? x and empty? y => intToPAdic n
      empty? x => intMinusPAdic(n,y)
      empty? y => intPlusPAdic(n,x)
      modp := modPInfo(n + frst x - frst y)
      concat(modp.digit,minusAux(modp.carry,rst x,rst y))

    x + y == padic plusAux(0,stream x,stream y)
    x - y == padic minusAux(0,stream x,stream y)
    - y   == padic intMinusPAdic(0,stream y)
    coerce(n:I) == padic intToPAdic n

    intMult:(I,ST I) -> ST I
    intMult(n,x) == delay
      empty? x => empty()
      modp := modPInfo(n * frst x)
      concat(modp.digit,intPlusPAdic(modp.carry,intMult(n,rst x)))

    (n:I) * (x:%) ==
      padic intMult(n,stream x)

    timesAux:(ST I,ST I) -> ST I
    timesAux(x,y) == delay
      empty? x or empty? y => empty()
      modp := modPInfo(frst x * frst y)
      car := modp.digit
      cdr : ST I --!!
      cdr := plusAux(modp.carry,intMult(frst x,rst y),timesAux(rst x,y))
      concat(car,cdr)

    (x:%) * (y:%) == padic timesAux(stream x,stream y)

    quotientAux:(ST I,ST I) -> ST I
    quotientAux(x,y) == delay
      empty? x => error "quotientAux: first argument"
      empty? y => empty()
      modP frst x = 0 =>
        modP frst y = 0 => quotientAux(rst x,rst y)
        error "quotient: quotient not integral"
      z0 := modP(invModP frst x * frst y)
      yy : ST I --!!
      yy := rest minusAux(0,y,intMult(z0,x))
      concat(z0,quotientAux(x,yy))

    recip x ==
      empty? x or modP frst x = 0 => "failed"
      padic quotientAux(stream x,concat(1,empty()))

    iExquo: (%,%,I) -> Union(%,"failed")
    iExquo(xx,yy,n) ==
      n > 1000 =>
        error "exquo: quotient by series with many leading zero coefs"
      empty? yy => "failed"
      empty? xx => 0
      zero? frst yy =>
        zero? frst xx => iExquo(rst xx,rst yy,n + 1)
        "failed"
      (rec := recip yy) case "failed" => "failed"
      xx * (rec :: %)

    x exquo y == iExquo(stream x,stream y,0)

    divide(x,y) ==
      (z:=x exquo y) case "failed" => [0,x]
      [z, 0]

    iSqrt: (I,I,I,%) -> %
    iSqrt(pn,an,bn,bSt) == delay
      bn1 := (empty? bSt => bn; bn + pn * frst(bSt))
      c := (bn1 - an*an) quo pn
      aa := modP(c * invmod(2*an,p))
      nSt := (empty? bSt => bSt; rst bSt)
      concat(aa,iSqrt(pn*p,an + pn*aa,bn1,nSt))

    sqrt(b,a) ==
      p = 2 =>
        error "sqrt: no square roots in Z2 yet"
      not zero? modP(a*a - (bb := moduloP b)) =>
        error "sqrt: not a square root (mod p)"
      b := (empty? b => b; rst b)
      a := modP a
      concat(a,iSqrt(p,a,bb,b))

    iRoot: (SUP I,I,I,I) -> ST I
    iRoot(f,alpha,invFpx0,pPow) == delay
      num := -((f(alpha) exquo pPow) :: I)
      digit := modP(num * invFpx0)
      concat(digit,iRoot(f,alpha + digit * pPow,invFpx0,p * pPow))

    root(f,x0) ==
      x0 := modP x0
      not zero? modP f(x0) =>
        error "root: not a root (mod p)"
      fpx0 := modP (differentiate f)(x0)
      zero? fpx0 =>
        error "root: approximate root must be a simple root (mod p)"
      invFpx0 := modP invModP fpx0
      padic concat(x0,iRoot(f,x0,invFpx0,p))

    termOutput:(I,I) -> OUT
    termOutput(k,c) ==
      k = 0 => c :: OUT
      mon := (k = 1 => PEXPR; PEXPR ** (k :: OUT))
      c = 1 => mon
      c = -1 => -mon
      (c :: OUT) * mon

    showAll?:() -> Boolean
    -- check a global Lisp variable
    showAll?() == true

    coerce(x:%):OUT ==
      empty?(st := stream x) => 0 :: OUT
      n : NNI ; count : NNI := _$streamCount$Lisp
      l : L OUT := empty()
      for n in 0..count while not empty? st repeat
        if frst(st) ~= 0 then
          l := concat(termOutput(n :: I,frst st),l)
        st := rst st
      if showAll?() then
        for n in (count + 1).. while explicitEntries? st and _
               not eq?(st,rst st) repeat
          if frst(st) ~= 0 then
            l := concat(termOutput(n pretend I,frst st),l)
          st := rst st
      l :=
        explicitlyEmpty? st => l
        eq?(st,rst st) and frst st = 0 => l
        concat(prefix("O" :: OUT,[PEXPR ** (n :: OUT)]),l)
      empty? l => 0 :: OUT
      reduce("+",reverse_! l)

@
\section{domain PADIC PAdicInteger}
<<domain PADIC PAdicInteger>>=
)abbrev domain PADIC PAdicInteger
++ Author: Clifton J. Williamson
++ Date Created: 20 August 1989
++ Date Last Updated: 15 May 1990
++ Basic Operations:
++ Related Domains:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ Keywords: p-adic, completion
++ Examples:
++ References:
++ Description:
++   Stream-based implementation of Zp: p-adic numbers are represented as
++   sum(i = 0.., a[i] * p^i), where the a[i] lie in 0,1,...,(p - 1).
PAdicInteger(p:Integer) == InnerPAdicInteger(p,true$Boolean)

@
\section{domain BPADIC BalancedPAdicInteger}
<<domain BPADIC BalancedPAdicInteger>>=
)abbrev domain BPADIC BalancedPAdicInteger
++ Author: Clifton J. Williamson
++ Date Created: 15 May 1990
++ Date Last Updated: 15 May 1990
++ Basic Operations:
++ Related Domains:
++ Also See:
++ AMS Classifications:
++ Keywords: p-adic, complementation
++ Examples:
++ References:
++ Description:
++   Stream-based implementation of Zp: p-adic numbers are represented as
++   sum(i = 0.., a[i] * p^i), where the a[i] lie in -(p - 1)/2,...,(p - 1)/2.
BalancedPAdicInteger(p:Integer) == InnerPAdicInteger(p,false$Boolean)

@
\section{domain PADICRC PAdicRationalConstructor}
<<domain PADICRC PAdicRationalConstructor>>=
)abbrev domain PADICRC PAdicRationalConstructor
++ Author: Clifton J. Williamson
++ Date Created: 10 May 1990
++ Date Last Updated: 10 May 1990
++ Basic Operations:
++ Related Domains:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ Keywords: p-adic, completion
++ Examples:
++ References:
++ Description: This is the category of stream-based representations of Qp.
PAdicRationalConstructor(p,PADIC): Exports == Implementation where
  p     :   Integer
  PADIC :   PAdicIntegerCategory p
  CF    ==> ContinuedFraction
  I     ==> Integer
  NNI   ==> NonNegativeInteger
  OUT   ==> OutputForm
  L     ==> List
  RN    ==> Fraction Integer
  ST    ==> Stream

  Exports ==> QuotientFieldCategory(PADIC) with
    approximate: (%,I) -> RN
      ++ \spad{approximate(x,n)} returns a rational number y such that
      ++ \spad{y = x (mod p^n)}.
    continuedFraction: % -> CF RN
      ++ \spad{continuedFraction(x)} converts the p-adic rational number x
      ++ to a continued fraction.
    removeZeroes: % -> %
      ++ \spad{removeZeroes(x)} removes leading zeroes from the
      ++ representation of the p-adic rational \spad{x}.
      ++ A p-adic rational is represented by (1) an exponent and
      ++ (2) a p-adic integer which may have leading zero digits.
      ++ When the p-adic integer has a leading zero digit, a 'leading zero'
      ++ is removed from the p-adic rational as follows:
      ++ the number is rewritten by increasing the exponent by 1 and
      ++ dividing the p-adic integer by p.
      ++ Note: \spad{removeZeroes(f)} removes all leading zeroes from f.
    removeZeroes: (I,%) -> %
      ++ \spad{removeZeroes(n,x)} removes up to n leading zeroes from
      ++ the p-adic rational \spad{x}.

  Implementation ==> add

    PEXPR := p :: OUT

--% representation

    Rep := Record(expon:I,pint:PADIC)

    getExpon: % -> I
    getZp   : % -> PADIC
    makeQp  : (I,PADIC) -> %

    getExpon x    == x.expon
    getZp x       == x.pint
    makeQp(r,int) == [r,int]

--% creation

    0 == makeQp(0,0)
    1 == makeQp(0,1)

    coerce(x:I)     == x :: PADIC :: %
    coerce(r:RN)    == (numer(r) :: %)/(denom(r) :: %)
    coerce(x:PADIC) == makeQp(0,x)

--% normalizations

    removeZeroes x ==
      empty? digits(xx := getZp x) => 0
      zero? moduloP xx =>
        removeZeroes makeQp(getExpon x + 1,quotientByP xx)
      x

    removeZeroes(n,x) ==
      n <= 0 => x
      empty? digits(xx := getZp x) => 0
      zero? moduloP xx =>
        removeZeroes(n - 1,makeQp(getExpon x + 1,quotientByP xx))
      x

--% arithmetic

    x = y ==
      EQ(x,y)$Lisp => true
      n := getExpon(x) - getExpon(y)
      n >= 0 =>
        (p**(n :: NNI) * getZp(x)) = getZp(y)
      (p**((- n) :: NNI) * getZp(y)) = getZp(x)

    x + y ==
      n := getExpon(x) - getExpon(y)
      n >= 0 =>
        makeQp(getExpon y,getZp(y) + p**(n :: NNI) * getZp(x))
      makeQp(getExpon x,getZp(x) + p**((-n) :: NNI) * getZp(y))

    -x == makeQp(getExpon x,-getZp(x))

    x - y ==
      n := getExpon(x) - getExpon(y)
      n >= 0 =>
        makeQp(getExpon y,p**(n :: NNI) * getZp(x) - getZp(y))
      makeQp(getExpon x,getZp(x) - p**((-n) :: NNI) * getZp(y))

    n:I * x:% == makeQp(getExpon x,n * getZp x)
    x:% * y:% == makeQp(getExpon x + getExpon y,getZp x * getZp y)

    x:% ** n:I ==
      zero? n => 1
      positive? n => expt(x,n :: PositiveInteger)$RepeatedSquaring(%)
      inv expt(x,(-n) :: PositiveInteger)$RepeatedSquaring(%)

    recip x ==
      x := removeZeroes(1000,x)
      zero? moduloP(xx := getZp x) => "failed"
      (inv := recip xx) case "failed" => "failed"
      makeQp(- getExpon x,inv :: PADIC)

    inv x ==
      (inv := recip x) case "failed" => error "inv: no inverse"
      inv :: %

    x:% / y:% == x * inv y
    x:PADIC / y:PADIC == (x :: %) / (y :: %)
    x:PADIC * y:% == makeQp(getExpon y,x * getZp y)

    approximate(x,n) ==
      k := getExpon x
      (p :: RN) ** k * approximate(getZp x,n - k)

    cfStream: % -> Stream RN
    cfStream x == delay
--    zero? x => empty()
      invx := inv x; x0 := approximate(invx,1)
      concat(x0,cfStream(invx - (x0 :: %)))

    continuedFraction x ==
      x0 := approximate(x,1)
      reducedContinuedFraction(x0,cfStream(x - (x0 :: %)))

    termOutput:(I,I) -> OUT
    termOutput(k,c) ==
      k = 0 => c :: OUT
      mon := (k = 1 => PEXPR; PEXPR ** (k :: OUT))
      c = 1 => mon
      c = -1 => -mon
      (c :: OUT) * mon

    showAll?:() -> Boolean
    -- check a global Lisp variable
    showAll?() == true

    coerce(x:%):OUT ==
      x := removeZeroes(_$streamCount$Lisp,x)
      m := getExpon x; zp := getZp x
      uu := digits zp
      l : L OUT := empty()
      empty? uu => 0 :: OUT
      n : NNI ; count : NNI := _$streamCount$Lisp
      for n in 0..count while not empty? uu repeat
        if frst(uu) ~= 0 then
          l := concat(termOutput((n :: I) + m,frst(uu)),l)
        uu := rst uu
      if showAll?() then
        for n in (count + 1).. while explicitEntries? uu and _
               not eq?(uu,rst uu) repeat
          if frst(uu) ~= 0 then
            l := concat(termOutput((n::I) + m,frst(uu)),l)
          uu := rst uu
      l :=
        explicitlyEmpty? uu => l
        eq?(uu,rst uu) and frst uu = 0 => l
        concat(prefix("O" :: OUT,[PEXPR ** ((n :: I) + m) :: OUT]),l)
      empty? l => 0 :: OUT
      reduce("+",reverse_! l)

@
\section{domain PADICRAT PAdicRational}
<<domain PADICRAT PAdicRational>>=
)abbrev domain PADICRAT PAdicRational
++ Author: Clifton J. Williamson
++ Date Created: 15 May 1990
++ Date Last Updated: 15 May 1990
++ Keywords: p-adic, complementation
++ Basic Operations:
++ Related Domains:
++ Also See:
++ AMS Classifications:
++ Keywords: p-adic, completion
++ Examples:
++ References:
++ Description:
++   Stream-based implementation of Qp: numbers are represented as
++   sum(i = k.., a[i] * p^i) where the a[i] lie in 0,1,...,(p - 1).
PAdicRational(p:Integer) == PAdicRationalConstructor(p,PAdicInteger p)

@
\section{domain BPADICRT BalancedPAdicRational}
<<domain BPADICRT BalancedPAdicRational>>=
)abbrev domain BPADICRT BalancedPAdicRational
++ Author: Clifton J. Williamson
++ Date Created: 15 May 1990
++ Date Last Updated: 15 May 1990
++ Keywords: p-adic, complementation
++ Basic Operations:
++ Related Domains:
++ Also See:
++ AMS Classifications:
++ Keywords: p-adic, completion
++ Examples:
++ References:
++ Description:
++   Stream-based implementation of Qp: numbers are represented as
++   sum(i = k.., a[i] * p^i), where the a[i] lie in -(p - 1)/2,...,(p - 1)/2.
BalancedPAdicRational(p:Integer) ==
  PAdicRationalConstructor(p,BalancedPAdicInteger p)

@
\section{License}
<<license>>=
--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
--All rights reserved.
--
--Redistribution and use in source and binary forms, with or without
--modification, are permitted provided that the following conditions are
--met:
--
--    - Redistributions of source code must retain the above copyright
--      notice, this list of conditions and the following disclaimer.
--
--    - Redistributions in binary form must reproduce the above copyright
--      notice, this list of conditions and the following disclaimer in
--      the documentation and/or other materials provided with the
--      distribution.
--
--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
--      names of its contributors may be used to endorse or promote products
--      derived from this software without specific prior written permission.
--
--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
@
<<*>>=
<<license>>

<<category PADICCT PAdicIntegerCategory>>
<<domain IPADIC InnerPAdicInteger>>
<<domain PADIC PAdicInteger>>
<<domain BPADIC BalancedPAdicInteger>>
<<domain PADICRC PAdicRationalConstructor>>
<<domain PADICRAT PAdicRational>>
<<domain BPADICRT BalancedPAdicRational>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}