Rosetta Code

Untouchable numbers: Difference between revisions

Page Discussion

Untouchable numbers (view source)

Revision as of 18:03, 26 March 2024

63,768 bytes added ,  1 month ago
Added FreeBASIC
(Added FreeBASIC)
 
(46 intermediate revisions by 13 users not shown)
Line 1:
{{draft task|Prime Numbers}}
 
;Definitions:
Line 16:
 
;Observations and conjectures:
All untouchable numbers &nbsp; <big>&gt;</big>&nbsp; '''5'''&nbsp; are composite numbers.
 
No untouchable number is perfect.
Line 61:
:* &nbsp; Wikipedia: [https://en.wikipedia.org/wiki/Goldbach%27s_conjecture Goldbach's conjecture].
<br><br>
 
=={{header|ALGOL 68}}==
Generates the proper divisor sums with a sieve-like process.
<br>
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.
<br>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 :).
<br>
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).
<br><br>
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]].
{{libheader|ALGOL 68-primes}}
<syntaxhighlight 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</syntaxhighlight>
{{out}}
<pre>
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
</pre>
 
=={{header|C}}==
Run time around 14 seconds on my Core i7 machine.
<syntaxhighlight lang="c">#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <locale.h>
 
bool *primeSieve(int limit) {
int i, p;
limit++;
// True denotes composite, false denotes prime.
bool *c = calloc(limit, sizeof(bool)); // all false by default
c[0] = true;
c[1] = true;
for (i = 4; i < limit; i += 2) c[i] = true;
p = 3; // Start from 3.
while (true) {
int p2 = p * p;
if (p2 >= limit) break;
for (i = p2; i < limit; i += 2 * p) c[i] = true;
while (true) {
p += 2;
if (!c[p]) break;
}
}
return c;
}
 
int main() {
const int limit = 1000000;
int i, j, n, uc = 2, p = 10, m = 63, ul = 151000;
bool *c = primeSieve(limit);
n = m * limit + 1;
int *sumDivs = (int *)calloc(n, sizeof(int));
for (i = 1; i < n; ++i) {
for (j = i; j < n; j += i) sumDivs[j] += i;
}
bool *s = (bool *)calloc(n, sizeof(bool)); // all false
for (i = 1; i < n; ++i) {
int sum = sumDivs[i] - i; // proper divs sum
if (sum <= n) s[sum] = true;
}
free(sumDivs);
int *untouchable = (int *)malloc(ul * sizeof(int));
untouchable[0] = 2;
untouchable[1] = 5;
for (n = 6; n <= limit; n += 2) {
if (!s[n] && c[n-1] && c[n-3]) untouchable[uc++] = n;
}
setlocale(LC_NUMERIC, "");
printf("List of untouchable numbers <= 2,000:\n");
for (i = 0; i < uc; ++i) {
j = untouchable[i];
if (j > 2000) break;
printf("%'6d ", j);
if (!((i+1) % 10)) printf("\n");
}
printf("\n\n%'7d untouchable numbers were found <= 2,000\n", i);
for (i = 0; i < uc; ++i) {
j = untouchable[i];
if (j > p) {
printf("%'7d untouchable numbers were found <= %'9d\n", i, p);
p *= 10;
if (p == limit) break;
}
}
printf("%'7d untouchable numbers were found <= %'d\n", uc, limit);
free(c);
free(s);
free(untouchable);
return 0;
}</syntaxhighlight>
 
{{out}}
<pre>
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
</pre>
 
