Square form factorization: Difference between revisions

=={{header|Nim}}==

{{trans|Phix}}

<lang 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}"</lang>

{{out}}

<pre>N f N/f

======================================

2501 61 41

12851 71 181

13289 97 137

75301 293 257

120787 43 2809

967009 1609 601

997417 257 3881

7091569 2663 2663

13290059 3119 4261

42854447 4423 9689

223553581 11213 19937

2027651281 44021 46061

11111111111 21649 513239

100895598169 112303 898423

1002742628021 fail

60012462237239 6862753 8744663

287129523414791 6059887 47381993

9007199254740931 10624181 847801751

11111111111111111 2071723 5363222357

314159265358979323 317213509 990371647

384307168202281507 415718707 924440401

419244183493398773 48009977 8732438749

658812288346769681 62222119 10588072199

922337203685477563 110075821 8379108103

1000000000000000127 111756107 8948056861

1152921505680588799 139001459 8294312261

1537228672809128917 26675843 57626245319

4611686018427387877 343242169 13435662733</pre>