modified: Math/nth_composite.sf

modified:   Math/nth_prime_power.sf
modified:   Math/nth_squarefree.sf

Author: trizen

Math/nth_composite.sf

@@ -10,25 +10,53 @@
 # 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))
+    #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 = bsearch_le(min, max, {|k|
-        #~ (k - prime_count(k) - 1) <=> n
-    #~ })
-
     var k = 0
     var c = 0
 
@@ -88,3 +116,7 @@ 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

Math/nth_prime_power.sf

@@ -76,11 +76,11 @@ func nth_prime_power(n) {
 }
 
 for n in (1..10) {
-    var p = nth_prime_power(10**n)
-    assert(p.is_prime_power)
-    assert_eq(p.prime_power_count, 10**n)
-    #assert_eq(10**n -> nth_prime_power, p)
-    say "P(10^#{n}) = #{p}"
+    var pp = nth_prime_power(10**n)
+    assert(pp.is_prime_power)
+    assert_eq(pp.prime_power_count, 10**n)
+    assert_eq(10**n -> nth_prime_power, pp)
+    say "PP(10^#{n}) = #{pp}"
 }
 
 assert_eq(
@@ -89,15 +89,18 @@ assert_eq(
 )
 
 __END__
-P(10^1) = 16
-P(10^2) = 419
-P(10^3) = 7517
-P(10^4) = 103511
-P(10^5) = 1295953
-P(10^6) = 15474787
-P(10^7) = 179390821
-P(10^8) = 2037968761
-P(10^9) = 22801415981
-P(10^10) = 252096675073
-P(10^11) = 2760723662941
-P(10^12) = 29996212395727
+PP(10^1) = 16
+PP(10^2) = 419
+PP(10^3) = 7517
+PP(10^4) = 103511
+PP(10^5) = 1295953
+PP(10^6) = 15474787
+PP(10^7) = 179390821
+PP(10^8) = 2037968761
+PP(10^9) = 22801415981
+PP(10^10) = 252096675073
+PP(10^11) = 2760723662941
+PP(10^12) = 29996212395727
+PP(10^13) = 323780470283789
+PP(10^14) = 3475385632514321
+PP(10^15) = 37124507635789309

Math/nth_squarefree.sf

@@ -94,3 +94,7 @@ 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