padic.spad

)abbrev category PADICCT PAdicIntegerCategory
++ Author: Clifton J. Williamson
++ Date Created: 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)}.

)abbrev domain IPADIC InnerPAdicInteger
++ Author: Clifton J. Williamson
++ Date Created: 20 August 1989
++ 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?() == _$streamsShowAll$Lisp

    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
        st1 := st
        for n in (count + 1).. while explicitEntries? st and _
               not eq?(st1, rst st) repeat
          if frst(st) ~= 0 then
            l := concat(termOutput(n pretend I, frst st), l)
          if odd?(n) then st1 := rst(st1)
          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)

)abbrev domain PADIC PAdicInteger
++ Author: Clifton J. Williamson
++ Date Created: 20 August 1989
++ 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)

)abbrev domain BPADIC BalancedPAdicInteger
++ Author: Clifton J. Williamson
++ Date Created: 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)

)abbrev domain PADICRC PAdicRationalConstructor
++ Author: Clifton J. Williamson
++ Date Created: 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)

    retractIfCan(x : %) : Union(PADIC, "failed") ==
        n := getExpon(x)
        if n < 0 then
            x := removeZeroes(-n, x)
            n := getExpon(x)
        n < 0 => "failed"
        p^(qcoerce(n)@NNI)* getZp(x)

    numer(x : %) : PADIC ==
        n := getExpon(x)
        if n < 0 then
            x := removeZeroes(-n, x)
            n := getExpon(x)
        n <= 0 => getZp(x)
        p^(qcoerce(n)@NNI)* getZp(x)

    denom(x : %) : PADIC ==
        n := getExpon(x)
        if n < 0 then
            x := removeZeroes(-n, x)
            n := getExpon(x)
        n <= 0 => (p^(qcoerce(-n)@NNI))::PADIC
        1

    wholePart(x : %) : PADIC ==
        n := getExpon(x)
        if n < 0 then
            x := removeZeroes(-n, x)
            n := getExpon(x)
        n >= 0 => p^(qcoerce(n)@NNI)* getZp(x)
        ix := getZp(x)
        for i in 1..-n repeat
            ix := quotientByP(ix)
        ix

--% 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?() == _$streamsShowAll$Lisp

    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
        uu1 := uu
        for n in (count + 1).. while explicitEntries? uu and _
               not eq?(uu1, rst uu) repeat
          if frst(uu) ~= 0 then
            l := concat(termOutput((n::I) + m, frst(uu)), l)
          if odd?(n) then uu1 := rst(uu1)
          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)

)abbrev domain PADICRAT PAdicRational
++ Author: Clifton J. Williamson
++ Date Created: 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)

)abbrev domain BPADICRT BalancedPAdicRational
++ Author: Clifton J. Williamson
++ Date Created: 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)

--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.