/usr/share/axiom-20170501/src/algebra/RADIX.spad is in axiom-source 20170501-3.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 | )abbrev domain RADIX RadixExpansion
++ Author: Stephen M. Watt
++ Date Created: October 1986
++ Date Last Updated: May 15, 1991
++ Description:
++ This domain allows rational numbers to be presented as repeating
++ decimal expansions or more generally as repeating expansions in any base.
RadixExpansion(bb) : SIG == CODE where
bb : Integer
I ==> Integer
NNI ==> NonNegativeInteger
OUT ==> OutputForm
RN ==> Fraction Integer
ST ==> Stream Integer
QuoRem ==> Record(quotient: Integer, remainder: Integer)
SIG ==> 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...}.
CODE ==> 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, nil(), nil(), nil()]
1 == [1, [1], nil(), nil()]
- 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] => nil(); 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 ==
n0 := 0
for r in a.pfx repeat n0 := bb*n0 + r
null 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, nil(), nil()]
fractRadix(lpfx, lcyc) ==
checkRagits lpfx; checkRagits lcyc
[1, nil(), 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 => intToExpr first(li)
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 := nil()
if not null a.cyc then le := concat(overBar a.cyc,le)
if not null a.pfx then le := concat(intgroup a.pfx,le)
if not null le then le := concat("." :: OUT,le)
if not null 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 := nil()
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
-- 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 := nil()
ritscyc: List I := nil()
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]
|