Greatest prime dividing the n-th cubefree number: Difference between revisions
Greatest prime dividing the n-th cubefree number (view source)
Revision as of 12:44, 6 March 2024
, 3 months ago→{{header|Wren}}: Added second far quicker version (translation of Phix).
(→{{header|Wren}}: Added a simple optimization & rewrote preamble.) |
(→{{header|Wren}}: Added second far quicker version (translation of Phix).) |
||
Line 771:
{{libheader|Wren-math}}
{{libheader|Wren-fmt}}
===Version 1===
Simple brute force approach though skipping numbers which are multiples of 8 or 27 as these can't be cubefree.
Line 835 ⟶ 837:
The 10,000,000th term of a[n] is 1,202,057.
</pre>
===Version 2===
{{trans|Phix}}
Astonishing speed-up, now taking only 0.04 seconds to reach 10 million and 38 seconds to reach 10 billion.
<syntaxhighlight lang="wren">import "./math" for Int
import "./fmt" for Conv, Fmt
import "./iterate" for Stepped
var primes = Int.primeSieve(15 * 1e6)
var countBits = Fn.new { |n| Conv.bin(n).count { |c| c == "1" } }
var cubesBefore = Fn.new { |n|
var r = 0
var xpm = true
var pm = []
for (i in 1..n) {
var p3 = primes[i-1].pow(3)
if (p3 > n) break
var k = (n/p3).floor
for (mask in Stepped.ascend(1..2.pow(pm.count)-1)) {
var m = mask
var mi = 0
var kpm = p3
while (m != 0) {
mi = mi + 1
if (m % 2 == 1) kpm = kpm * pm[mi-1]
m = (m/2).floor
}
if (kpm > n) {
if (countBits.call(mask) == 1) {
xpm = false
pm = pm[0...mi-1]
break
}
} else {
var l = (n/kpm).floor
if (countBits.call(mask) % 2 == 1) {
k = k - l
} else {
k = k + l
}
}
}
r = r + k
if (xpm) pm.add(p3)
}
return r
}
var cubeFree = Fn.new { |nth|
var lo = nth
var hi = lo * 2
var mid
var cb
var k
while (hi - cubesBefore.call(hi) < nth) {
lo = hi
hi = lo * 2
}
while (true) {
mid = ((lo + hi)/2).floor
cb = cubesBefore.call(mid)
k = mid - cb
if (k == nth) {
while (cubesBefore.call(mid-1) != cb) {
mid = mid - 1
cb = cb - 1
}
break
} else if (k < nth) {
lo = mid
} else {
hi = mid
}
}
return mid
}
var a370833 = Fn.new { |n|
if (n == 1) return 1
var nth = cubeFree.call(n)
return Int.primeFactors(nth)[-1]
}
var start = System.clock
System.print("First 100 terms of a[n]:")
Fmt.tprint("$3d", (1..100).map { |i| a370833.call(i) }, 10)
System.print()
var n = 1000
while (n <= 1e10) {
Fmt.print("The $,r term of a[n] is $,d.", n, a370833.call(n))
n = n * 10
}
System.print("\nTook %(System.clock - start) seconds.")></syntaxhighlight>
{{out}}
<pre>
First 100 terms of a[n]:
1 2 3 2 5 3 7 3 5 11
3 13 7 5 17 3 19 5 7 11
23 5 13 7 29 5 31 11 17 7
3 37 19 13 41 7 43 11 5 23
47 7 5 17 13 53 11 19 29 59
5 61 31 7 13 11 67 17 23 7
71 73 37 5 19 11 13 79 41 83
7 17 43 29 89 5 13 23 31 47
19 97 7 11 5 101 17 103 7 53
107 109 11 37 113 19 23 29 13 59
The 1,000th term of a[n] is 109.
The 10,000th term of a[n] is 101.
The 100,000th term of a[n] is 1,693.
The 1,000,000th term of a[n] is 1,202,057.
The 10,000,000th term of a[n] is 1,202,057.
The 100,000,000th term of a[n] is 20,743.
The 1,000,000,000th term of a[n] is 215,461.
The 10,000,000,000th term of a[n] is 1,322,977.
Took 37.888876 seconds.
</pre>
| |||