60 lines
1.4 KiB
Zig
60 lines
1.4 KiB
Zig
const std = @import("std");
|
|
|
|
pub fn main() !void {
|
|
const writer = std.io.getStdOut().writer();
|
|
|
|
try printSolvedPell(61, writer);
|
|
try printSolvedPell(109, writer);
|
|
try printSolvedPell(181, writer);
|
|
try printSolvedPell(277, writer);
|
|
}
|
|
|
|
const Pair = struct {
|
|
v1: u256,
|
|
v2: u256,
|
|
|
|
fn init(a: u256, b: u256) Pair {
|
|
return Pair{
|
|
.v1 = a,
|
|
.v2 = b,
|
|
};
|
|
}
|
|
};
|
|
|
|
fn solvePell(n: u256) Pair {
|
|
const x: u256 = std.math.sqrt(n);
|
|
|
|
// n is a perfect square - no solution other than 1,0
|
|
if (x * x == n)
|
|
return Pair.init(1, 0);
|
|
|
|
// there are non-trivial solutions
|
|
var y = x;
|
|
var z: u256 = 1;
|
|
var r = 2 * x;
|
|
var e = Pair.init(1, 0);
|
|
var f = Pair.init(0, 1);
|
|
var a: u256 = 0;
|
|
var b: u256 = 0;
|
|
|
|
while (true) {
|
|
y = r * z - y;
|
|
z = (n - y * y) / z;
|
|
r = (x + y) / z;
|
|
e = Pair.init(e.v2, r * e.v2 + e.v1);
|
|
f = Pair.init(f.v2, r * f.v2 + f.v1);
|
|
a = e.v2 + x * f.v2;
|
|
b = f.v2;
|
|
const ov = @subWithOverflow(a * a, n * b * b);
|
|
if (ov[1] != 0)
|
|
continue;
|
|
if (ov[0] == 1) // a * a, n * b * b == 1
|
|
break;
|
|
}
|
|
return Pair.init(a, b);
|
|
}
|
|
|
|
fn printSolvedPell(n: u256, writer: anytype) !void {
|
|
const r = solvePell(n);
|
|
try writer.print("x^2 - {d:3} * y^2 = 1 for x = {d:21} and y = {d:19}\n", .{ n, r.v1, r.v2 });
|
|
}
|