sidef-scripts/Math/sum_of_number_of_unitary_divisors.sf

#!/usr/bin/ruby

# Daniel "Trizen" Șuteu
# Date: 10 January 2019
# https://github.com/trizen

# Two fast algorithms for computing the sum of number of unitary divisors from 1 to n.
#   a(n) = Sum_{k=1..n} usigma_0(k)

# Based on the formula:
#   a(n) = Sum_{k=1..n} moebius(k)^2 * floor(n/k)

# See also:
#   https://oeis.org/A034444    -- Partial sums of A034444: sum of number of unitary divisors from 1 to n.
#   https://oeis.org/A180361    -- Sum of number of unitary divisors (A034444) from 1 to 10^n
#   https://oeis.org/A268732    -- Sum of the numbers of divisors of gcd(x,y) with x*y <= n.

# Asymptotic formula:
#   a(n) ~ n*log(n)/zeta(2) + O(n)

# Better asymptotic formula:
#   a(n) ~ (n/zeta(2))*(log(n) + 2*γ - 1 - c) + O(sqrt(n) * log(n))
#
# where γ is the Euler-Mascheroni constant and c = 2*zeta'(2)/zeta(2) = -1.1399219861890656127997287200...

func asymptotic_formula(n) {

    # c = 2*Zeta'(2)/Zeta(2) = (12 * Zeta'(2))/π^2 = 2 (-12 log(A) + γ + log(2) + log(π))
    const c = -1.13992198618906561279972872003946000480696456161386195911639472087583455473348121357

    # Asymptotic formula based on Merten's theorem (1874) (see: https://oeis.org/A064608)
    (n/zeta(2)) * (log(n) + 2*Num.EulerGamma - 1 - c)
}

func usigma0_partial_sum_1 (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 usigma0_partial_sum_2 (n) {      # based on formula by Jerome Raulin (https://oeis.org/A064608)

    var total = 0

    n.isqrt.each_squarefree {|k|

        var v = moebius(k)
        var t = 2*sum(1..isqrt(idiv(n, k*k)), {|j|
            idiv(n, j*k*k)
        })

        total += v*(t - isqrt(idiv(n, k*k))**2)
    }

    return total
}

func usigma0_partial_sum_3(n) {     # O(sqrt(n)) complexity, using Dirichlet's hyperbola method
    n.dirichlet_sum({1}, {.mu.abs}, {_}, {.squarefree_count})
}

say 20.of(usigma0_partial_sum_1)
say 20.of(usigma0_partial_sum_2)
say 20.of(usigma0_partial_sum_3)

for k in (0..7) {

    var n = 10**k
    var t = usigma0_partial_sum_1(n)
    var u = asymptotic_formula(n)

    printf("a(10^%s) = %10s ~ %-15s -> %s\n", k, t, u.round(-2), t/u)
}

__END__
[0, 1, 3, 5, 7, 9, 13, 15, 17, 19, 23, 25, 29, 31, 35, 39, 41, 43, 47, 49]
[0, 1, 3, 5, 7, 9, 13, 15, 17, 19, 23, 25, 29, 31, 35, 39, 41, 43, 47, 49]

a(10^0)  =            1 ~ 0.79            -> 1.27085398285349342897812915198984638968899591751
a(10^1)  =           23 ~ 21.87           -> 1.05182461403816051734935994402113331145060974294
a(10^2)  =          359 ~ 358.65          -> 1.00098140095602073835866744824992972185806123685
a(10^3)  =         4987 ~ 4986.28         -> 1.00014357239778054254970740667091143421188177813
a(10^4)  =        63869 ~ 63860.88        -> 1.00012715302552355451250212258735392366329621935
a(10^5)  =       778581 ~ 778589.19       -> 0.999989484576929013867264739526374966823956960403
a(10^6)  =      9185685 ~ 9185695.75      -> 0.99999882923368455522780513812504287278271814501
a(10^7)  =    105854997 ~ 105854996.37    -> 1.00000000598372061072117962943109677794267023891
a(10^8)  =   1198530315 ~ 1198530351.90   -> 0.999999969211002320383540850995519903094748492418
a(10^9)  =  13385107495 ~ 13385107401.37  -> 1.00000000699496540213133746406895764726726792391
a(10^10) = 147849112851 ~ 147849112837.28 -> 1.00000000009281141854332921757852421030396550125