-
Notifications
You must be signed in to change notification settings - Fork 6
/
675 2 to omega of n.jl
58 lines (41 loc) · 1.14 KB
/
675 2 to omega of n.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
#!/usr/bin/julia
# Daniel "Trizen" Șuteu
# Date: 27 June 2019
# https://github.com/trizen
# Based on the identity:
# Sum_{d|n} 2^omega(d) = sigma_0(n^2)
# The algorithm iteraters over each number in k = 1..N,
# and computes sigma_0(k^2) = Prod_{p^e | k} (2*e + 1).
# By keeping track of the partial products, we find sigma_0(k!^2).
# Runtime: 17.373s
# https://projecteuler.net/problem=675
using Primes
using Printf
function F(N::Int64, MOD::Int64)
table = Dict{Int64,Int64}()
total = 0
partial = 1
function mulmod(a,b)
(a * b) % MOD
end
for k in 2:N
old = 1 # old product
for (p,e) in factor(k)
if haskey(table, p)
old = mulmod(old, 2*table[p] + 1)
else
table[p] = 0
end
table[p] += e
partial = mulmod(partial, 2*table[p] + 1)
end
partial = mulmod(partial, invmod(old, MOD)) # remove the old product
total += partial
total %= MOD
end
return total
end
const MOD = 1000000087
for n in 1:7
@printf("F(10^%d) = %10d (mod %d)\n", n, F(10^n, MOD), MOD)
end