123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101 |
- #!/usr/bin/ruby
- # Daniel "Trizen" Șuteu
- # Date: 07 July 2022
- # https://github.com/trizen
- # Fast algorithm for computing the n-th squarefree number.
- # See also:
- # https://oeis.org/A005117
- # PARI/GP program:
- # S(n) = my(s); forsquarefree(k=1,sqrtint(n),s+=n\k[1]^2*moebius(k)); s;
- # a(n) = my(min=1, max=231, k=0, sc=0); if(n >= 144, min=floor(zeta(2)*n - 5*sqrt(n)); max=ceil(zeta(2)*n + 5*sqrt(n))); while(min <= max, k=(min+max)\2; sc=S(k); if(abs(sc-n) <= sqrtint(n), break); if(sc > n, max=k-1, if(sc < n, min=k+1, break))); while(!issquarefree(k), k-=1); while(sc != n, my(j=1); if(sc > n, j=-1); k += j; sc += j; while(!issquarefree(k), k+=j)); k;
- func nth_squarefree(n) {
- n == 0 && return 0 # not squarefree, but...
- n <= 0 && return NaN
- n == 1 && return 1
- var min = 1
- var max = 231
- # Bounds on squarefree numbers:
- # https://mathoverflow.net/questions/66701/bounds-on-squarefree-numbers
- if (n >= 268293) {
- min = int(zeta(2)*n - 0.058377*sqrt(n))
- max = int(zeta(2)*n + 0.058377*sqrt(n))
- }
- elsif (n >= 144) {
- min = int(zeta(2)*n - 5*sqrt(n))
- max = int(zeta(2)*n + 5*sqrt(n))
- }
- var k = 0
- var c = 0
- loop {
- k = (min + max)>>1
- c = k.squarefree_count
- if (abs(c - n) <= k.isqrt) {
- break
- }
- given (c <=> n) {
- when (+1) { max = k-1 }
- when (-1) { min = k+1 }
- else { break }
- }
- }
- while (!is_squarefree(k)) {
- --k
- }
- while (c != n) {
- var j = (n <=> c)
- k += j
- c += j
- k += j while !k.is_squarefree
- }
- return k
- }
- for n in (1..10) {
- var s = nth_squarefree(10**n)
- assert(s.is_squarefree)
- assert_eq(s.squarefree_count, 10**n)
- assert_eq(10**n -> nth_squarefree, s)
- say "S(10^#{n}) = #{s}"
- }
- assert_eq(
- nth_squarefree.map(1..100),
- 100.by { .is_squarefree },
- )
- __END__
- S(10^1) = 14
- S(10^2) = 163
- S(10^3) = 1637
- S(10^4) = 16446
- S(10^5) = 164498
- S(10^6) = 1644918
- S(10^7) = 16449369
- S(10^8) = 164493390
- S(10^9) = 1644934081
- S(10^10) = 16449340709
- S(10^11) = 164493406178
- S(10^12) = 1644934067511
- S(10^13) = 16449340668746
- S(10^14) = 164493406685659
- S(10^15) = 1644934066850410
- S(10^16) = 16449340668485215
- S(10^17) = 164493406684817902
- S(10^18) = 1644934066848209910
|