123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110 |
- #!/usr/bin/ruby
- # Author: Daniel "Trizen" Șuteu
- # Date: 20 November 2018
- # Edit: 06 June 2021
- # https://github.com/trizen
- # A new algorithm for computing the partial-sums of `ϕ(k)`, for `1 <= k <= n`:
- #
- # Sum_{k=1..n} phi(k)
- #
- # where phi(k) is the Euler totient function.
- # Based on the formula:
- # Sum_{k=1..n} phi(k) = (1/2)*Sum_{k=1..n} moebius(k) * floor(n/k) * floor(1+n/k)
- # Example:
- # a(10^1) = 32
- # a(10^2) = 3044
- # a(10^3) = 304192
- # a(10^4) = 30397486
- # a(10^5) = 3039650754
- # a(10^6) = 303963552392
- # a(10^7) = 30396356427242
- # a(10^8) = 3039635516365908
- # a(10^9) = 303963551173008414
- # This algorithm can be improved.
- # See also:
- # https://oeis.org/A002088
- # https://oeis.org/A064018
- # https://en.wikipedia.org/wiki/Mertens_function
- # https://en.wikipedia.org/wiki/M%C3%B6bius_function
- # https://en.wikipedia.org/wiki/Euler%27s_totient_function
- # https://en.wikipedia.org/wiki/Dirichlet_hyperbola_method
- # https://trizenx.blogspot.com/2018/08/interesting-formulas-and-exercises-in.html
- func euler_totient_partial_sum (n) {
- var total = 0
- var s = n.isqrt
- var u = idiv(n, s+1)
- var prev = mertens(n)
- for k in (1..s) {
- with (mertens(idiv(n, k+1))) {|curr|
- total += ((prev - curr) * faulhaber(k, 1))
- prev = curr
- }
- }
- u.each_squarefree {|k|
- total += (moebius(k) * faulhaber(idiv(n,k), 1))
- }
- return total
- }
- func euler_totient_partial_sum_2 (n) { # using Dirichlet's hyperbola method
- var total = 0
- var s = n.isqrt
- for k in (1..s) {
- total += k*mertens(idiv(n,k))
- }
- s.each_squarefree {|k|
- total += moebius(k)*faulhaber(idiv(n,k), 1)
- }
- total -= mertens(s)*faulhaber(s, 1)
- return total
- }
- func euler_totient_partial_sum_test (n) { # just for testing
- 1..n -> sum { .euler_phi }
- }
- for m in (0 .. 10) {
- var n = 10000.irand
- var t1 = euler_totient_partial_sum(n)
- var t2 = euler_totient_partial_sum_2(n)
- var t3 = euler_totient_partial_sum_test(n)
- assert_eq(t1, t2)
- assert_eq(t1, t3)
- say "Sum_{k=1..#{n}} phi(k) = #{t1}"
- }
- __END__
- Sum_{k=1..6000} phi(k) = 10943164
- Sum_{k=1..9647} phi(k) = 28292296
- Sum_{k=1..1767} phi(k) = 949684
- Sum_{k=1..4643} phi(k) = 6554816
- Sum_{k=1..5814} phi(k) = 10276132
- Sum_{k=1..3120} phi(k) = 2959328
- Sum_{k=1..1998} phi(k) = 1213790
- Sum_{k=1..1269} phi(k) = 489854
- Sum_{k=1..1298} phi(k) = 512500
- Sum_{k=1..6555} phi(k) = 13062444
- Sum_{k=1..1101} phi(k) = 368850
|