Numbers which are the cube roots of the product of their proper divisors: Difference between revisions
Content added Content deleted
(added julia entry) |
(Added Lua version) |
||
| Line 771: | Line 771: | ||
50000th: 223735 |
50000th: 223735 |
||
</pre> |
</pre> |
||
=={{header|Lua}}== |
|||
The OEIS page gives a formula of "1 together with numbers with 8 divisors", so that's what we test. |
|||
<syntaxhighlight lang="lua">function is_1_or_has_eight_divisors (n) |
|||
if n == 1 then return true end |
|||
local divCount, sqr = 2, math.sqrt(n) |
|||
for d = 2, sqr do |
|||
if n % d == 0 then |
|||
divCount = d == sqr and divCount + 1 or divCount + 2 |
|||
end |
|||
if divCount > 8 then return false end |
|||
end |
|||
return divCount == 8 |
|||
end |
|||
-- First 50 |
|||
local count, x = 1, 0 |
|||
while count <= 50 do |
|||
x = x + 1 |
|||
if is_1_or_has_eight_divisors(x) then |
|||
io.write(x .. " ") |
|||
count = count + 1 |
|||
end |
|||
end |
|||
-- 500th, 5,000th and 50,000th |
|||
while count <= 50000 do |
|||
x = x + 1 |
|||
if is_1_or_has_eight_divisors(x) then |
|||
if count == 500 then print("\n\n500th: " .. x) end |
|||
if count == 5000 then print("5,000th: " .. x) end |
|||
count = count + 1 |
|||
end |
|||
end |
|||
print("50,000th: " .. x)</syntaxhighlight> |
|||
{{out}} |
|||
<pre>1 24 30 40 42 54 56 66 70 78 88 102 104 105 110 114 128 130 135 136 138 152 154 165 170 174 182 184 186 189 190 195 222 230 231 232 238 246 248 250 255 258 266 273 282 285 286 290 296 297 |
|||
500th: 2526 |
|||
5,000th: 23118 |
|||
50,000th: 223735</pre> |
|||
=={{header|Pascal}}== |
=={{header|Pascal}}== |
||
Revision as of 10:46, 6 April 2023
You are encouraged to solve this task according to the task description, using any language you may know.
- Example
Consider the number 24. Its proper divisors are: 1, 2, 3, 4, 6, 8 and 12. Their product is 13,824 and the cube root of this is 24. So 24 satisfies the definition in the task title.
- Task
Compute and show here the first 50 positive integers which are the cube roots of the product of their proper divisors.
Also show the 500th and 5,000th such numbers.
- Stretch
Compute and show the 50,000th such number.
- Reference
- Note
OEIS considers 1 to be the first number in this sequence even though, strictly speaking, it has no proper divisors. Please therefore do likewise.
11l
F product_of_proper_divisors(n)
V prod = Int64(1)
L(d) 2 .< Int(sqrt(n) + 1)
I n % d == 0
prod *= d
V otherD = n I/ d
I otherD != d
prod *= otherD
R prod
print(‘First 50 numbers which are the cube roots of the products of their proper divisors:’)
V found = 0
L(num) 1..
I Int64(num) ^ 3 == product_of_proper_divisors(num)
found++
I found <= 50
print(f:‘{num:3}’, end' I found % 10 == 0 {"\n"} E ‘ ’)
E I found C (500, 5000, 50000)
print(f:‘{commatize(found):6}th: {commatize(num)}’)
I found == 50000
L.break- Output:
First 50 numbers which are the cube roots of the products of their proper divisors: 1 24 30 40 42 54 56 66 70 78 88 102 104 105 110 114 128 130 135 136 138 152 154 165 170 174 182 184 186 189 190 195 222 230 231 232 238 246 248 250 255 258 266 273 282 285 286 290 296 297 500th: 2,526 5,000th: 23,118 50,000th: 223,735
ALGOL 68
As with the second Wren sample, uses the observation on the OEIS page to reduce the problem to finding numbers that are 1 or have 8 divisors (or 7 proper divisors).
BEGIN # find some numbers which are the cube roots of the product of their #
# proper divisors #
# the Online Encyclopedia of Integer Sequences states that these #
# numbers are 1 and those with eight divisors #
# NB: numbers with 8 divisors have 7 proper divisors #
INT max number = 500 000; # maximum number we will consider #
# form a table of proper divisor counts - assume the pdc of 1 is 7 #
[ 1 : max number ]INT pdc; FOR i TO UPB pdc DO pdc[ i ] := 1 OD;
pdc[ 1 ] := 7;
FOR i FROM 2 TO UPB pdc DO
FOR j FROM i + i BY i TO UPB pdc DO pdc[ j ] +:= 1 OD
OD;
# show the numbers which are the cube root of their proper divisor #
# product - equivalent to finding the numbers with a proper divisor #
# count of 7 ( we have "cheated" and set the pdc of 1 to 7 ) #
INT next show := 500;
INT c count := 0;
FOR n TO UPB pdc DO
IF pdc[ n ] = 7 THEN
# found a suitable number #
IF ( c count +:= 1 ) <= 50 THEN
print( ( " ", whole( n, -3 ) ) );
IF c count MOD 10 = 0 THEN print( ( newline ) ) FI
ELIF c count = next show THEN
print( ( whole( c count, -9 ), "th: ", whole( n, 0 ), newline ) );
next show *:= 10
FI
FI
OD
END- Output:
1 24 30 40 42 54 56 66 70 78
88 102 104 105 110 114 128 130 135 136
138 152 154 165 170 174 182 184 186 189
190 195 222 230 231 232 238 246 248 250
255 258 266 273 282 285 286 290 296 297
500th: 2526
5000th: 23118
50000th: 223735
AppleScript
Like other solutions here, this checks for numbers having seven proper divisors rather than doing the multiplications, which saves time and avoids products that are too large for AppleScript numbers.
on properDivisors(n)
set output to {}
if (n > 1) then
set sqrt to n ^ 0.5
set limit to sqrt div 1
if (limit = sqrt) then
set end of output to limit
set limit to limit - 1
end if
repeat with i from limit to 2 by -1
if (n mod i is 0) then
set beginning of output to i
set end of output to n div i
end if
end repeat
set beginning of output to 1
end if
return output
end properDivisors
on join(lst, delim)
set astid to AppleScript's text item delimiters
set AppleScript's text item delimiters to delim
set txt to lst as text
set AppleScript's text item delimiters to astid
return txt
end join
on task()
set output to {"First 50 numbers whose cubes are the products of their proper divisors", ¬
"(and of course whose fourth powers are the products of ALL their positive divisors):"}
set pad to " "
set n to 1
set first50 to {" 1"}
repeat 49 times
set n to n + 1
repeat until ((count properDivisors(n)) = 7)
set n to n + 1
end repeat
set end of first50 to text -5 thru -1 of (pad & n)
end repeat
repeat with i from 1 to 41 by 10
set end of output to join(first50's items i thru (i + 9), "")
end repeat
set |count| to 50
repeat with target in {500, 5000, 50000}
repeat with |count| from (|count| + 1) to target
set n to n + 1
repeat until ((count properDivisors(n)) = 7)
set n to n + 1
end repeat
end repeat
set end of output to text -6 thru -1 of (pad & |count|) & "th: " & text -6 thru -1 of (pad & n)
end repeat
return join(output, linefeed)
end task
task()
- Output:
"First 50 numbers whose cubes are the products of their proper divisors
(and of course whose fourth powers are the products of ALL their positive divisors):
1 24 30 40 42 54 56 66 70 78
88 102 104 105 110 114 128 130 135 136
138 152 154 165 170 174 182 184 186 189
190 195 222 230 231 232 238 246 248 250
255 258 266 273 282 285 286 290 296 297
500th: 2526
5000th: 23118
50000th: 223735"
BASIC
BASIC256
arraybase 1
limite = 500000
dim pdc(limite) fill 1
pdc[1] = 7
for i = 2 to pdc[?]
for j = i + i to pdc[?] step i
pdc[j] += 1
next j
next i
n5 = 500
cont = 0
print "First 50 numbers which are the cube roots"
print "of the products of their proper divisors:"
for i = 1 to pdc[?]
if pdc[i] = 7 then
cont += 1
if cont <= 50 then
print rjust(string(i),5);
if cont mod 10 = 0 then print
else
if cont = n5 then
print
print rjust(string(cont),9); "th: "; i;
n5 *= 10
end if
end if
end if
next i- Output:
Same as FreeBASIC entry.
True BASIC
LET limite = 500000
DIM pdc(1 to 500000)
FOR i = 1 to ubound(pdc)
LET pdc(i) = 1
NEXT i
LET pdc(1) = 7
FOR i = 2 to ubound(pdc)
FOR j = i+i to ubound(pdc) step i
LET pdc(j) = pdc(j)+1
NEXT j
NEXT i
LET n5 = 500
LET count = 0
PRINT "First 50 numbers which are the cube roots"
PRINT "of the products of their proper divisors:"
FOR i = 1 to ubound(pdc)
IF pdc(i) = 7 then
LET count = count + 1
IF count <= 50 THEN
PRINT using "####": i;
IF remainder(count, 10) = 0 THEN PRINT
ELSE
IF count = n5 THEN
PRINT
PRINT USING "#########th:": count;
PRINT i;
LET n5 = n5*10
END IF
END IF
END IF
NEXT i
END
- Output:
Same as FreeBASIC entry.
XBasic
PROGRAM "progname"
VERSION "0.0000"
DECLARE FUNCTION Entry ()
FUNCTION Entry ()
limite = 500000
DIM pdc[limite] '(1 TO limite)
FOR i = 1 TO UBOUND(pdc[])
pdc[i] = 1
NEXT i
pdc[1] = 7
FOR i = 2 TO UBOUND(pdc[])
FOR j = i + i TO UBOUND(pdc[]) STEP i
INC pdc[j]
NEXT j
NEXT i
n5 = 500
cont = 0
PRINT "First 50 numbers which are the cube roots"
PRINT "of the products of their proper divisors:"
FOR i = 1 TO UBOUND(pdc[])
IF pdc[i] = 7 THEN
INC cont
IF cont <= 50 THEN
PRINT RJUST$ (STRING$(i), 4);
IF cont MOD 10 = 0 THEN PRINT
ELSE
IF cont = n5 THEN
PRINT "\n"; FORMAT$("#########", cont); "th: "; i;
n5 = n5 * 10
END IF
END IF
END IF
NEXT i
END FUNCTION
END PROGRAM- Output:
Same as FreeBASIC entry.
Yabasic
limite = 500000
dim pdc(limite)
for i = 1 to arraysize(pdc(), 1)
pdc(i) = 1
next i
pdc(1) = 7
for i = 2 to arraysize(pdc(), 1)
for j = i + i to arraysize(pdc(), 1) step i
pdc(j) = pdc(j) + 1
next j
next i
n5 = 500
cont = 0
print "First 50 numbers which are the cube roots"
print "of the products of their proper divisors:"
for i = 1 to arraysize(pdc(), 1)
if pdc(i) = 7 then
cont = cont + 1
if cont <= 50 then
print i using("###");
if mod(cont, 10) = 0 print
else
if cont = n5 then
print "\n", cont using("#########"), "th: ", i;
n5 = n5 * 10
end if
end if
end if
next i- Output:
Same as FreeBASIC entry.
C
The faster version.
#include <stdio.h>
#include <locale.h>
int divisorCount(int n) {
int i, j, count = 0, k = !(n%2) ? 1 : 2;
for (i = 1; i*i <= n; i += k) {
if (!(n%i)) {
++count;
j = n/i;
if (j != i) ++count;
}
}
return count;
}
int main() {
int i, numbers50[50], count = 0, n = 1, dc;
printf("First 50 numbers which are the cube roots of the products of their proper divisors:\n");
setlocale(LC_NUMERIC, "");
while (1) {
dc = divisorCount(n);
if (n == 1|| dc == 8) {
++count;
if (count <= 50) {
numbers50[count-1] = n;
if (count == 50) {
for (i = 0; i < 50; ++i) {
printf("%3d ", numbers50[i]);
if (!((i+1) % 10)) printf("\n");
}
}
} else if (count == 500) {
printf("\n500th : %'d\n", n);
} else if (count == 5000) {
printf("5,000th : %'d\n", n);
} else if (count == 50000) {
printf("50,000th: %'d\n", n);
break;
}
}
++n;
}
return 0;
}
- Output:
First 50 numbers which are the cube roots of the products of their proper divisors: 1 24 30 40 42 54 56 66 70 78 88 102 104 105 110 114 128 130 135 136 138 152 154 165 170 174 182 184 186 189 190 195 222 230 231 232 238 246 248 250 255 258 266 273 282 285 286 290 296 297 500th : 2,526 5,000th : 23,118 50,000th: 223,735
C++
#include <iomanip>
#include <iostream>
unsigned int divisor_count(unsigned int n) {
unsigned int total = 1;
for (; (n & 1) == 0; n >>= 1)
++total;
for (unsigned int p = 3; p * p <= n; p += 2) {
unsigned int count = 1;
for (; n % p == 0; n /= p)
++count;
total *= count;
}
if (n > 1)
total *= 2;
return total;
}
int main() {
std::cout.imbue(std::locale(""));
std::cout << "First 50 numbers which are the cube roots of the products of "
"their proper divisors:\n";
for (unsigned int n = 1, count = 0; count < 50000; ++n) {
if (n == 1 || divisor_count(n) == 8) {
++count;
if (count <= 50)
std::cout << std::setw(3) << n
<< (count % 10 == 0 ? '\n' : ' ');
else if (count == 500 || count == 5000 || count == 50000)
std::cout << std::setw(6) << count << "th: " << n << '\n';
}
}
}
- Output:
First 50 numbers which are the cube roots of the products of their proper divisors: 1 24 30 40 42 54 56 66 70 78 88 102 104 105 110 114 128 130 135 136 138 152 154 165 170 174 182 184 186 189 190 195 222 230 231 232 238 246 248 250 255 258 266 273 282 285 286 290 296 297 500th: 2,526 5,000th: 23,118 50,000th: 223,735
Factor
USING: formatting grouping io kernel lists lists.lazy math
prettyprint project-euler.common ;
: A111398 ( -- list )
L{ 1 } 2 lfrom [ tau 8 = ] lfilter lappend-lazy ;
50 A111398 ltake list>array 10 group simple-table. nl
499 4999 49999
[ [ 1 + ] keep A111398 lnth "%5dth: %d\n" printf ] tri@
- Output:
1 24 30 40 42 54 56 66 70 78 88 102 104 105 110 114 128 130 135 136 138 152 154 165 170 174 182 184 186 189 190 195 222 230 231 232 238 246 248 250 255 258 266 273 282 285 286 290 296 297 500th: 2526 5000th: 23118 50000th: 223735
FreeBASIC
Dim As Single limite = 500000
Dim As Integer pdc(1 To limite)
Dim As Integer i, j
For i = 1 To Ubound(pdc)
pdc(i) = 1
Next i
pdc(1) = 7
For i = 2 To Ubound(pdc)
For j = i + i To Ubound(pdc) Step i
pdc(j) += 1
Next j
Next i
Dim As Integer n5 = 500, cont = 0
Print "First 50 numbers which are the cube roots"
Print "of the products of their proper divisors:"
For i = 1 To Ubound(pdc)
If pdc(i) = 7 Then
cont += 1
If cont <= 50 Then
Print Using "####"; i;
If cont Mod 10 = 0 Then Print
Elseif cont = n5 Then
Print Using !"\n#########th: &"; cont; i;
n5 *= 10
End If
End If
Next i
Sleep- Output:
First 50 numbers which are the cube roots
of the products of their proper divisors:
1 24 30 40 42 54 56 66 70 78
88 102 104 105 110 114 128 130 135 136
138 152 154 165 170 174 182 184 186 189
190 195 222 230 231 232 238 246 248 250
255 258 266 273 282 285 286 290 296 297
500th: 2526
5000th: 23118
50000th: 223735
Forth
500000 constant limit
variable pdc limit cells allot
: main
limit 0 do
1 pdc i cells + !
loop
7 pdc !
limit 2 +do
limit i 2* 1- +do
1 pdc i cells + +!
j +loop
loop
." First 50 numbers which are the cube roots" cr
." of the products of their proper divisors:" cr
500 0
limit 0 do
pdc i cells + @ 7 = if
1+
dup 50 <= if
i 1+ 3 .r
dup 10 mod 0= if cr else space then
else
2dup = if
cr over 5 .r ." th: " i 1+ .
swap 10 * swap
then
then
then
loop
2drop cr ;
main
bye
- Output:
First 50 numbers which are the cube roots of the products of their proper divisors: 1 24 30 40 42 54 56 66 70 78 88 102 104 105 110 114 128 130 135 136 138 152 154 165 170 174 182 184 186 189 190 195 222 230 231 232 238 246 248 250 255 258 266 273 282 285 286 290 296 297 500th: 2526 5000th: 23118 50000th: 223735
Go
The faster version.
package main
import (
"fmt"
"math"
"rcu"
)
func divisorCount(n int) int {
k := 1
if n%2 == 1 {
k = 2
}
count := 0
sqrt := int(math.Sqrt(float64(n)))
for i := 1; i <= sqrt; i += k {
if n%i == 0 {
count++
j := n / i
if j != i {
count++
}
}
}
return count
}
func main() {
var numbers50 []int
count := 0
for n := 1; count < 50000; n++ {
dc := divisorCount(n)
if n == 1 || dc == 8 {
count++
if count <= 50 {
numbers50 = append(numbers50, n)
if count == 50 {
rcu.PrintTable(numbers50, 10, 3, false)
}
} else if count == 500 {
fmt.Printf("\n500th : %s", rcu.Commatize(n))
} else if count == 5000 {
fmt.Printf("\n5,000th : %s", rcu.Commatize(n))
} else if count == 50000 {
fmt.Printf("\n50,000th: %s\n", rcu.Commatize(n))
}
}
}
}
- Output:
1 24 30 40 42 54 56 66 70 78 88 102 104 105 110 114 128 130 135 136 138 152 154 165 170 174 182 184 186 189 190 195 222 230 231 232 238 246 248 250 255 258 266 273 282 285 286 290 296 297 500th : 2,526 5,000th : 23,118 50,000th: 223,735
J
Note that the cube root of the product of the proper divisors is the fourth root of the product of all divisors of a positive integer. That said, we do not need to find roots here -- we only need to inspect the powers of the prime factors of the number:
F=: 1 8 e.~_ */@:>:@q:"0 ]
Task examples:
N=: 1+I.F 1+i.2^18
5 10$N
1 24 30 40 42 54 56 66 70 78
88 102 104 105 110 114 128 130 135 136
138 152 154 165 170 174 182 184 186 189
190 195 222 230 231 232 238 246 248 250
255 258 266 273 282 285 286 290 296 297
499{N
2526
4999{N
23118
49999{N
223735
Julia
using Printf
function proper_divisors(n::Integer)
uptosqr = 1:isqrt(n)
divs = Iterators.filter(uptosqr) do m
n % m == 0
end
pd_pairs = Iterators.map(divs) do d1
d2 = div(n, d1)
(d1 == d2 || d1 == 1) ? (d1,) : (d1, d2)
end
return Iterators.flatten(pd_pairs)
end
function show_divisors_print(n::Integer, found::Integer)
if found <= 50
@printf "%5i" n
if found % 10 == 0
println()
end
elseif found in (500, 5_000, 50_000)
th = "$(found)th: "
@printf "%10s%i\n" th n
end
end
function show_divisors()
found = 0
n = 1
while found <= 50_000
pds = proper_divisors(n)
if n^3 == prod(pds)
found += 1
show_divisors_print(n, found)
end
n += 1
end
end
show_divisors()
- Output:
1 24 30 40 42 54 56 66 70 78 88 102 104 105 110 114 128 130 135 136 138 152 154 165 170 174 182 184 186 189 190 195 222 230 231 232 238 246 248 250 255 258 266 273 282 285 286 290 296 297 500th: 2526 5000th: 23118 50000th: 223735
jq
(subject to IEEE 754 limitations)
Also works with gojq, the Go implementation of jq (without such limitations)
Generic utilities
# Notice that `prod(empty)` evaluates to 1.
def prod(s): reduce s as $x (1; . * $x);
# Output: the unordered stream of proper divisors of .
def proper_divisors:
. as $n
| if $n > 1 then 1,
( range(2; 1 + (sqrt|floor)) as $i
| if ($n % $i) == 0 then $i,
(($n / $i) | if . == $i then empty else . end)
else empty
end)
else empty
end;The Task
# Emit a stream beginning with 1 and followed by the integers that are
# cube-roots of their proper divisors
def numbers_being_cube_roots_of_their_proper_divisors:
range(1; infinite)
| select(prod(proper_divisors) == .*.*.);
# print first 50 and then the 500th, 5000th, and $limit-th
def harness(generator; $limit):
label $out
| foreach generator as $n (
{ numbers50: [],
count: 0 };
.emit = null
| .count += 1
| if .count > $limit
then break $out
else if .count <= 50
then .numbers50 += [$n]
else .
end
| if .count == 50
then .emit = .numbers50
elif .count | IN(500, 5000, $limit)
then .emit = "\(.count)th: \($n)"
else .
end
end )
| .emit // empty ;
"First 50 numbers which are the cube roots of the products of their proper divisors:",
harness(numbers_being_cube_roots_of_their_proper_divisors; 50000)- Output:
First 50 numbers which are the cube roots of the products of their proper divisors: [1,24,30,40,42,54,56,66,70,78,88,102,104,105,110,114,128,130,135,136,138,152,154,165,170,174,182,184,186,189,190,195,222,230,231,232,238,246,248,250,255,258,266,273,282,285,286,290,296,297] 500th: 2526 5000th: 23118 50000th: 223735
Lua
The OEIS page gives a formula of "1 together with numbers with 8 divisors", so that's what we test.
function is_1_or_has_eight_divisors (n)
if n == 1 then return true end
local divCount, sqr = 2, math.sqrt(n)
for d = 2, sqr do
if n % d == 0 then
divCount = d == sqr and divCount + 1 or divCount + 2
end
if divCount > 8 then return false end
end
return divCount == 8
end
-- First 50
local count, x = 1, 0
while count <= 50 do
x = x + 1
if is_1_or_has_eight_divisors(x) then
io.write(x .. " ")
count = count + 1
end
end
-- 500th, 5,000th and 50,000th
while count <= 50000 do
x = x + 1
if is_1_or_has_eight_divisors(x) then
if count == 500 then print("\n\n500th: " .. x) end
if count == 5000 then print("5,000th: " .. x) end
count = count + 1
end
end
print("50,000th: " .. x)
- Output:
1 24 30 40 42 54 56 66 70 78 88 102 104 105 110 114 128 130 135 136 138 152 154 165 170 174 182 184 186 189 190 195 222 230 231 232 238 246 248 250 255 258 266 273 282 285 286 290 296 297 500th: 2526 5,000th: 23118 50,000th: 223735
Pascal
Free Pascal
As stated, the result are the numbers with 8 = 2^3 divisors.Therefor only numbers with prime decomposition of the form:
8 = 2^3 ( all powers+1 must be a power of 2 )
a^7 , a^3*b ( a <> b) and a*b*c (a>b>c ( oBdA ) ), of cause all prime
Avoid sorting by using an array of limit size for only marking those numbers.
program Root3rd_divs_n.pas;
{$IFDEF FPC}
{$MODE DELPHI} {$OPTIMIZATION ON,ALL} {$COPERATORS ON}
{$ENDIF}
{$IFDEF WINDOWS}
{$APPTYPE CONSOLE}
{$ENDIF}
uses
sysutils
{$IFDEF WINDOWS},Windows{$ENDIF}
;
const
limit = 110*1000 *1000;
var
sol : array [0..limit] of byte;
primes : array of Uint32;
gblCount: Uint64;
procedure SievePrimes(lmt:Uint32);
var
sieve : array of byte;
p,i,delta : NativeInt;
Begin
setlength(sieve,lmt DIV 2);
//estimate count of prime
i := trunc(lmt/(ln(lmt)-1.1));
setlength(primes,i);
p := 1;
repeat
delta := 2*p+1;
// ((2*p+1)^2 ) -1)/ 2 = ((4*p*p+4*p+1) -1)/2 = 2*p*(p+1)
i := 2*p*(p+1);
if i>High(sieve) then
BREAK;
while i <= High(sieve) do
begin
sieve[i] := 1;
i += delta;
end;
repeat
inc(p);
until sieve[p] = 0;
until false;
primes[0] := 2;
i := 1;
For p := 1 to High(sieve) do
if sieve[p] = 0 then
begin
primes[i] := 2*p+1;
inc(i);
end;
setlength(primes,i);
end;
procedure Get_a7;
var
q3,n : UInt64;
i : nativeInt;
begin
sol[1] := 1;
gblCount +=1;
For i := 0 to High(primes) do
begin
q3 := primes[i];
n := sqr(sqr(sqr(q3))) DIV q3;
if n > limit then
break;
sol[n] := 1;
gblCount +=1;
end;
end;
procedure Get_a3_b;
var
i,j,q3,n : nativeInt;
begin
For i := 0 to High(primes) do
begin
q3 := primes[i];
q3 := q3*q3*q3;
if q3 > limit then
BREAK;
For j := 0 to High(primes) do
begin
if j = i then
continue;
n := Primes[j]*q3;
if n > limit then
BREAK;
sol[n] := 1;
gblCount +=1;
end;
end;
end;
procedure Get_a_b_c;
var
i,j,k,q1,q2,n : nativeInt;
begin
For i := 0 to High(primes)-2 do
begin
q1 := primes[i];
For j := i+1 to High(primes)-1 do
Begin
q2:= q1*Primes[j];
if q2 > limit then
BREAK;
For k := j+1 to High(primes) do
Begin
n:= q2*Primes[k];
if n > limit then
BREAK;
sol[n] := 1;
gblCount +=1;
end;
end;
end;
end;
var
i,cnt,lmt : Int32;
begin
SievePrimes(limit DIV 8);// 2^3 * (prime > 2)
gblCount := 0;
Get_a7;
Get_a3_b;
Get_a_b_c;
Writeln('First 50 numbers which are the cube roots of the products of their proper divisors:');
cnt := 0;
i := 1;
while cnt < 50 do
begin
if sol[i] <> 0 then
begin
write(i:5);
cnt +=1;
if cnt mod 10 = 0 then writeln;
end;
inc(i);
end;
dec(i);
lmt := 500;
repeat
while cnt < lmt do
begin
inc(i);
if sol[i] <> 0 then
cnt +=1;
if i > limit then
break;
end;
if i > limit then
break;
writeln(lmt:8,'.th:',i:12);
lmt *= 10;
until lmt> limit;
writeln('Total found: ', gblCount, ' til ',limit);
end.
- @TIO.RUN:
First 50 numbers which are the cube roots of the products of their proper divisors:
1 24 30 40 42 54 56 66 70 78
88 102 104 105 110 114 128 130 135 136
138 152 154 165 170 174 182 184 186 189
190 195 222 230 231 232 238 246 248 250
255 258 266 273 282 285 286 290 296 297
500.th: 2526
5000.th: 23118
50000.th: 223735
Total found: 243069 til 1100000
Real time: 0.144 s CPU share: 99.00 %
..
500000.th: 2229229
5000000.th: 22553794
Total found: 23797493 til 110000000
Real time: 1.452 s CPU share: 99.05 %
Perl
use v5.36;
use ntheory 'divisors';
use List::Util <max product>;
sub table ($c, @V) { my $t = $c * (my $w = 2 + length max @V); ( sprintf( ('%'.$w.'d')x@V, @V) ) =~ s/.{1,$t}\K/\n/gr }
sub proper_divisors ($n) { my @d = divisors($n); pop @d; @d }
sub is_N ($n) {
state @N = 1;
state $p = 1;
do { push @N, $p if ++$p**3 == product proper_divisors($p); } until $N[$n];
$N[$n-1]
}
say table 10, map { is_N $_ } 1..50;
printf "%5d %d\n", $_, is_N $_ for 500, 5000, 50000;
- Output:
1 24 30 40 42 54 56 66 70 78 88 102 104 105 110 114 128 130 135 136 138 152 154 165 170 174 182 184 186 189 190 195 222 230 231 232 238 246 248 250 255 258 266 273 282 285 286 290 296 297 500 2526 5000 23118 50000 223735
Phix
with javascript_semantics sequence n50 = {} integer count = 0, n = 1, n5 = 500 atom t0 = time() printf(1,"First 50 numbers which are the cube roots\n") printf(1," of the products of their proper divisors:\n") while count<500000 do -- if product(factors(n))=n*n*n then if n=1 or length(factors(n))=6 then -- safer/smidge faster count += 1 if count<=50 then n50 &= n if count=50 then printf(1,"%s\n",join_by(n50,1,10,"",fmt:="%4d")) end if elsif count=n5 then printf(1,"%,8dth: %,d (%s)\n",{n5,n,elapsed(time()-t0)}) n5 *= 10 end if end if n += 1 end while
By default factors() does not include 1 and n, or I could use length(factors(n,1))=8, both 25% faster than using product(), which exceeds precision limits on 32-bit for n=180, and on 64bit for n=240, though since you'll get exactly the same precision error on the n*n*n it kinda "worked by chance".
- Output:
First 50 numbers which are the cube roots
of the products of their proper divisors:
1 24 30 40 42 54 56 66 70 78
88 102 104 105 110 114 128 130 135 136
138 152 154 165 170 174 182 184 186 189
190 195 222 230 231 232 238 246 248 250
255 258 266 273 282 285 286 290 296 297
500th: 2,526 (0.0s)
5,000th: 23,118 (0.0s)
50,000th: 223,735 (0.6s)
500,000th: 2,229,229 (14.1s)
For comparison, the gcc/C++ entry gets the 500kth about 8* faster, roughly about what I'd expect... 🤥
PL/M
Solves the basic task by counting the proper divisors as per the OEIS page (the 50 000th number is too large for 16 bits).
... under CP/M (or an emulator)
100H: /* FIND NUMBERS THAT ARE THE CUBE ROOT OF THEIR PROPER DIVISORS */
DECLARE FALSE LITERALLY '0', TRUE LITERALLY '0FFH';
/* CP/M SYSTEM CALL AND I/O ROUTINES */
BDOS: PROCEDURE( FN, ARG ); DECLARE FN BYTE, ARG ADDRESS; GOTO 5; END;
PR$CHAR: PROCEDURE( C ); DECLARE C BYTE; CALL BDOS( 2, C ); END;
PR$STRING: PROCEDURE( S ); DECLARE S ADDRESS; CALL BDOS( 9, S ); END;
PR$NL: PROCEDURE; CALL PR$CHAR( 0DH ); CALL PR$CHAR( 0AH ); END;
PR$NUMBER: PROCEDURE( N ); /* PRINTS A NUMBER IN THE MINIMUN FIELD WIDTH */
DECLARE N ADDRESS;
DECLARE V ADDRESS, N$STR ( 6 )BYTE, W BYTE;
V = N;
W = LAST( N$STR );
N$STR( W ) = '$';
N$STR( W := W - 1 ) = '0' + ( V MOD 10 );
DO WHILE( ( V := V / 10 ) > 0 );
N$STR( W := W - 1 ) = '0' + ( V MOD 10 );
END;
CALL PR$STRING( .N$STR( W ) );
END PR$NUMBER;
/* END SYSTEM CALL AND I/O ROUTINES */
DECLARE PDC ( 5000 )ADDRESS;
DECLARE ( I, I2, J, COUNT ) ADDRESS;
DO I = 1 TO LAST( PDC ); PDC( I ) = 1; END;
DO I = 2 TO LAST( PDC );
I2 = I + I;
DO J = I2 TO LAST( PDC ) BY I;
PDC( J ) = PDC( J ) + 1;
END;
END;
PDC( 1 ) = 7;
COUNT, I = 0;
DO WHILE COUNT < 500 AND I < LAST( PDC );
I = I + 1;
IF PDC( I ) = 7 THEN DO;
IF ( COUNT := COUNT + 1 ) < 51 THEN DO;
CALL PR$CHAR( ' ' );
IF I < 10 THEN CALL PR$CHAR( ' ' );
IF I < 100 THEN CALL PR$CHAR( ' ' );
IF I < 1000 THEN CALL PR$CHAR( ' ' );
CALL PR$NUMBER( I );
IF COUNT MOD 10 = 0 THEN CALL PR$NL;
END;
ELSE IF COUNT = 500 THEN DO;
CALL PR$NUMBER( COUNT );
CALL PR$STRING( .'TH: $' );
CALL PR$NUMBER( I );
CALL PR$NL;
END;
END;
END;
EOF- Output:
1 24 30 40 42 54 56 66 70 78 88 102 104 105 110 114 128 130 135 136 138 152 154 165 170 174 182 184 186 189 190 195 222 230 231 232 238 246 248 250 255 258 266 273 282 285 286 290 296 297 500TH: 2526
Alternative version, calculating the proper divisor products and cubes modulo 65536 (as PL/M uses unsigned 16 bit arithmetic and doesn't check for overflow, all calculations are modulo 65536). This is sufficient to detect the numbers apart from those where the product/cube is 0 mod 65536. To handle the zero cases, it uses Rdm's hints (see J sample and Discussion page) that if x = n^3 then the prime factors of x must be the same as the prime factors of n and the prime factors of x must have powers three times those of n - additionally, we don't have to calclate the product of the proper divisors, we only need to factorise them and aggregate their powers.
Using this technique, the first 50 numbers can be found in a few seconds but to find the 5000th takes several minutes. As the candidates increase, the proportion that have cubes that are 0 mod 65536 increases and the factorisation and aggregation is quite expensive (the code could doubtless be improved).
... under CP/M (or an emulator)
100H: /* FIND NUMBERS THAT ARE THE CUBE ROOT OF THEIR PROPER DIVISORS */
DECLARE FALSE LITERALLY '0', TRUE LITERALLY '0FFH';
/* CP/M SYSTEM CALL AND I/O ROUTINES */
BDOS: PROCEDURE( FN, ARG ); DECLARE FN BYTE, ARG ADDRESS; GOTO 5; END;
PR$CHAR: PROCEDURE( C ); DECLARE C BYTE; CALL BDOS( 2, C ); END;
PR$STRING: PROCEDURE( S ); DECLARE S ADDRESS; CALL BDOS( 9, S ); END;
PR$NL: PROCEDURE; CALL PR$CHAR( 0DH ); CALL PR$CHAR( 0AH ); END;
PR$NUMBER: PROCEDURE( N ); /* PRINTS A NUMBER IN THE MINIMUN FIELD WIDTH */
DECLARE N ADDRESS;
DECLARE V ADDRESS, N$STR ( 6 )BYTE, W BYTE;
V = N;
W = LAST( N$STR );
N$STR( W ) = '$';
N$STR( W := W - 1 ) = '0' + ( V MOD 10 );
DO WHILE( ( V := V / 10 ) > 0 );
N$STR( W := W - 1 ) = '0' + ( V MOD 10 );
END;
CALL PR$STRING( .N$STR( W ) );
END PR$NUMBER;
/* END SYSTEM CALL AND I/O ROUTINES */
DECLARE MAX$PF LITERALLY '200';
/* SETS PF$A AND PFC$A TO THE PRIME FACTORS AND COUNTS OF F, THE NUMBER */
/* NUMBER OF FACTORS IS RETURNED IN PF$POS$PTR */
/* PF$POS$PTR MUST BE INITIALISED BEFORE THE CALL */
FACTORISE: PROCEDURE( F, PF$POS$PTR, PF$A, PFC$A );
DECLARE ( F, PF$POS$PTR, PF$A, PFC$A ) ADDRESS;
DECLARE PF$POS BASED PF$POS$PTR ADDRESS;
DECLARE PF BASED PF$A ( 0 )ADDRESS;
DECLARE PFC BASED PFC$A ( 0 )ADDRESS;
DECLARE ( FF, V, POWER ) ADDRESS;
/* START WITH 2 */
V = F;
FF = 2;
DO WHILE V > 1;
IF V MOD FF = 0 THEN DO;
/* FF IS A PRIME FACTOR OF V */
DECLARE P ADDRESS;
POWER = 0;
DO WHILE V MOD FF = 0;
POWER = POWER + 1;
V = V / FF;
END;
P = 0;
DO WHILE P < PF$POS AND PF( P ) <> FF;
P = P + 1;
END;
IF P >= PF$POS THEN DO;
/* FIRST TIME FF HAS APPEARED AS A PRIME FACTOR */
P = PF$POS;
PFC( P ) = 0;
PF$POS = PF$POS + 1;
END;
PF( P ) = FF;
PFC( P ) = PFC( P ) + POWER;
END;
IF FF = 2 THEN FF = 3; ELSE FF = FF + 2;
END;
END FACTORISE;
/* RETURNS TRUE THE PRODUCT OF THE PROPER DIVISORS OF N IS THE CUBE OF N */
/* FALSE OTHERWISE */
PD$PRODUCT$IS$CUBE: PROCEDURE( N )ADDRESS;
DECLARE N ADDRESS;
DECLARE IS$CUBE BYTE;
IF N < 2
THEN IS$CUBE = TRUE;
ELSE DO;
DECLARE ( I, PF$POS, NF$POS ) ADDRESS;
DECLARE ( PF, PFC, NF, NFC ) ( MAX$PF ) ADDRESS;
PFC( 0 ), PF( 0 ), PF$POS, NFC( 0 ), NF( 0 ), NF$POS = 0;
/* FACTORISE N */
CALL FACTORISE( N, .NF$POS, .NF, .NFC );
/* COPY FACTORS BUT ZERO THE COUNTS SO WE CAN EASILY CHECK THE */
/* FACTORS OF N ARE THE SAME AS THOSE OF THE PROPER DIVISOR PRODUCT */
DO I = 0 TO NF$POS - 1;
PF( I ) = NF( I );
PFC( I ) = 0;
END;
/* FIND THE PROPER DIVISORS AND FACTORISE THEM, ACCUMULATING THE */
/* PRIME FACTOR COUNTS */
I = 2;
DO WHILE I * I <= N;
IF N MOD I = 0 THEN DO;
/* I IS A DIVISOR OF N */
DECLARE ( F, G ) ADDRESS;
F = I; /* FIRST FACTOR */
G = N / F; /* SECOND FACTOR */
/* FACTORISE F, COUNTING THE PRIME FACTORS */
CALL FACTORISE( F, .PF$POS, .PF, .PFC );
/* FACTORISE G, IF IT IS NOT THE SAME AS F */
IF F <> G THEN CALL FACTORISE( G, .PF$POS, .PF, .PFC );
END;
I = I + 1;
END;
IS$CUBE = PF$POS = NF$POS;
IF IS$CUBE THEN DO;
/* N AND ITS PROPER DIVISOR PRODUCT HAVE THE SAME PRIME FACTOR */
/* COUNT - CHECK THE PRIME FACTLORS ARE THE SAME AND THAT THE */
/* PRODUCTS FACTORS APPEAR 3 TIMEs THOSE OF N */
I = 0;
DO WHILE I < PF$POS AND IS$CUBE;
IS$CUBE = ( PF( I ) = NF( I ) )
AND ( PFC( I ) = NFC( I ) * 3 );
I = I + 1;
END;
END;
END;
RETURN IS$CUBE;
END;
/* RETURNS THE PROPER DIVISOR PRODUCT OF N, MOD 65536 */
PDP: PROCEDURE( N )ADDRESS;
DECLARE N ADDRESS;
DECLARE ( I, I2, PRODUCT ) ADDRESS;
PRODUCT = 1;
I = 2;
DO WHILE ( I2 := I * I ) <= N;
IF N MOD I = 0 THEN DO;
PRODUCT = PRODUCT * I;
IF I2 <> N THEN DO;
PRODUCT = PRODUCT * ( N / I );
END;
END;
I = I + 1;
END;
RETURN PRODUCT;
END PDP;
DECLARE ( I, I3, J, COUNT ) ADDRESS;
COUNT, I = 0;
DO WHILE COUNT < 5$000;
I = I + 1;
I3 = I * I * I;
IF PDP( I ) = I3 THEN DO;
/* THE PROPER DIVISOR PRODUCT MOD 65536 IS THE SAME AS N CUBED ALSO */
/* MOD 65536, IF THE CUBE IS 0 MOD 65536, WE NEED TO CHECK THE */
/* PRIME FACTORS */
DECLARE IS$NUMBER BYTE;
IF I3 <> 0 THEN IS$NUMBER = TRUE;
ELSE IS$NUMBER = PD$PRODUCT$IS$CUBE( I );
IF IS$NUMBER THEN DO;
IF ( COUNT := COUNT + 1 ) < 51 THEN DO;
CALL PR$CHAR( ' ' );
IF I < 10 THEN CALL PR$CHAR( ' ' );
IF I < 100 THEN CALL PR$CHAR( ' ' );
IF I < 1000 THEN CALL PR$CHAR( ' ' );
CALL PR$NUMBER( I );
IF COUNT MOD 10 = 0 THEN CALL PR$NL;
END;
ELSE IF COUNT = 500 OR COUNT = 5000 THEN DO;
IF COUNT < 1000 THEN CALL PR$CHAR( ' ' );
CALL PR$STRING( .' $' );
CALL PR$NUMBER( COUNT );
CALL PR$STRING( .'TH: $' );
CALL PR$NUMBER( I );
CALL PR$NL;
END;
END;
END;
END;
EOF- Output:
1 24 30 40 42 54 56 66 70 78
88 102 104 105 110 114 128 130 135 136
138 152 154 165 170 174 182 184 186 189
190 195 222 230 231 232 238 246 248 250
255 258 266 273 282 285 286 290 296 297
500TH: 2526
5000TH: 23118
Python
''' Rosetta code rosettacode.org/wiki/Numbers_which_are_the_cube_roots_of_the_product_of_their_proper_divisors '''
from functools import reduce
from sympy import divisors
FOUND = 0
for num in range(1, 1_000_000):
divprod = reduce(lambda x, y: x * y, divisors(num)[:-1])if num > 1 else 1
if num * num * num == divprod:
FOUND += 1
if FOUND <= 50:
print(f'{num:5}', end='\n' if FOUND % 10 == 0 else '')
if FOUND == 500:
print(f'\nFive hundreth: {num:,}')
if FOUND == 5000:
print(f'\nFive thousandth: {num:,}')
if FOUND == 50000:
print(f'\nFifty thousandth: {num:,}')
break
- Output:
1 24 30 40 42 54 56 66 70 78 88 102 104 105 110 114 128 130 135 136 138 152 154 165 170 174 182 184 186 189 190 195 222 230 231 232 238 246 248 250 255 258 266 273 282 285 286 290 296 297 Five hundreth: 2,526 Five thousandth: 23,118 Fifty thousandth: 223,735
OEIS algorithm (see talk pages)
from sympy import divisors
numfound = 0
for num in range(1, 1_000_000):
if num == 1 or len(divisors(num)) == 8:
numfound += 1
if numfound <= 50:
print(f'{num:5}', end='\n' if numfound % 10 == 0 else '')
if numfound == 500:
print(f'\nFive hundreth: {num:,}')
if numfound == 5000:
print(f'\nFive thousandth: {num:,}')
if numfound == 50000:
print(f'\nFifty thousandth: {num:,}')
break
Output same as first algorithm.
Quackery
factors is defined at Factors of an integer#Quackery.
' [ 1 ] 1
[ 1+ dup
factors size 8 = until
tuck join swap
over size 50000 = until ]
drop
dup 50 split drop echo cr cr
dup 499 peek echo cr cr
dup 4999 peek echo cr cr
49999 peek echo- Output:
[ 1 24 30 40 42 54 56 66 70 78 88 102 104 105 110 114 128 130 135 136 138 152 154 165 170 174 182 184 186 189 190 195 222 230 231 232 238 246 248 250 255 258 266 273 282 285 286 290 296 297 ] 2526 23118 223735
Raku
use Prime::Factor;
use Lingua::EN::Numbers;
my @cube-div = lazy 1, |(2..∞).hyper.grep: { .³ == [×] .&proper-divisors }
put "First 50 numbers which are the cube roots of the products of their proper divisors:\n" ~
@cube-div[^50]».fmt("%3d").batch(10).join: "\n";
printf "\n%16s: %s\n", .Int.&ordinal.tc, comma @cube-div[$_ - 1] for 5e2, 5e3, 5e4;
- Output:
First 50 numbers which are the cube roots of the products of their proper divisors: 1 24 30 40 42 54 56 66 70 78 88 102 104 105 110 114 128 130 135 136 138 152 154 165 170 174 182 184 186 189 190 195 222 230 231 232 238 246 248 250 255 258 266 273 282 285 286 290 296 297 Five hundredth: 2,526 Five thousandth: 23,118 Fifty thousandth: 223,735
Ruby
require 'prime'
def tau(n) = n.prime_division.inject(1){|res, (d, exp)| res *= exp+1}
a111398 = [1].chain (1..).lazy.select{|n| tau(n) == 8}
puts "The first 50 numbers which are the cube roots of the products of their proper divisors:"
p a111398.first(50)
[500, 5000, 50000].each{|n| puts "#{n}th: #{a111398.drop(n-1).next}" }
- Output:
The first 50 numbers which are the cube roots of the products of their proper divisors: [1, 24, 30, 40, 42, 54, 56, 66, 70, 78, 88, 102, 104, 105, 110, 114, 128, 130, 135, 136, 138, 152, 154, 165, 170, 174, 182, 184, 186, 189, 190, 195, 222, 230, 231, 232, 238, 246, 248, 250, 255, 258, 266, 273, 282, 285, 286, 290, 296, 297] 500th: 2526 5000th: 23118 50000th: 223735
Wren
import "./math" for Int, Nums
import "./long" for ULong, ULongs
import "./fmt" for Fmt
var numbers50 = []
var count = 0
var n = 1
var ln
var maxSafe = Num.maxSafeInteger.cbrt.floor
System.print("First 50 numbers which are the cube roots of the products of their proper divisors:")
while (true) {
var pd = Int.properDivisors(n)
if ((n <= maxSafe && Nums.prod(pd) == n * n * n) ||
(ULongs.prod(pd.map { |f| ULong.new(f) }) == (ln = ULong.new(n)) * ln * ln )) {
count = count + 1
if (count <= 50) {
numbers50.add(n)
if (count == 50) Fmt.tprint("$3d", numbers50, 10)
} else if (count == 500) {
Fmt.print("\n500th : $,d", n)
} else if (count == 5000) {
Fmt.print("5,000th : $,d", n)
} else if (count == 50000) {
Fmt.print("50,000th: $,d", n)
break
}
}
n = n + 1
}- Output:
First 50 numbers which are the cube roots of the products of their proper divisors: 1 24 30 40 42 54 56 66 70 78 88 102 104 105 110 114 128 130 135 136 138 152 154 165 170 174 182 184 186 189 190 195 222 230 231 232 238 246 248 250 255 258 266 273 282 285 286 290 296 297 500th : 2,526 5,000th : 23,118 50,000th: 223,735
Alternatively and a bit quicker, inspired by the C++ entry and the OEIS comment that (apart from 1) n must have exactly 8 divisors:
import "./fmt" for Fmt
var divisorCount = Fn.new { |n|
var i = 1
var k = (n%2 == 0) ? 1 : 2
var count = 0
while (i <= n.sqrt) {
if (n%i == 0) {
count = count + 1
var j = (n/i).floor
if (j != i) count = count + 1
}
i = i + k
}
return count
}
var numbers50 = []
var count = 0
var n = 1
System.print("First 50 numbers which are the cube roots of the products of their proper divisors:")
while (true) {
var dc = divisorCount.call(n)
if (n == 1 || dc == 8) {
count = count + 1
if (count <= 50) {
numbers50.add(n)
if (count == 50) Fmt.tprint("$3d", numbers50, 10)
} else if (count == 500) {
Fmt.print("\n500th : $,d", n)
} else if (count == 5000) {
Fmt.print("5,000th : $,d", n)
} else if (count == 50000) {
Fmt.print("50,000th: $,d", n)
break
}
}
n = n + 1
}- Output:
Same as first version.
XPL0
func DivisorCount(N); \Return count of divisors
int N, Total, P, Count;
[Total:= 1;
while (N&1) = 0 do
[Total:= Total+1;
N:= N>>1;
];
P:= 3;
while P*P <= N do
[Count:= 1;
while rem(N/P) = 0 do
[Count:= Count+1;
N:= N/P;
];
Total:= Total*Count;
P:= P+2;
];
if N > 1 then
Total:= Total*2;
return Total;
];
int N, Count;
[Text(0, "First 50 numbers which are the cube roots of the products of ");
Text(0, "their proper divisors:^m^j");
N:= 1; Count:= 0;
repeat if N = 1 or DivisorCount(N) = 8 then
[Count:= Count+1;
if Count <= 50 then
[Format(4, 0);
RlOut(0, float(N));
if rem(Count/10) = 0 then CrLf(0);
]
else if Count = 500 or Count = 5000 or Count = 50000 then
[Format(6, 0);
RlOut(0, float(Count));
Text(0, "th: ");
IntOut(0, N);
CrLf(0);
];
];
N:= N+1;
until Count >= 50000;
]- Output:
First 50 numbers which are the cube roots of the products of their proper divisors: 1 24 30 40 42 54 56 66 70 78 88 102 104 105 110 114 128 130 135 136 138 152 154 165 170 174 182 184 186 189 190 195 222 230 231 232 238 246 248 250 255 258 266 273 282 285 286 290 296 297 500th: 2526 5000th: 23118 50000th: 223735