12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273 |
- #!/usr/bin/ruby
- # Let f(n) be the number of couples (x,y) with x and y positive integers, x ≤ y and the least common multiple of x and y equal to n.
- # Let a(n) = A007875(n), with a(1) = 1, for n > 1 (due to Vladeta Jovovic, Jan 25 2002):
- # a(n) = (1/2)*Sum_{d|n} abs(mu(d))
- # = 2^(omega(n)-1)
- # = usigma_0(n)/2
- # This gives us f(n) as:
- # f(n) = Sum_{d|n} a(d)
- # This script implements a sub-linear formula for computing partial sums of f(n):
- # S(n) = Sum_{k=1..n} f(k)
- # = Sum_{k=1..n} Sum_{d|k} a(d)
- # = Sum_{k=1..n} a(k) * floor(n/k)
- # See also:
- # https://oeis.org/A007875
- # https://oeis.org/A064608
- # https://oeis.org/A182082
- # Problem from:
- # https://projecteuler.net/problem=379
- # Several values for S(10^n):
- # S(10^1) = 29
- # S(10^2) = 647
- # S(10^3) = 11751
- # S(10^4) = 186991
- # S(10^5) = 2725630
- # S(10^6) = 37429395
- # S(10^7) = 492143953
- # S(10^8) = 6261116500
- # S(10^9) = 77619512018
- # S(10^10) = 942394656385
- # S(10^11) = 11247100884096
- func usigma0_partial_sum (n) { # O(sqrt(n)) complexity
- var total = 0
- var s = n.isqrt
- var u = idiv(n, s+1)
- var prev = squarefree_count(n)
- for k in (1..s) {
- var curr = squarefree_count(idiv(n, k+1))
- total += (prev - curr)*k
- prev = curr
- }
- u.each_squarefree {|k|
- total += idiv(n, k)
- }
- return total
- }
- func S(n) {
- dirichlet_sum(n,
- {|k| (k == 1) ? 1 : usigma0(k)/2 }, # a(n)
- { 1 },
- {|k| (usigma0_partial_sum(k)-1)/2 }, # partial sums of a(n)
- { _ }
- )
- }
- for k in (1..5) {
- say "S(10^#{k}) = #{S(10**k)}"
- }
|