47 lines
1.4 KiB
Nim
47 lines
1.4 KiB
Nim
import strutils, math, sequtils, times
|
|
|
|
proc getStep(n: int64) : int64 {.inline.} =
|
|
result = 1 + n*4 - int64(n /% 2)*2
|
|
|
|
proc primeFac(n: int64): seq[int64] =
|
|
var res: seq[int64] = @[]
|
|
var maxq = int64(floor(sqrt(float(n))))
|
|
var d = 1
|
|
var q: int64 = (n %% 2) and 2 or 3 # either 2 or 3, alternating
|
|
while (q <= maxq) and ((n %% q) != 0):
|
|
q = getStep(d)
|
|
d += 1
|
|
if q <= maxq:
|
|
var q1: seq[int64] = primeFac(n /% q)
|
|
var q2: seq[int64] = primeFac(q)
|
|
res = concat(q2, q1, res)
|
|
else:
|
|
res.add(n)
|
|
result = res
|
|
|
|
var is_prime: seq[Bool] = @[]
|
|
is_prime.add(False)
|
|
is_prime.add(False)
|
|
|
|
iterator primes(limit: int): int =
|
|
for n in high(is_prime) .. limit+2: is_prime.add(True)
|
|
for n in 2 .. limit + 1:
|
|
if is_prime[n]:
|
|
yield n
|
|
for i in countup((n *% n), limit+1, n): # start at ``n`` squared
|
|
try:
|
|
is_prime[i] = False
|
|
except EInvalidIndex: break
|
|
|
|
# Example: calculate factors of Mersenne numbers to M59 #
|
|
|
|
for m in primes(59):
|
|
var p = int64(pow(2.0,float(m)) - 1)
|
|
write(stdout,"2**$1-1 = $2, with factors: " % [$m, $p] )
|
|
var start = cpuTime()
|
|
var f = primeFac(p)
|
|
for factor in f:
|
|
write(stdout, factor)
|
|
write(stdout, ", ")
|
|
FlushFile(stdout)
|
|
writeln(stdout, "=> $#ms" % $int(1000*(cpuTime()-start)) )
|