radix.spad

)abbrev domain RADIX RadixExpansion
++ Author: Stephen M. Watt
++ Date Created: October 1986
++ Basic Operations: wholeRadix, fractRadix, wholeRagits, fractRagits
++ Related Domains: BinaryExpansion, DecimalExpansion, HexadecimalExpansion,
++    RadixUtilities
++ Also See:
++ AMS Classifications:
++ Keywords: radix, base, repeating decimal
++ Examples:
++ References:
++ Description:
++   This domain allows rational numbers to be presented as repeating
++   decimal expansions or more generally as repeating expansions in any base.

RadixExpansion(bb) : Exports == Implementation where
  bb   :  Integer
  I   ==> Integer
  NNI ==> NonNegativeInteger
  OUT ==> OutputForm
  RN  ==> Fraction Integer
  ST  ==> Stream Integer
  QuoRem ==> Record(quotient : Integer, remainder : Integer)

  Exports ==> QuotientFieldCategory(Integer) with
    coerce : % -> Fraction Integer
      ++ coerce(rx) converts a radix expansion to a rational number.
    fractionPart : % -> Fraction Integer
      ++ fractionPart(rx) returns the fractional part of a radix expansion.
    wholeRagits : % -> List Integer
      ++ wholeRagits(rx) returns the ragits of the integer part
      ++ of a radix expansion.
    fractRagits : % -> Stream Integer
      ++ fractRagits(rx) returns the ragits of the fractional part
      ++ of a radix expansion.
    prefixRagits : % -> List Integer
      ++ prefixRagits(rx) returns the non-cyclic part of the ragits
      ++ of the fractional part of a radix expansion.
      ++ For example, if \spad{x = 3/28 = 0.10 714285 714285 ...},
      ++ then \spad{prefixRagits(x)=[1, 0]}.
    cycleRagits : % -> List Integer
      ++ cycleRagits(rx) returns the cyclic part of the ragits of the
      ++ fractional part of a radix expansion.
      ++ For example, if \spad{x = 3/28 = 0.10 714285 714285 ...},
      ++ then \spad{cycleRagits(x) = [7, 1, 4, 2, 8, 5]}.
    wholeRadix : List Integer -> %
      ++ wholeRadix(l) creates an integral radix expansion from a list
      ++ of ragits.
      ++ For example, \spad{wholeRadix([1, 3, 4])} will return \spad{134}.
    fractRadix : (List Integer, List Integer) -> %
      ++ fractRadix(pre, cyc) creates a fractional radix expansion
      ++ from a list of prefix ragits and a list of cyclic ragits.
      ++ For example, \spad{fractRadix([1], [6])} will return \spad{0.16666666...}.

  Implementation ==> add
    -- The efficiency of arithmetic operations is poor.
    -- Could use a lazy eval where either rational rep
    -- or list of ragit rep (the current) or both are kept
    -- as demanded.

    bb < 2 => error "Radix base must be at least 2"
    Rep := Record(sgn : Integer,      int : List Integer,
                  pfx : List Integer, cyc : List Integer)

    q :     RN
    qr :    QuoRem
    a, b :   %
    n :     I

    radixInt :    (I, I)    -> List I
    radixFrac :   (I, I, I) -> Record(pfx : List I, cyc : List I)
    checkRagits : List I    -> Boolean

    -- Arithmetic operations
    characteristic() == 0
    differentiate a == 0

    0     == [1, [],  [], []]
    1     == [1, [1], [], []]
    - a   == (a = 0 => 0; [-a.sgn, a.int, a.pfx, a.cyc])
    a + b == (a::RN + b::RN)::%
    a - b == (a::RN - b::RN)@RN::%
    n * a == (n     * a::RN)::%
    a * b == (a::RN * b::RN)::%
    a / b == (a::RN / b::RN)::%
    (i : I) / (j : I) == (i/j)@RN :: %
    a < b == a::RN < b::RN
    a = b == a.sgn = b.sgn and a.int = b.int and
             a.pfx = b.pfx and a.cyc = b.cyc
    numer a == numer(a::RN)
    denom a == denom(a::RN)

    -- Algebraic coercions
    coerce(a) : RN == (wholePart a) :: RN + fractionPart a
    coerce(n) : %  == n :: RN :: %
    coerce(q) : %  ==
      s := 1; if q < 0 then (s := -1; q := -q)
      qr      := divide(numer q, denom q)
      whole   := radixInt (qr.quotient, bb)
      fractn  := radixFrac(qr.remainder, denom q, bb)
      cycle   := (fractn.cyc = [0] => []; fractn.cyc)
      [s, whole, fractn.pfx, cycle]

    retractIfCan(a):Union(RN,"failed") == a::RN
    retractIfCan(a):Union(I,"failed") ==
      empty?(a.pfx) and empty?(a.cyc) => wholePart a
      "failed"

    -- Exported constructor/destructors
    ceiling a == ceiling(a::RN)
    floor a == floor(a::RN)

    wholePart a ==
      n0 := 0
      for r in a.int repeat n0 := bb*n0 + r
      a.sgn*n0

    fractionPart(a : %) : RN ==
      n0 := 0
      for r in a.pfx repeat n0 := bb*n0 + r
      empty?(a.cyc) =>
          a.sgn*n0/bb^((#a.pfx)::NNI)
      n1 := n0
      for r in a.cyc repeat n1 := bb*n1 + r
      n := n1 - n0
      d := (bb^((#a.cyc)::NNI) - 1) * bb^((#a.pfx)::NNI)
      a.sgn*n/d

    wholeRagits  a == a.int
    fractRagits  a == concat(construct(a.pfx)@ST, repeating a.cyc)
    prefixRagits a == a.pfx
    cycleRagits  a == a.cyc

    wholeRadix li ==
      checkRagits li
      [1, li, [], []]
    fractRadix(lpfx, lcyc) ==
      checkRagits lpfx; checkRagits lcyc
      [1, [], lpfx, lcyc]

    -- Output

    ALPHAS : String := "ABCDEFGHIJKLMNOPQRSTUVWXYZ"

    intToExpr(i : I) : OUT ==
      -- computes a digit for bases between 11 and 36
      i < 10 => i :: OUT
      elt(ALPHAS, (i-10) + minIndex(ALPHAS)) :: OUT

    exprgroup(le : List OUT) : OUT ==
      empty? le      => error "exprgroup needs non-null list"
      empty? rest le => first le
      abs bb <= 36 => hconcat le
      blankSeparate le

    intgroup(li : List I) : OUT ==
      empty? li      => error "intgroup needs non-null list"
      empty? rest li =>
         abs bb <= 36 => intToExpr first(li)
         first(li) :: OUT
      abs bb <= 10 => hconcat [i :: OUT for i in li]
      abs bb <= 36 => hconcat [intToExpr(i) for i in li]
      blankSeparate [i :: OUT for i in li]

    overBar(li : List I) : OUT == overbar intgroup li

    coerce(a) : OUT ==
      le : List OUT := []
      if not(empty?(a.cyc)) then le := concat(overBar(a.cyc), le)
      if not(empty?(a.pfx)) then le := concat(intgroup(a.pfx), le)
      if not(empty?(le))    then le := concat(message("."), le)
      if not(empty?(a.int)) then le := concat(intgroup(a.int), le)
      else le := concat(0 :: OUT, le)
      rex := exprgroup le
      if a.sgn < 0 then -rex else rex

    -- Construction utilities
    checkRagits li ==
      for i in li repeat if i < 0 or i >= bb then
        error "Each ragit (digit) must be between 0 and base-1"
      true

    radixInt(n, bas) ==
      rits : List I := []
      while abs n ~= 0 repeat
        qr   := divide(n, bas)
        n    := qr.quotient
        rits := concat(qr.remainder, rits)
      rits

    radixFrac(num, den, bas) ==
      -- Rits is the sequence of quotient/remainder pairs
      -- in calculating the radix expansion of the rational number.
      -- We wish to find p and c such that
      --    rits.i are distinct    for 0<=i<=p+c-1
      --    rits.i = rits.(i+p)    for i>p
      -- I.e. p is the length of the non-periodic prefix and c is
      -- the length of the cycle.

      -- Compute p and c using Floyd's algorithm.
      -- 1. Find smallest n s.t. rits.n = rits.(2*n)
      qr    := divide(bas * num, den)
      i : I := 0
      qr1i  := qr2i := qr
      rits : List QuoRem := [qr]
      until qr1i = qr2i repeat
        qr1i := divide(bas * qr1i.remainder, den)
        qrt  := divide(bas * qr2i.remainder, den)
        qr2i := divide(bas * qrt.remainder, den)
        rits := concat(qr2i, concat(qrt, rits))
        i    := i + 1
      rits := reverse! rits
      n    := i
      -- 2. Find p = first i such that rits.i = rits.(i+n)
      ritsi := rits
      ritsn := rits; for i in 1..n repeat ritsn := rest ritsn
      i := 0
      while first(ritsi) ~= first(ritsn) repeat
        ritsi := rest ritsi
        ritsn := rest ritsn
        i     := i + 1
      p := i
      -- 3. Find c = first i such that rits.p = rits.(p+i)
      ritsn := rits; for i in 1..n repeat ritsn := rest ritsn
      rn    := first ritsn
      cfound:= false
      c : I := 0
      for i in 1..p while not cfound repeat
        ritsn := rest ritsn
        if rn = first(ritsn) then
          c := i
          cfound := true
      if not cfound then c := n
      -- 4. Now produce the lists of ragits.
      ritspfx : List I := []
      ritscyc : List I := []
      for i in 1..p repeat
        ritspfx := concat(first(rits).quotient, ritspfx)
        rits    := rest rits
      for i in 1..c repeat
        ritscyc := concat(first(rits).quotient, ritscyc)
        rits    := rest rits
      [reverse! ritspfx, reverse! ritscyc]

)abbrev domain BINARY BinaryExpansion
++ Author: Clifton J. Williamson
++ Date Created: April 26, 1990
++ Basic Operations:
++ Related Domains: RadixExpansion
++ Also See:
++ AMS Classifications:
++ Keywords: radix, base, binary
++ Examples:
++ References:
++ Description:
++   This domain allows rational numbers to be presented as repeating
++   binary expansions.

BinaryExpansion() : Exports == Implementation where
  Exports ==> QuotientFieldCategory(Integer) with
    coerce : % -> Fraction Integer
      ++ coerce(b) converts a binary expansion to a rational number.
    coerce : % -> RadixExpansion(2)
      ++ coerce(b) converts a binary expansion to a radix expansion with base 2.
    fractionPart : % -> Fraction Integer
      ++ fractionPart(b) returns the fractional part of a binary expansion.
    binary : Fraction Integer -> %
      ++ binary(r) converts a rational number to a binary expansion.

  Implementation ==> RadixExpansion(2) add
    binary r == r :: %
    coerce(x : %) : RadixExpansion(2) == x pretend RadixExpansion(2)

)abbrev domain DECIMAL DecimalExpansion
++ Author: Stephen M. Watt
++ Date Created: October, 1986
++ Basic Operations:
++ Related Domains: RadixExpansion
++ Also See:
++ AMS Classifications:
++ Keywords: radix, base, repeating decimal
++ Examples:
++ References:
++ Description:
++   This domain allows rational numbers to be presented as repeating
++   decimal expansions.
DecimalExpansion() : Exports == Implementation where
  Exports ==> QuotientFieldCategory(Integer) with
    coerce : % -> Fraction Integer
      ++ coerce(d) converts a decimal expansion to a rational number.
    coerce : % -> RadixExpansion(10)
      ++ coerce(d) converts a decimal expansion to a radix expansion
      ++ with base 10.
    fractionPart : % -> Fraction Integer
      ++ fractionPart(d) returns the fractional part of a decimal expansion.
    decimal : Fraction Integer -> %
      ++ decimal(r) converts a rational number to a decimal expansion.

  Implementation ==> RadixExpansion(10) add
    decimal r == r :: %
    coerce(x : %) : RadixExpansion(10) == x pretend RadixExpansion(10)

)abbrev domain HEXADEC HexadecimalExpansion
++ Author: Clifton J. Williamson
++ Date Created: April 26, 1990
++ Basic Operations:
++ Related Domains: RadixExpansion
++ Also See:
++ AMS Classifications:
++ Keywords: radix, base, hexadecimal
++ Examples:
++ References:
++ Description:
++   This domain allows rational numbers to be presented as repeating
++   hexadecimal expansions.

HexadecimalExpansion() : Exports == Implementation where
  Exports ==> QuotientFieldCategory(Integer) with
    coerce : % -> Fraction Integer
      ++ coerce(h) converts a hexadecimal expansion to a rational number.
    coerce : % -> RadixExpansion(16)
      ++ coerce(h) converts a hexadecimal expansion to a radix expansion
      ++ with base 16.
    fractionPart : % -> Fraction Integer
      ++ fractionPart(h) returns the fractional part of a hexadecimal expansion.
    hex : Fraction Integer -> %
      ++ hex(r) converts a rational number to a hexadecimal expansion.

  Implementation ==> RadixExpansion(16) add
    hex r == r :: %
    coerce(x : %) : RadixExpansion(16) == x pretend RadixExpansion(16)

)abbrev package RADUTIL RadixUtilities
++ Author: Stephen M. Watt
++ Date Created: October 1986
++ Basic Operations:
++ Related Domains: RadixExpansion
++ Also See:
++ AMS Classifications:
++ Keywords: radix, base, repeading decimal
++ Examples:
++ References:
++ Description:
++   This package provides tools for creating radix expansions.
RadixUtilities : Exports == Implementation where
  Exports ==> with
    radix : (Fraction Integer, Integer) -> Any
      ++ radix(x, b) converts x to a radix expansion in base b.
  Implementation ==> add
    radix(q, b) ==
      coerce(q :: RadixExpansion(b))$AnyFunctions1(RadixExpansion b)

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