Untouchable numbers: Difference between revisions
Content added Content deleted
(→{{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}} |