Anonymous user
Untouchable numbers: Difference between revisions
added pascal fast but with problems after 40,000,000 count -> 6,430,223 < wrong by -1
(→{{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:
There are 1212 untouchable numbers ≤ 10000.
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}}==
{{libheader|ntheory}}
| |||