123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293 |
- #!/usr/bin/ruby
- # Daniel "Trizen" Șuteu
- # Date: 10 January 2019
- # https://github.com/trizen
- # Fast formula for computing:
- # a(n) = sum of the number of divisors of gcd(x,y) with x*y <= n.
- # 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.
- # Adrian Dudek, on the Success of Mishandling Euclid's Lemma:
- # https://arxiv.org/abs/1602.03555
- # Asymptotic formula:
- # a(n) ~ 1/6 * π^2 * n * (2 * (-12 * log(A) + γ + log(2) + log(π)) + log(n) + 2*γ - 1) + O(sqrt(n)*log(n))
- #
- # where γ is the Euler-Mascheroni constant and "A" is the Glaisher-Kinkelin constant.
- # Alternative asymptotic formula:
- # a(n) ~ (n*zeta(2) * (log(n) + 2*γ - 1 + c)) + O(sqrt(n)*log(n))
- #
- # where γ is the Euler-Mascheroni 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 asymptotic_formula2(n) {
- # The Glaisher-Kinkelin constant
- const A = 1.28242712910062263687534256886979172776768892732500119206374002174040630885882646112973649195820237439420646120399
- # Asymptotic formula in terms of the Glaisher-Kinkelin constant
- zeta(2) * n * (2 * (-12 * log(A) + Num.EulerGamma + log(Num.tau)) + log(n) + 2*Num.EulerGamma - 1)
- }
- func sum_of_number_of_divisors_of_gcd (n) {
- var total = 0
- for k in (1..n.isqrt) {
- var t = 2*sum(1..isqrt(n/(k*k)), {|j|
- int(n / (j*k*k))
- })
- total += (t - isqrt(n/(k*k))**2)
- }
- return total
- }
- say 20.of(sum_of_number_of_divisors_of_gcd)
- for k in (1..6) {
- var n = 10**k
- var t = sum_of_number_of_divisors_of_gcd(n)
- var u = asymptotic_formula(n)
- printf("a(10^%s) = %10s ~ %-15s -> %s\n", k, t, u.round(-2), t/u)
- assert((t - u) < (sqrt(n) * log(n))) if (n > 100)
- }
- __END__
- [0, 1, 3, 5, 9, 11, 15, 17, 23, 27, 31, 33, 41, 43, 47, 51, 60, 62, 70, 72]
- a(10^1) = 31 ~ 21.66 -> 1.43085716814724731567697388362262512087796085132
- a(10^2) = 629 ~ 595.41 -> 1.05640884486870073219427770179934635115325018838
- a(10^3) = 9823 ~ 9741.73 -> 1.00834196073027036381381492602216565392721426965
- a(10^4) = 135568 ~ 135293.35 -> 1.00202999682691312076763529313619057317755443274
- a(10^5) = 1732437 ~ 1731693.62 -> 1.00042928187585922855456102384841804396816671626
- a(10^6) = 21107131 ~ 21104536.81 -> 1.00012292075282086768302929969619689628662614091
- a(10^7) = 248928748 ~ 248921374.75 -> 1.00002962076965424120327637576433900637540389794
- a(10^8) = 2867996696 ~ 2867973813.70 -> 1.00000797855916535575678041071686222727851109258
- a(10^9) = 32467409097 ~ 32467338798.29 -> 1.00000216521302261846703873643427029189711363986
- a(10^10) = 362549612240 ~ 362549394595.78 -> 1.00000060031604804834071744691960444929352043683
- a(10^11) = 4004254692640 ~ 4004254012086.08 -> 1.00000016995772897612356184672572401706556174343
- a(10^12) = 43830142939380 ~ 43830140782143.61 -> 1.00000004921810301432497497420018745129768545890
- a(10^13) = 476177421208658 ~ 476177414434264.13 -> 1.00000001422661735697513455710167585383336332082
- a(10^14) = 5140534231877816 ~ 5140534210470921.03 -> 1.00000000416433275074901946766616776434113033877
- a(10^15) = 55192942833495679 ~ 55192942765992007.53 -> 1.00000000122304896383936291361582837193299642341
|