331 lines
6.2 KiB
Go
331 lines
6.2 KiB
Go
package main
|
|
|
|
import "fmt"
|
|
|
|
// constants
|
|
const EMX = 64 // exponent maximum (if indexing starts at -EMX)
|
|
const DMX = 100000 // approximation loop maximum
|
|
const AMX = 1048576 // argument maximum
|
|
const PMAX = 32749 // prime maximum
|
|
|
|
// global variables
|
|
var p1 = 0
|
|
var p = 7 // default prime
|
|
var k = 11 // precision
|
|
|
|
func abs(a int) int {
|
|
if a >= 0 {
|
|
return a
|
|
}
|
|
return -a
|
|
}
|
|
|
|
func min(a, b int) int {
|
|
if a < b {
|
|
return a
|
|
}
|
|
return b
|
|
}
|
|
|
|
type Ratio struct {
|
|
a, b int
|
|
}
|
|
|
|
type Padic struct {
|
|
v int
|
|
d [2 * EMX]int // add EMX to index to be consistent wih FB
|
|
}
|
|
|
|
// (re)initialize receiver from Ratio, set 'sw' to print
|
|
func (pa *Padic) r2pa(q Ratio, sw int) int {
|
|
a := q.a
|
|
b := q.b
|
|
if b == 0 {
|
|
return 1
|
|
}
|
|
if b < 0 {
|
|
b = -b
|
|
a = -a
|
|
}
|
|
if abs(a) > AMX || b > AMX {
|
|
return -1
|
|
}
|
|
if p < 2 || k < 1 {
|
|
return 1
|
|
}
|
|
p = min(p, PMAX) // maximum short prime
|
|
k = min(k, EMX-1) // maxumum array length
|
|
if sw != 0 {
|
|
fmt.Printf("%d/%d + ", a, b) // numerator, denominator
|
|
fmt.Printf("0(%d^%d)\n", p, k) // prime, precision
|
|
}
|
|
|
|
// (re)initialize
|
|
pa.v = 0
|
|
p1 = p - 1
|
|
pa.d = [2 * EMX]int{}
|
|
if a == 0 {
|
|
return 0
|
|
}
|
|
i := 0
|
|
|
|
// find -exponent of p in b
|
|
for b%p == 0 {
|
|
b = b / p
|
|
i--
|
|
}
|
|
s := 0
|
|
r := b % p
|
|
|
|
// modular inverse for small p
|
|
b1 := 1
|
|
for b1 <= p1 {
|
|
s += r
|
|
if s > p1 {
|
|
s -= p
|
|
}
|
|
if s == 1 {
|
|
break
|
|
}
|
|
b1++
|
|
}
|
|
if b1 == p {
|
|
fmt.Println("r2pa: impossible inverse mod")
|
|
return -1
|
|
}
|
|
pa.v = EMX
|
|
for {
|
|
// find exponent of P in a
|
|
for a%p == 0 {
|
|
a = a / p
|
|
i++
|
|
}
|
|
|
|
// valuation
|
|
if pa.v == EMX {
|
|
pa.v = i
|
|
}
|
|
|
|
// upper bound
|
|
if i >= EMX {
|
|
break
|
|
}
|
|
|
|
// check precision
|
|
if (i - pa.v) > k {
|
|
break
|
|
}
|
|
|
|
// next digit
|
|
pa.d[i+EMX] = a * b1 % p
|
|
if pa.d[i+EMX] < 0 {
|
|
pa.d[i+EMX] += p
|
|
}
|
|
|
|
// remainder - digit * divisor
|
|
a -= pa.d[i+EMX] * b
|
|
if a == 0 {
|
|
break
|
|
}
|
|
}
|
|
return 0
|
|
}
|
|
|
|
// Horner's rule
|
|
func (pa *Padic) dsum() int {
|
|
t := min(pa.v, 0)
|
|
s := 0
|
|
for i := k - 1 + t; i >= t; i-- {
|
|
r := s
|
|
s *= p
|
|
if r != 0 && (s/r-p != 0) {
|
|
// overflow
|
|
s = -1
|
|
break
|
|
}
|
|
s += pa.d[i+EMX]
|
|
}
|
|
return s
|
|
}
|
|
|
|
// add b to receiver
|
|
func (pa *Padic) add(b Padic) *Padic {
|
|
c := 0
|
|
r := Padic{}
|
|
r.v = min(pa.v, b.v)
|
|
for i := r.v; i <= k+r.v; i++ {
|
|
c += pa.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 of receiver
|
|
func (pa *Padic) cmpt() *Padic {
|
|
c := 1
|
|
r := Padic{}
|
|
r.v = pa.v
|
|
for i := pa.v; i <= k+pa.v; i++ {
|
|
c += p1 - pa.d[i+EMX]
|
|
if c > p1 {
|
|
r.d[i+EMX] = c - p
|
|
c = 1
|
|
} else {
|
|
r.d[i+EMX] = c
|
|
c = 0
|
|
}
|
|
}
|
|
return &r
|
|
}
|
|
|
|
// rational reconstruction
|
|
func (pa *Padic) crat() {
|
|
fl := false
|
|
s := pa
|
|
j := 0
|
|
i := 1
|
|
|
|
// denominator count
|
|
for i <= DMX {
|
|
// check for integer
|
|
j = k - 1 + pa.v
|
|
for j >= pa.v {
|
|
if s.d[j+EMX] != 0 {
|
|
break
|
|
}
|
|
j--
|
|
}
|
|
fl = ((j - pa.v) * 2) < k
|
|
if fl {
|
|
fl = false
|
|
break
|
|
}
|
|
|
|
// check negative integer
|
|
j = k - 1 + pa.v
|
|
for j >= pa.v {
|
|
if p1-s.d[j+EMX] != 0 {
|
|
break
|
|
}
|
|
j--
|
|
}
|
|
fl = ((j - pa.v) * 2) < k
|
|
if fl {
|
|
break
|
|
}
|
|
|
|
// repeatedly add self to s
|
|
s = s.add(*pa)
|
|
i++
|
|
}
|
|
if fl {
|
|
s = s.cmpt()
|
|
}
|
|
|
|
// numerator: weighted digit sum
|
|
x := s.dsum()
|
|
y := i
|
|
if x < 0 || y > DMX {
|
|
fmt.Println(x, y)
|
|
fmt.Println("crat: fail")
|
|
} else {
|
|
// negative powers
|
|
i = pa.v
|
|
for i <= -1 {
|
|
y *= p
|
|
i++
|
|
}
|
|
|
|
// negative rational
|
|
if fl {
|
|
x = -x
|
|
}
|
|
fmt.Print(x)
|
|
if y > 1 {
|
|
fmt.Printf("/%d", y)
|
|
}
|
|
fmt.Println()
|
|
}
|
|
}
|
|
|
|
// print expansion
|
|
func (pa *Padic) printf(sw int) {
|
|
t := min(pa.v, 0)
|
|
for i := k - 1 + t; i >= t; i-- {
|
|
fmt.Print(pa.d[i+EMX])
|
|
if i == 0 && pa.v < 0 {
|
|
fmt.Print(".")
|
|
}
|
|
fmt.Print(" ")
|
|
}
|
|
fmt.Println()
|
|
// rational approximation
|
|
if sw != 0 {
|
|
pa.crat()
|
|
}
|
|
}
|
|
|
|
func main() {
|
|
data := [][]int{
|
|
/* 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},
|
|
}
|
|
|
|
sw := 0
|
|
a := Padic{}
|
|
b := Padic{}
|
|
|
|
for _, d := range data {
|
|
q := Ratio{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)
|
|
c := a.add(b)
|
|
fmt.Println("+ =")
|
|
c.printf(1)
|
|
}
|
|
fmt.Println()
|
|
}
|
|
}
|