123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123 |
- #!/usr/bin/ruby
- # Daniel "Trizen" Șuteu
- # Date: 06 June 2021
- # Edit: 07 July 2022
- # https://github.com/trizen
- # Return the n-th composite number, using binary search and the prime counting function.
- # See also:
- # https://oeis.org/A002808
- func composite_count_lower(n) {
- n - n.prime_count_upper - 1
- }
- func composite_count_upper(n) {
- n - n.prime_count_lower - 1
- }
- func simple_composite_lower(n) {
- int(n + n/log(n) + n/(log(n)**2))
- }
- func simple_composite_upper(n) {
- int(n + n/log(n) + (3*n)/(log(n)**2))
- }
- func nth_composite_lower(n) {
- bsearch_min(simple_composite_lower(n), simple_composite_upper(n), {|k|
- composite_count_upper(k) <=> n
- })
- }
- func nth_composite_upper(n) {
- bsearch_max(simple_composite_lower(n), simple_composite_upper(n), {|k|
- composite_count_lower(k) <=> n
- })
- }
- func nth_composite(n) {
- n == 0 && return 1 # not composite, but...
- n <= 0 && return NaN
- n == 1 && return 4
- # Lower and upper bounds from A002808 (for n >= 4).
- #var min = int(n + n/log(n) + n/(log(n)**2))
- #var max = int(n + n/log(n) + (3*n)/(log(n)**2))
- # Better bounds for the n-th composite number
- var min = nth_composite_lower(n)
- var max = nth_composite_upper(n)
- if (n < 4) {
- min = 4;
- max = 8;
- }
- var k = 0
- var c = 0
- loop {
- k = (min + max)>>1
- c = (k - prime_count(k) - 1)
- if (abs(c - n) <= k.iroot(3)) {
- break
- }
- given (c <=> n) {
- when (+1) { max = k-1 }
- when (-1) { min = k+1 }
- else { break }
- }
- }
- if (!k.is_composite) {
- --k
- }
- while (c != n) {
- var j = (n <=> c)
- k += j
- c += j
- k += j while !k.is_composite
- }
- return k
- }
- for n in (1..10) {
- var c = nth_composite(10**n)
- assert(c.is_composite)
- assert_eq(c.composite_count, 10**n)
- assert_eq(10**n -> nth_composite, c)
- say "C(10^#{n}) = #{c}"
- }
- assert_eq(
- nth_composite.map(1..100),
- 100.by { .is_composite }
- )
- __END__
- C(10^1) = 18
- C(10^2) = 133
- C(10^3) = 1197
- C(10^4) = 11374
- C(10^5) = 110487
- C(10^6) = 1084605
- C(10^7) = 10708555
- C(10^8) = 106091745
- C(10^9) = 1053422339
- C(10^10) = 10475688327
- C(10^11) = 104287176419
- C(10^12) = 1039019056246
- C(10^13) = 10358018863853
- C(10^14) = 103307491450820
- C(10^15) = 1030734020030318
- C(10^16) = 10287026204717358
|