=={{header|C++}}==
This solution implements [[Talk:Untouchable_numbers#Nice_recursive_solution]]
===The Function===
<syntaxhighlight 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 Z2{{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 Z2{{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);
}
};
</syntaxhighlight>
===The Task===
;Less than 2000
<syntaxhighlight 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");
}
</syntaxhighlight>
{{out}}
<pre>
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
</pre>
;Count less than or equal 100000
<syntaxhighlight lang="cpp">
int main(int argc, char *argv[]) {
int z{100000}; auto n{uT{z}}; cout<<"untouchables below "<<z<<"->"<<n.count()<<endl;
}
</syntaxhighlight>
{{out}}
<pre>
untouchables below 100000->13863
real 0m2.928s
</pre>
;Count less than or equal 1000000
<syntaxhighlight lang="cpp">
int main(int argc, char *argv[]) {
int z{1000000}; auto n{uT{z}}; cout<<"untouchables below "<<z<<"->"<<n.count()<<endl;
}
</syntaxhighlight>
{{out}}
<pre>
untouchables below 1000000->150232
real 3m6.909s
</pre>
;Count less than or equal 2000000
<syntaxhighlight lang="cpp">
int main(int argc, char *argv[]) {
int z{2000000}; auto n{uT{z}}; cout<<"untouchables below "<<z<<"->"<<n.count()<<endl;
}
</syntaxhighlight>
{{out}}
<pre>
untouchables below 2000000->305290
real 11m28.031s
</pre>
 
=={{header|Delphi}}==
 
{{output?}}
 
{{libheader| System.SysUtils}}
{{Trans|Go}}
<syntaxhighlight 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.</syntaxhighlight>
 
=={{header|FreeBASIC}}==
<syntaxhighlight lang="vbnet">Const max_untouchable = 1e6
 
Dim Shared untouchable(1 To max_untouchable) As Uinteger
For i As Uinteger = 1 To Ubound(untouchable)
untouchable(i) = True
Next i
 
Sub show_untouchable_statistics()
Dim As Uinteger i, cnt = 0
For i = 1 To Ubound(untouchable)
If untouchable(i) Then cnt += 1
If i = 10 Orelse i = 100 Orelse i = 1000 Orelse i = 10000 Orelse i = 1e5 Then
Print Using "###,### untouchable numbers were found <= #,###,###"; cnt; i
End If
Next i
End Sub
 
Sub print_untouchables(n As Uinteger)
Print Using "List of untouchable numbers <= #,###:"; n
Dim As Uinteger i, cnt = 0
For i = 1 To n
If untouchable(i) Then
Print Using "##,###"; i;
cnt += 1
Print Iif(cnt Mod 16 = 0, !"\n", " ");
End If
Next i
Print: Print
Print Using "###,### untouchable numbers were found <= #,###,###"; cnt; n
End Sub
 
Dim As Uinteger i, j
untouchable(1) = False
untouchable(3) = False
For i = 7 To Ubound(untouchable) Step 2
untouchable(i) = False
Next i
 
Dim Shared spd(1 To max_untouchable * 64) As Uinteger
Dim As Uinteger ub = Ubound(spd)
For i = 1 To ub
spd(i) = 1
Next i
For i = 2 To ub
For j = i + i To ub Step i
spd(j) += i
Next j
Next i
For i = 1 To ub
If spd(i) <= Ubound(untouchable) Then untouchable(spd(i)) = False
Next i
 
' Show the untouchable numbers up to 2000
print_untouchables(2000)
' Show the counts of untouchable numbers
show_untouchable_statistics()
 
Sleep</syntaxhighlight>
 
=={{header|F_Sharp|F#}}==
===The Function===
This task uses [[Extensible_prime_generator#The_functions|Extensible Prime Generator (F#)]].<br>
It implements [[Talk:Untouchable_numbers#Nice_recursive_solution]]
<syntaxhighlight 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]
</syntaxhighlight>
===The Task===
;Less than 2000
<syntaxhighlight 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 "")
</syntaxhighlight>
{{out}}
<pre>
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
</pre>
;Count less than or equal 100000
<syntaxhighlight lang="fsharp">
printfn "%d" (uT 100000|>Array.filter id|>Array.length)
</syntaxhighlight>
{{out}}
<pre>
13863
Real: 00:00:02.784, CPU: 00:00:02.750
</pre>
;Count less than or equal 1000000
<syntaxhighlight lang="fsharp">
printfn "%d" (uT 1000000|>Array.filter id|>Array.length)
</syntaxhighlight>
{{out}}
<pre>
150232
Real: 00:03:08.929, CPU: 00:03:08.859
</pre>
;Count less than or equal 2000000
<syntaxhighlight lang="fsharp">
printfn "%d" (uT 2000000|>Array.filter id|>Array.length)
</syntaxhighlight>
{{out}}
<pre>
305290
Real: 00:11:22.325, CPU: 00:11:21.828
</pre>
;Count less than or equal 3000000
<syntaxhighlight lang="fsharp">
printfn "%d" (uT 3000000|>Array.filter id|>Array.length)
</syntaxhighlight>
{{out}}
<pre>
462110
Real: 00:24:00.126, CPU: 00:23:59.203
</pre>
 
=={{header|Go}}==
===Version1 ===
{{trans|Wren}}
This was originally based on the Wren Version 1 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.
... but now goes up to 1,000,000 as I'm reasonably sure (see talk page) that the figure obtained is correct.
<langsyntaxhighlight lang="go">package main
 
import "fmt"
 
func sumDivisors(n int) int {
sum := 01
i := 1
k := 2
if n%2 == 0 {
k = 1
}
for i := 1 + k; i*i <= n; i += k {
if n%i == 0 {
sum += i
Line 84 ⟶ 649:
}
}
i += k
}
return sum - n
}
 
func sieve(n int) []bool {
n++
s := make([]bool, n*4+1) // all false by default
for i := 06; i <= n; i++ {
s[sd := sumDivisors(i)] = true
if sd <= n {
s[sd] = true
}
}
return s
}
 
func isPrimeprimeSieve(nlimit int) []bool {
if n < 2 {limit++
// True denotes returncomposite, false denotes prime.
c := make([]bool, limit) // all false by default
}
if n%2c[0] == 0 {true
c[1] return n == 2true
// no need to bother with even numbers over 2 for this task
}
ifp := n%3 ==// 0Start {from 3.
for return n == 3{
} p2 := p * p
d : if p2 >= 5limit {
for d*d <= n { break
if n%d == 0 {
return false
}
dfor i := p2; i < limit; i += 2 * p {
if n%d == 0 {c[i] = true
return false}
for {
p += 2
if !c[p] {
break
}
}
d += 4
}
return truec
}
 
Line 137 ⟶ 706:
}
 
func main() {
limit := 1000000
sc := sieveprimeSieve(20 * limit)
s := sieve(63 * limit)
untouchable := []int{2, 5}
for n := 6; n <= limit; n += 2 {
if !s[n] && !isPrime(c[n-1)] && !isPrime(c[n-3)] {
untouchable = append(untouchable, n)
}
Line 173 ⟶ 743:
cl := commatize(limit)
fmt.Printf("%7s untouchable numbers were found <= %s\n", cu, cl)
}</langsyntaxhighlight>
 
{{out}}
Line 205 ⟶ 775:
1,212 untouchable numbers were found <= 10,000
13,863 untouchable numbers were found <= 100,000
150,252232 untouchable numbers were found <= 1,000,000
</pre>
 
===Version 2===
{{trans|C}}
{{libheader|Go-rcu}}
This version uses a much more efficient algorithm for calculating the sums of divisors in bulk.
 
Run time for finding untouchable numbers up to 1 million is now only 11 seconds which (surprisingly) is faster than C itself on the same machine.
 
Also, since the first version was written, some of the functions used have now been incorporated into the above library.
<syntaxhighlight lang="go">package main
 
import (
"fmt"
"rcu"
)
 
func main() {
limit := 1000000
m := 63
c := rcu.PrimeSieve(limit, false)
n := m*limit + 1
sumDivs := make([]int, n)
for i := 1; i < n; i++ {
for j := i; j < n; j += i {
sumDivs[j] += i
}
}
s := make([]bool, n) // all false
for i := 1; i < n; i++ {
sum := sumDivs[i] - i // proper divs sum
if sum <= n {
s[sum] = true
}
}
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", rcu.Commatize(untouchable[i]))
if (i+1)%10 == 0 {
fmt.Println()
}
count++
}
fmt.Printf("\n\n%7s untouchable numbers were found <= 2,000\n", rcu.Commatize(count))
p := 10
count = 0
for _, n := range untouchable {
count++
if n > p {
cc := rcu.Commatize(count - 1)
cp := rcu.Commatize(p)
fmt.Printf("%7s untouchable numbers were found <= %9s\n", cc, cp)
p = p * 10
if p == limit {
break
}
}
}
cu := rcu.Commatize(len(untouchable))
cl := rcu.Commatize(limit)
fmt.Printf("%7s untouchable numbers were found <= %s\n", cu, cl)
}</syntaxhighlight>
 
{{out}}
<pre>
Same as Version 1.
</pre>
 
=={{header|J}}==
<syntaxhighlight 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
</syntaxhighlight>
We may revisit to correct the larger untouchable tallies. The straightforward algorithm presented is a little bit efficient, and, I claim, the verb <tt>(candidates-.spds)</tt> 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. [https://www.eecg.utoronto.ca/~jzhu/csc326/readings/iverson.pdf]
 
Time (seconds) and space (bytes) on a high end six year old (new in 2015) laptop computer.
<pre>
timespacex'(-.candidates@#)/:~~.spds 10000'
600.285 4.29497e9
</pre>
<pre>
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 │
└──────┴──────────────────────────────────────────┘
</pre>
 
Line 212 ⟶ 922:
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.
<langsyntaxhighlight lang="julia">using Primes
 
function properfactorsum(n)
Line 245 ⟶ 955:
println("The count of untouchable numbers ≤ $N is: ", count(x -> untouchables[x], 1:N))
end
</langsyntaxhighlight>{{out}}
<pre>
The untouchable numbers ≤ 2000 are:
Line 277 ⟶ 987:
</pre>
 
=={{header|PhixMathematica}}/{{header|Wolfram Language}}==
<syntaxhighlight lang="mathematica">f = DivisorSigma[1, #] - # &;
Note: the limit of 18*n is unsound, albeit matching the above b005114.txt list, see talk page (1464)
limit = 10^5;
<lang Phix>procedure untouchable(integer n, cols=0, tens=0)
c = Not /@ PrimeQ[Range[limit]];
bool tell = n>0
slimit = 15 limit;
n = abs(n)
s = ConstantArray[False, slimit + 1];
sequence sums = repeat(0,n+3)
untouchable = {2, 5};
for i=1 to n do
Do[
integer p = get_prime(i)
val = f[i];
if p>n then exit end if
If[val <= slimit,
sums[p+1] = 1
sumss[p+3[val]] = 1True
]
end for
,
sums[5] = 0
{i, 6, slimit}
]
Do[
If[! s[[n]],
If[c[[n - 1]],
If[c[[n - 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])]</syntaxhighlight>
{{out}}
<pre>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</pre>
 
=={{header|Nim}}==
-- for j=2 to 2*n do
I borrowed some ideas from Go version 1, but keep the limit to 100_000 as in Wren version 1.
for j=2 to 18*n do -- see talk page (1464)
<syntaxhighlight lang="nim">import math, strutils
integer y = sum(factors(j,-1))
 
if y<=n then
const
sums[y] = 1
Lim1 = 100_000 # Limit for untouchable numbers.
end if
Lim2 = 14 * Lim1 # Limit for computation of sum of divisors.
end for
 
if tell then
proc sumdiv(n: uint): uint =
printf(1,"The list of all untouchable numbers <= %d:\n",{n})
## Return the sum of the strict divisors of "n".
end if
result = 1
string line = " 2 5"
let r = sqrt(n.float).uint
integer cnt = 2
let k for t=6 toif (n byand 1) == 0: 1u 2else: do2u
for d in countup(k + 1, ifr, sums[t]=0 thenk):
if n mod d cnt +== 10:
result += if tell thend
let q = n div line &= sprintf("%,8d",t)d
if q != d: result if remainder(cnt,cols)+=0 thenq
 
printf(1,"%s\n",line)
var
line = ""
isSumDiv: array[1..Lim2, bool]
end if
isPrime: array[1..Lim1, bool]
end if
 
end if
# Fill both sieves in a single pass.
end for
for n in 1u..Lim2:
if tell then
let s = sumdiv(n)
if line!="" then
if s <= Lim2:
printf(1,"%s\n",line)
isSumDiv[s] = end iftrue
if s == printf(1,"\ and n") <= Lim1:
isPrime[n] = true
end if
isPrime[1] = false
printf(1,"%,20d untouchable numbers were found <= %,d\n",{cnt,n})
 
for p=1 to tens do
# Build list of untouchable numbers.
untouchable(-power(10,p))
var list = @[2, 5]
end for
for n in countup(6, Lim1, 2):
end procedure
if not (isSumDiv[n] or isPrime[n - 1] or isPrime[n - 3]):
list.add n
untouchable(2000, 10, 5)</lang>
 
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)</syntaxhighlight>
 
{{out}}
<pre>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.</pre>
 
=={{header|Pascal}}==
modified [[Factors_of_an_integer#using_Prime_decomposition]] to calculate only sum of divisors<BR>
Appended a list of count of untouchable numbers out of math.dartmouth.edu/~carlp/uupaper3.pdf<BR>
Nigel is still right, that this procedure is not the way to get proven results.
<syntaxhighlight lang="pascal">program UntouchableNumbers;
program UntouchableNumbers;
{$IFDEF FPC}
{$MODE DELPHI} {$OPTIMIZATION ON,ALL} {$COPERATORS ON}
{$CODEALIGN proc=16,loop=4}
{$ELSE}
{$APPTYPE CONSOLE}
{$ENDIF}
uses
sysutils,strutils
{$IFDEF WINDOWS},Windows{$ENDIF}
;
const
MAXPRIME = 1742537;
//sqr(MaxPrime) = 3e12
LIMIT = 5*1000*1000;
LIMIT_mul = trunc(exp(ln(LIMIT)/3))+1;
 
const
SizePrDeFe = 16*8192;//*size of(tprimeFac) =16 byte 2 Mb ~ level 3 cache
type
tdigits = array [0..31] of Uint32;
tprimeFac = packed record
pfSumOfDivs,
pfRemain : Uint64;
end;
tpPrimeFac = ^tprimeFac;
 
tPrimeDecompField = array[0..SizePrDeFe-1] of tprimeFac;
 
tPrimes = array[0..1 shl 17-1] of Uint32;
 
var
{$ALIGN 16}
PrimeDecompField :tPrimeDecompField;
{$ALIGN 16}
SmallPrimes: tPrimes;
pdfIDX,pdfOfs: NativeInt;
TD : Int64;
 
procedure OutCounts(pUntouch:pByte);
var
n,cnt,lim,deltaLim : NativeInt;
Begin
n := 0;
cnt := 0;
deltaLim := 100;
lim := deltaLim;
repeat
cnt += 1-pUntouch[n];
if n = lim then
Begin
writeln(Numb2USA(IntToStr(lim)):13,' ',Numb2USA(IntToStr(cnt)):12);
lim += deltaLim;
if lim = 10*deltaLim then
begin
deltaLim *=10;
lim := deltaLim;
writeln;
end;
end;
 
inc(n);
until n > LIMIT;
end;
 
function OutN(n:UInt64):UInt64;
begin
write(Numb2USA(IntToStr(n)):15,' dt ',(GettickCount64-TD)/1000:5:3,' s'#13);
TD := GettickCount64;
result := n+LIMIT;
end;
 
//######################################################################
//gets sum of divisors of consecutive integers fast
procedure InitSmallPrimes;
//get primes. Sieving only odd numbers
var
pr : array[0..MAXPRIME] 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>MAXPRIME then
BREAK;
d := 2*p+1;
repeat
pr[j] := 1;
j += d;
until j>MAXPRIME;
until false;
 
SmallPrimes[1] := 3;
SmallPrimes[2] := 5;
j := 3;
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] := 2*p+1;
inc(j);
end;
p+=flipflop;
flipflop := 3-flipflop;
until (p > MAXPRIME) OR (j>High(SmallPrimes));
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 IncByOneInBase(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;
 
procedure CalcSumOfDivs(var pdf:tPrimeDecompField;var dgt:tDigits;n,k,pr:Uint64);
var
fac,s :Uint64;
j : Int32;
Begin
//j is power of prime
j := CnvtoBASE(dgt,n+k,pr);
repeat
fac := 1;
s := 1;
repeat
fac *= pr;
dec(j);
s += fac;
until j<= 0;
with pdf[k] do
Begin
pfSumOfDivs *= s;
pfRemain := pfRemain DIV fac;
end;
j := IncByOneInBase(dgt,pr);
k += pr;
until k >= SizePrDeFe;
end;
 
function SieveOneSieve(var pdf:tPrimeDecompField):boolean;
var
dgt:tDigits;
i,j,k,pr,n,MaxP : Uint64;
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
pfSumOfDivs := 1;
pfRemain := n+i;
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);
pfRemain := (n+i) shr j;
pfSumOfDivs := (Uint64(1) shl (j+1))-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 > maxP then
BREAK;
 
CalcSumOfDivs(pdf,dgt,n,k,pr);
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
pfSumOFDivs *= (j+1);
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
Begin
writeln('of limits ');
EXIT(NIL);
end;
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;
//gets sum of divisors of consecutive integers fast
//######################################################################
 
procedure CheckRest(n: Uint64;pUntouch:pByte);
var
k,lim : Uint64;
begin
lim := 2*LIMIT;
repeat
k := GetNextPrimeDecomp^.pfSumOfDivs-n;
inc(n);
if Not(ODD(k)) AND (k<=LIMIT) then
pUntouch[k] := 1;
// showing still alive not for TIO.RUN
// if n >= lim then lim := OutN(n);
until n >LIMIT_mul*LIMIT;
end;
 
var
Untouch : array of byte;
pUntouch: pByte;
puQW : pQword;
T0:Int64;
n,k : NativeInt;
Begin
if sqrt(LIMIT_mul*LIMIT) >=MAXPRIME then
Begin
writeln('Need to extend count of primes > ',
trunc(sqrt(LIMIT_mul*LIMIT))+1);
HALT(0);
end;
 
setlength(untouch,LIMIT+8+1);
pUntouch := @untouch[0];
//Mark all odd as touchable
puQW := @pUntouch[0];
For n := 0 to LIMIT DIV 8 do puQW[n] := $0100010001000100;
 
InitSmallPrimes;
T0 := GetTickCount64;
writeln('LIMIT = ',Numb2USA(IntToStr(LIMIT)));
writeln('factor beyond LIMIT ',LIMIT_mul);
 
n := 0;
Init_Sieve(n);
 
pUntouch[1] := 1;//all primes
repeat
k := GetNextPrimeDecomp^.pfSumOfDivs-n;
inc(n);//n-> n+1
if k <= LIMIT then
begin
If k <> 1 then
pUntouch[k] := 1
else
begin
//n-1 is prime p
//mark p*p
pUntouch[n] := 1;
//mark 2*p
//5 marked by prime 2 but that is p*p, but 4 has factor sum = 3
pUntouch[n+2] := 1;
end;
end;
until n > LIMIT-2;
//unmark 5 and mark 0
puntouch[5] := 0;
pUntouch[0] := 1;
 
//n=limit-1
k := GetNextPrimeDecomp^.pfSumOfDivs-n;
inc(n);
If (k <> 1) AND (k<=LIMIT) then
pUntouch[k] := 1
else
pUntouch[n] := 1;
//n=limit
k := GetNextPrimeDecomp^.pfSumOfDivs-n;
If Not(odd(k)) AND (k<=LIMIT) then
pUntouch[k] := 1;
 
 
n:= limit+1;
writeln('runtime for n<= LIMIT ',(GetTickCount64-T0)/1000:0:3,' s');
writeln('Check the rest ',Numb2USA(IntToStr((LIMIT_mul-1)*Limit)));
TD := GettickCount64;
CheckRest(n,pUntouch);
writeln('runtime ',(GetTickCount64-T0)/1000:0:3,' s');
T0 := GetTickCount64-T0;
 
OutCounts(pUntouch);
end.</syntaxhighlight>
{{out}}
<pre>
TIO.RUN
LIMIT = 5,000,000
factor beyond LIMIT 171
runtime for n smaller than LIMIT 0.204 s
Check the rest 850,000,000
runtime 32.643 s
100 5
200 10
300 22
400 30
500 38
600 48
700 55
800 69
900 80
 
1,000 89
2,000 196
3,000 318
4,000 443
5,000 570
6,000 689
7,000 801
8,000 936
9,000 1,082
 
10,000 1,212
20,000 2,566
30,000 3,923
40,000 5,324
50,000 6,705
60,000 8,153
70,000 9,586
80,000 10,994
90,000 12,429
 
100,000 13,863
200,000 28,572
300,000 43,515
400,000 58,459
500,000 73,565
600,000 88,828
700,000 104,062
800,000 119,302
900,000 134,758
 
1,000,000 150,232
2,000,000 305,290
3,000,000 462,110
4,000,000 619,638
5,000,000 777,672
 
Real time: 32.827 s CPU share: 99.30 %
//url=https://math.dartmouth.edu/~carlp/uupaper3.pdf
100000 13863
200000 28572
300000 43515
400000 58459
500000 73565
600000 88828
700000 104062
800000 119302
900000 134758
1000000 150232
2000000 305290
3000000 462110
4000000 619638
5000000 777672
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
 
... at home up to 1E8
50,000,000 8,060,163
60,000,000 9,694,467
70,000,000 11,330,312
80,000,000 12,967,239
90,000,000 14,606,549
 
100,000,000 16,246,940
 
real 18m51,234s
</pre>
 
=={{header|Perl}}==
{{libheader|ntheory}}
<syntaxhighlight 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)
}
}</syntaxhighlight>
{{out}}
<pre>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</pre>
 
=={{header|Phix}}==
<!--<syntaxhighlight lang="phix">(phixonline)-->
<span style="color: #008080;">constant</span> <span style="color: #000000;">limz</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #000000;">8</span><span style="color: #0000FF;">,</span><span style="color: #000000;">9</span><span style="color: #0000FF;">,</span><span style="color: #000000;">18</span><span style="color: #0000FF;">,</span><span style="color: #000000;">64</span><span style="color: #0000FF;">}</span> <span style="color: #000080;font-style:italic;">-- found by experiment</span>
<span style="color: #008080;">procedure</span> <span style="color: #000000;">untouchable</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">cols</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">tens</span><span style="color: #0000FF;">=</span><span style="color: #000000;">0</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">atom</span> <span style="color: #000000;">t0</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">time</span><span style="color: #0000FF;">(),</span> <span style="color: #000000;">t1</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">t0</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span>
<span style="color: #004080;">bool</span> <span style="color: #000000;">tell</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">n</span><span style="color: #0000FF;">></span><span style="color: #000000;">0</span>
<span style="color: #000000;">n</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">abs</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">)</span>
<span style="color: #004080;">sequence</span> <span style="color: #000000;">sums</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">repeat</span><span style="color: #0000FF;">(</span><span style="color: #000000;">0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">n</span><span style="color: #0000FF;">+</span><span style="color: #000000;">3</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">i</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">n</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">p</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">get_prime</span><span style="color: #0000FF;">(</span><span style="color: #000000;">i</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">></span><span style="color: #000000;">n</span> <span style="color: #008080;">then</span> <span style="color: #008080;">exit</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #000000;">sums</span><span style="color: #0000FF;">[</span><span style="color: #000000;">p</span><span style="color: #0000FF;">+</span><span style="color: #000000;">1</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span>
<span style="color: #000000;">sums</span><span style="color: #0000FF;">[</span><span style="color: #000000;">p</span><span style="color: #0000FF;">+</span><span style="color: #000000;">3</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #000000;">sums</span><span style="color: #0000FF;">[</span><span style="color: #000000;">5</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">0</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">m</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">floor</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">log10</span><span style="color: #0000FF;">(</span><span style="color: #000000;">n</span><span style="color: #0000FF;">))</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">lim</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">limz</span><span style="color: #0000FF;">[</span><span style="color: #000000;">m</span><span style="color: #0000FF;">]*</span><span style="color: #000000;">n</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">j</span><span style="color: #0000FF;">=</span><span style="color: #000000;">2</span> <span style="color: #008080;">to</span> <span style="color: #000000;">lim</span> <span style="color: #008080;">do</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">y</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">sum</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">factors</span><span style="color: #0000FF;">(</span><span style="color: #000000;">j</span><span style="color: #0000FF;">,-</span><span style="color: #000000;">1</span><span style="color: #0000FF;">))</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">y</span><span style="color: #0000FF;"><=</span><span style="color: #000000;">n</span> <span style="color: #008080;">then</span>
<span style="color: #000000;">sums</span><span style="color: #0000FF;">[</span><span style="color: #000000;">y</span><span style="color: #0000FF;">]</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">1</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">if</span> <span style="color: #7060A8;">platform</span><span style="color: #0000FF;">()!=</span><span style="color: #004600;">JS</span> <span style="color: #008080;">and</span> <span style="color: #7060A8;">time</span><span style="color: #0000FF;">()></span><span style="color: #000000;">t1</span> <span style="color: #008080;">then</span>
<span style="color: #7060A8;">progress</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"j:%,d/%,d (%3.2f%%)\r"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">j</span><span style="color: #0000FF;">,</span><span style="color: #000000;">lim</span><span style="color: #0000FF;">,(</span><span style="color: #000000;">j</span><span style="color: #0000FF;">/</span><span style="color: #000000;">lim</span><span style="color: #0000FF;">)*</span><span style="color: #000000;">100</span><span style="color: #0000FF;">})</span>
<span style="color: #000000;">t1</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">time</span><span style="color: #0000FF;">()+</span><span style="color: #000000;">1</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">if</span> <span style="color: #7060A8;">platform</span><span style="color: #0000FF;">()!=</span><span style="color: #004600;">JS</span> <span style="color: #008080;">then</span> <span style="color: #7060A8;">progress</span><span style="color: #0000FF;">(</span><span style="color: #008000;">""</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">tell</span> <span style="color: #008080;">then</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"The list of all untouchable numbers &lt;= %d:\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">n</span><span style="color: #0000FF;">})</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #004080;">string</span> <span style="color: #000000;">line</span> <span style="color: #0000FF;">=</span> <span style="color: #008000;">" 2 5"</span>
<span style="color: #004080;">integer</span> <span style="color: #000000;">cnt</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">2</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">t</span><span style="color: #0000FF;">=</span><span style="color: #000000;">6</span> <span style="color: #008080;">to</span> <span style="color: #000000;">n</span> <span style="color: #008080;">by</span> <span style="color: #000000;">2</span> <span style="color: #008080;">do</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">sums</span><span style="color: #0000FF;">[</span><span style="color: #000000;">t</span><span style="color: #0000FF;">]=</span><span style="color: #000000;">0</span> <span style="color: #008080;">then</span>
<span style="color: #000000;">cnt</span> <span style="color: #0000FF;">+=</span> <span style="color: #000000;">1</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">tell</span> <span style="color: #008080;">then</span>
<span style="color: #000000;">line</span> <span style="color: #0000FF;">&=</span> <span style="color: #7060A8;">sprintf</span><span style="color: #0000FF;">(</span><span style="color: #008000;">"%,8d"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">t</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">if</span> <span style="color: #7060A8;">remainder</span><span style="color: #0000FF;">(</span><span style="color: #000000;">cnt</span><span style="color: #0000FF;">,</span><span style="color: #000000;">cols</span><span style="color: #0000FF;">)=</span><span style="color: #000000;">0</span> <span style="color: #008080;">then</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%s\n"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">line</span><span style="color: #0000FF;">)</span>
<span style="color: #000000;">line</span> <span style="color: #0000FF;">=</span> <span style="color: #008000;">""</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">tell</span> <span style="color: #008080;">then</span>
<span style="color: #008080;">if</span> <span style="color: #000000;">line</span><span style="color: #0000FF;">!=</span><span style="color: #008000;">""</span> <span style="color: #008080;">then</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%s\n"</span><span style="color: #0000FF;">,</span><span style="color: #000000;">line</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"\n"</span><span style="color: #0000FF;">)</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
<span style="color: #004080;">string</span> <span style="color: #000000;">t</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">elapsed</span><span style="color: #0000FF;">(</span><span style="color: #7060A8;">time</span><span style="color: #0000FF;">()-</span><span style="color: #000000;">t0</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">" (%s)"</span><span style="color: #0000FF;">)</span>
<span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"%,20d untouchable numbers were found &lt;= %,d%s\n"</span><span style="color: #0000FF;">,{</span><span style="color: #000000;">cnt</span><span style="color: #0000FF;">,</span><span style="color: #000000;">n</span><span style="color: #0000FF;">,</span><span style="color: #000000;">t</span><span style="color: #0000FF;">})</span>
<span style="color: #008080;">for</span> <span style="color: #000000;">p</span><span style="color: #0000FF;">=</span><span style="color: #000000;">1</span> <span style="color: #008080;">to</span> <span style="color: #000000;">tens</span> <span style="color: #008080;">do</span>
<span style="color: #000000;">untouchable</span><span style="color: #0000FF;">(-</span><span style="color: #7060A8;">power</span><span style="color: #0000FF;">(</span><span style="color: #000000;">10</span><span style="color: #0000FF;">,</span><span style="color: #000000;">p</span><span style="color: #0000FF;">))</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">for</span>
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
<span style="color: #000000;">untouchable</span><span style="color: #0000FF;">(</span><span style="color: #000000;">2000</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">10</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">6</span><span style="color: #0000FF;">-(</span><span style="color: #7060A8;">platform</span><span style="color: #0000FF;">()==</span><span style="color: #004600;">JS</span><span style="color: #0000FF;">))</span>
<!--</syntaxhighlight>-->
{{out}}
<pre>
Line 354 ⟶ 1,757:
196 untouchable numbers were found <= 2,000
2 untouchable numbers were found <= 10
56 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)
</pre>
for comparison, on the same box, the Julia entry took 58 minutes and 40s.
 
=={{header|Raku}}==
Borrow the proper divisors routine from [https://rosettacode.org/wiki/Proper_divisors#Raku here].
{{trans|Wren}}
<syntaxhighlight lang="raku" line># 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</syntaxhighlight>
{{out}}
<pre>
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
</pre>
 
=={{header|REXX}}==
Some optimization was done to the code to produce prime numbers, &nbsp; since a simple test could be made to exclude
<br>the calculation of touchability for primes as the aliquot sum of a prime is always unity.
<br>This saved around &nbsp; '''15%''' &nbsp; of the running time.
 
A fair amount of code was put into the generation of primes, &nbsp; but the computation of the aliquot sum wasis the area
<br>that consumedconsumes the most CPU time.
 
<lang rexx>/*REXX pgm finds N untouchable numbers (numbers that can't be equal to any aliquot sum).*/
The REXX code below will accurately calculate &nbsp; ''untouchable numbers'' &nbsp; up to and including &nbsp; 100,000. &nbsp; Beyond that,
<br>a higher overreach &nbsp; (the '''over''' option) &nbsp; 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.
<syntaxhighlight 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.*/
Line 427 ⟶ 1,916:
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. */</langsyntaxhighlight>
{{out|output|text=&nbsp; when using the default inputs:}}
<pre>
Line 461 ⟶ 1,950:
1,212 untouchable numbers were found ≤ 10,000
13,863 untouchable numbers were found ≤ 100,000
150,252 untouchable numbers were found ≤ 1,000,000
</pre>
 
=={{header|Wren}}==
{{libheader|Wren-seq}}
{{libheader|Wren-math}}
{{libheader|Wren-fmt}}
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.
===Version 1===
<lang ecmascript>import "/math" for Int, Nums
Run time about 70 seconds on my Core i7 machine.
import "/seq" for Lst
<syntaxhighlight lang="wren">import "./fmtmath" for FmtInt, Nums
import "./fmt" for Fmt
 
var sieve = Fn.new { |n|
n = n + 1
var s = List.filled(n*4+1, false)
for (i in 0..n) {
var sum = Nums.sum(Int.properDivisors(i))
if (sum <= n) s[sum] = true
}
return s
Line 484 ⟶ 1,972:
 
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] && !Int.isPrime(c[n-1)] && !Int.isPrime(c[n-3)]) untouchable.add(n)
n = n + 2
}
 
System.print("List of untouchable numbers <= 2,000:")
for Fmt.tprint(chunk"$,6d", 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 })
Line 508 ⟶ 1,995:
}
}
Fmt.print("$,6d untouchable numbers were found <= $,d", untouchable.count, limit)</langsyntaxhighlight>
 
{{out}}
Line 540 ⟶ 2,027:
1,212 untouchable numbers were found <= 10,000
13,863 untouchable numbers were found <= 100,000
</pre>
 
===Version 2===
{{trans|C}}
This version uses a much more efficient algorithm for calculating the sums of divisors in bulk.
 
Run time for untouchable numbers up to 100,000 (m = 14) is now only 1.4 seconds and 1,000,000 (m = 63) is reached in 132 seconds.
<syntaxhighlight lang="wren">import "./math" for Int, Nums
import "./fmt" for Fmt
 
var limit = 1e6
var m = 63
var c = Int.primeSieve(limit, false)
var n = m * limit + 1
var sumDivs = List.filled(n, 0)
for (i in 1...n) {
var j = i
while (j < n) {
sumDivs[j] = sumDivs[j] + i
j = j + i
}
}
var s = List.filled(n, false)
for (i in 1...n) {
var sum = sumDivs[i] - i // proper divs sum
if (sum <= n) s[sum] = true
}
var untouchable = [2, 5]
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:")
Fmt.tprint("$,6d", untouchable.where { |n| n <= 2000 }, 10)
System.print()
Fmt.print("$,7d 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("$,7d untouchable numbers were found <= $,9d", count-1, p)
p = p * 10
if (p == limit) break
}
}
Fmt.print("$,7d untouchable numbers were found <= $,d", untouchable.count, limit)</syntaxhighlight>
 
{{out}}
As Version 1, plus:
<pre>
150,232 untouchable numbers were found <= 1,000,000
</pre>
Jjuanhdez
2,122

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/307811...362949"