Untouchable numbers: Difference between revisions
(→{{header|ALGOL 68}}: Use ALGOL 68-primes) |
(added pascal fast but with problems after 40,000,000 count -> 6,430,223 < wrong by -1) |
||
Line 888: | Line 888: | ||
There are 1212 untouchable numbers ≤ 10000. |
There are 1212 untouchable numbers ≤ 10000. |
||
There are 13863 untouchable numbers ≤ 100000.</pre> |
There are 13863 untouchable numbers ≤ 100000.</pre> |
||
=={{header|Pascal}}== |
|||
modified [[Factors_of_an_integer#using_Prime_decomposition]]<BR> |
|||
see also math.dartmouth.edu/~carlp/uupaper3.pdf with list count up to 1e9 |
|||
<lang pascal>program UntouchableNumbers; |
|||
//limited to 6.75E11 (Max smallprimes^2) |
|||
{$IFDEF FPC} |
|||
{$MODE DELPHI} {$OPTIMIZATION ON,ALL} {$COPERATORS ON} |
|||
{$ELSE} |
|||
{$APPTYPE CONSOLE} |
|||
{$ENDIF} |
|||
uses |
|||
sysutils |
|||
{$IFDEF WINDOWS},Windows{$ENDIF} |
|||
; |
|||
//###################################################################### |
|||
//prime decomposition |
|||
//http://rosettacode.org/wiki/Factors_of_an_integer#using_Prime_decomposition |
|||
const |
|||
SizePrDeFe = 6*8192;//*size of(tprimeFac) =40 byte level I or 2 Mb ~ level 2 cache |
|||
type |
|||
tdigits = array [0..31] of Uint32; |
|||
//the first number with 12 different prime factors = |
|||
//2*3*5*7*11*13*17*19*23*29*31*37 = 7,42E12 |
|||
//40 byte |
|||
tprimeFac = packed record |
|||
pfSumOfDivs, |
|||
pfRemain : Uint64; |
|||
pfpotPrimIdx : array[0..11] of word; |
|||
pfDivCnt : Uint16; |
|||
pfpotMax : array[0..12] of byte; |
|||
pfMaxIdx : byte; |
|||
end; |
|||
tpPrimeFac = ^tprimeFac; |
|||
tPrimeDecompField = array[0..SizePrDeFe-1] of tprimeFac; |
|||
tPrimes = array[0..65535] of Uint32; |
|||
var |
|||
{$ALIGN 8} |
|||
SmallPrimes: tPrimes; |
|||
{$ALIGN 32} |
|||
PrimeDecompField :tPrimeDecompField; |
|||
pdfIDX,pdfOfs: NativeInt; |
|||
procedure InitSmallPrimes; |
|||
//get primes. #0..65535.Sieving only odd numbers |
|||
const |
|||
MAXLIMIT = (821641-1) shr 1; |
|||
var |
|||
pr : array[0..MAXLIMIT] of byte; |
|||
p,j,d,flipflop :NativeUInt; |
|||
Begin |
|||
SmallPrimes[0] := 2; |
|||
fillchar(pr[0],SizeOf(pr),#0); |
|||
p := 0; |
|||
repeat |
|||
repeat |
|||
p +=1 |
|||
until pr[p]= 0; |
|||
j := (p+1)*p*2; |
|||
if j>MAXLIMIT then |
|||
BREAK; |
|||
d := 2*p+1; |
|||
repeat |
|||
pr[j] := 1; |
|||
j += d; |
|||
until j>MAXLIMIT; |
|||
until false; |
|||
SmallPrimes[1] := 3; |
|||
SmallPrimes[2] := 5; |
|||
j := 3; |
|||
d := 7; |
|||
flipflop := (2+1)-1;//7+2*2->11+2*1->13 ,17 ,19 , 23 |
|||
p := 3; |
|||
repeat |
|||
if pr[p] = 0 then |
|||
begin |
|||
SmallPrimes[j] := d; |
|||
inc(j); |
|||
end; |
|||
d += 2*flipflop; |
|||
p+=flipflop; |
|||
flipflop := 3-flipflop; |
|||
until (p > MAXLIMIT) OR (j>High(SmallPrimes)); |
|||
end; |
|||
function smplPrimeDecomp(n:Uint64):tprimeFac; |
|||
var |
|||
pr,i,pot,fac,q :NativeUInt; |
|||
Begin |
|||
with result do |
|||
Begin |
|||
pfDivCnt := 1; |
|||
pfSumOfDivs := 1; |
|||
pfRemain := n; |
|||
pfMaxIdx := 0; |
|||
pfpotPrimIdx[0] := 1; |
|||
pfpotMax[0] := 0; |
|||
i := 0; |
|||
while i < High(SmallPrimes) do |
|||
begin |
|||
pr := SmallPrimes[i]; |
|||
q := n DIV pr; |
|||
//if n < pr*pr |
|||
if pr > q then |
|||
BREAK; |
|||
if n = pr*q then |
|||
Begin |
|||
pfpotPrimIdx[pfMaxIdx] := i; |
|||
pot := 0; |
|||
fac := pr; |
|||
repeat |
|||
n := q; |
|||
q := n div pr; |
|||
pot+=1; |
|||
fac *= pr; |
|||
until n <> pr*q; |
|||
pfpotMax[pfMaxIdx] := pot; |
|||
pfDivCnt *= pot+1; |
|||
pfSumOfDivs *= (fac-1)DIV(pr-1); |
|||
inc(pfMaxIdx); |
|||
end; |
|||
inc(i); |
|||
end; |
|||
pfRemain := n; |
|||
if n > 1 then |
|||
Begin |
|||
pfDivCnt *= 2; |
|||
pfSumOfDivs *= n+1 |
|||
end; |
|||
end; |
|||
end; |
|||
function CnvtoBASE(var dgt:tDigits;n:Uint64;base:NativeUint):NativeInt; |
|||
//n must be multiple of base aka n mod base must be 0 |
|||
var |
|||
q,r: Uint64; |
|||
i : NativeInt; |
|||
Begin |
|||
fillchar(dgt,SizeOf(dgt),#0); |
|||
i := 0; |
|||
n := n div base; |
|||
result := 0; |
|||
repeat |
|||
r := n; |
|||
q := n div base; |
|||
r -= q*base; |
|||
n := q; |
|||
dgt[i] := r; |
|||
inc(i); |
|||
until (q = 0); |
|||
//searching lowest pot in base |
|||
result := 0; |
|||
while (result<i) AND (dgt[result] = 0) do |
|||
inc(result); |
|||
inc(result); |
|||
end; |
|||
function IncByBaseInBase(var dgt:tDigits;base:NativeInt):NativeInt; |
|||
var |
|||
q :NativeInt; |
|||
Begin |
|||
result := 0; |
|||
q := dgt[result]+1; |
|||
if q = base then |
|||
repeat |
|||
dgt[result] := 0; |
|||
inc(result); |
|||
q := dgt[result]+1; |
|||
until q <> base; |
|||
dgt[result] := q; |
|||
result +=1; |
|||
end; |
|||
function SieveOneSieve(var pdf:tPrimeDecompField):boolean; |
|||
var |
|||
dgt:tDigits; |
|||
i,j,k,pr,fac,n,MaxP : NativeUInt; |
|||
begin |
|||
n := pdfOfs; |
|||
if n+SizePrDeFe >= sqr(SmallPrimes[High(SmallPrimes)]) then |
|||
EXIT(FALSE); |
|||
//init |
|||
for i := 0 to SizePrDeFe-1 do |
|||
begin |
|||
with pdf[i] do |
|||
Begin |
|||
pfDivCnt := 1; |
|||
pfSumOfDivs := 1; |
|||
pfRemain := n+i; |
|||
pfMaxIdx := 0; |
|||
pfpotPrimIdx[0] := 0; |
|||
pfpotMax[0] := 0; |
|||
end; |
|||
end; |
|||
//first factor 2. Make n+i even |
|||
i := (pdfIdx+n) AND 1; |
|||
IF (n = 0) AND (pdfIdx<2) then |
|||
i := 2; |
|||
repeat |
|||
with pdf[i] do |
|||
begin |
|||
j := BsfQWord(n+i); |
|||
pfMaxIdx := 1; |
|||
pfpotPrimIdx[0] := 0; |
|||
pfpotMax[0] := j; |
|||
pfRemain := (n+i) shr j; |
|||
pfSumOfDivs := (1 shl (j+1))-1; |
|||
pfDivCnt := j+1; |
|||
end; |
|||
i += 2; |
|||
until i >=SizePrDeFe; |
|||
//i now index in SmallPrimes |
|||
i := 0; |
|||
maxP := trunc(sqrt(n+SizePrDeFe))+1; |
|||
repeat |
|||
//search next prime that is in bounds of sieve |
|||
if n = 0 then |
|||
begin |
|||
repeat |
|||
inc(i); |
|||
pr := SmallPrimes[i]; |
|||
k := pr-n MOD pr; |
|||
if k < SizePrDeFe then |
|||
break; |
|||
until pr > MaxP; |
|||
end |
|||
else |
|||
begin |
|||
repeat |
|||
inc(i); |
|||
pr := SmallPrimes[i]; |
|||
k := pr-n MOD pr; |
|||
if (k = pr) AND (n>0) then |
|||
k:= 0; |
|||
if k < SizePrDeFe then |
|||
break; |
|||
until pr > MaxP; |
|||
end; |
|||
//no need to use higher primes |
|||
if pr*pr > n+SizePrDeFe then |
|||
BREAK; |
|||
//j is power of prime |
|||
j := CnvtoBASE(dgt,n+k,pr); |
|||
repeat |
|||
with pdf[k] do |
|||
Begin |
|||
pfpotPrimIdx[pfMaxIdx] := i; |
|||
pfpotMax[pfMaxIdx] := j; |
|||
pfDivCnt *= j+1; |
|||
fac := pr; |
|||
repeat |
|||
pfRemain := pfRemain DIV pr; |
|||
dec(j); |
|||
fac *= pr; |
|||
until j<= 0; |
|||
pfSumOfDivs *= (fac-1)DIV(pr-1); |
|||
inc(pfMaxIdx); |
|||
k += pr; |
|||
j := IncByBaseInBase(dgt,pr); |
|||
end; |
|||
until k >= SizePrDeFe; |
|||
until false; |
|||
//correct sum of & count of divisors |
|||
for i := 0 to High(pdf) do |
|||
Begin |
|||
with pdf[i] do |
|||
begin |
|||
j := pfRemain; |
|||
if j <> 1 then |
|||
begin |
|||
pfSumOFDivs *= (j+1); |
|||
pfDivCnt *=2; |
|||
end; |
|||
end; |
|||
end; |
|||
result := true; |
|||
end; |
|||
function NextSieve:boolean; |
|||
begin |
|||
dec(pdfIDX,SizePrDeFe); |
|||
inc(pdfOfs,SizePrDeFe); |
|||
result := SieveOneSieve(PrimeDecompField); |
|||
end; |
|||
function GetNextPrimeDecomp:tpPrimeFac; |
|||
begin |
|||
if pdfIDX >= SizePrDeFe then |
|||
if Not(NextSieve) then |
|||
EXIT(NIL); |
|||
result := @PrimeDecompField[pdfIDX]; |
|||
inc(pdfIDX); |
|||
end; |
|||
function Init_Sieve(n:NativeUint):boolean; |
|||
//Init Sieve pdfIdx,pdfOfs are Global |
|||
begin |
|||
pdfIdx := n MOD SizePrDeFe; |
|||
pdfOfs := n-pdfIdx; |
|||
result := SieveOneSieve(PrimeDecompField); |
|||
end; |
|||
const |
|||
LIMIT = 100*1000*1000; |
|||
LIMIT_mul = 6 * (1 shl (trunc(ln(Limit)/ln(10))-2)); |
|||
var |
|||
Untouch : array of byte; |
|||
pUntouch: pByte; |
|||
pPrimeDecomp :tpPrimeFac; |
|||
T0:Int64; |
|||
n,k,lim,prmIdx : NativeInt; |
|||
Begin |
|||
setlength(untouch,LIMIT+1); |
|||
pUntouch := @untouch[0]; |
|||
InitSmallPrimes; |
|||
T0 := GetTickCount64; |
|||
prmIdx := 0; |
|||
REPEAT |
|||
inc(prmIdx); |
|||
until smallprimes[prmIdx]> 2*sqrt(LIMIT); |
|||
writeln(prmIdx:10,smallprimes[prmIdx]:10); |
|||
writeln(LIMIT_mul); |
|||
n := 0; |
|||
Init_Sieve(0); |
|||
repeat |
|||
pPrimeDecomp:= GetNextPrimeDecomp; |
|||
k := pPrimeDecomp^.pfSumOfDivs-n; |
|||
If k <> 1 then |
|||
Begin |
|||
if k > 1 then |
|||
if k <= LIMIT then |
|||
pUntouch[k] := 1; |
|||
end |
|||
else |
|||
begin |
|||
//n= prime, would be marked by n*n with proper factors 1,n |
|||
if n<= LIMIT+1 then |
|||
Begin |
|||
pUntouch[n+1] := 1; |
|||
//marked by prime*n with proper factors 1,(prime),n |
|||
For lim := 0 to prmIdx do |
|||
begin |
|||
k := smallprimes[lim]+n+1; |
|||
If k > LIMIT then |
|||
begin |
|||
prmIdx:= lim-1; |
|||
BREAK; |
|||
end; |
|||
pUntouch[k] := 1; |
|||
end; |
|||
end; |
|||
end; |
|||
inc(n); |
|||
until n >LIMIT_mul*LIMIT; |
|||
T0 := GetTickCount64-T0; |
|||
writeln('runtime ',T0/1000:0:3,' s'); |
|||
k := 0; |
|||
lim :=10; |
|||
for n := 0 to LIMIT DIV 10 do |
|||
Begin |
|||
if n = lim then |
|||
Begin |
|||
writeln(lim:10,k:10); |
|||
lim *= 10; |
|||
end; |
|||
k += 1-pUntouch[n]; |
|||
end; |
|||
lim := 2*LIMIT DIV 10; |
|||
for n := LIMIT DIV 10+1 to LIMIT do |
|||
Begin |
|||
if n = lim then |
|||
Begin |
|||
writeln(lim:10,k:10); |
|||
lim += LIMIT DIV 10; |
|||
end; |
|||
k += 1-pUntouch[n]; |
|||
end; |
|||
end. |
|||
</lang> |
|||
{{out}} |
|||
<pre> |
|||
runtime 2.486 s |
|||
10 2 |
|||
100 5 |
|||
1000 89 |
|||
10000 1212 |
|||
100000 13863 |
|||
200000 28572 |
|||
300000 43514 |
|||
400000 58459 |
|||
500000 73565 |
|||
600000 88828 |
|||
700000 104061 |
|||
800000 119302 |
|||
900000 134757 |
|||
1000000 150232 |
|||
#### |
|||
823 6329 |
|||
192 -> test til 192x10,000,000= 1.92e9 |
|||
runtime 75.065 s |
|||
10 2 |
|||
100 5 |
|||
1000 89 |
|||
10000 1212 |
|||
100000 13863 |
|||
1000000 150232 |
|||
2000000 305290 |
|||
3000000 462110 |
|||
4000000 619638 |
|||
5000000 777672 |
|||
6000000 936243 |
|||
7000000 1095710 |
|||
8000000 1255015 |
|||
9000000 1414783 |
|||
10000000 1574973 |
|||
#### |
|||
2262 20011 |
|||
384 -> test til 384x100,000,000= 38.4e9 |
|||
runtime 1750.728 s |
|||
10 2 |
|||
100 5 |
|||
1000 89 |
|||
10000 1212 |
|||
100000 13863 |
|||
1000000 150232 |
|||
10000000 1574973 |
|||
20000000 3184111 |
|||
30000000 4804331 |
|||
40000000 6430223 < wrong -1 |
|||
50000000 8060162 < wrong -1 |
|||
60000000 9694467 |
|||
70000000 11330312 |
|||
80000000 12967238 < wrong -1 |
|||
90000000 14606549 |
|||
100000000 16246941 < wrong +1 |
|||
real 29m10,824s |
|||
//url=https://math.dartmouth.edu/~carlp/uupaper3.pdf |
|||
6000000 936244 |
|||
7000000 1095710 |
|||
8000000 1255016 |
|||
9000000 1414783 |
|||
10000000 1574973 |
|||
20000000 3184111 |
|||
30000000 4804331 |
|||
40000000 6430224 |
|||
50000000 8060163 |
|||
60000000 9694467 |
|||
70000000 11330312 |
|||
80000000 12967239 |
|||
90000000 14606549 |
|||
100000000 16246940 |
|||
</pre> |
|||
=={{header|Perl}}== |
=={{header|Perl}}== |
||
{{libheader|ntheory}} |
{{libheader|ntheory}} |
Revision as of 20:12, 1 September 2021
You are encouraged to solve this task according to the task description, using any language you may know.
- Definitions
-
- Untouchable numbers are also known as nonaliquot numbers.
- An untouchable number is a positive integer that cannot be expressed as the sum of all the proper divisors of any positive integer. (From Wikipedia)
- The sum of all the proper divisors is also known as the aliquot sum.
- An untouchable are those numbers that are not in the image of the aliquot sum function. (From Wikipedia)
- Untouchable numbers: impossible values for the sum of all aliquot parts function. (From OEIS: The On-line Encyclopedia of Integer Sequences®)
- An untouchable number is a positive integer that is not the sum of the proper divisors of any number. (From MathWorld™)
- Observations and conjectures
All untouchable numbers > 5 are composite numbers.
No untouchable number is perfect.
No untouchable number is sociable.
No untouchable number is a Mersenne prime.
No untouchable number is one more than a prime number, since if p is prime, then the sum of the proper divisors of p2 is p + 1.
No untouchable number is three more than an odd prime number, since if p is an odd prime, then the sum of the proper divisors of 2p is p + 3.
The number 5 is believed to be the only odd untouchable number, but this has not been proven: it would follow from a slightly stronger version of the Goldbach's conjecture, since the sum of the proper divisors of pq (with p, q being distinct primes) is 1 + p + q.
There are infinitely many untouchable numbers, a fact that was proven by Paul Erdős.
According to Chen & Zhao, their natural density is at least d > 0.06.
- Task
-
- show (in a grid format) all untouchable numbers ≤ 2,000.
- show (for the above) the count of untouchable numbers.
- show the count of untouchable numbers from unity up to (inclusive):
- 10
- 100
- 1,000
- 10,000
- 100,000
- ... or as high as is you think is practical.
- all output is to be shown here, on this page.
- See also
-
- Wolfram MathWorld: untouchable number.
- OEIS: A005114 untouchable numbers.
- OEIS: a list of all untouchable numbers below 100,000 (inclusive).
- Wikipedia: untouchable number.
- Wikipedia: Goldbach's conjecture.
ALGOL 68
Generates the proper divisor sums with a sieve-like process.
As noted by the Wren, Go, etc. samples, it is hard to determine what upper limit to use to eliminate non-untouchable numbers. It seems that using the proper divisor sums only, to find untouchables up to limit, limit^2 must be considered.
However, using the observations of the Wren, Go etc. solutions that by using the facts prime + 1 and prime + 3 are not untouchable, and that (probably) 5 is the only odd untouchable, the untouchables up to 1 000 000 can be found by considering numbers up to 64 000 000 (possibly less could be used). At least, this sample finds the same values as the other ones :).
Possibly, this works because the prime + 1 value eliminates the need to calculate the proper divisor sum of p^2 which appears to help a lot when p is close to the upper limit. Why the 64 * limit (or 63 * limit) is valid is unclear (to me, anyway).
Note that under Windows (and possibly under Linux), Algol 68G requires that the heap size be increased in order to allow arrays big enough to handle 100 000 and 1 000 000 untouchable numbers. See ALGOL_68_Genie#Using_a_Large_Heap.
<lang algol68>BEGIN # find some untouchable numbers - numbers not equal to the sum of the #
# proper divisors of any +ve integer # INT max untouchable = 1 000 000; # a table of the untouchable numbers # [ 1 : max untouchable ]BOOL untouchable; FOR i TO UPB untouchable DO untouchable[ i ] := TRUE OD; # show the counts of untouchable numbers found # PROC show untouchable statistics = VOID: BEGIN print( ( "Untouchable numbers:", newline ) ); INT u count := 0; FOR i TO UPB untouchable DO IF untouchable[ i ] THEN u count +:= 1 FI; IF i = 10 OR i = 100 OR i = 1 000 OR i = 10 000 OR i = 100 000 OR i = 1 000 000 THEN print( ( whole( u count, -7 ), " to ", whole( i, -8 ), newline ) ) FI OD END; # show untouchable counts # # prints the untouchable numbers up to n # PROC print untouchables = ( INT n )VOID: BEGIN print( ( "Untouchable numbers up to ", whole( n, 0 ), newline ) ); INT u count := 0; FOR i TO n DO IF untouchable[ i ] THEN print( ( whole( i, -4 ) ) ); IF u count +:= 1; u count MOD 16 = 0 THEN print( ( newline ) ) ELSE print( ( " " ) ) FI FI OD; print( ( newline ) ); print( ( whole( u count, -7 ), " to ", whole( n, -8 ), newline ) ) END; # print untouchables # # find the untouchable numbers # # to find untouchable numbers up to e.g.: 10 000, we need to sieve up to # # 10 000 ^2 i.e. 100 000 000 # # however if we also use the facts that no untouchable = prime + 1 # # and no untouchable = odd prime + 3 and 5 is (very probably) the only # # odd untouchable, other samples suggest we can use limit * 64 to find # # untlouchables up to 1 000 000 - experimentation reveals this to be true # # assume the conjecture that there are no odd untouchables except 5 # BEGIN untouchable[ 1 ] := FALSE; untouchable[ 3 ] := FALSE; FOR i FROM 7 BY 2 TO UPB untouchable DO untouchable[ i ] := FALSE OD END; # sieve the primes to max untouchable and flag the non untouchables # BEGIN PR read "primes.incl.a68" PR []BOOL prime = PRIMESIEVE max untouchable; FOR i FROM 3 BY 2 TO UPB prime DO IF prime[ i ] THEN IF i < max untouchable THEN untouchable[ i + 1 ] := FALSE; IF i < ( max untouchable - 2 ) THEN untouchable[ i + 3 ] := FALSE FI FI FI OD; untouchable[ 2 + 1 ] := FALSE # special case for the only even prime # END; # construct the proper divisor sums and flag the non untouchables # BEGIN [ 1 : max untouchable * 64 ]INT spd; FOR i TO UPB spd DO spd[ i ] := 1 OD; FOR i FROM 2 TO UPB spd DO FOR j FROM i + i BY i TO UPB spd DO spd[ j ] +:= i OD OD; FOR i TO UPB spd DO IF spd[ i ] <= UPB untouchable THEN untouchable[ spd[ i ] ] := FALSE FI OD END; # show the untouchable numbers up to 2000 # print untouchables( 2 000 ); # show the counts of untouchable numbers # show untouchable statistics
END</lang>
- Output:
2 5 52 88 96 120 124 146 162 188 206 210 216 238 246 248 262 268 276 288 290 292 304 306 322 324 326 336 342 372 406 408 426 430 448 472 474 498 516 518 520 530 540 552 556 562 576 584 612 624 626 628 658 668 670 708 714 718 726 732 738 748 750 756 766 768 782 784 792 802 804 818 836 848 852 872 892 894 896 898 902 926 934 936 964 966 976 982 996 1002 1028 1044 1046 1060 1068 1074 1078 1080 1102 1116 1128 1134 1146 1148 1150 1160 1162 1168 1180 1186 1192 1200 1212 1222 1236 1246 1248 1254 1256 1258 1266 1272 1288 1296 1312 1314 1316 1318 1326 1332 1342 1346 1348 1360 1380 1388 1398 1404 1406 1418 1420 1422 1438 1476 1506 1508 1510 1522 1528 1538 1542 1566 1578 1588 1596 1632 1642 1650 1680 1682 1692 1716 1718 1728 1732 1746 1758 1766 1774 1776 1806 1816 1820 1822 1830 1838 1840 1842 1844 1852 1860 1866 1884 1888 1894 1896 1920 1922 1944 1956 1958 1960 1962 1972 1986 1992 196 to 2000 Untouchable numbers: 2 to 10 5 to 100 89 to 1000 1212 to 10000 13863 to 100000 150232 to 1000000
C++
This solution implements Talk:Untouchable_numbers#Nice_recursive_solution
The Function
<lang cpp> // Untouchable Numbers : Nigel Galloway - March 4th., 2021;
- include <functional>
- include <bitset>
- include <iostream>
- include <cmath>
using namespace std; using Z0=long long; using Z1=optional<Z0>; using Z2=optional<array<int,3>>; using Z3=function<Z2()>; const int maxUT{3000000}, dL{(int)log2(maxUT)}; struct uT{
bitset<maxUT+1>N; vector<int> G{}; array<Z3,int(dL+1)>L{Z3{}}; int sG{0},mUT{}; void _g(int n,int g){if(g<=mUT){N[g]=false; return _g(n,n+g);}} Z1 nxt(const int n){if(n>mUT) return Z1{}; if(N[n]) return Z1(n); return nxt(n+1);} Z3 fN(const Z0 n,const Z0 i,int g){return [=]()mutable{if(g<sG && ((n+i)*(1+G[g])-n*G[g]<=mUT)) return Z2Template:N,i,g++; return Z2{};};} Z3 fG(Z0 n,Z0 i,const int g){Z0 e{n+i},l{1},p{1}; return [=]()mutable{n=n*G[g]; p=p*G[g]; l=l+p; i=e*l-n; if(i<=mUT) return Z2Template:N,i,g; return Z2{};};} void fL(Z3 n, int g){for(;;){ if(auto i=n()){N[(*i)[1]]=false; L[g+1]=fN((*i)[0],(*i)[1],(*i)[2]+1); g=g+1; continue;} if(auto i=L[g]()){n=fG((*i)[0],(*i)[1],(*i)[2]); continue;} if(g>0) if(auto i=L[g-1]()){ g=g-1; n=fG((*i)[0],(*i)[1],(*i)[2]); continue;} if(g>0){ n=[](){return Z2{};}; g=g-1; continue;} break;} } int count(){int g{0}; for(auto n=nxt(0); n; n=nxt(*n+1)) ++g; return g;} uT(const int n):mUT{n}{ N.set(); N[0]=false; N[1]=false; for(auto n=nxt(0);*n<=sqrt(mUT);n=nxt(*n+1)) _g(*n,*n+*n); for(auto n=nxt(0); n; n=nxt(*n+1)) G.push_back(*n); sG=G.size(); N.set(); N[0]=false; L[0]=fN(1,0,0); fL([](){return Z2{};},0); }
}; </lang>
The Task
- Less than 2000
<lang cpp> int main(int argc, char *argv[]) {
int c{0}; auto n{uT{2000}}; for(auto g=n.nxt(0); g; g=n.nxt(*g+1)){if(c++==30){c=1; printf("\n");} printf("%4d ",*g);} printf("\n");
} </lang>
- Output:
2 5 52 88 96 120 124 146 162 188 206 210 216 238 246 248 262 268 276 288 290 292 304 306 322 324 326 336 342 372 406 408 426 430 448 472 474 498 516 518 520 530 540 552 556 562 576 584 612 624 626 628 658 668 670 708 714 718 726 732 738 748 750 756 766 768 782 784 792 802 804 818 836 848 852 872 892 894 896 898 902 926 934 936 964 966 976 982 996 1002 1028 1044 1046 1060 1068 1074 1078 1080 1102 1116 1128 1134 1146 1148 1150 1160 1162 1168 1180 1186 1192 1200 1212 1222 1236 1246 1248 1254 1256 1258 1266 1272 1288 1296 1312 1314 1316 1318 1326 1332 1342 1346 1348 1360 1380 1388 1398 1404 1406 1418 1420 1422 1438 1476 1506 1508 1510 1522 1528 1538 1542 1566 1578 1588 1596 1632 1642 1650 1680 1682 1692 1716 1718 1728 1732 1746 1758 1766 1774 1776 1806 1816 1820 1822 1830 1838 1840 1842 1844 1852 1860 1866 1884 1888 1894 1896 1920 1922 1944 1956 1958 1960 1962 1972 1986 1992
- Count less than or equal 100000
<lang cpp> int main(int argc, char *argv[]) {
int z{100000}; auto n{uT{z}}; cout<<"untouchables below "<<z<<"->"<<n.count()<<endl;
} </lang>
- Output:
untouchables below 100000->13863 real 0m2.928s
- Count less than or equal 1000000
<lang cpp> int main(int argc, char *argv[]) {
int z{1000000}; auto n{uT{z}}; cout<<"untouchables below "<<z<<"->"<<n.count()<<endl;
} </lang>
- Output:
untouchables below 1000000->150232 real 3m6.909s
- Count less than or equal 2000000
<lang cpp> int main(int argc, char *argv[]) {
int z{2000000}; auto n{uT{z}}; cout<<"untouchables below "<<z<<"->"<<n.count()<<endl;
} </lang>
- Output:
untouchables below 2000000->305290 real 11m28.031s
Delphi
<lang Delphi> program Untouchable_numbers;
{$APPTYPE CONSOLE}
uses
System.SysUtils;
function SumDivisors(n: Integer): Integer; begin
Result := 1; var k := 2; if not odd(n) then k := 1; var i := 1 + k; while i * i <= n do begin if (n mod i) = 0 then begin inc(Result, i); var j := n div i; if j <> i then inc(Result, j); end; inc(i, k); end;
end;
function Sieve(n: Integer): TArray<Boolean>; begin
inc(n); SetLength(result, n + 1); for var i := 6 to n do begin var sd := SumDivisors(i); if sd <= n then result[sd] := True; end;
end;
function PrimeSieve(limit: Integer): TArray<Boolean>; begin
inc(limit); SetLength(result, limit); Result[0] := True; Result[1] := True;
var p := 3; repeat var p2 := p * p; if p2 >= limit then Break; var i := p2; while i < limit do begin
Result[i] := True; inc(i, 2 * p); end;
repeat inc(p, 2); until not Result[p];
until (False);
end;
function Commatize(n: Double): string; var
fmt: TFormatSettings;
begin
fmt := TFormatSettings.Create('en-US'); Result := n.ToString(ffNumber, 64, 0, fmt);
end;
begin
var limit := 1000000; var c := primeSieve(limit); var s := sieve(63 * limit); var untouchable: TArray<Integer> := [2, 5]; var n := 6; while n <= limit do begin if not s[n] and c[n - 1] and c[n - 3] then begin SetLength(untouchable, Length(untouchable) + 1); untouchable[High(untouchable)] := n; end; inc(n, 2); end; writeln('List of untouchable numbers <= 2,000:'); var count := 0; var i := 0; while untouchable[i] <= 2000 do begin write(commatize(untouchable[i]): 6); if ((i + 1) mod 10) = 0 then writeln; inc(i); end; writeln(#10#10, commatize(count): 7, ' untouchable numbers were found <= 2,000');
var p := 10; count := 0; for n in untouchable do begin inc(count); if n > p then begin var cc := commatize(count - 1); var cp := commatize(p); writeln(cc, ' untouchable numbers were found <= ', cp); p := p * 10; if p = limit then Break; end; end;
var cu := commatize(Length(untouchable)); var cl := commatize(limit); writeln(cu:7, ' untouchable numbers were found <= ', cl); readln;
end.</lang>
F#
The Function
This task uses Extensible Prime Generator (F#).
It implements Talk:Untouchable_numbers#Nice_recursive_solution
<lang fsharp>
// Applied dendrology. Nigel Galloway: February 15., 2021
let uT a=let N,G=Array.create(a+1) true, [|yield! primes64()|>Seq.takeWhile((>)(int64 a))|]
let fN n i e=let mutable p=e-1 in (fun()->p<-p+1; if p<G.Length && (n+i)*(1L+G.[p])-n*G.[p]<=(int64 a) then Some(n,i,p) else None) let fG n i e=let g=n+i in let mutable n,l,p=n,1L,1L (fun()->n<-n*G.[e]; p<-p*G.[e]; l<-l+p; let i=g*l-n in if i<=(int64 a) then Some(n,i,e) else None) let rec fL n g=match n() with Some(f,i,e)->N.[(int i)]<-false; fL n ((fN f i (e+1))::g) |_->match g with n::t->match n() with Some (n,i,e)->fL (fG n i e) g |_->fL n t |_->N.[0]<-false; N fL (fG 1L 0L 0) [fN 1L 0L 1]
</lang>
The Task
- Less than 2000
<lang fsharp> uT 2000|>Array.mapi(fun n g->(n,g))|>Array.filter(fun(_,n)->n)|>Array.chunkBySize 30|>Array.iter(fun n->n|>Array.iter(fst>>printf "%5d");printfn "") </lang>
- Output:
2 5 52 88 96 120 124 146 162 188 206 210 216 238 246 248 262 268 276 288 290 292 304 306 322 324 326 336 342 372 406 408 426 430 448 472 474 498 516 518 520 530 540 552 556 562 576 584 612 624 626 628 658 668 670 708 714 718 726 732 738 748 750 756 766 768 782 784 792 802 804 818 836 848 852 872 892 894 896 898 902 926 934 936 964 966 976 982 996 1002 1028 1044 1046 1060 1068 1074 1078 1080 1102 1116 1128 1134 1146 1148 1150 1160 1162 1168 1180 1186 1192 1200 1212 1222 1236 1246 1248 1254 1256 1258 1266 1272 1288 1296 1312 1314 1316 1318 1326 1332 1342 1346 1348 1360 1380 1388 1398 1404 1406 1418 1420 1422 1438 1476 1506 1508 1510 1522 1528 1538 1542 1566 1578 1588 1596 1632 1642 1650 1680 1682 1692 1716 1718 1728 1732 1746 1758 1766 1774 1776 1806 1816 1820 1822 1830 1838 1840 1842 1844 1852 1860 1866 1884 1888 1894 1896 1920 1922 1944 1956 1958 1960 1962 1972 1986 1992
- Count less than or equal 100000
<lang fsharp> printfn "%d" (uT 100000|>Array.filter id|>Array.length) </lang>
- Output:
13863 Real: 00:00:02.784, CPU: 00:00:02.750
- Count less than or equal 1000000
<lang fsharp> printfn "%d" (uT 1000000|>Array.filter id|>Array.length) </lang>
- Output:
150232 Real: 00:03:08.929, CPU: 00:03:08.859
- Count less than or equal 2000000
<lang fsharp> printfn "%d" (uT 2000000|>Array.filter id|>Array.length) </lang>
- Output:
305290 Real: 00:11:22.325, CPU: 00:11:21.828
- Count less than or equal 3000000
<lang fsharp> <lang fsharp> printfn "%d" (uT 3000000|>Array.filter id|>Array.length) </lang>
- Output:
462110 Real: 00:24:00.126, CPU: 00:23:59.203
Go
This was originally based on the Wren example but has been modified somewhat to find untouchable numbers up to 1 million rather than 100,000 in a reasonable time. On my machine, the former (with a sieve factor of 63) took 31 minutes 9 seconds and the latter (with a sieve factor of 14) took 6.2 seconds. <lang go>package main
import "fmt"
func sumDivisors(n int) int {
sum := 1 k := 2 if n%2 == 0 { k = 1 } for i := 1 + k; i*i <= n; i += k { if n%i == 0 { sum += i j := n / i if j != i { sum += j } } } return sum
}
func sieve(n int) []bool {
n++ s := make([]bool, n+1) // all false by default for i := 6; i <= n; i++ { sd := sumDivisors(i) if sd <= n { s[sd] = true } } return s
}
func primeSieve(limit int) []bool {
limit++ // True denotes composite, false denotes prime. c := make([]bool, limit) // all false by default c[0] = true c[1] = true // no need to bother with even numbers over 2 for this task p := 3 // Start from 3. for { p2 := p * p if p2 >= limit { break } for i := p2; i < limit; i += 2 * p { c[i] = true } for { p += 2 if !c[p] { break } } } return c
}
func commatize(n int) string {
s := fmt.Sprintf("%d", n) if n < 0 { s = s[1:] } le := len(s) for i := le - 3; i >= 1; i -= 3 { s = s[0:i] + "," + s[i:] } if n >= 0 { return s } return "-" + s
}
func main() {
limit := 1000000 c := primeSieve(limit) s := sieve(63 * limit) untouchable := []int{2, 5} for n := 6; n <= limit; n += 2 { if !s[n] && c[n-1] && c[n-3] { untouchable = append(untouchable, n) } } fmt.Println("List of untouchable numbers <= 2,000:") count := 0 for i := 0; untouchable[i] <= 2000; i++ { fmt.Printf("%6s", commatize(untouchable[i])) if (i+1)%10 == 0 { fmt.Println() } count++ } fmt.Printf("\n\n%7s untouchable numbers were found <= 2,000\n", commatize(count)) p := 10 count = 0 for _, n := range untouchable { count++ if n > p { cc := commatize(count - 1) cp := commatize(p) fmt.Printf("%7s untouchable numbers were found <= %9s\n", cc, cp) p = p * 10 if p == limit { break } } } cu := commatize(len(untouchable)) cl := commatize(limit) fmt.Printf("%7s untouchable numbers were found <= %s\n", cu, cl)
}</lang>
- Output:
List of untouchable numbers <= 2,000: 2 5 52 88 96 120 124 146 162 188 206 210 216 238 246 248 262 268 276 288 290 292 304 306 322 324 326 336 342 372 406 408 426 430 448 472 474 498 516 518 520 530 540 552 556 562 576 584 612 624 626 628 658 668 670 708 714 718 726 732 738 748 750 756 766 768 782 784 792 802 804 818 836 848 852 872 892 894 896 898 902 926 934 936 964 966 976 982 996 1,002 1,028 1,044 1,046 1,060 1,068 1,074 1,078 1,080 1,102 1,116 1,128 1,134 1,146 1,148 1,150 1,160 1,162 1,168 1,180 1,186 1,192 1,200 1,212 1,222 1,236 1,246 1,248 1,254 1,256 1,258 1,266 1,272 1,288 1,296 1,312 1,314 1,316 1,318 1,326 1,332 1,342 1,346 1,348 1,360 1,380 1,388 1,398 1,404 1,406 1,418 1,420 1,422 1,438 1,476 1,506 1,508 1,510 1,522 1,528 1,538 1,542 1,566 1,578 1,588 1,596 1,632 1,642 1,650 1,680 1,682 1,692 1,716 1,718 1,728 1,732 1,746 1,758 1,766 1,774 1,776 1,806 1,816 1,820 1,822 1,830 1,838 1,840 1,842 1,844 1,852 1,860 1,866 1,884 1,888 1,894 1,896 1,920 1,922 1,944 1,956 1,958 1,960 1,962 1,972 1,986 1,992 196 untouchable numbers were found <= 2,000 2 untouchable numbers were found <= 10 5 untouchable numbers were found <= 100 89 untouchable numbers were found <= 1,000 1,212 untouchable numbers were found <= 10,000 13,863 untouchable numbers were found <= 100,000 150,232 untouchable numbers were found <= 1,000,000
J
<lang J> factor=: 3 : 0 NB. explicit
'primes powers'=. __&q: y input_to_cartesian_product=. primes ^&.> i.&.> >: powers cartesian_product=. , { input_to_cartesian_product , */&> cartesian_product
)
factor=: [: , [: */&> [: { [: (^&.> i.&.>@>:)/ __&q: NB. tacit
proper_divisors=: [: }: factor
sum_of_proper_divisors=: +/@proper_divisors
candidates=: 5 , [: +: [: #\@i. >.@-: NB. within considered range, all but one candidate are even. spds=:([:sum_of_proper_divisors"0(#\@i.-.i.&.:(p:inv))@*:)f. NB. remove primes which contribute 1 </lang> We may revisit to correct the larger untouchable tallies. The straightforward algorithm presented is a little bit efficient, and, I claim, the verb (candidates-.spds) produces the full list of untouchable numbers up to its argument. It considers the sum of proper divisors through the argument squared, less primes. Since J is an algorithm description language, it may be "fairer" to state in J that "more resources required" than in some other language. [1]
Time (seconds) and space (bytes) on a high end six year old (new in 2015) laptop computer.
timespacex'(-.candidates@#)/:~~.spds 10000' 600.285 4.29497e9
UNTOUCHABLES=:(candidates-.spds) 2000 NB. compute untouchable numbers /:~factor#UNTOUCHABLES NB. look for nice display size 1 2 4 7 14 28 49 98 196 _14[\UNTOUCHABLES 5 2 52 88 96 120 124 146 162 188 206 210 216 238 246 248 262 268 276 288 290 292 304 306 322 324 326 336 342 372 406 408 426 430 448 472 474 498 516 518 520 530 540 552 556 562 576 584 612 624 626 628 658 668 670 708 714 718 726 732 738 748 750 756 766 768 782 784 792 802 804 818 836 848 852 872 892 894 896 898 902 926 934 936 964 966 976 982 996 1002 1028 1044 1046 1060 1068 1074 1078 1080 1102 1116 1128 1134 1146 1148 1150 1160 1162 1168 1180 1186 1192 1200 1212 1222 1236 1246 1248 1254 1256 1258 1266 1272 1288 1296 1312 1314 1316 1318 1326 1332 1342 1346 1348 1360 1380 1388 1398 1404 1406 1418 1420 1422 1438 1476 1506 1508 1510 1522 1528 1538 1542 1566 1578 1588 1596 1632 1642 1650 1680 1682 1692 1716 1718 1728 1732 1746 1758 1766 1774 1776 1806 1816 1820 1822 1830 1838 1840 1842 1844 1852 1860 1866 1884 1888 1894 1896 1920 1922 1944 1956 1958 1960 1962 1972 1986 1992 SBD=:spds 10000 NB. sums of proper divisors of many integers T=:/:~~.SBD NB. sort the nub #T NB. leaving this many touchable numbers, which are dense through at least 10000 47269787 U=:(-.~candidates@#)T NB. possible untouchable numbers head=: 'n';'tally of possible untouchable numbers to n' head,:,.L:0;/|:(,.U&I.)10^#\i.5 NB. testing to ten thousand squared ┌──────┬──────────────────────────────────────────┐ │n │tally of possible untouchable numbers to n│ ├──────┼──────────────────────────────────────────┤ │ 10│ 2 │ │ 100│ 5 │ │ 1000│ 89 │ │ 10000│ 1212 │ │100000│17538 │ └──────┴──────────────────────────────────────────┘
Julia
I can prove that the number to required to sieve to assure only untouchables for the interval 1:N is less than (N/2 - 1)^2 for larger N, but the 512,000,000 sieved below is just from doubling 1,000,000 and running the sieve until we get 150232 for the number of untouchables under 1,000,000. <lang julia>using Primes
function properfactorsum(n)
f = [one(n)] for (p,e) in factor(n) f = reduce(vcat, [f*p^j for j in 1:e], init=f) end pop!(f) return sum(f)
end
const maxtarget, sievelimit = 1_000_000, 512_000_000 const untouchables = ones(Bool, maxtarget)
for i in 2:sievelimit
n = properfactorsum(i) if n <= maxtarget untouchables[n] = false end
end for i in 6:maxtarget
if untouchables[i] && (isprime(i - 1) || isprime(i - 3)) untouchables[i] = false end
end
println("The untouchable numbers ≤ 2000 are: ") for (i, n) in enumerate(filter(x -> untouchables[x], 1:2000))
print(rpad(n, 5), i % 10 == 0 || i == 196 ? "\n" : "")
end for N in [2000, 10, 100, 1000, 10_000, 100_000, 1_000_000]
println("The count of untouchable numbers ≤ $N is: ", count(x -> untouchables[x], 1:N))
end
</lang>
- Output:
The untouchable numbers ≤ 2000 are: 2 5 52 88 96 120 124 146 162 188 206 210 216 238 246 248 262 268 276 288 290 292 304 306 322 324 326 336 342 372 406 408 426 430 448 472 474 498 516 518 520 530 540 552 556 562 576 584 612 624 626 628 658 668 670 708 714 718 726 732 738 748 750 756 766 768 782 784 792 802 804 818 836 848 852 872 892 894 896 898 902 926 934 936 964 966 976 982 996 1002 1028 1044 1046 1060 1068 1074 1078 1080 1102 1116 1128 1134 1146 1148 1150 1160 1162 1168 1180 1186 1192 1200 1212 1222 1236 1246 1248 1254 1256 1258 1266 1272 1288 1296 1312 1314 1316 1318 1326 1332 1342 1346 1348 1360 1380 1388 1398 1404 1406 1418 1420 1422 1438 1476 1506 1508 1510 1522 1528 1538 1542 1566 1578 1588 1596 1632 1642 1650 1680 1682 1692 1716 1718 1728 1732 1746 1758 1766 1774 1776 1806 1816 1820 1822 1830 1838 1840 1842 1844 1852 1860 1866 1884 1888 1894 1896 1920 1922 1944 1956 1958 1960 1962 1972 1986 1992 The count of untouchable numbers ≤ 2000 is: 196 The count of untouchable numbers ≤ 10 is: 2 The count of untouchable numbers ≤ 100 is: 5 The count of untouchable numbers ≤ 1000 is: 89 The count of untouchable numbers ≤ 10000 is: 1212 The count of untouchable numbers ≤ 100000 is: 13863 The count of untouchable numbers ≤ 1000000 is: 150232
Mathematica/Wolfram Language
<lang Mathematica>f = DivisorSigma[1, #] - # &; limit = 10^5; c = Not /@ PrimeQ[Range[limit]]; slimit = 15 limit; s = ConstantArray[False, slimit + 1]; untouchable = {2, 5}; Do[
val = f[i]; If[val <= slimit, sval = True ] , {i, 6, slimit}
] Do[
If[! sn, If[cn - 1, If[cn - 3, AppendTo[untouchable, n] ] ] ] , {n, 6, limit, 2}
] Multicolumn[Select[untouchable, LessEqualThan[2000]]] Count[untouchable, _?(LessEqualThan[2000])] Count[untouchable, _?(LessEqualThan[10])] Count[untouchable, _?(LessEqualThan[100])] Count[untouchable, _?(LessEqualThan[1000])] Count[untouchable, _?(LessEqualThan[10000])] Count[untouchable, _?(LessEqualThan[100000])]</lang>
- Output:
2 246 342 540 714 804 964 1102 1212 1316 1420 1596 1774 1884 5 248 372 552 718 818 966 1116 1222 1318 1422 1632 1776 1888 52 262 406 556 726 836 976 1128 1236 1326 1438 1642 1806 1894 88 268 408 562 732 848 982 1134 1246 1332 1476 1650 1816 1896 96 276 426 576 738 852 996 1146 1248 1342 1506 1680 1820 1920 120 288 430 584 748 872 1002 1148 1254 1346 1508 1682 1822 1922 124 290 448 612 750 892 1028 1150 1256 1348 1510 1692 1830 1944 146 292 472 624 756 894 1044 1160 1258 1360 1522 1716 1838 1956 162 304 474 626 766 896 1046 1162 1266 1380 1528 1718 1840 1958 188 306 498 628 768 898 1060 1168 1272 1388 1538 1728 1842 1960 206 322 516 658 782 902 1068 1180 1288 1398 1542 1732 1844 1962 210 324 518 668 784 926 1074 1186 1296 1404 1566 1746 1852 1972 216 326 520 670 792 934 1078 1192 1312 1406 1578 1758 1860 1986 238 336 530 708 802 936 1080 1200 1314 1418 1588 1766 1866 1992 196 2 5 89 1212 13863
Nim
I borrowed some ideas from Go version, but keep the limit to 100_000 as in Wren version. <lang Nim>import math, strutils
const
Lim1 = 100_000 # Limit for untouchable numbers. Lim2 = 14 * Lim1 # Limit for computation of sum of divisors.
proc sumdiv(n: uint): uint =
## Return the sum of the strict divisors of "n". result = 1 let r = sqrt(n.float).uint let k = if (n and 1) == 0: 1u else: 2u for d in countup(k + 1, r, k): if n mod d == 0: result += d let q = n div d if q != d: result += q
var
isSumDiv: array[1..Lim2, bool] isPrime: array[1..Lim1, bool]
- Fill both sieves in a single pass.
for n in 1u..Lim2:
let s = sumdiv(n) if s <= Lim2: isSumDiv[s] = true if s == 1 and n <= Lim1: isPrime[n] = true
isPrime[1] = false
- Build list of untouchable numbers.
var list = @[2, 5] for n in countup(6, Lim1, 2):
if not (isSumDiv[n] or isPrime[n - 1] or isPrime[n - 3]): list.add n
echo "Untouchable numbers ≤ 2000:" var count, lcount = 0 for n in list:
if n <= 2000: stdout.write ($n).align(5) inc count inc lcount if lcount == 20: echo() lcount = 0 else: if lcount > 0: echo() break
const CountMessage = "There are $1 untouchable numbers ≤ $2." echo CountMessage.format(count, 2000), '\n'
count = 0 var lim = 10 for n in list:
if n > lim: echo CountMessage.format(count, lim) lim *= 10 inc count
if lim == Lim1:
# Emit last message. echo CountMessage.format(count, lim)</lang>
- Output:
Untouchable numbers ≤ 2000: 2 5 52 88 96 120 124 146 162 188 206 210 216 238 246 248 262 268 276 288 290 292 304 306 322 324 326 336 342 372 406 408 426 430 448 472 474 498 516 518 520 530 540 552 556 562 576 584 612 624 626 628 658 668 670 708 714 718 726 732 738 748 750 756 766 768 782 784 792 802 804 818 836 848 852 872 892 894 896 898 902 926 934 936 964 966 976 982 996 1002 1028 1044 1046 1060 1068 1074 1078 1080 1102 1116 1128 1134 1146 1148 1150 1160 1162 1168 1180 1186 1192 1200 1212 1222 1236 1246 1248 1254 1256 1258 1266 1272 1288 1296 1312 1314 1316 1318 1326 1332 1342 1346 1348 1360 1380 1388 1398 1404 1406 1418 1420 1422 1438 1476 1506 1508 1510 1522 1528 1538 1542 1566 1578 1588 1596 1632 1642 1650 1680 1682 1692 1716 1718 1728 1732 1746 1758 1766 1774 1776 1806 1816 1820 1822 1830 1838 1840 1842 1844 1852 1860 1866 1884 1888 1894 1896 1920 1922 1944 1956 1958 1960 1962 1972 1986 1992 There are 196 untouchable numbers ≤ 2000. There are 2 untouchable numbers ≤ 10. There are 5 untouchable numbers ≤ 100. There are 89 untouchable numbers ≤ 1000. There are 1212 untouchable numbers ≤ 10000. There are 13863 untouchable numbers ≤ 100000.
Pascal
modified Factors_of_an_integer#using_Prime_decomposition
see also math.dartmouth.edu/~carlp/uupaper3.pdf with list count up to 1e9
<lang pascal>program UntouchableNumbers;
//limited to 6.75E11 (Max smallprimes^2)
{$IFDEF FPC}
{$MODE DELPHI} {$OPTIMIZATION ON,ALL} {$COPERATORS ON}
{$ELSE}
{$APPTYPE CONSOLE}
{$ENDIF} uses
sysutils
{$IFDEF WINDOWS},Windows{$ENDIF}
;
//###################################################################### //prime decomposition //http://rosettacode.org/wiki/Factors_of_an_integer#using_Prime_decomposition const
SizePrDeFe = 6*8192;//*size of(tprimeFac) =40 byte level I or 2 Mb ~ level 2 cache
type
tdigits = array [0..31] of Uint32; //the first number with 12 different prime factors = //2*3*5*7*11*13*17*19*23*29*31*37 = 7,42E12 //40 byte tprimeFac = packed record pfSumOfDivs, pfRemain : Uint64; pfpotPrimIdx : array[0..11] of word; pfDivCnt : Uint16; pfpotMax : array[0..12] of byte; pfMaxIdx : byte; end; tpPrimeFac = ^tprimeFac;
tPrimeDecompField = array[0..SizePrDeFe-1] of tprimeFac; tPrimes = array[0..65535] of Uint32;
var
{$ALIGN 8} SmallPrimes: tPrimes; {$ALIGN 32} PrimeDecompField :tPrimeDecompField; pdfIDX,pdfOfs: NativeInt;
procedure InitSmallPrimes; //get primes. #0..65535.Sieving only odd numbers const
MAXLIMIT = (821641-1) shr 1;
var
pr : array[0..MAXLIMIT] of byte; p,j,d,flipflop :NativeUInt;
Begin
SmallPrimes[0] := 2; fillchar(pr[0],SizeOf(pr),#0); p := 0; repeat repeat p +=1 until pr[p]= 0; j := (p+1)*p*2; if j>MAXLIMIT then BREAK; d := 2*p+1; repeat pr[j] := 1; j += d; until j>MAXLIMIT; until false;
SmallPrimes[1] := 3; SmallPrimes[2] := 5; j := 3; d := 7; flipflop := (2+1)-1;//7+2*2->11+2*1->13 ,17 ,19 , 23 p := 3; repeat if pr[p] = 0 then begin SmallPrimes[j] := d; inc(j); end; d += 2*flipflop; p+=flipflop; flipflop := 3-flipflop; until (p > MAXLIMIT) OR (j>High(SmallPrimes));
end;
function smplPrimeDecomp(n:Uint64):tprimeFac; var
pr,i,pot,fac,q :NativeUInt;
Begin
with result do Begin pfDivCnt := 1; pfSumOfDivs := 1; pfRemain := n; pfMaxIdx := 0; pfpotPrimIdx[0] := 1; pfpotMax[0] := 0;
i := 0; while i < High(SmallPrimes) do begin pr := SmallPrimes[i]; q := n DIV pr; //if n < pr*pr if pr > q then BREAK; if n = pr*q then Begin pfpotPrimIdx[pfMaxIdx] := i; pot := 0; fac := pr; repeat n := q; q := n div pr; pot+=1; fac *= pr; until n <> pr*q; pfpotMax[pfMaxIdx] := pot; pfDivCnt *= pot+1; pfSumOfDivs *= (fac-1)DIV(pr-1); inc(pfMaxIdx); end; inc(i); end; pfRemain := n; if n > 1 then Begin pfDivCnt *= 2; pfSumOfDivs *= n+1 end; end;
end;
function CnvtoBASE(var dgt:tDigits;n:Uint64;base:NativeUint):NativeInt; //n must be multiple of base aka n mod base must be 0 var
q,r: Uint64; i : NativeInt;
Begin
fillchar(dgt,SizeOf(dgt),#0); i := 0; n := n div base; result := 0; repeat r := n; q := n div base; r -= q*base; n := q; dgt[i] := r; inc(i); until (q = 0); //searching lowest pot in base result := 0; while (result<i) AND (dgt[result] = 0) do inc(result); inc(result);
end;
function IncByBaseInBase(var dgt:tDigits;base:NativeInt):NativeInt; var
q :NativeInt;
Begin
result := 0; q := dgt[result]+1; if q = base then repeat dgt[result] := 0; inc(result); q := dgt[result]+1; until q <> base; dgt[result] := q; result +=1;
end;
function SieveOneSieve(var pdf:tPrimeDecompField):boolean; var
dgt:tDigits; i,j,k,pr,fac,n,MaxP : NativeUInt;
begin
n := pdfOfs; if n+SizePrDeFe >= sqr(SmallPrimes[High(SmallPrimes)]) then EXIT(FALSE); //init for i := 0 to SizePrDeFe-1 do begin with pdf[i] do Begin pfDivCnt := 1; pfSumOfDivs := 1; pfRemain := n+i; pfMaxIdx := 0; pfpotPrimIdx[0] := 0; pfpotMax[0] := 0; end; end; //first factor 2. Make n+i even i := (pdfIdx+n) AND 1; IF (n = 0) AND (pdfIdx<2) then i := 2;
repeat with pdf[i] do begin j := BsfQWord(n+i); pfMaxIdx := 1; pfpotPrimIdx[0] := 0; pfpotMax[0] := j; pfRemain := (n+i) shr j; pfSumOfDivs := (1 shl (j+1))-1; pfDivCnt := j+1; end; i += 2; until i >=SizePrDeFe; //i now index in SmallPrimes i := 0; maxP := trunc(sqrt(n+SizePrDeFe))+1; repeat //search next prime that is in bounds of sieve if n = 0 then begin repeat inc(i); pr := SmallPrimes[i]; k := pr-n MOD pr; if k < SizePrDeFe then break; until pr > MaxP; end else begin repeat inc(i); pr := SmallPrimes[i]; k := pr-n MOD pr; if (k = pr) AND (n>0) then k:= 0; if k < SizePrDeFe then break; until pr > MaxP; end;
//no need to use higher primes if pr*pr > n+SizePrDeFe then BREAK;
//j is power of prime j := CnvtoBASE(dgt,n+k,pr); repeat with pdf[k] do Begin pfpotPrimIdx[pfMaxIdx] := i; pfpotMax[pfMaxIdx] := j; pfDivCnt *= j+1; fac := pr; repeat pfRemain := pfRemain DIV pr; dec(j); fac *= pr; until j<= 0; pfSumOfDivs *= (fac-1)DIV(pr-1); inc(pfMaxIdx); k += pr; j := IncByBaseInBase(dgt,pr); end; until k >= SizePrDeFe; until false;
//correct sum of & count of divisors for i := 0 to High(pdf) do Begin with pdf[i] do begin j := pfRemain; if j <> 1 then begin pfSumOFDivs *= (j+1); pfDivCnt *=2; end; end; end; result := true;
end;
function NextSieve:boolean; begin
dec(pdfIDX,SizePrDeFe); inc(pdfOfs,SizePrDeFe); result := SieveOneSieve(PrimeDecompField);
end;
function GetNextPrimeDecomp:tpPrimeFac; begin
if pdfIDX >= SizePrDeFe then if Not(NextSieve) then EXIT(NIL); result := @PrimeDecompField[pdfIDX]; inc(pdfIDX);
end;
function Init_Sieve(n:NativeUint):boolean; //Init Sieve pdfIdx,pdfOfs are Global begin
pdfIdx := n MOD SizePrDeFe; pdfOfs := n-pdfIdx; result := SieveOneSieve(PrimeDecompField);
end;
const
LIMIT = 100*1000*1000; LIMIT_mul = 6 * (1 shl (trunc(ln(Limit)/ln(10))-2));
var
Untouch : array of byte; pUntouch: pByte; pPrimeDecomp :tpPrimeFac; T0:Int64; n,k,lim,prmIdx : NativeInt;
Begin
setlength(untouch,LIMIT+1); pUntouch := @untouch[0]; InitSmallPrimes; T0 := GetTickCount64; prmIdx := 0; REPEAT inc(prmIdx); until smallprimes[prmIdx]> 2*sqrt(LIMIT); writeln(prmIdx:10,smallprimes[prmIdx]:10);
writeln(LIMIT_mul); n := 0; Init_Sieve(0); repeat pPrimeDecomp:= GetNextPrimeDecomp; k := pPrimeDecomp^.pfSumOfDivs-n; If k <> 1 then Begin if k > 1 then if k <= LIMIT then pUntouch[k] := 1; end else begin //n= prime, would be marked by n*n with proper factors 1,n if n<= LIMIT+1 then Begin pUntouch[n+1] := 1; //marked by prime*n with proper factors 1,(prime),n For lim := 0 to prmIdx do begin k := smallprimes[lim]+n+1; If k > LIMIT then begin prmIdx:= lim-1; BREAK; end; pUntouch[k] := 1; end; end; end;
inc(n); until n >LIMIT_mul*LIMIT; T0 := GetTickCount64-T0; writeln('runtime ',T0/1000:0:3,' s'); k := 0; lim :=10; for n := 0 to LIMIT DIV 10 do Begin if n = lim then Begin writeln(lim:10,k:10); lim *= 10; end; k += 1-pUntouch[n]; end; lim := 2*LIMIT DIV 10; for n := LIMIT DIV 10+1 to LIMIT do Begin if n = lim then Begin writeln(lim:10,k:10); lim += LIMIT DIV 10; end; k += 1-pUntouch[n]; end;
end. </lang>
- Output:
runtime 2.486 s 10 2 100 5 1000 89 10000 1212 100000 13863 200000 28572 300000 43514 400000 58459 500000 73565 600000 88828 700000 104061 800000 119302 900000 134757 1000000 150232 #### 823 6329 192 -> test til 192x10,000,000= 1.92e9 runtime 75.065 s 10 2 100 5 1000 89 10000 1212 100000 13863 1000000 150232 2000000 305290 3000000 462110 4000000 619638 5000000 777672 6000000 936243 7000000 1095710 8000000 1255015 9000000 1414783 10000000 1574973 #### 2262 20011 384 -> test til 384x100,000,000= 38.4e9 runtime 1750.728 s 10 2 100 5 1000 89 10000 1212 100000 13863 1000000 150232 10000000 1574973 20000000 3184111 30000000 4804331 40000000 6430223 < wrong -1 50000000 8060162 < wrong -1 60000000 9694467 70000000 11330312 80000000 12967238 < wrong -1 90000000 14606549 100000000 16246941 < wrong +1 real 29m10,824s //url=https://math.dartmouth.edu/~carlp/uupaper3.pdf 6000000 936244 7000000 1095710 8000000 1255016 9000000 1414783 10000000 1574973 20000000 3184111 30000000 4804331 40000000 6430224 50000000 8060163 60000000 9694467 70000000 11330312 80000000 12967239 90000000 14606549 100000000 16246940
Perl
<lang perl>use strict; use warnings; use enum qw(False True); use ntheory qw/divisor_sum is_prime/;
sub sieve {
my($n) = @_; my %s; for my $k (0 .. $n+1) { my $sum = divisor_sum($k) - $k; $s{$sum} = True if $sum <= $n+1; } %s
}
my(%s,%c); my($max, $limit, $cnt) = (2000, 1e5, 0);
%s = sieve 14 * $limit; !is_prime($_) and $c{$_} = True for 1..$limit; my @untouchable = (2, 5); for ( my $n = 6; $n <= $limit; $n += 2 ) {
push @untouchable, $n if !$s{$n} and $c{$n-1} and $c{$n-3};
} map { $cnt++ if $_ <= $max } @untouchable; print "Number of untouchable numbers ≤ $max : $cnt \n\n" .
(sprintf "@{['%6d' x $cnt]}", @untouchable[0..$cnt-1]) =~ s/(.{84})/$1\n/gr . "\n";
my($p, $count) = (10, 0); my $fmt = "%6d untouchable numbers were found ≤ %7d\n"; for my $n (@untouchable) {
$count++; if ($n > $p) { printf $fmt, $count-1, $p; printf($fmt, scalar @untouchable, $limit) and last if $limit == ($p *= 10) }
}</lang>
- Output:
Number of untouchable numbers ≤ 2000 : 196 2 5 52 88 96 120 124 146 162 188 206 210 216 238 246 248 262 268 276 288 290 292 304 306 322 324 326 336 342 372 406 408 426 430 448 472 474 498 516 518 520 530 540 552 556 562 576 584 612 624 626 628 658 668 670 708 714 718 726 732 738 748 750 756 766 768 782 784 792 802 804 818 836 848 852 872 892 894 896 898 902 926 934 936 964 966 976 982 996 1002 1028 1044 1046 1060 1068 1074 1078 1080 1102 1116 1128 1134 1146 1148 1150 1160 1162 1168 1180 1186 1192 1200 1212 1222 1236 1246 1248 1254 1256 1258 1266 1272 1288 1296 1312 1314 1316 1318 1326 1332 1342 1346 1348 1360 1380 1388 1398 1404 1406 1418 1420 1422 1438 1476 1506 1508 1510 1522 1528 1538 1542 1566 1578 1588 1596 1632 1642 1650 1680 1682 1692 1716 1718 1728 1732 1746 1758 1766 1774 1776 1806 1816 1820 1822 1830 1838 1840 1842 1844 1852 1860 1866 1884 1888 1894 1896 1920 1922 1944 1956 1958 1960 1962 1972 1986 1992 2 untouchable numbers were found ≤ 10 5 untouchable numbers were found ≤ 100 89 untouchable numbers were found ≤ 1000 1212 untouchable numbers were found ≤ 10000 13863 untouchable numbers were found ≤ 100000
Phix
constant limz = {1,1,8,9,18,64} -- found by experiment procedure untouchable(integer n, cols=0, tens=0) atom t0 = time(), t1 = t0+1 bool tell = n>0 n = abs(n) sequence sums = repeat(0,n+3) for i=1 to n do integer p = get_prime(i) if p>n then exit end if sums[p+1] = 1 sums[p+3] = 1 end for sums[5] = 0 integer m = floor(log10(n)) integer lim = limz[m]*n for j=2 to lim do integer y = sum(factors(j,-1)) if y<=n then sums[y] = 1 end if if platform()!=JS and time()>t1 then progress("j:%,d/%,d (%3.2f%%)\r",{j,lim,(j/lim)*100}) t1 = time()+1 end if end for if platform()!=JS then progress("") end if if tell then printf(1,"The list of all untouchable numbers <= %d:\n",{n}) end if string line = " 2 5" integer cnt = 2 for t=6 to n by 2 do if sums[t]=0 then cnt += 1 if tell then line &= sprintf("%,8d",t) if remainder(cnt,cols)=0 then printf(1,"%s\n",line) line = "" end if end if end if end for if tell then if line!="" then printf(1,"%s\n",line) end if printf(1,"\n") end if string t = elapsed(time()-t0,1," (%s)") printf(1,"%,20d untouchable numbers were found <= %,d%s\n",{cnt,n,t}) for p=1 to tens do untouchable(-power(10,p)) end for end procedure untouchable(2000, 10, 6-(platform()==JS))
- Output:
The list of all untouchable numbers <= 2000: 2 5 52 88 96 120 124 146 162 188 206 210 216 238 246 248 262 268 276 288 290 292 304 306 322 324 326 336 342 372 406 408 426 430 448 472 474 498 516 518 520 530 540 552 556 562 576 584 612 624 626 628 658 668 670 708 714 718 726 732 738 748 750 756 766 768 782 784 792 802 804 818 836 848 852 872 892 894 896 898 902 926 934 936 964 966 976 982 996 1,002 1,028 1,044 1,046 1,060 1,068 1,074 1,078 1,080 1,102 1,116 1,128 1,134 1,146 1,148 1,150 1,160 1,162 1,168 1,180 1,186 1,192 1,200 1,212 1,222 1,236 1,246 1,248 1,254 1,256 1,258 1,266 1,272 1,288 1,296 1,312 1,314 1,316 1,318 1,326 1,332 1,342 1,346 1,348 1,360 1,380 1,388 1,398 1,404 1,406 1,418 1,420 1,422 1,438 1,476 1,506 1,508 1,510 1,522 1,528 1,538 1,542 1,566 1,578 1,588 1,596 1,632 1,642 1,650 1,680 1,682 1,692 1,716 1,718 1,728 1,732 1,746 1,758 1,766 1,774 1,776 1,806 1,816 1,820 1,822 1,830 1,838 1,840 1,842 1,844 1,852 1,860 1,866 1,884 1,888 1,894 1,896 1,920 1,922 1,944 1,956 1,958 1,960 1,962 1,972 1,986 1,992 196 untouchable numbers were found <= 2,000 2 untouchable numbers were found <= 10 6 untouchable numbers were found <= 100 89 untouchable numbers were found <= 1,000 1,212 untouchable numbers were found <= 10,000 13,863 untouchable numbers were found <= 100,000 (12.4s) 150,232 untouchable numbers were found <= 1,000,000 (33 minutes and 55s)
for comparison, on the same box, the Julia entry took 58 minutes and 40s.
Raku
Borrow the proper divisors routine from here.
<lang perl6># 20210220 Raku programming solution
sub propdiv (\x) {
my @l = 1 if x > 1; (2 .. x.sqrt.floor).map: -> \d { unless x % d { @l.push: d; my \y = x div d; @l.push: y if y != d } } @l
}
sub sieve (\n) {
my %s; for (0..(n+1)) -> \k { given ( [+] propdiv k ) { %s{$_} = True if $_ ≤ (n+1) } } %s;
}
my \limit = 1e5; my %c = ( grep { !.is-prime }, 1..limit ).Set; # store composites my %s = sieve(14 * limit); my @untouchable = 2, 5;
loop ( my \n = $ = 6 ; n ≤ limit ; n += 2 ) {
@untouchable.append(n) if (!%s{n} && %c{n-1} && %c{n-3})
}
my ($c, $last) = 0, False; for @untouchable.rotor(10) {
say [~] @_».&{$c++ ; $_ > 2000 ?? ( $last = True and last ) !! .fmt: "%6d "} $c-- and last if $last
}
say "\nList of untouchable numbers ≤ 2,000 : $c \n";
my ($p, $count) = 10,0; BREAK: for @untouchable -> \n {
$count++; if (n > $p) { printf "%6d untouchable numbers were found ≤ %7d\n", $count-1, $p; last BREAK if limit == ($p *= 10) }
} printf "%6d untouchable numbers were found ≤ %7d\n", +@untouchable, limit</lang>
- Output:
2 5 52 88 96 120 124 146 162 188 206 210 216 238 246 248 262 268 276 288 290 292 304 306 322 324 326 336 342 372 406 408 426 430 448 472 474 498 516 518 520 530 540 552 556 562 576 584 612 624 626 628 658 668 670 708 714 718 726 732 738 748 750 756 766 768 782 784 792 802 804 818 836 848 852 872 892 894 896 898 902 926 934 936 964 966 976 982 996 1002 1028 1044 1046 1060 1068 1074 1078 1080 1102 1116 1128 1134 1146 1148 1150 1160 1162 1168 1180 1186 1192 1200 1212 1222 1236 1246 1248 1254 1256 1258 1266 1272 1288 1296 1312 1314 1316 1318 1326 1332 1342 1346 1348 1360 1380 1388 1398 1404 1406 1418 1420 1422 1438 1476 1506 1508 1510 1522 1528 1538 1542 1566 1578 1588 1596 1632 1642 1650 1680 1682 1692 1716 1718 1728 1732 1746 1758 1766 1774 1776 1806 1816 1820 1822 1830 1838 1840 1842 1844 1852 1860 1866 1884 1888 1894 1896 1920 1922 1944 1956 1958 1960 1962 1972 1986 1992 List of untouchable numbers ≤ 2,000 : 196 2 untouchable numbers were found ≤ 10 5 untouchable numbers were found ≤ 100 89 untouchable numbers were found ≤ 1000 1212 untouchable numbers were found ≤ 10000 13863 untouchable numbers were found ≤ 100000
REXX
Some optimization was done to the code to produce prime numbers, since a simple test could be made to exclude
the calculation of touchability for primes as the aliquot sum of a prime is always unity.
This saved around 15% of the running time.
A fair amount of code was put into the generation of primes, but the computation of the aliquot sum is the area
that consumes the most CPU time.
The REXX code below will accurately calculate untouchable numbers up to and including 100,000. Beyond that,
a higher overreach (the over option) would be need to specified.
This version of REXX would need a 64-bit version to calculate the number of untouchable numbers for one million. <lang rexx>/*REXX pgm finds N untouchable numbers (numbers that can't be equal to any aliquot sum).*/ parse arg n cols tens over . /*obtain optional arguments from the CL*/ if n= | n=="," then n=2000 /*Not specified? Then use the default.*/ if cols= | cols=="," | cols==0 then cols= 10 /* " " " " " " */ if tens= | tens=="," then tens= 0 /* " " " " " " */ if over= | over=="," then over= 20 /* " " " " " " */ tell= n>0; n= abs(n) /*N>0? Then display the untouchable #s*/ call genP n * over /*call routine to generate some primes.*/ u.= 0 /*define all possible aliquot sums ≡ 0.*/
do p=1 for #; _= @.p + 1; u._= 1 /*any prime+1 is not an untouchable.*/ _= @.p + 3; u._= 1 /* " prime+3 " " " " */ end /*p*/ /* [↑] this will also rule out 5. */
u.5= 0 /*special case as prime 2 + 3 sum to 5.*/
do j=2 for lim; if !.j then iterate /*Is J a prime? Yes, then skip it. */ y= sigmaP() /*compute: aliquot sum (sigma P) of J.*/ if y<=n then u.y= 1 /*mark Y as a touchable if in range. */ end /*j*/
call show /*maybe show untouchable #s and a count*/ if tens>0 then call powers /*Any "tens" specified? Calculate 'em.*/ exit cnt /*stick a fork in it, we're all done. */ /*──────────────────────────────────────────────────────────────────────────────────────*/ commas: parse arg ?; do jc=length(?)-3 to 1 by -3; ?=insert(',', ?, jc); end; return ? genSq: do _=1 until _*_>lim; q._= _*_; end; q._= _*_; _= _+1; q._= _*_; return grid: $= $ right( commas(t), w); if cnt//cols==0 then do; say $; $=; end; return powers: do pr=1 for tens; call 'UNTOUCHA' -(10**pr); end /*recurse*/; return /*──────────────────────────────────────────────────────────────────────────────────────*/ genP: #= 9; @.1=2; @.2=3; @.3=5; @.4=7; @.5=11; @.6=13; @.7=17; @.8=19; @.9=23 /*a list*/
!.=0; !.2=1; !.3=1; !.5=1; !.7=1; !.11=1; !.13=1; !.17=1; !.19=1 !.23=1 /*primes*/ parse arg lim; call genSq /*define the (high) limit for searching*/ qq.10= 100 /*define square of the 10th prime index*/ do j=@.#+6 by 2 to lim /*find odd primes from here on forward.*/ parse var j -1 _; if _==5 then iterate; if j// 3==0 then iterate if j// 7==0 then iterate; if j//11==0 then iterate; if j//13==0 then iterate if j//17==0 then iterate; if j//19==0 then iterate; if j//23==0 then iterate /*start dividing by the tenth prime: 29*/ do k=10 while qq.k <= j /* [↓] divide J by known odd primes.*/ if j//@.k==0 then iterate j /*J ÷ by a prime? Then ¬prime. ___ */ end /*k*/ /* [↑] only process numbers ≤ √ J */ #= #+1; @.#= j /*bump prime count; assign a new prime.*/ !.j= 1; qq.#= j*j /*mark prime; compute square of prime.*/ end /*j*/; return /*#: is the number of primes generated*/
/*──────────────────────────────────────────────────────────────────────────────────────*/ show: w=7; $= right(2, w+1) right(5, w) /*start the list of an even prime and 5*/
cnt= 2 /*count of the only two primes in list.*/ do t=6 by 2 to n; if u.t then iterate /*Is T touchable? Then skip it. */ cnt= cnt + 1; if tell then call grid /*bump count; maybe show a grid line. */ end /*t*/ if tell & $\== then say $ /*display a residual grid line, if any.*/ if tell then say /*show a spacing blank line for output.*/ if n>0 then say right( commas(cnt), 20) , /*indent the output a bit.*/ ' untouchable numbers were found ≤ ' commas(n); return
/*──────────────────────────────────────────────────────────────────────────────────────*/ sigmaP: s= 1 /*set initial sigma sum (S) to 1. ___*/
if j//2 then do m=3 by 2 while q.m<j /*divide by odd integers up to the √ J */ if j//m==0 then s=s+m+j%m /*add the two divisors to the sum. */ end /*m*/ /* [↑] process an odd integer. ___*/ else do m=2 while q.m<j /*divide by all integers up to the √ J */ if j//m==0 then s=s+m+j%m /*add the two divisors to the sum. */ end /*m*/ /* [↑] process an even integer. ___*/ if q.m==j then return s + m /*Was J a square? If so, add √ J */ return s /* No, just return. */</lang>
- output when using the default inputs:
2 5 52 88 96 120 124 146 162 188 206 210 216 238 246 248 262 268 276 288 290 292 304 306 322 324 326 336 342 372 406 408 426 430 448 472 474 498 516 518 520 530 540 552 556 562 576 584 612 624 626 628 658 668 670 708 714 718 726 732 738 748 750 756 766 768 782 784 792 802 804 818 836 848 852 872 892 894 896 898 902 926 934 936 964 966 976 982 996 1,002 1,028 1,044 1,046 1,060 1,068 1,074 1,078 1,080 1,102 1,116 1,128 1,134 1,146 1,148 1,150 1,160 1,162 1,168 1,180 1,186 1,192 1,200 1,212 1,222 1,236 1,246 1,248 1,254 1,256 1,258 1,266 1,272 1,288 1,296 1,312 1,314 1,316 1,318 1,326 1,332 1,342 1,346 1,348 1,360 1,380 1,388 1,398 1,404 1,406 1,418 1,420 1,422 1,438 1,476 1,506 1,508 1,510 1,522 1,528 1,538 1,542 1,566 1,578 1,588 1,596 1,632 1,642 1,650 1,680 1,682 1,692 1,716 1,718 1,728 1,732 1,746 1,758 1,766 1,774 1,776 1,806 1,816 1,820 1,822 1,830 1,838 1,840 1,842 1,844 1,852 1,860 1,866 1,884 1,888 1,894 1,896 1,920 1,922 1,944 1,956 1,958 1,960 1,962 1,972 1,986 1,992 196 untouchable numbers were found ≤ 2,000
- output when using the inputs: 0 , 5
2 untouchable numbers were found ≤ 10 5 untouchable numbers were found ≤ 100 89 untouchable numbers were found ≤ 1,000 1,212 untouchable numbers were found ≤ 10,000 13,863 untouchable numbers were found ≤ 100,000
Wren
Not an easy task as, even allowing for the prime tests, it's difficult to know how far you need to sieve to get the right answers. The parameters here were found by trial and error. <lang ecmascript>import "/math" for Int, Nums import "/seq" for Lst import "/fmt" for Fmt
var sieve = Fn.new { |n|
n = n + 1 var s = List.filled(n+1, false) for (i in 0..n) { var sum = Nums.sum(Int.properDivisors(i)) if (sum <= n) s[sum] = true } return s
}
var limit = 1e5 var c = Int.primeSieve(limit, false) var s = sieve.call(14 * limit) var untouchable = [2, 5] var n = 6 while (n <= limit) {
if (!s[n] && c[n-1] && c[n-3]) untouchable.add(n) n = n + 2
}
System.print("List of untouchable numbers <= 2,000:") for (chunk in Lst.chunks(untouchable.where { |n| n <= 2000 }.toList, 10)) {
Fmt.print("$,6d", chunk)
} System.print() Fmt.print("$,6d untouchable numbers were found <= 2,000", untouchable.count { |n| n <= 2000 }) var p = 10 var count = 0 for (n in untouchable) {
count = count + 1 if (n > p) { Fmt.print("$,6d untouchable numbers were found <= $,7d", count-1, p) p = p * 10 if (p == limit) break }
} Fmt.print("$,6d untouchable numbers were found <= $,d", untouchable.count, limit)</lang>
- Output:
List of untouchable numbers <= 2,000: 2 5 52 88 96 120 124 146 162 188 206 210 216 238 246 248 262 268 276 288 290 292 304 306 322 324 326 336 342 372 406 408 426 430 448 472 474 498 516 518 520 530 540 552 556 562 576 584 612 624 626 628 658 668 670 708 714 718 726 732 738 748 750 756 766 768 782 784 792 802 804 818 836 848 852 872 892 894 896 898 902 926 934 936 964 966 976 982 996 1,002 1,028 1,044 1,046 1,060 1,068 1,074 1,078 1,080 1,102 1,116 1,128 1,134 1,146 1,148 1,150 1,160 1,162 1,168 1,180 1,186 1,192 1,200 1,212 1,222 1,236 1,246 1,248 1,254 1,256 1,258 1,266 1,272 1,288 1,296 1,312 1,314 1,316 1,318 1,326 1,332 1,342 1,346 1,348 1,360 1,380 1,388 1,398 1,404 1,406 1,418 1,420 1,422 1,438 1,476 1,506 1,508 1,510 1,522 1,528 1,538 1,542 1,566 1,578 1,588 1,596 1,632 1,642 1,650 1,680 1,682 1,692 1,716 1,718 1,728 1,732 1,746 1,758 1,766 1,774 1,776 1,806 1,816 1,820 1,822 1,830 1,838 1,840 1,842 1,844 1,852 1,860 1,866 1,884 1,888 1,894 1,896 1,920 1,922 1,944 1,956 1,958 1,960 1,962 1,972 1,986 1,992 196 untouchable numbers were found <= 2,000 2 untouchable numbers were found <= 10 5 untouchable numbers were found <= 100 89 untouchable numbers were found <= 1,000 1,212 untouchable numbers were found <= 10,000 13,863 untouchable numbers were found <= 100,000