Multiplicatively perfect numbers: Difference between revisions
Content added Content deleted
(→{{header|C}}: Optimized and extended to 5 million,) |
(→{{header|C}}: Oops, even better :)) |
||
| Line 101: | Line 101: | ||
=={{header|C}}== |
=={{header|C}}== |
||
This includes '1' as an MPN. Run time around |
This includes '1' as an MPN. Run time around 2.3 seconds. |
||
<syntaxhighlight lang="c">#include <stdio.h> |
<syntaxhighlight lang="c">#include <stdio.h> |
||
#include <stdbool.h> |
#include <stdbool.h> |
||
| Line 127: | Line 127: | ||
if (!(n%i)) { |
if (!(n%i)) { |
||
divs[c++] = i; |
divs[c++] = i; |
||
if (c > |
if (c > 2) break; // not eligible if has > 2 divisors |
||
j = n / i; |
j = n / i; |
||
if (j != i) divs[c++] = j; |
if (j != i) divs[c++] = j; |
||
| Line 137: | Line 137: | ||
int main() { |
int main() { |
||
int i, d, j, k, t, length, prod; |
int i, d, j, k, t, length, prod; |
||
int divs[ |
int divs[4], count = 0, limit = 500, s = 3, c = 3, squares = 1, cubes = 1; |
||
printf("Multiplicatively perfect numbers under %d:\n", limit); |
printf("Multiplicatively perfect numbers under %d:\n", limit); |
||
setlocale(LC_NUMERIC, ""); |
setlocale(LC_NUMERIC, ""); |
||
Revision as of 18:04, 12 May 2023
Multiplicatively perfect numbers is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.
If the product of the divisors of an integer n (including n itself) is equal to n^2, then n is a multiplicatively perfect number. Alternatively: the product of the proper divisors of n (i.e. excluding n) is equal to n.
Note that the integer '1' qualifies under the first definition (1 = 1 x 1), but not under the second since '1' has no proper divisors. Given this ambiguity, it is optional whether '1' is included or not in solutions to this task.
- Task
Find and show on this page the multiplicatively perfect numbers below 500.
- Stretch
Find and show the number of multiplicatively perfect numbers under 500, 5,000, 50,000 and 500,000 and for each of these limits deduce (avoid counting separately) and show the number of semi-primes (numbers which are the product of exactly two primes) under that limit.
- Related (and near duplicate) task
- See also
- The OEIS sequence: A007422: Multiplicatively perfect numbers
ALGOL 68
Uses sieves, counts 1 as a multiplicatively perfect number.
As with the Phix and probably other samples, uses the fact that the multiplicatively perfect numbers (other than 1) must have three proper divisors - see OIES A007422.
BEGIN # find multiplicatively perfect numbers - numbers whose proper #
# divisor product is the number itself ( includes 1 ) #
PR read "primes.incl.a68" PR # include prime utilties #
INT max number = 500 000; # largest number we wlil consider #
[]BOOL prime = PRIMESIEVE max number; # sieve the primes to max number #
[ 1 : max number ]INT pdc; # table of proper divisor counts #
[ 1 : max number ]INT pfc; # table of prime factor counts #
# count the proper divisors and prime factors #
FOR n TO UPB pdc DO pdc[ n ] := 1; pfc[ n ] := 0 OD;
FOR i FROM 2 TO UPB pdc DO
BOOL is prime = prime[ i ];
INT m := 1; # j will be the mth multiple of i #
FOR j FROM i + i BY i TO UPB pdc DO
pdc[ j ] +:= 1;
IF prime[ m +:= 1 ] THEN # j is a prime multiple of i #
pfc[ j ] +:= 1;
IF i = m THEN # j is i squared #
pfc[ j ] +:= 1
FI
ELSE # j is not a prime multiple of i #
pfc[ j ] +:= pfc[ m ] # add the prime factor count of m #
FI
OD
OD;
# find the mutiplicatively perfect numbers #
INT mp count := 0;
INT sp count := 0;
INT next to show := 500;
FOR n TO UPB pdc DO
IF n = 1 OR pdc[ n ] = 3 THEN
# n is 1 or has 3 proper divisors so is multiplicatively perfect #
# number - see OEIS A007422 #
mp count +:= 1;
IF n < 500 THEN
print( ( " ", whole( n, -3 ) ) );
IF mp count MOD 10 = 0 THEN print( ( newline ) ) FI
FI
FI;
IF pfc[ n ] = 2 THEN
# 2 prime factors, n is semi-prime #
sp count +:= 1
FI;
IF n = next to show THEN
print( ( "Up to ", whole( next to show, -8 )
, " there are ", whole( mp count, -6 )
, " multiplicatively perfect numbers and ", whole( sp count, -6 )
, " semi-primes", newline
)
);
next to show *:= 10
FI
OD
END- Output:
1 6 8 10 14 15 21 22 26 27 33 34 35 38 39 46 51 55 57 58 62 65 69 74 77 82 85 86 87 91 93 94 95 106 111 115 118 119 122 123 125 129 133 134 141 142 143 145 146 155 158 159 161 166 177 178 183 185 187 194 201 202 203 205 206 209 213 214 215 217 218 219 221 226 235 237 247 249 253 254 259 262 265 267 274 278 287 291 295 298 299 301 302 303 305 309 314 319 321 323 326 327 329 334 335 339 341 343 346 355 358 362 365 371 377 381 382 386 391 393 394 395 398 403 407 411 413 415 417 422 427 437 445 446 447 451 453 454 458 466 469 471 473 478 481 482 485 489 493 497 Up to 500 there are 150 multiplicatively perfect numbers and 153 semi-primes Up to 5000 there are 1354 multiplicatively perfect numbers and 1365 semi-primes Up to 50000 there are 12074 multiplicatively perfect numbers and 12110 semi-primes Up to 500000 there are 108223 multiplicatively perfect numbers and 108326 semi-primes
C
This includes '1' as an MPN. Run time around 2.3 seconds.
#include <stdio.h>
#include <stdbool.h>
#include <locale.h>
bool isPrime(int n) {
if (n < 2) return false;
if (n%2 == 0) return n == 2;
if (n%3 == 0) return n == 3;
int d = 5;
while (d*d <= n) {
if (n%d == 0) return false;
d += 2;
if (n%d == 0) return false;
d += 4;
}
return true;
}
void divisors(int n, int *divs, int *length) {
int i, j, k = 1, c = 0;
if (n%2) k = 2;
for (i = 1; i*i <= n; i += k) {
if (i == 1) continue; // exclude 1 and n
if (!(n%i)) {
divs[c++] = i;
if (c > 2) break; // not eligible if has > 2 divisors
j = n / i;
if (j != i) divs[c++] = j;
}
}
*length = c;
}
int main() {
int i, d, j, k, t, length, prod;
int divs[4], count = 0, limit = 500, s = 3, c = 3, squares = 1, cubes = 1;
printf("Multiplicatively perfect numbers under %d:\n", limit);
setlocale(LC_NUMERIC, "");
for (i = 1; ; ++i) {
if (i != 1) {
divisors(i, divs, &length);
} else {
divs[1] = divs[0] = 1;
length = 2;
}
if (length == 2 && divs[0] * divs[1] == i) {
++count;
if (i < 500) {
printf("%3d ", i);
if (!(count%10)) printf("\n");
}
}
if (i == 499) printf("\n");
if (i >= limit - 1) {
for (j = s; j * j < limit; j += 2) if (isPrime(j)) ++squares;
for (k = c; k * k * k < limit; k +=2 ) if (isPrime(k)) ++cubes;
t = count + squares - cubes - 1;
printf("Counts under %'9d: MPNs = %'7d Semi-primes = %'7d\n", limit, count, t);
if (limit == 5000000) break;
s = j;
c = k;
limit *= 10;
}
}
return 0;
}
- Output:
Multiplicatively perfect numbers under 500: 1 6 8 10 14 15 21 22 26 27 33 34 35 38 39 46 51 55 57 58 62 65 69 74 77 82 85 86 87 91 93 94 95 106 111 115 118 119 122 123 125 129 133 134 141 142 143 145 146 155 158 159 161 166 177 178 183 185 187 194 201 202 203 205 206 209 213 214 215 217 218 219 221 226 235 237 247 249 253 254 259 262 265 267 274 278 287 291 295 298 299 301 302 303 305 309 314 319 321 323 326 327 329 334 335 339 341 343 346 355 358 362 365 371 377 381 382 386 391 393 394 395 398 403 407 411 413 415 417 422 427 437 445 446 447 451 453 454 458 466 469 471 473 478 481 482 485 489 493 497 Counts under 500: MPNs = 150 Semi-primes = 153 Counts under 5,000: MPNs = 1,354 Semi-primes = 1,365 Counts under 50,000: MPNs = 12,074 Semi-primes = 12,110 Counts under 500,000: MPNs = 108,223 Semi-primes = 108,326 Counts under 5,000,000: MPNs = 978,983 Semi-primes = 979,274
FreeBASIC
#define ceil(x) (-((-x*2.0-0.5) Shr 1))
Dim As Integer limit = 500
Dim As Integer n, pro, Divisors(), m, c = 0, ub
Print "Special numbers under"; limit; ":"
For n = 1 To limit
pro = 1
For m = 2 To ceil(n / 2)
If n Mod m = 0 Then
pro *= m
Redim Preserve Divisors(c) : Divisors(c) = m
c += 1
End If
Next m
ub = Ubound(Divisors)
If n = pro And ub > 1 Then
Print Using "### = ## x ###"; n; Divisors(ub-1); Divisors(ub)
End If
Next n
Sleep- Output:
Similar to Ring entry.
Go
package main
import (
"fmt"
"rcu"
)
func main() {
count := 0
limit := 500
s := 3
c := 3
squares := 1
cubes := 1
fmt.Printf("Multiplicatively perfect numbers under %d:\n", limit)
var divs []int
for i := 0; ; i++ {
if i != 1 {
divs = rcu.ProperDivisors(i)
} else {
divs = []int{1, 1}
}
length := len(divs)
if length > 1 {
prod := 1
for d := 1; d < length; d++ {
prod *= divs[d]
}
if prod == i {
count++
if i < 500 {
fmt.Printf("%3d ", i)
if count%10 == 0 {
fmt.Println()
}
}
}
}
if i == 499 {
fmt.Println()
}
if i >= limit-1 {
var j, k int
for j = s; j*j < limit; j += 2 {
if rcu.IsPrime(j) {
squares++
}
}
for k = c; k*k*k < limit; k += 2 {
if rcu.IsPrime(k) {
cubes++
}
}
t := count + squares - cubes - 1
slimit := rcu.Commatize(limit)
scount := rcu.Commatize(count)
st := rcu.Commatize(t)
fmt.Printf("Counts under %7s: MPNs = %7s Semi-primes = %7s\n", slimit, scount, st)
if limit == 500000 {
break
}
s, c = j, k
limit *= 10
}
}
}
- Output:
Same as C example.
Haskell
divisors :: Int -> [Int]
divisors n = [d | d <- [1..n] , mod n d == 0 ]
isMultiplicativelyPerfect :: Int -> Bool
isMultiplicativelyPerfect n = (product $ divisors n) == n ^ 2
solution :: [Int]
solution = filter isMultiplicativelyPerfect [1..500]
- Output:
[1,6,8,10,14,15,21,22,26,27,33,34,35,38,39,46,51,55,57,58,62,65,69,74,77,82,85,86,87,91,93,94,95,106,111,115,118,119,122,123,125,129,133,134,141,142,143,145,146,155,158,159,161,166,177,178,183,185,187,194,201,202,203,205,206,209,213,214,215,217,218,219,221,226,235,237,247,249,253,254,259,262,265,267,274,278,287,291,295,298,299,301,302,303,305,309,314,319,321,323,326,327,329,334,335,339,341,343,346,355,358,362,365,371,377,381,382,386,391,393,394,395,398,403,407,411,413,415,417,422,427,437,445,446,447,451,453,454,458,466,469,471,473,478,481,482,485,489,493,497]
J
Implementation:
factors=: [: /:~@, */&>@{@((^ i.@>:)&.>/)@q:~&__
isMPerfect=: *: = */@(factors ::_:)"0
Task example:
#(#~ isMPerfect) i.500
150
10 15$(#~ isMPerfect) i.500
1 6 8 10 14 15 21 22 26 27 33 34 35 38 39
46 51 55 57 58 62 65 69 74 77 82 85 86 87 91
93 94 95 106 111 115 118 119 122 123 125 129 133 134 141
142 143 145 146 155 158 159 161 166 177 178 183 185 187 194
201 202 203 205 206 209 213 214 215 217 218 219 221 226 235
237 247 249 253 254 259 262 265 267 274 278 287 291 295 298
299 301 302 303 305 309 314 319 321 323 326 327 329 334 335
339 341 343 346 355 358 362 365 371 377 381 382 386 391 393
394 395 398 403 407 411 413 415 417 422 427 437 445 446 447
451 453 454 458 466 469 471 473 478 481 482 485 489 493 497
For the stretch goal, we need to determine the number of semi-primes, given the number of multiplicatively perfect numbers less than N:
adjSemiPrime=: + _1 + %: -&(p:inv) 3&%:
Thus (first number in following results is count of multiplicatively perfect numbers, second is count of semiprimes):
{{ (, adjSemiPrime&y) +/isMPerfect i.y}} 500
150 153
{{ (, adjSemiPrime&y) +/isMPerfect i.y}} 5000
1354 1365
{{ (, adjSemiPrime&y) +/isMPerfect i.y}} 50000
12074 12110
{{ (, adjSemiPrime&y) +/isMPerfect i.y}} 500000
108223 108326
Julia
using Printf
using Primes
""" Find and count multiplicatively perfect integers up to thresholds """
function testmultiplicativelyperfects(thresholds = [500, 5000, 50_000, 500_000])
mpcount, scount = 0, 0
pmask = primesmask(thresholds[end])
println("Multiplicatively perfect numbers under $(thresholds[begin]):")
for n in 1:thresholds[end]
f = factor(n).pe
flen = length(f)
if flen == 2 && f[1][2] == 1 == f[2][2] || flen == 1 && f[1][2] == 3
mpcount += 1
if n < thresholds[begin]
@printf("%3d * %3d = %3d ", f[1][1], n ÷ f[1][1], n)
mpcount % 5 == 0 && println()
end
end
if n in thresholds
cbsum, sqsum = sum(pmask[1:Int(floor(n^(1/3)))]), sum(pmask[1:isqrt(n)])
scount = mpcount - cbsum + sqsum
@printf("\nCounts under %d: MPNs = %7d Semi-primes = %7d\n", n, mpcount, scount)
end
end
end
testmultiplicativelyperfects()
- Output:
Multiplicatively perfect numbers under 500: 2 * 3 = 6 2 * 4 = 8 2 * 5 = 10 2 * 7 = 14 3 * 5 = 15 3 * 7 = 21 2 * 11 = 22 2 * 13 = 26 3 * 9 = 27 3 * 11 = 33 2 * 17 = 34 5 * 7 = 35 2 * 19 = 38 3 * 13 = 39 2 * 23 = 46 3 * 17 = 51 5 * 11 = 55 3 * 19 = 57 2 * 29 = 58 2 * 31 = 62 5 * 13 = 65 3 * 23 = 69 2 * 37 = 74 7 * 11 = 77 2 * 41 = 82 5 * 17 = 85 2 * 43 = 86 3 * 29 = 87 7 * 13 = 91 3 * 31 = 93 2 * 47 = 94 5 * 19 = 95 2 * 53 = 106 3 * 37 = 111 5 * 23 = 115 2 * 59 = 118 7 * 17 = 119 2 * 61 = 122 3 * 41 = 123 5 * 25 = 125 3 * 43 = 129 7 * 19 = 133 2 * 67 = 134 3 * 47 = 141 2 * 71 = 142 11 * 13 = 143 5 * 29 = 145 2 * 73 = 146 5 * 31 = 155 2 * 79 = 158 3 * 53 = 159 7 * 23 = 161 2 * 83 = 166 3 * 59 = 177 2 * 89 = 178 3 * 61 = 183 5 * 37 = 185 11 * 17 = 187 2 * 97 = 194 3 * 67 = 201 2 * 101 = 202 7 * 29 = 203 5 * 41 = 205 2 * 103 = 206 11 * 19 = 209 3 * 71 = 213 2 * 107 = 214 5 * 43 = 215 7 * 31 = 217 2 * 109 = 218 3 * 73 = 219 13 * 17 = 221 2 * 113 = 226 5 * 47 = 235 3 * 79 = 237 13 * 19 = 247 3 * 83 = 249 11 * 23 = 253 2 * 127 = 254 7 * 37 = 259 2 * 131 = 262 5 * 53 = 265 3 * 89 = 267 2 * 137 = 274 2 * 139 = 278 7 * 41 = 287 3 * 97 = 291 5 * 59 = 295 2 * 149 = 298 13 * 23 = 299 7 * 43 = 301 2 * 151 = 302 3 * 101 = 303 5 * 61 = 305 3 * 103 = 309 2 * 157 = 314 11 * 29 = 319 3 * 107 = 321 17 * 19 = 323 2 * 163 = 326 3 * 109 = 327 7 * 47 = 329 2 * 167 = 334 5 * 67 = 335 3 * 113 = 339 11 * 31 = 341 7 * 49 = 343 2 * 173 = 346 5 * 71 = 355 2 * 179 = 358 2 * 181 = 362 5 * 73 = 365 7 * 53 = 371 13 * 29 = 377 3 * 127 = 381 2 * 191 = 382 2 * 193 = 386 17 * 23 = 391 3 * 131 = 393 2 * 197 = 394 5 * 79 = 395 2 * 199 = 398 13 * 31 = 403 11 * 37 = 407 3 * 137 = 411 7 * 59 = 413 5 * 83 = 415 3 * 139 = 417 2 * 211 = 422 7 * 61 = 427 19 * 23 = 437 5 * 89 = 445 2 * 223 = 446 3 * 149 = 447 11 * 41 = 451 3 * 151 = 453 2 * 227 = 454 2 * 229 = 458 2 * 233 = 466 7 * 67 = 469 3 * 157 = 471 11 * 43 = 473 2 * 239 = 478 13 * 37 = 481 2 * 241 = 482 5 * 97 = 485 3 * 163 = 489 17 * 29 = 493 7 * 71 = 497 Counts under 500: MPNs = 149 Semi-primes = 153 Counts under 5000: MPNs = 1353 Semi-primes = 1365 Counts under 50000: MPNs = 12073 Semi-primes = 12110 Counts under 500000: MPNs = 108222 Semi-primes = 108326
Phix
Includes 1, to exclude just start with multiplicatively_perfect_numbers = 0 and r = {}.
with javascript_semantics integer multiplicatively_perfect_numbers = 1, semiprime_numbers = 0, five_e_n = 5e2 sequence r = {1} for n=1 to 5e6 do sequence pn = vslice(prime_powers(n),2) if product(sq_add(pn,1))=4 then multiplicatively_perfect_numbers += 1 if n<=500 then r &= n end if end if if sum(pn)=2 then semiprime_numbers += 1 end if if n=five_e_n then if n=5e2 then printf(1,"%d multiplicatively perfect numbers under 500: %s\n", {length(r),join(shorten(r,"",5,"%d"),",")}) end if printf(1,"Counts under %,9d: MPNs = %,7d Semi-primes = %,7d\n", {five_e_n,multiplicatively_perfect_numbers,semiprime_numbers}) five_e_n *= 10 end if end for
- Output:
150 multiplicatively perfect numbers under 500: 1,6,8,10,14,...,482,485,489,493,497 Counts under 500: MPNs = 150 Semi-primes = 153 Counts under 5,000: MPNs = 1,354 Semi-primes = 1,365 Counts under 50,000: MPNs = 12,074 Semi-primes = 12,110 Counts under 500,000: MPNs = 108,223 Semi-primes = 108,326 Counts under 5,000,000: MPNs = 978,983 Semi-primes = 979,274
PL/M
... under CP/M (or an emulator)
100H: /* FIND MUKTIPLICATIVELY PERFECT NUMBERS - NUMBERS WHOSE PROPER */
/* DIVISOR PRODUCT IS THE NUMBER ITSELF */
/* NOTE THIS IS EQUIVALENT TO FINDING NUMBERS WITH 3 PROPER DIVISORS */
/* SEE OEIS A007422 */
/* CP/M BDOS SYSTEM CALL */
BDOS: PROCEDURE( FN, ARG ); DECLARE FN BYTE, ARG ADDRESS; GOTO 5; END;
/* I/O ROUTINES */
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;
/* FIND THE NUMBERS UP TO 500 */
DECLARE PDC ( 501 )ADDRESS; /* TABLE OF PROPER DIVISOR COUNTS */
DECLARE ( I, J, COUNT ) ADDRESS;
/* COUNT THE PROPER DIVISORS OF NUMBERS TO 500 */
DO I = 0 TO LAST( PDC ); PDC( I ) = 1; END;
DO I = 2 TO LAST( PDC );
DO J = I + I TO LAST( PDC ) BY I;
PDC( J ) = PDC( J ) + 1;
END;
END;
PDC( 1 ) = 3; /* PRETEND 1 HAS 3 PROPER DIVISORS SO IT IS INCLUDED */
/* SHOW YHE MULTIPLICATIVELY PERFECT NUMBERS */
COUNT = 0;
DO I = 1 TO LAST( PDC );
IF PDC( I ) = 3 THEN DO;
CALL PR$CHAR( ' ' );
IF I < 100 THEN DO;
IF I < 10 THEN CALL PR$CHAR( ' ' );
CALL PR$CHAR( ' ' );
END;
CALL PR$NUMBER( I );
IF ( COUNT := COUNT + 1 ) MOD 10 = 0 THEN CALL PR$NL;
END;
END;
EOF- Output:
1 6 8 10 14 15 21 22 26 27 33 34 35 38 39 46 51 55 57 58 62 65 69 74 77 82 85 86 87 91 93 94 95 106 111 115 118 119 122 123 125 129 133 134 141 142 143 145 146 155 158 159 161 166 177 178 183 185 187 194 201 202 203 205 206 209 213 214 215 217 218 219 221 226 235 237 247 249 253 254 259 262 265 267 274 278 287 291 295 298 299 301 302 303 305 309 314 319 321 323 326 327 329 334 335 339 341 343 346 355 358 362 365 371 377 381 382 386 391 393 394 395 398 403 407 411 413 415 417 422 427 437 445 446 447 451 453 454 458 466 469 471 473 478 481 482 485 489 493 497
Raku
use List::Divvy;
use Lingua::EN::Numbers;
constant @primes = (^∞).grep: &is-prime;
constant @cubes = @primes.map: *³;
state $cube = 0;
sub is-mpn(Int $n ) {
return False if $n.is-prime;
if $n == @cubes[$cube] {
++$cube;
return True
}
my $factor = find-factor($n);
return True if ($factor.is-prime && ( my $div = $n div $factor ).is-prime && ($div != $factor));
False;
}
sub find-factor ( Int $n, $constant = 1 ) {
my $x = 2;
my $rho = 1;
my $factor = 1;
while $factor == 1 {
$rho *= 2;
my $fixed = $x;
for ^$rho {
$x = ( $x * $x + $constant ) % $n;
$factor = ( $x - $fixed ) gcd $n;
last if 1 < $factor;
}
}
$factor = find-factor( $n, $constant + 1 ) if $n == $factor;
$factor;
}
constant @mpn = lazy 1, |(2..*).grep: &is-mpn;
say 'Multiplicatively perfect numbers less than 500:';
put @mpn.&upto(500).batch(10)».fmt("%3d").join: "\n";
put "\nThere are:";
for 5e2, 5e3, 5e4, 5e5, 5e6 {
printf "%8s MPNs less than %9s, %7s semiprimes.\n",
comma(my $count = +@mpn.&upto($_)), .Int.&comma,
comma $count + @primes.map(*²).&upto($_) - @cubes.&upto($_) - 1;
}
- Output:
Multiplicatively perfect numbers less than 500:
1 6 8 10 14 15 21 22 26 27
33 34 35 38 39 46 51 55 57 58
62 65 69 74 77 82 85 86 87 91
93 94 95 106 111 115 118 119 122 123
125 129 133 134 141 142 143 145 146 155
158 159 161 166 177 178 183 185 187 194
201 202 203 205 206 209 213 214 215 217
218 219 221 226 235 237 247 249 253 254
259 262 265 267 274 278 287 291 295 298
299 301 302 303 305 309 314 319 321 323
326 327 329 334 335 339 341 343 346 355
358 362 365 371 377 381 382 386 391 393
394 395 398 403 407 411 413 415 417 422
427 437 445 446 447 451 453 454 458 466
469 471 473 478 481 482 485 489 493 497
There are:
150 MPNs less than 500, 153 semiprimes.
1,354 MPNs less than 5,000, 1,365 semiprimes.
12,074 MPNs less than 50,000, 12,110 semiprimes.
108,223 MPNs less than 500,000, 108,326 semiprimes.
978,983 MPNs less than 5,000,000, 979,274 semiprimes.
Ring
see "working..." + nl
see "Special numbers under 500:" + nl
limit = 500
Divisors = []
for n = 1 to limit
pro = 1
Divisors = []
for m = 2 to ceil(n/2)
if n % m = 0
pro = pro * m
add(Divisors,m)
ok
next
str = ""
if n = pro and len(Divisors) > 1
for m = 1 to len(Divisors)
str = str + Divisors[m] + " * "
if m = len(Divisors)
str = left(str,len(str)-2)
ok
next
see "" + n + " = " + str + nl
ok
next
see "done..." + nl- Output:
working... Special numbers under 500: 6 = 2 x 3 8 = 2 x 4 10 = 2 x 5 14 = 2 x 7 15 = 3 x 5 21 = 3 x 7 22 = 2 x 11 26 = 2 x 13 27 = 3 x 9 33 = 3 x 11 34 = 2 x 17 35 = 5 x 7 38 = 2 x 19 39 = 3 x 13 46 = 2 x 23 51 = 3 x 17 55 = 5 x 11 57 = 3 x 19 58 = 2 x 29 62 = 2 x 31 65 = 5 x 13 69 = 3 x 23 74 = 2 x 37 77 = 7 x 11 82 = 2 x 41 85 = 5 x 17 86 = 2 x 43 87 = 3 x 29 91 = 7 x 13 93 = 3 x 31 94 = 2 x 47 95 = 5 x 19 106 = 2 x 53 111 = 3 x 37 115 = 5 x 23 118 = 2 x 59 119 = 7 x 17 122 = 2 x 61 123 = 3 x 41 125 = 5 x 25 129 = 3 x 43 133 = 7 x 19 134 = 2 x 67 141 = 3 x 47 142 = 2 x 71 143 = 11 x 13 145 = 5 x 29 146 = 2 x 73 155 = 5 x 31 158 = 2 x 79 159 = 3 x 53 161 = 7 x 23 166 = 2 x 83 177 = 3 x 59 178 = 2 x 89 183 = 3 x 61 185 = 5 x 37 187 = 11 x 17 194 = 2 x 97 201 = 3 x 67 202 = 2 x 101 203 = 7 x 29 205 = 5 x 41 206 = 2 x 103 209 = 11 x 19 213 = 3 x 71 214 = 2 x 107 215 = 5 x 43 217 = 7 x 31 218 = 2 x 109 219 = 3 x 73 221 = 13 x 17 226 = 2 x 113 235 = 5 x 47 237 = 3 x 79 247 = 13 x 19 249 = 3 x 83 253 = 11 x 23 254 = 2 x 127 259 = 7 x 37 262 = 2 x 131 265 = 5 x 53 267 = 3 x 89 274 = 2 x 137 278 = 2 x 139 287 = 7 x 41 291 = 3 x 97 295 = 5 x 59 298 = 2 x 149 299 = 13 x 23 301 = 7 x 43 302 = 2 x 151 303 = 3 x 101 305 = 5 x 61 309 = 3 x 103 314 = 2 x 157 319 = 11 x 29 321 = 3 x 107 323 = 17 x 19 326 = 2 x 163 327 = 3 x 109 329 = 7 x 47 334 = 2 x 167 335 = 5 x 67 339 = 3 x 113 341 = 11 x 31 343 = 7 x 49 346 = 2 x 173 355 = 5 x 71 358 = 2 x 179 362 = 2 x 181 365 = 5 x 73 371 = 7 x 53 377 = 13 x 29 381 = 3 x 127 382 = 2 x 191 386 = 2 x 193 391 = 17 x 23 393 = 3 x 131 394 = 2 x 197 395 = 5 x 79 398 = 2 x 199 403 = 13 x 31 407 = 11 x 37 411 = 3 x 137 413 = 7 x 59 415 = 5 x 83 417 = 3 x 139 422 = 2 x 211 427 = 7 x 61 437 = 19 x 23 445 = 5 x 89 446 = 2 x 223 447 = 3 x 149 451 = 11 x 41 453 = 3 x 151 454 = 2 x 227 458 = 2 x 229 466 = 2 x 233 469 = 7 x 67 471 = 3 x 157 473 = 11 x 43 478 = 2 x 239 481 = 13 x 37 482 = 2 x 241 485 = 5 x 97 489 = 3 x 163 493 = 17 x 29 497 = 7 x 71 done...
Wren
This includes '1' as an MPN. Not very quick at around 112 seconds but reasonable for Wren.
import "./math" for Int, Nums
import "./fmt" for Fmt
// library method customized for this task
var divisors = Fn.new { |n|
var divisors = []
var i = 1
var k = (n%2 == 0) ? 1 : 2
while (i <= n.sqrt) {
if (i > 1 && n%i == 0) { // exclude 1 and n
divisors.add(i)
if (divisors.count > 2) break // not eligible if has > 2 divisors
var j = (n/i).floor
if (j != i) divisors.add(j)
}
i = i + k
}
return divisors
}
var limit = 500
var count = 0
var i = 1
System.print("Multiplicatively perfect numbers under %(limit):")
while (true) {
var pd = (i != 1) ? divisors.call(i) : [1, 1]
if (pd.count == 2 && pd[0] * pd[1] == i) {
count = count + 1
if (i < 500) {
var pds = Fmt.swrite("$3d x $3d", pd[0], pd[1])
Fmt.write("$3d = $s ", i, pds)
if (count % 5 == 0) System.print()
}
}
if (i == 499) System.print()
if (i >= limit - 1) {
var squares = Int.primeCount((limit - 1).sqrt.floor)
var cubes = Int.primeCount((limit - 1).cbrt.floor)
var count2 = count + squares - cubes - 1
Fmt.print("Counts under $,9d: MPNs = $,7d Semi-primes = $,7d", limit, count, count2)
if (limit == 5000000) return
limit = limit * 10
}
i = i + 1
}- Output:
Multiplicatively perfect numbers under 500: 1 = 1 x 1 6 = 2 x 3 8 = 2 x 4 10 = 2 x 5 14 = 2 x 7 15 = 3 x 5 21 = 3 x 7 22 = 2 x 11 26 = 2 x 13 27 = 3 x 9 33 = 3 x 11 34 = 2 x 17 35 = 5 x 7 38 = 2 x 19 39 = 3 x 13 46 = 2 x 23 51 = 3 x 17 55 = 5 x 11 57 = 3 x 19 58 = 2 x 29 62 = 2 x 31 65 = 5 x 13 69 = 3 x 23 74 = 2 x 37 77 = 7 x 11 82 = 2 x 41 85 = 5 x 17 86 = 2 x 43 87 = 3 x 29 91 = 7 x 13 93 = 3 x 31 94 = 2 x 47 95 = 5 x 19 106 = 2 x 53 111 = 3 x 37 115 = 5 x 23 118 = 2 x 59 119 = 7 x 17 122 = 2 x 61 123 = 3 x 41 125 = 5 x 25 129 = 3 x 43 133 = 7 x 19 134 = 2 x 67 141 = 3 x 47 142 = 2 x 71 143 = 11 x 13 145 = 5 x 29 146 = 2 x 73 155 = 5 x 31 158 = 2 x 79 159 = 3 x 53 161 = 7 x 23 166 = 2 x 83 177 = 3 x 59 178 = 2 x 89 183 = 3 x 61 185 = 5 x 37 187 = 11 x 17 194 = 2 x 97 201 = 3 x 67 202 = 2 x 101 203 = 7 x 29 205 = 5 x 41 206 = 2 x 103 209 = 11 x 19 213 = 3 x 71 214 = 2 x 107 215 = 5 x 43 217 = 7 x 31 218 = 2 x 109 219 = 3 x 73 221 = 13 x 17 226 = 2 x 113 235 = 5 x 47 237 = 3 x 79 247 = 13 x 19 249 = 3 x 83 253 = 11 x 23 254 = 2 x 127 259 = 7 x 37 262 = 2 x 131 265 = 5 x 53 267 = 3 x 89 274 = 2 x 137 278 = 2 x 139 287 = 7 x 41 291 = 3 x 97 295 = 5 x 59 298 = 2 x 149 299 = 13 x 23 301 = 7 x 43 302 = 2 x 151 303 = 3 x 101 305 = 5 x 61 309 = 3 x 103 314 = 2 x 157 319 = 11 x 29 321 = 3 x 107 323 = 17 x 19 326 = 2 x 163 327 = 3 x 109 329 = 7 x 47 334 = 2 x 167 335 = 5 x 67 339 = 3 x 113 341 = 11 x 31 343 = 7 x 49 346 = 2 x 173 355 = 5 x 71 358 = 2 x 179 362 = 2 x 181 365 = 5 x 73 371 = 7 x 53 377 = 13 x 29 381 = 3 x 127 382 = 2 x 191 386 = 2 x 193 391 = 17 x 23 393 = 3 x 131 394 = 2 x 197 395 = 5 x 79 398 = 2 x 199 403 = 13 x 31 407 = 11 x 37 411 = 3 x 137 413 = 7 x 59 415 = 5 x 83 417 = 3 x 139 422 = 2 x 211 427 = 7 x 61 437 = 19 x 23 445 = 5 x 89 446 = 2 x 223 447 = 3 x 149 451 = 11 x 41 453 = 3 x 151 454 = 2 x 227 458 = 2 x 229 466 = 2 x 233 469 = 7 x 67 471 = 3 x 157 473 = 11 x 43 478 = 2 x 239 481 = 13 x 37 482 = 2 x 241 485 = 5 x 97 489 = 3 x 163 493 = 17 x 29 497 = 7 x 71 Counts under 500: MPNs = 150 Semi-primes = 153 Counts under 5,000: MPNs = 1,354 Semi-primes = 1,365 Counts under 50,000: MPNs = 12,074 Semi-primes = 12,110 Counts under 500,000: MPNs = 108,223 Semi-primes = 108,326 Counts under 5,000,000: MPNs = 978,983 Semi-primes = 979,274
XPL0
func Special(N);
int N, D, P;
[D:= 2; P:= 1;
while D < N do
[if rem(N/D) = 0 then P:= P*D;
D:= D+1;
];
return P = N;
];
int N, C;
[C:= 0;
Format(4, 0);
for N:= 2 to 500-1 do
if Special(N) then
[RlOut(0, float(N));
C:= C+1;
if rem(C/20) = 0 then CrLf(0);
];
]- Output:
6 8 10 14 15 21 22 26 27 33 34 35 38 39 46 51 55 57 58 62 65 69 74 77 82 85 86 87 91 93 94 95 106 111 115 118 119 122 123 125 129 133 134 141 142 143 145 146 155 158 159 161 166 177 178 183 185 187 194 201 202 203 205 206 209 213 214 215 217 218 219 221 226 235 237 247 249 253 254 259 262 265 267 274 278 287 291 295 298 299 301 302 303 305 309 314 319 321 323 326 327 329 334 335 339 341 343 346 355 358 362 365 371 377 381 382 386 391 393 394 395 398 403 407 411 413 415 417 422 427 437 445 446 447 451 453 454 458 466 469 471 473 478 481 482 485 489 493 497