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.


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.