97 lines
2.4 KiB
Nim
97 lines
2.4 KiB
Nim
import math, strformat
|
|
|
|
const M = [uint64 1, 3, 5, 7, 11]
|
|
|
|
template isqrt(n: uint64): uint64 = uint64(sqrt(float(n)))
|
|
template isEven(n: uint64): bool = (n and 1) == 0
|
|
|
|
proc squfof(n: uint64): uint64 =
|
|
|
|
if n.isEven: return 2
|
|
var h = uint64(sqrt(float(n)) + 0.5)
|
|
if h * h == n: return h
|
|
|
|
for m in M:
|
|
if m > 1 and (n mod m == 0): return m
|
|
# Check overflow m * n.
|
|
if n > uint64.high div m: break
|
|
let mn = m * n
|
|
var r = isqrt(mn)
|
|
if r * r > mn: dec r
|
|
let rn = r
|
|
|
|
# Principal form.
|
|
var b = r
|
|
var a = 1u64
|
|
h = (rn + b) div a * a - b
|
|
var c = (mn - h * h) div a
|
|
|
|
for i in 2..<(4 * isqrt(2 * r)):
|
|
# Search principal cycle.
|
|
swap a, c
|
|
var q = (rn + b) div a
|
|
let t = b
|
|
b = q * a - b
|
|
c += q * (t - b)
|
|
|
|
if i.isEven:
|
|
r = uint64(sqrt(float(c)) + 0.5)
|
|
if r * r == c: # Square form found?
|
|
|
|
# Inverse square root.
|
|
q = (rn - b) div r
|
|
var v = q * r + b
|
|
var w = (mn - v * v) div r
|
|
|
|
# Search ambiguous cycle.
|
|
var u = r
|
|
while true:
|
|
swap w, u
|
|
r = v
|
|
q = (rn + v) div u
|
|
v = q * u - v
|
|
if v == r: break
|
|
w += q * (r - v)
|
|
|
|
# Symmetry point.
|
|
h = gcd(n, u)
|
|
if h != 1: return h
|
|
|
|
result = 1
|
|
|
|
const Data = [2501u64,
|
|
12851u64,
|
|
13289u64,
|
|
75301u64,
|
|
120787u64,
|
|
967009u64,
|
|
997417u64,
|
|
7091569u64,
|
|
13290059u64,
|
|
42854447u64,
|
|
223553581u64,
|
|
2027651281u64,
|
|
11111111111u64,
|
|
100895598169u64,
|
|
1002742628021u64,
|
|
60012462237239u64,
|
|
287129523414791u64,
|
|
9007199254740931u64,
|
|
11111111111111111u64,
|
|
314159265358979323u64,
|
|
384307168202281507u64,
|
|
419244183493398773u64,
|
|
658812288346769681u64,
|
|
922337203685477563u64,
|
|
1000000000000000127u64,
|
|
1152921505680588799u64,
|
|
1537228672809128917u64,
|
|
4611686018427387877u64]
|
|
|
|
echo "N f N/f"
|
|
echo "======================================"
|
|
for n in Data:
|
|
let f = squfof(n)
|
|
let res = if f == 1: "fail" else: &"{f:<10} {n div f}"
|
|
echo &"{n:<22} {res}"
|