275 lines
6.3 KiB
Plaintext
275 lines
6.3 KiB
Plaintext
import "./dynamic" for Struct
|
|
|
|
// constants
|
|
var EMX = 64 // exponent maximum (if indexing starts at -EMX)
|
|
var DMX = 1e5 // approximation loop maximum
|
|
var AMX = 1048576 // argument maximum
|
|
var PMAX = 32749 // prime maximum
|
|
|
|
// global variables
|
|
var P1 = 0
|
|
var P = 7 // default prime
|
|
var K = 11 // precision
|
|
|
|
var Ratio = Struct.create("Ratio", ["a", "b"])
|
|
|
|
class Padic {
|
|
// uninitialized
|
|
construct new() {
|
|
_v = 0
|
|
_d = List.filled(2 * EMX, 0) // add EMX to index to be consistent wih FB
|
|
}
|
|
|
|
// properties
|
|
v { _v }
|
|
v=(o) { _v = o }
|
|
d { _d }
|
|
|
|
// (re)initialize 'this' from Ratio, set 'sw' to print
|
|
r2pa(q, sw) {
|
|
var a = q.a
|
|
var b = q.b
|
|
if (b == 0) return 1
|
|
if (b < 0) {
|
|
b = -b
|
|
a = -a
|
|
}
|
|
if (a.abs > AMX || b > AMX) return -1
|
|
if (P < 2 || K < 1) return 1
|
|
P = P.min(PMAX) // maximum short prime
|
|
K = K.min(EMX-1) // maximum array length
|
|
if (sw != 0) {
|
|
System.write("%(a)/%(b) + ") // numerator, denominator
|
|
System.print("0(%(P)^%(K))") // prime, precision
|
|
}
|
|
|
|
// (re)initialize
|
|
_v = 0
|
|
P1 = P - 1
|
|
_d = List.filled(2 * EMX, 0)
|
|
if (a == 0) return 0
|
|
var i = 0
|
|
// find -exponent of P in b
|
|
while (b%P == 0) {
|
|
b = (b/P).truncate
|
|
i = i - 1
|
|
}
|
|
var s = 0
|
|
var r = b % P
|
|
|
|
// modular inverse for small P
|
|
var b1 = 1
|
|
while (b1 <= P1) {
|
|
s = s + r
|
|
if (s > P1) s = s - P
|
|
if (s == 1) break
|
|
b1 = b1 + 1
|
|
}
|
|
if (b1 == P) {
|
|
System.print("r2pa: impossible inverse mod")
|
|
return -1
|
|
}
|
|
_v = EMX
|
|
while (true) {
|
|
// find exponent of P in a
|
|
while (a%P == 0) {
|
|
a = (a/P).truncate
|
|
i = i + 1
|
|
}
|
|
|
|
// valuation
|
|
if (_v == EMX) _v = i
|
|
|
|
// upper bound
|
|
if (i >= EMX) break
|
|
|
|
// check precision
|
|
if ((i - _v) > K) break
|
|
|
|
// next digit
|
|
_d[i+EMX] = a * b1 % P
|
|
if (_d[i+EMX] < 0) _d[i+EMX] = _d[i+EMX] + P
|
|
|
|
// remainder - digit * divisor
|
|
a = a - _d[i+EMX]*b
|
|
if (a == 0) break
|
|
}
|
|
return 0
|
|
}
|
|
|
|
// Horner's rule
|
|
dsum() {
|
|
var t = _v.min(0)
|
|
var s = 0
|
|
for (i in K - 1 + t..t) {
|
|
var r = s
|
|
s = s * P
|
|
if (r != 0 && ((s/r).truncate - P) != 0) {
|
|
// overflow
|
|
s = -1
|
|
break
|
|
}
|
|
s = s + _d[i+EMX]
|
|
}
|
|
return s
|
|
}
|
|
|
|
// rational reconstruction
|
|
crat() {
|
|
var fl = false
|
|
var s = this
|
|
var j = 0
|
|
var i = 1
|
|
|
|
// denominator count
|
|
while (i <= DMX) {
|
|
// check for integer
|
|
j = K - 1 + _v
|
|
while (j >= _v) {
|
|
if (s.d[j+EMX] != 0) break
|
|
j = j - 1
|
|
}
|
|
fl = ((j - _v) * 2) < K
|
|
if (fl) {
|
|
fl = false
|
|
break
|
|
}
|
|
|
|
// check negative integer
|
|
j = K - 1 + _v
|
|
while (j >= _v) {
|
|
if (P1 - s.d[j+EMX] != 0) break
|
|
j = j - 1
|
|
}
|
|
fl = ((j - _v) * 2) < K
|
|
if (fl) break
|
|
|
|
// repeatedly add self to s
|
|
s = s + this
|
|
i = i + 1
|
|
}
|
|
if (fl) s = s.cmpt
|
|
|
|
// numerator: weighted digit sum
|
|
var x = s.dsum()
|
|
var y = i
|
|
if (x < 0 || y > DMX) {
|
|
System.print("crat: fail")
|
|
} else {
|
|
// negative powers
|
|
i = _v
|
|
while (i <= -1) {
|
|
y = y * P
|
|
i = i + 1
|
|
}
|
|
|
|
// negative rational
|
|
if (fl) x = -x
|
|
System.write(x)
|
|
if (y > 1) System.write("/%(y)")
|
|
System.print()
|
|
}
|
|
}
|
|
|
|
// print expansion
|
|
printf(sw) {
|
|
var t = _v.min(0)
|
|
for (i in K - 1 + t..t) {
|
|
System.write(_d[i + EMX])
|
|
if (i == 0 && _v < 0) System.write(".")
|
|
System.write(" ")
|
|
}
|
|
System.print()
|
|
// rational approximation
|
|
if (sw != 0) crat()
|
|
}
|
|
|
|
// add b to 'this'
|
|
+(b) {
|
|
var c = 0
|
|
var r = Padic.new()
|
|
r.v = _v.min(b.v)
|
|
for (i in r.v..K + r.v) {
|
|
c = c + _d[i+EMX] + b.d[i+EMX]
|
|
if (c > P1) {
|
|
r.d[i+EMX] = c - P
|
|
c = 1
|
|
} else {
|
|
r.d[i+EMX] = c
|
|
c = 0
|
|
}
|
|
}
|
|
return r
|
|
}
|
|
|
|
// complement
|
|
cmpt {
|
|
var c = 1
|
|
var r = Padic.new()
|
|
r.v = _v
|
|
for (i in _v..K + _v) {
|
|
c = c + P1 - _d[i+EMX]
|
|
if (c > P1) {
|
|
r.d[i+EMX] = c - P
|
|
c = 1
|
|
} else {
|
|
r.d[i+EMX] = c
|
|
c = 0
|
|
}
|
|
}
|
|
return r
|
|
}
|
|
}
|
|
|
|
var data = [
|
|
/* rational reconstruction depends on the precision
|
|
until the dsum-loop overflows */
|
|
[2, 1, 2, 4, 1, 1],
|
|
[4, 1, 2, 4, 3, 1],
|
|
[4, 1, 2, 5, 3, 1],
|
|
[4, 9, 5, 4, 8, 9],
|
|
[26, 25, 5, 4, -109, 125],
|
|
[49, 2, 7, 6, -4851, 2],
|
|
[-9, 5, 3, 8, 27, 7],
|
|
[5, 19, 2, 12, -101, 384],
|
|
/* two decadic pairs */
|
|
[2, 7, 10, 7, -1, 7],
|
|
[34, 21, 10, 9, -39034, 791],
|
|
/* familiar digits */
|
|
[11, 4, 2, 43, 679001, 207],
|
|
[-8, 9, 23, 9, 302113, 92],
|
|
[-22, 7, 3, 23, 46071, 379],
|
|
[-22, 7, 32749, 3, 46071, 379],
|
|
[35, 61, 5, 20, 9400, 109],
|
|
[-101, 109, 61, 7, 583376, 6649],
|
|
[-25, 26, 7, 13, 5571, 137],
|
|
[1, 4, 7, 11, 9263, 2837],
|
|
[122, 407, 7, 11, -517, 1477],
|
|
/* more subtle */
|
|
[5, 8, 7, 11, 353, 30809]
|
|
]
|
|
|
|
var sw = 0
|
|
var a = Padic.new()
|
|
var b = Padic.new()
|
|
|
|
for (d in data) {
|
|
var q = Ratio.new(d[0], d[1])
|
|
P = d[2]
|
|
K = d[3]
|
|
sw = a.r2pa(q, 1)
|
|
if (sw == 1) break
|
|
a.printf(0)
|
|
q.a = d[4]
|
|
q.b = d[5]
|
|
sw = sw | b.r2pa(q, 1)
|
|
if (sw == 1) break
|
|
if (sw == 0) {
|
|
b.printf(0)
|
|
var c = a + b
|
|
System.print("+ =")
|
|
c.printf(1)
|
|
}
|
|
System.print()
|
|
}
|