358 lines
8.4 KiB
Zig
358 lines
8.4 KiB
Zig
const std = @import("std");
|
|
const math = std.math;
|
|
const stdout = std.io.getStdOut().writer();
|
|
const print = std.debug.print;
|
|
|
|
// Constants
|
|
const EMX: i32 = 64; // exponent maximum (if indexing starts at -EMX)
|
|
const DMX: i32 = 100000; // approximation loop maximum
|
|
const AMX: i32 = 1048576; // argument maximum
|
|
const PMAX: i32 = 32749; // prime maximum
|
|
|
|
// Global variables
|
|
var P1: i32 = 0;
|
|
var P: i32 = 7; // default prime
|
|
var K: i32 = 11; // precision
|
|
|
|
// Helper functions
|
|
fn abs(a: i32) i32 {
|
|
return if (a >= 0) a else -a;
|
|
}
|
|
|
|
fn min(a: i32, b: i32) i32 {
|
|
return if (a < b) a else b;
|
|
}
|
|
|
|
const Ratio = struct {
|
|
a: i32,
|
|
b: i32,
|
|
};
|
|
|
|
const Padic = struct {
|
|
v: i32,
|
|
d: [2 * EMX]i32, // add EMX to index to be consistent
|
|
|
|
// Create a new Padic with default values
|
|
fn new() Padic {
|
|
return Padic{
|
|
.v = 0,
|
|
.d = [_]i32{0} ** (2 * EMX),
|
|
};
|
|
}
|
|
|
|
// (re)initialize receiver from Ratio, set 'sw' to print
|
|
fn r2pa(self: *Padic, q: Ratio, sw: i32) i32 {
|
|
var a = q.a;
|
|
var b = q.b;
|
|
if (b == 0) {
|
|
return 1;
|
|
}
|
|
if (b < 0) {
|
|
b = -b;
|
|
a = -a;
|
|
}
|
|
if (abs(a) > AMX or b > AMX) {
|
|
return -1;
|
|
}
|
|
|
|
if (P < 2 or K < 1) {
|
|
return 1;
|
|
}
|
|
|
|
// Set P to minimum of P and PMAX
|
|
P = @min(P, PMAX);
|
|
|
|
// Set K to minimum of K and EMX-1
|
|
K = @min(K, EMX - 1);
|
|
|
|
if (sw != 0) {
|
|
print("{d}/{d} + \n", .{a, b}); // numerator, denominator
|
|
print("0({d}^{d})\n", .{P, K}); // prime, precision
|
|
}
|
|
|
|
// (re)initialize
|
|
self.v = 0;
|
|
P1 = P - 1;
|
|
// std.mem.set(i32, &self.d, 0);
|
|
|
|
if (a == 0) {
|
|
return 0;
|
|
}
|
|
|
|
var i: i32 = 0;
|
|
|
|
// find -exponent of p in b
|
|
while ( @rem(b, P) == 0) {
|
|
b = @divTrunc(b, P);
|
|
i -= 1;
|
|
}
|
|
|
|
var s: i32 = 0;
|
|
const r = @rem(b, P);
|
|
|
|
// modular inverse for small p
|
|
var b1: i32 = 1;
|
|
while (b1 <= P1) {
|
|
s += r;
|
|
if (s > P1) {
|
|
s -= P;
|
|
}
|
|
if (s == 1) {
|
|
break;
|
|
}
|
|
b1 += 1;
|
|
}
|
|
|
|
if (b1 == P) {
|
|
print("r2pa: impossible inverse mod\n", .{});
|
|
return -1;
|
|
}
|
|
|
|
self.v = EMX;
|
|
|
|
while (true) {
|
|
// find exponent of P in a
|
|
while ( @rem(a, P)==0) {
|
|
a = @divTrunc(a, P);
|
|
i += 1;
|
|
}
|
|
|
|
// valuation
|
|
if (self.v == EMX) {
|
|
self.v = i;
|
|
}
|
|
|
|
// upper bound
|
|
if (i >= EMX) {
|
|
break;
|
|
}
|
|
|
|
// check precision
|
|
if ((i - self.v) > K) {
|
|
break;
|
|
}
|
|
|
|
// next digit
|
|
self.d[@intCast(i + EMX)] = @mod(a * b1, P);
|
|
if (self.d[@intCast(i + EMX)] < 0) {
|
|
self.d[@intCast(i + EMX)] += P;
|
|
}
|
|
|
|
// remainder - digit * divisor
|
|
a -= self.d[@intCast(i + EMX)] * b;
|
|
if (a == 0) {
|
|
break;
|
|
}
|
|
}
|
|
|
|
return 0;
|
|
}
|
|
|
|
// Horner's rule
|
|
fn dsum(self: *const Padic) i32 {
|
|
const t = @min(self.v, 0);
|
|
var s: i32 = 0;
|
|
|
|
var i = K - 1 + t;
|
|
while (i >= t) : (i -= 1) {
|
|
const r = s;
|
|
s *= P;
|
|
if (r != 0 and (@divTrunc(s, r) - P != 0)) {
|
|
// overflow
|
|
s = -1;
|
|
break;
|
|
}
|
|
s += self.d[@intCast(i + EMX)];
|
|
}
|
|
return s;
|
|
}
|
|
|
|
// add b to receiver
|
|
fn add(self: *const Padic, b: *const Padic) Padic {
|
|
var c: i32 = 0;
|
|
var r = Padic.new();
|
|
r.v = @min(self.v, b.v);
|
|
|
|
var i: i32 = r.v;
|
|
while (i <= K + r.v) : (i += 1) {
|
|
c += self.d[@intCast(i + EMX)] + b.d[@intCast(i + EMX)];
|
|
if (c > P1) {
|
|
r.d[@intCast(i + EMX)] = c - P;
|
|
c = 1;
|
|
} else {
|
|
r.d[@intCast(i + EMX)] = c;
|
|
c = 0;
|
|
}
|
|
}
|
|
|
|
return r;
|
|
}
|
|
|
|
// complement of receiver
|
|
fn cmpt(self: *const Padic) Padic {
|
|
var c: i32 = 1;
|
|
var r = Padic.new();
|
|
r.v = self.v;
|
|
|
|
var i: i32 = self.v;
|
|
while (i <= K + self.v) : (i += 1) {
|
|
c += P1 - self.d[@intCast(i + EMX)];
|
|
if (c > P1) {
|
|
r.d[@intCast(i + EMX)] = c - P;
|
|
c = 1;
|
|
} else {
|
|
r.d[@intCast(i + EMX)] = c;
|
|
c = 0;
|
|
}
|
|
}
|
|
|
|
return r;
|
|
}
|
|
|
|
// rational reconstruction
|
|
fn crat(self: *const Padic) void {
|
|
var fl = false;
|
|
var s = self.*;
|
|
var j: i32 = 0;
|
|
var i: i32 = 1;
|
|
|
|
// denominator count
|
|
while (i <= DMX) : (i += 1) {
|
|
// check for integer
|
|
j = K - 1 + self.v;
|
|
while (j >= self.v) : (j -= 1) {
|
|
if (s.d[@intCast(j + EMX)] != 0) {
|
|
break;
|
|
}
|
|
}
|
|
fl = ((j - self.v) * 2) < K;
|
|
if (fl) {
|
|
fl = false;
|
|
break;
|
|
}
|
|
|
|
// check negative integer
|
|
j = K - 1 + self.v;
|
|
while (j >= self.v) : (j -= 1) {
|
|
if (P1 - s.d[@intCast(j + EMX)] != 0) {
|
|
break;
|
|
}
|
|
}
|
|
fl = ((j - self.v) * 2) < K;
|
|
if (fl) {
|
|
break;
|
|
}
|
|
|
|
// repeatedly add self to s
|
|
s = s.add(self);
|
|
}
|
|
|
|
if (fl) {
|
|
s = s.cmpt();
|
|
}
|
|
|
|
// numerator: weighted digit sum
|
|
const x = s.dsum();
|
|
var y = i;
|
|
|
|
if (x < 0 or y > DMX) {
|
|
print("{d} {d}\n", .{x, y});
|
|
print("crat: fail\n", .{});
|
|
} else {
|
|
// negative powers
|
|
var i_pow = self.v;
|
|
while (i_pow <= -1) : (i_pow += 1) {
|
|
y *= P;
|
|
}
|
|
|
|
// negative rational
|
|
const final_x = if (fl) -x else x;
|
|
print("{d}", .{final_x});
|
|
if (y > 1) {
|
|
print("/{d}", .{y});
|
|
}
|
|
print("\n", .{});
|
|
}
|
|
}
|
|
|
|
// print expansion
|
|
fn printf(self: *const Padic, sw: i32) void {
|
|
const t = @min(self.v, 0);
|
|
|
|
var i = K - 1 + t;
|
|
while (i >= t) : (i -= 1) {
|
|
print("{d}", .{self.d[@intCast(i + EMX)]});
|
|
if (i == 0 and self.v < 0) {
|
|
print(".", .{});
|
|
}
|
|
print(" ", .{});
|
|
}
|
|
print("\n", .{});
|
|
// rational approximation
|
|
if (sw != 0) {
|
|
self.crat();
|
|
}
|
|
}
|
|
};
|
|
|
|
pub fn main() !void {
|
|
const data = [_][6]i32{
|
|
// rational reconstruction depends on the precision
|
|
// until the dsum-loop overflows
|
|
[_]i32{ 2, 1, 2, 4, 1, 1 },
|
|
[_]i32{ 4, 1, 2, 4, 3, 1 },
|
|
[_]i32{ 4, 1, 2, 5, 3, 1 },
|
|
[_]i32{ 4, 9, 5, 4, 8, 9 },
|
|
[_]i32{ 26, 25, 5, 4, -109, 125 },
|
|
[_]i32{ 49, 2, 7, 6, -4851, 2 },
|
|
[_]i32{ -9, 5, 3, 8, 27, 7 },
|
|
[_]i32{ 5, 19, 2, 12, -101, 384 },
|
|
// two decadic pairs
|
|
[_]i32{ 2, 7, 10, 7, -1, 7 },
|
|
[_]i32{ 34, 21, 10, 9, -39034, 791 },
|
|
// familiar digits
|
|
[_]i32{ 11, 4, 2, 43, 679001, 207 },
|
|
[_]i32{ -8, 9, 23, 9, 302113, 92 },
|
|
[_]i32{ -22, 7, 3, 23, 46071, 379 },
|
|
[_]i32{ -22, 7, 32749, 3, 46071, 379 },
|
|
[_]i32{ 35, 61, 5, 20, 9400, 109 },
|
|
[_]i32{ -101, 109, 61, 7, 583376, 6649 },
|
|
[_]i32{ -25, 26, 7, 13, 5571, 137 },
|
|
[_]i32{ 1, 4, 7, 11, 9263, 2837 },
|
|
[_]i32{ 122, 407, 7, 11, -517, 1477 },
|
|
// more subtle
|
|
[_]i32{ 5, 8, 7, 11, 353, 30809 },
|
|
};
|
|
|
|
var sw: i32 = 0;
|
|
var a = Padic.new();
|
|
var b = Padic.new();
|
|
|
|
for (data) |d| {
|
|
const q = Ratio{ .a = d[0], .b = d[1] };
|
|
|
|
P = d[2];
|
|
K = d[3];
|
|
|
|
sw = a.r2pa(q, 1);
|
|
if (sw == 1) {
|
|
break;
|
|
}
|
|
a.printf(0);
|
|
|
|
const q_b = Ratio{ .a = d[4], .b = d[5] };
|
|
sw = sw | b.r2pa(q_b, 1);
|
|
if (sw == 1) {
|
|
break;
|
|
}
|
|
|
|
if (sw == 0) {
|
|
b.printf(0);
|
|
const c = a.add(&b);
|
|
print("+ =\n", .{});
|
|
c.printf(1);
|
|
}
|
|
print("\n", .{});
|
|
}
|
|
}
|