Erdős-Nicolas numbers

From Rosetta Code
Task
Erdős-Nicolas numbers
You are encouraged to solve this task according to the task description, using any language you may know.
Definition

An Erdős–Nicolas number is a positive integer which is not perfect but is equal to the sum of its first k divisors (arranged in ascending order and including one) for some value of k greater than one.

Examples

24 is an Erdős–Nicolas number because the sum of its first 6 divisors (1, 2, 3, 4, 6 and 8) is equal to 24 and it is not perfect because 12 is also a divisor.

6 is not an Erdős–Nicolas number because it is perfect (1 + 2 + 3 = 6).

48 is not an Erdős–Nicolas number because its divisors are: 1, 2, 3, 4, 6, 8, 12, 16, 24 and 48. The first seven of these add up to 36, but the first eight add up to 52 which is more than 48.

Task

Find and show here the first 8 Erdős–Nicolas numbers and the number of divisors needed (i.e. the value of 'k') to satisfy the definition.

Stretch

Do the same for any further Erdős–Nicolas numbers which you have the patience for.

Note

As all known Erdős–Nicolas numbers are even you may assume this to be generally true in order to quicken up the search. However, it is not obvious (to me at least) why this should necessarily be the case.

Reference



ALGOL 68

Builds tables of proper divisor counts and sums and finds the numbers whilst doing it. This means that the numbers are not found in numerical order.

BEGIN # find some Erdos-Nicolas numbers: numbers equal to the sum of their    #
      # first k proper divisors but k is not the count of all their proper    #
      # divisors ( so the numbers aren't perfect )                            #
    INT max number  = 2 000 000;      # largest number we will consider       #
    # construct tables of the divisor counts and divisor sums and check for   #
    # the numbers as we do it - note they will not necessarily be found in    #
    #                                order                                    #
    [ 1 : max number ]INT dsum;   FOR i TO UPB dsum   DO dsum[   i ] := 1 OD;
    [ 1 : max number ]INT dcount; FOR i TO UPB dcount DO dcount[ i ] := 1 OD;
    FOR i FROM 2 TO UPB dsum
        DO FOR j FROM i + i BY i TO UPB dsum DO
            # have another proper divisor                                     #
            IF dsum[ j ] = j THEN
                # the divisor sum is currently equal to the number but is     #
                # about to increase, so we have an Erdos-Nicolas number       #
                print( ( whole( j, -10 ), " equals the sum of its first "
                       , whole( dcount[ j ], 0 ), " divisors"
                       , newline
                       )
                     )
            FI;
            dsum[   j ] +:= i;
            dcount[ j ] +:= 1
        OD
    OD
END
Output:
        24 equals the sum of its first 6 divisors
      2016 equals the sum of its first 31 divisors
      8190 equals the sum of its first 43 divisors
     42336 equals the sum of its first 66 divisors
     45864 equals the sum of its first 66 divisors
    714240 equals the sum of its first 113 divisors
    392448 equals the sum of its first 68 divisors
   1571328 equals the sum of its first 115 divisors

With max number set to 200 000 000, another two numbers can be found. Note that that would probably be too large for ALGOL 68G under Windows, so another compiler would need to be used:

        24 equals the sum of its first 6 divisors
      2016 equals the sum of its first 31 divisors
      8190 equals the sum of its first 43 divisors
     42336 equals the sum of its first 66 divisors
     45864 equals the sum of its first 66 divisors
    714240 equals the sum of its first 113 divisors
    392448 equals the sum of its first 68 divisors
   1571328 equals the sum of its first 115 divisors
  61900800 equals the sum of its first 280 divisors
  91963648 equals the sum of its first 142 divisors


Arturo

erdosNicolas: function [n][
    facts: factors n
    if facts > 2 [
        loop 1..(size facts)-2 'k [
            if n = sum first.n:k facts -> return @[n, k]
        ]
    ]

    return ø
]

cnt: 0
i: 2
while [cnt < 8][
    if enNum: <= erdosNicolas i [
        print[enNum\0 "equals the sum of its first" enNum\1 "divisors"]
        cnt: cnt + 1
    ]
    i: i + 2
]
Output:
24 equals the sum of its first 6 divisors 
2016 equals the sum of its first 31 divisors 
8190 equals the sum of its first 43 divisors 
42336 equals the sum of its first 66 divisors 
45864 equals the sum of its first 66 divisors 
392448 equals the sum of its first 68 divisors 
714240 equals the sum of its first 113 divisors 
1571328 equals the sum of its first 115 divisors

BASIC

Rosetta Code problem: https://rosettacode.org/wiki/Erdős-Nicolas_numbers

by Jjuanhdez, 02/2023

BASIC256

Translation of: FreeBASIC
limite = 400000
dim DSum(limite+1) fill 1
dim DCount(limite+1) fill 1

for i = 2 to limite
	j = i + i
	while j <= limite
		if DSum[j] = j then
			print rjust(j,8); " equals the sum of its first "; rjust(DCount[j],3); " divisors"
		end if
		DSum[j] += i
		DCount[j] += 1
		j += i
	end while
next i
Output:
24 equals the sum of its first 6 divisors
2016 equals the sum of its first 31 divisors
8190 equals the sum of its first 43 divisors
42336 equals the sum of its first 66 divisors
45864 equals the sum of its first 66 divisors

FreeBASIC

Dim As Uinteger limite = 2e6
Dim As Uinteger DSum(limite+1), DCount(limite+1)
Dim As Integer i, j

For i = 0 To limite
    DSum(i) = 1
    DCount(i) = 1
Next i

For i = 2 To limite
    j = i + i
    While j <= limite
        If DSum(j) = j Then
            Print Using "######## equals the sum of its first ### divisors"; j; DCount(j)
        End If
        DSum(j) += i
        DCount(j) += 1
        j += i
    Wend
Next i
Output:
      24 equals the sum of its first   6 divisors
    2016 equals the sum of its first  31 divisors
    8190 equals the sum of its first  43 divisors
   42336 equals the sum of its first  66 divisors
   45864 equals the sum of its first  66 divisors
  714240 equals the sum of its first 113 divisors
  392448 equals the sum of its first  68 divisors
 1571328 equals the sum of its first 115 divisors

QBasic

Works with: QBasic version 1.1
Works with: QuickBasic version 4.5
limite = 50000
DIM DSum(limite + 1), DCount(limite + 1)

FOR i = 0 TO limite
    DSum(i) = 1
    DCount(i) = 1
NEXT i

FOR i = 2 TO limite
    j = i + i
    WHILE j <= limite
        IF DSum(j) = j THEN
            PRINT USING "######## equals the sum of its first ### divisors"; j; DCount(j)
        END IF
        DSum(j) = DSum(j) + i
        DCount(j) = DCount(j) + 1
        j = j + i
    WEND
NEXT i

Run BASIC

Works with: Just BASIC
Works with: Liberty BASIC
limite = 1000000-1
'Maximum array size is 1,000,001 elements
dim DSum(limite+1)
dim DCount(limite+1)

for i = 0 to limite
    DSum(i) = 1
    DCount(i) = 1
next i

for i = 2 to limite
    j = i + i
    while j <= limite
        if DSum(j) = j then
            print using("########", j); " equals the sum of its first"; using("###", DCount(j)); " divisors"
        end if
        DSum(j) = DSum(j) + i
        DCount(j) = DCount(j) + 1
        j = j + i
    wend
next i

PureBasic

Translation of: FreeBASIC
OpenConsole()

limite.l = 2e6
Dim DSum.l(limite+1)
Dim DCount.l(limite+1)

For i.l = 0 To limite
  DSum(i) = 1
  DCount(i) = 1
Next i

For i = 2 To limite
  j.l = i + i
  While j <= limite
    If DSum(j) = j
      PrintN(RSet(Str(j), 8) + " equals the sum of its first " + RSet(Str(DCount(j)), 3) + " divisors")
    EndIf
    DSum(j) = DSum(j) + i
    DCount(j) = DCount(j) + 1
    j = j + i
  Wend
Next i

PrintN(#CRLF$ + "--- terminado, pulsa RETURN---"): Input()
CloseConsole()
Output:
Same as FreeBASIC entry.

Yabasic

limite = 2e6
dim DSum(limite+1), DCount(limite+1)

for i = 0 to limite
    DSum(i) = 1
    DCount(i) = 1
next i

for i = 2 to limite
    j = i + i
    while j <= limite
        if DSum(j) = j  print j using ("########"), " equals the sum of its first ", DCount(j) using ("###"), " divisors"
        DSum(j) = DSum(j) + i
        DCount(j) = DCount(j) + 1
        j = j + i
    wend
next i
Output:
Same as FreeBASIC entry.

C

Translation of: C++

Run time about 46 seconds.

#include <stdio.h>
#include <stdlib.h>

int main() {
    const int maxNumber = 100000000;
    int *dsum = (int *)malloc((maxNumber + 1) * sizeof(int));
    int *dcount = (int *)malloc((maxNumber + 1) * sizeof(int));
    int i, j;
    for (i = 0; i <= maxNumber; ++i) {
        dsum[i] = 1;
        dcount[i] = 1;
    }
    for (i = 2; i <= maxNumber; ++i) {
        for (j = i + i; j <= maxNumber; j += i) {
            if (dsum[j] == j) {
                printf("%8d equals the sum of its first %d divisors\n", j, dcount[j]);
            }
            dsum[j] += i;
            ++dcount[j];
        }
    }
    free(dsum);
    free(dcount);
    return 0;
}
Output:
      24 equals the sum of its first 6 divisors
    2016 equals the sum of its first 31 divisors
    8190 equals the sum of its first 43 divisors
   42336 equals the sum of its first 66 divisors
   45864 equals the sum of its first 66 divisors
  714240 equals the sum of its first 113 divisors
  392448 equals the sum of its first 68 divisors
 1571328 equals the sum of its first 115 divisors
61900800 equals the sum of its first 280 divisors
91963648 equals the sum of its first 142 divisors

slight improvement

dsum and dcount are in "far" distance -> cache miss. Using an struct "typedef struct { int divsum, divcnt;} divs;" improves runtime to 32 s. Calculating the count of divisors only at output saves 50% space and runtime too.

#include <stdio.h>
#include <stdlib.h>

void get_div_cnt(int n){
	int lmt,f,divcnt,divsum;
	divsum = 1;
	divcnt = 1;
	lmt = n/2;
    f = 2;
	for (;;) {
	  if (f > lmt ) break;
	  if (!(n % f)){
		  divsum +=f;
		  divcnt++;
	  }
	  if (divsum == n) break;
      f++;  
	}
    printf("%8d equals the sum of its first %d divisors\n", n, divcnt);	
}

int main() {
    const int maxNumber = 100*1000*1000;
    int *dsum = (int *)malloc((maxNumber + 1) * sizeof(int));
    int i, j;
    for (i = 0; i <= maxNumber; ++i) {
        dsum[i] = 1;
    }
    for (i = 2; i <= maxNumber; ++i) {
        for (j = i + i; j <= maxNumber; j += i) {
        if (dsum[j] == j) get_div_cnt(j);
        dsum[j] += i;
        }
    }
    free(dsum);
    return 0;
}
TIO.RUN:
 the same.
compiler flags -march=native -O2
Real time: 23.893 s User time: 23.475 s Sys. time: 0.203 s CPU share: 99.10 %
//Sys. time: up to 8s for version before???
//with lmt 2,000,000 the results on TIO.RUN are more stable.
version before
User time: 0.304 s CPU share: 98.93 %
slight improvement version
User time: 0.139 s CPU share: 98.80 %

C#

Translation of: Java
using System;

class ErdosNicolasNumbers
{
    static void Main(string[] args)
    {
        const int limit = 100_000_000;
        
        int[] divisorSum = new int[limit + 1];
        int[] divisorCount = new int[limit + 1];
        for (int i = 0; i <= limit; i++)
        {
            divisorSum[i] = 1;
            divisorCount[i] = 1;
        }
        
        for (int index = 2; index <= limit / 2; index++)
        {
            for (int number = 2 * index; number <= limit; number += index)
            {
                if (divisorSum[number] == number)
                {
                    Console.WriteLine($"{number,8} equals the sum of its first {divisorCount[number],3} divisors");
                }
                
                divisorSum[number] += index;
                divisorCount[number]++;
            }
        }
    }
}
Output:
      24 equals the sum of its first   6 divisors
    2016 equals the sum of its first  31 divisors
    8190 equals the sum of its first  43 divisors
   42336 equals the sum of its first  66 divisors
   45864 equals the sum of its first  66 divisors
  714240 equals the sum of its first 113 divisors
  392448 equals the sum of its first  68 divisors
 1571328 equals the sum of its first 115 divisors
61900800 equals the sum of its first 280 divisors
91963648 equals the sum of its first 142 divisors

C++

Translation of: ALGOL 68
#include <iomanip>
#include <iostream>
#include <vector>

int main() {
    const int max_number = 100000000;
    std::vector<int> dsum(max_number + 1, 1);
    std::vector<int> dcount(max_number + 1, 1);
    for (int i = 2; i <= max_number; ++i) {
        for (int j = i + i; j <= max_number; j += i) {
            if (dsum[j] == j) {
                std::cout << std::setw(8) << j
                          << " equals the sum of its first " << dcount[j]
                          << " divisors\n";
            }
            dsum[j] += i;
            ++dcount[j];
        }
    }
}
Output:
      24 equals the sum of its first 6 divisors
    2016 equals the sum of its first 31 divisors
    8190 equals the sum of its first 43 divisors
   42336 equals the sum of its first 66 divisors
   45864 equals the sum of its first 66 divisors
  714240 equals the sum of its first 113 divisors
  392448 equals the sum of its first 68 divisors
 1571328 equals the sum of its first 115 divisors
61900800 equals the sum of its first 280 divisors
91963648 equals the sum of its first 142 divisors

Delphi

Works with: Delphi version 6.0

Translation from the C version. Runs in 38 seconds

const MaxNumber = 100000000;
var DSum: array [0..MaxNumber-1] of integer;
var DCount: array [0..MaxNumber-1] of integer;

procedure ShowErdosNicolasNumbers(Memo: TMemo);
var I,J: integer;
begin
for I:=0 to MaxNumber-1 do
	begin
	DSum[I]:=1;
	DCount[I]:=1;
	end;
for I:=2 to MaxNumber-1 do
	begin
	J:=I*2;
	while J<MaxNumber do
		begin
		if dsum[J] = j then
			begin
			Memo.Lines.Add(Format('%8d equals the sum of its first %d divisors', [j, dcount[j]]));
			end;
		Inc(dsum[J],I);
		Inc(DCount[J]);
		Inc(J,I);
		end;
        end;
end;
Output:
      24 equals the sum of its first 6 divisors
    2016 equals the sum of its first 31 divisors
    8190 equals the sum of its first 43 divisors
   42336 equals the sum of its first 66 divisors
   45864 equals the sum of its first 66 divisors
  714240 equals the sum of its first 113 divisors
  392448 equals the sum of its first 68 divisors
 1571328 equals the sum of its first 115 divisors
61900800 equals the sum of its first 280 divisors
91963648 equals the sum of its first 142 divisors

EasyLang

Translation of: FreeBASIC
limit = 2000000
for i to limit
   dsum[] &= 1
   dcnt[] &= 1
.
for i = 2 to limit
   j = i + i
   while j <= limit
      if dsum[j] = j
         print j & " equals the sum of its first " & dcnt[j] & " divisors"
      .
      dsum[j] += i
      dcnt[j] += 1
      j += i
   .
.

FutureBasic

_limit = 500000

void local fn ErdosNicolasNumbers
  long i, j, sum( _limit ), count( _limit )
  
  for i = 0 to _limit
    sum(i) = 1
    count(i) = 1
  next
  
  for i = 2 to _limit
    j = i + i
    while ( j <= _limit )
      if sum(j) == j then printf @"%8ld == sum of its first %3ld divisors", j, count(j)
      sum(j) = sum(j) + i
      count(j) = count(j) + 1
      j = j + i
    wend
  next
end fn

fn ErdosNicolasNumbers

HandleEvents
Output:
      24 == sum of its first   6 divisors
    2016 == sum of its first  31 divisors
    8190 == sum of its first  43 divisors
   42336 == sum of its first  66 divisors
   45864 == sum of its first  66 divisors
  392448 == sum of its first  68 divisors
 

Go

Translation of: C++

This runs in about 30 seconds which, surprisingly, is faster than C++ (43 seconds) - usually Go is a good bit slower.

package main

import "fmt"

func main() {
	const maxNumber = 100000000
	dsum := make([]int, maxNumber+1)
	dcount := make([]int, maxNumber+1)
	for i := 0; i <= maxNumber; i++ {
		dsum[i] = 1
		dcount[i] = 1
	}
	for i := 2; i <= maxNumber; i++ {
		for j := i + i; j <= maxNumber; j += i {
			if dsum[j] == j {
				fmt.Printf("%8d equals the sum of its first %d divisors\n", j, dcount[j])
			}
			dsum[j] += i
			dcount[j]++
		}
	}
}
Output:
      24 equals the sum of its first 6 divisors
    2016 equals the sum of its first 31 divisors
    8190 equals the sum of its first 43 divisors
   42336 equals the sum of its first 66 divisors
   45864 equals the sum of its first 66 divisors
  714240 equals the sum of its first 113 divisors
  392448 equals the sum of its first 68 divisors
 1571328 equals the sum of its first 115 divisors
61900800 equals the sum of its first 280 divisors
91963648 equals the sum of its first 142 divisors

J

Implementation:

divisors=: {{ /:~ ,*/@> { (^ i.@>:)&.>/__ q: y}} ::_:
erdosnicolas=: {{ y e. +/\ _2}. divisors y }}"0

Task example:

   I.erdosnicolas i.1e7
24 2016 8190 42336 45864 392448 714240 1571328
   (,. 1++/\@divisors i. ])@>24 2016 8190 42336 45864 392448 714240 1571328
     24   6
   2016  31
   8190  43
  42336  66
  45864  66
 392448  68
 714240 113
1571328 115

Java

import java.util.Arrays;

public final class ErdosNicolasNumbers {

	public static void main(String[] aArgs) {
		final int limit = 100_000_000;
		
	    int[] divisorSum = new int[limit + 1];
	    int[] divisorCount = new int[limit + 1];
	    Arrays.fill(divisorSum, 1);
	    Arrays.fill(divisorCount, 1);
	    
	    for ( int index = 2; index <= limit / 2; index++ ) {
	        for ( int number = 2 * index; number <= limit; number += index ) {
	            if ( divisorSum[number] == number ) {
	                System.out.println(String.format("%8d", number) + " equals the sum of its first "
	                	             + String.format("%3d", divisorCount[number]) + " divisors");
	            }
	            
	            divisorSum[number] += index;
	            divisorCount[number]++;
	        }
	    }
	}

}
Output:
      24 equals the sum of its first   6 divisors
    2016 equals the sum of its first  31 divisors
    8190 equals the sum of its first  43 divisors
   42336 equals the sum of its first  66 divisors
   45864 equals the sum of its first  66 divisors
  714240 equals the sum of its first 113 divisors
  392448 equals the sum of its first  68 divisors
 1571328 equals the sum of its first 115 divisors
61900800 equals the sum of its first 280 divisors
91963648 equals the sum of its first 142 divisors

jq

Adapted from #Wren

Works with jq and gojq, the C and Go implementations of jq

The following program will also work using jaq provided `sqrt` is defined appropriately and other minor adjustments are made.

# Output a stream of the (unsorted) proper divisors of . including 1
def proper_divisors:
  . as $n
  | if $n > 1 then 1,
      ( range(2; 1 + sqrt) as $i
        | if ($n % $i) == 0 then $i,
            (($n / $i) | if . == $i then empty else . end)
         else empty
         end)
    else empty
    end;

# Emit k if . is an Erdos-Nicolas number, otherwise emit 0
def erdosNicolas:
  . as $n
  | ([proper_divisors] | sort) as $divisors
  | ($divisors|length) as $dc
  | if $dc < 3 then 0
    else {sum: ($divisors[0] + $divisors[1])}
    # An Erdos-Nicolas is not perfect, and hence $dc-1 in the following line:
    | first(
        foreach range(2; $dc-1) as $i (.;
          .sum += $divisors[$i]
          | if .sum == $n then .emit = $i + 1
            elif .sum > $n then .emit = 0
            else .
            end )
          | select(.emit).emit  ) // 0
    end ;

limit(8;
  range(2; infinite)
  | . as $n 
  | erdosNicolas as $k
  | select($k > 0)
  | "\($n) from \($k)" )
Output:
24 from 6
2016 from 31
8190 from 43
42336 from 66
45864 from 66
392448 from 68
714240 from 113
1571328 from 115

Julia

using Primes

function isErdősNicolas_with_k(n)
    @assert n > 2
    d = [one(n)]
    for (p, e) in eachfactor(n)
        d = reduce(vcat, [d * p^j for j in 1:e], init=d)
    end
    sort!(d)
    pop!(d)
    len = length(d)
    (len < 2 || sum(d) <= n) && return false, 0
    for k in 2:len
        sum(@view d[1:k]) == n && return true, k
    end
    return false, 0
end

for n in 3:2_000_000
    isEN, k = isErdősNicolas_with_k(n)
    isEN && println(lpad(n, 8), " equals the sum of its first $k divisors.")
end
Output:
      24 equals the sum of its first 6 divisors.
    2016 equals the sum of its first 31 divisors.
    8190 equals the sum of its first 43 divisors.
   42336 equals the sum of its first 66 divisors.
   45864 equals the sum of its first 66 divisors.
  392448 equals the sum of its first 68 divisors.
  714240 equals the sum of its first 113 divisors.
 1571328 equals the sum of its first 115 divisors.

Nim

Translation of: Go

To get better performances, we use 32 bits integers rather than 64 bits integers. We got the best (and excellent) execution time with compilation command: nim c -d:release -d:lto --gc:boehm erdos_nicolas_numbers.nim

import std/[sequtils, strformat]

proc main() =
  const MaxNumber = 100_000_000i32
  var dsum, dcount = repeat(1'i32, MaxNumber + 1)
  for i in 2i32..MaxNumber:
    for j in countup(i + i, MaxNumber, i):
      if dsum[j] == j:
        echo &"{j:>8} equals the sum of its first {dcount[j]} divisors"
      inc dsum[j], i
      inc dcount[j]
main()
Output:
      24 equals the sum of its first 6 divisors
    2016 equals the sum of its first 31 divisors
    8190 equals the sum of its first 43 divisors
   42336 equals the sum of its first 66 divisors
   45864 equals the sum of its first 66 divisors
  714240 equals the sum of its first 113 divisors
  392448 equals the sum of its first 68 divisors
 1571328 equals the sum of its first 115 divisors
61900800 equals the sum of its first 280 divisors
91963648 equals the sum of its first 142 divisors

Pascal

Free Pascal

using Factors_of_an_integer#using_Prime_decomposition , using a sort of sieve.

program ErdoesNumb;

// gets factors of consecutive integers fast
// limited to 1.2e11
{$IFDEF FPC}
  {$MODE DELPHI}  {$OPTIMIZATION ON,ALL}  {$COPERATORS ON}
{$ELSE}
  {$APPTYPE CONSOLE}
{$ENDIF}
uses
  sysutils
{$IFDEF WINDOWS},Windows{$ENDIF}
  ;
//######################################################################
//prime decomposition
const
//HCN(86) > 1.2E11 = 128,501,493,120     count of divs = 4096   7 3 1 1 1 1 1 1 1
  HCN_DivCnt  = 4096;
type
  tItem     = Uint64;
  tDivisors = array [0..HCN_DivCnt] of tItem;
  tpDivisor = pUint64;
const
  SizePrDeFe = 32768;//*56 <= 64kb level I or 2 Mb ~ level 2 cache or more
type
  tdigits = array [0..31] of Uint32;
  //the first number with 11 different prime factors =
  //2*3*5*7*11*13*17*19*23*29*31 = 2E11
  //56 byte
  tprimeFac = packed record
                 pfSumOfDivs,
                 pfRemain  : Uint64;
                 pfDivCnt  : Uint32;
                 pfMaxIdx  : Uint32;
                 pfpotMax  : array[0..11] of byte;
                 pfpotPrimIdx : array[0..9] of word;
               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 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 : 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
      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 := (Uint64(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;

procedure InsertSort(pDiv:tpDivisor; Left, Right : NativeInt );
var
  I, J: NativeInt;
  Pivot : tItem;
begin
  for i:= 1 + Left to Right do
  begin
    Pivot:= pDiv[i];
    j:= i - 1;
    while (j >= Left) and (pDiv[j] > Pivot) do
    begin
      pDiv[j+1]:=pDiv[j];
      Dec(j);
    end;
    pDiv[j+1]:= pivot;
  end;
end;

procedure GetDivisors(pD:tpPrimeFac;var Divs:tDivisors);
var
  pDivs : tpDivisor;
  pPot : UInt64;
  i,len,j,l,p,k: Int32;
Begin
  pDivs := @Divs[0];
  pDivs[0] := 1;
  len := 1;
  l   := 1;
  with pD^ do
  Begin
    For i := 0 to pfMaxIdx-1 do
    begin
      //Multiply every divisor before with the new primefactors
      //and append them to the list
      k := pfpotMax[i];
      p := SmallPrimes[pfpotPrimIdx[i]];
      pPot :=1;
      repeat
        pPot *= p;
        For j := 0 to len-1 do
        Begin
          pDivs[l]:= pPot*pDivs[j];
          inc(l);
        end;
        dec(k);
      until k<=0;
      len := l;
    end;
    p := pfRemain;
    If p >1 then
    begin
      For j := 0 to len-1 do
      Begin
        pDivs[l]:= p*pDivs[j];
        inc(l);
      end;
      len := l;
    end;
  end;
  //Sort. Insertsort much faster than QuickSort in this special case
  InsertSort(pDivs,0,len-1);
  //end marker
  pDivs[len] :=0;
end;

var
  pPrimeDecomp :tpPrimeFac;
  Divs:tDivisors;
  T0:Int64;
  n,s : NativeUInt;
  i : Int32;
Begin
  T0 := GetTickCount64;
  InitSmallPrimes;
  Init_Sieve(0);
  //jump over 0
  pPrimeDecomp:= GetNextPrimeDecomp;
  n := 1;
  repeat
    pPrimeDecomp:= GetNextPrimeDecomp;
    s := pPrimeDecomp^.pfSumOfDivs;
    if  (s > 2*n) then
    begin
      s -= n;
      // 75% of runtime
      GetDivisors(pPrimeDecomp,Divs);
      //calculate downwards. Not really an impact
      For i := pPrimeDecomp^.pfDivCnt-2 downto 0 do
      Begin
        s -= Divs[i];
        if s<n then
          break;
        if s = n then
          writeln(Format('%8d equals the sum of its first %4d divisors',
                         [n,i]));
      end;
    end;
    n += 1;
  until n > 100*1000*1000+1;
  T0 := GetTickCount64-T0;
  writeln('runtime ',T0/1000:0:3,' s');
end.
@TIO.RUN:
      24 equals the sum of its first    6 divisors
    2016 equals the sum of its first   31 divisors
    8190 equals the sum of its first   43 divisors
   42336 equals the sum of its first   66 divisors
   45864 equals the sum of its first   66 divisors
  392448 equals the sum of its first   68 divisors
  714240 equals the sum of its first  113 divisors
 1571328 equals the sum of its first  115 divisors
61900800 equals the sum of its first  280 divisors
91963648 equals the sum of its first  142 divisors
runtime 15.760 s

Real time: 15.909 s User time: 15.721 s CPU share: 99.09 %

PascalABC.NET

Translation of: C++
const
  max_number = 100_000_000;

begin
  var dsum := arrfill(max_number + 1, 1);
  var dcount := arrfill(max_number + 1, 1);
  for var i := 2 to max_number do
    for var j := i + i to max_number step i do 
    begin
      if dsum[j] = j then 
        writeln(j:8, ' equals the sum of its first ', dcount[j], ' divisors');
      dsum[j] += i;
      dcount[j] += 1;
    end;
end.
Output:
    2016 equals the sum of its first 31 divisors
    8190 equals the sum of its first 43 divisors
   42336 equals the sum of its first 66 divisors
   45864 equals the sum of its first 66 divisors
  714240 equals the sum of its first 113 divisors
  392448 equals the sum of its first 68 divisors
 1571328 equals the sum of its first 115 divisors
61900800 equals the sum of its first 280 divisors
91963648 equals the sum of its first 142 divisors

Perl

Library: ntheory
use v5.36;
use enum    qw(False True);
use ntheory qw(divisors divisor_sum);

sub is_Erdos_Nicolas ($n) {

    divisor_sum($n) > 2*$n or return False;

    my $sum   = 0;
    my $count = 0;

    foreach my $d (divisors($n)) {
        ++$count;        $sum += $d;
        return $count if ($sum == $n);
        return False  if ($sum > $n);
    }
}

my ($n, $count) = (2, 0);
until ($count == 8) {
    next unless 0 == ++$n % 2;
    if (my $key = is_Erdos_Nicolas $n) {
        printf "%8d == sum of its first %3d divisors\n", $n, $key;
        $count++;
    }
}
Output:
      24 == sum of its first   6 divisors
    2016 == sum of its first  31 divisors
    8190 == sum of its first  43 divisors
   42336 == sum of its first  66 divisors
   45864 == sum of its first  66 divisors
  392448 == sum of its first  68 divisors
  714240 == sum of its first 113 divisors
 1571328 == sum of its first 115 divisors

Phix

Translation of: Wren
with javascript_semantics
function erdos_nicolas(integer n)
    sequence divisors = factors(n)
    integer tot = 1
    for i=1 to length(divisors)-1 do
        tot += divisors[i]
        if tot=n then return i+1 end if
        if tot>n then exit end if
    end for
    return 0
end function
 
constant limit = 8
integer n = 2, count = 0
while count<limit do
    integer k = erdos_nicolas(n)
    if k>0 then
        printf(1,"%8d equals the sum of its first %d divisors.\n",{n,k})
        count += 1
    end if
    n += 2
end while

Aside: The default for factors() is to return neither 1 nor n, though you can change that if you want, ie ",1" -> 1 and n; ",-1" -> 1 but not n.
Output same as Julia

Python

Translation of: C

This is a translation of the "slight improvement" C code. I ran this script using PyPy7.3.13 (Python3.10.13) and it completes in ~10.5s on a AMD Ryzen 7 7800X3D.

# erdos-nicolas.py by Xing216
from time import perf_counter
start = perf_counter()
def get_div_cnt(n: int) -> None:
    divcnt,divsum = 1,1
    lmt = n/2
    f = 2
    while True:
        if f > lmt: break
        if not (n % f):
            divsum += f
            divcnt += 1
        if divsum == n: break
        f+=1
    print(f"{n:>8} equals the sum of its first {divcnt} divisors")
max_number = 91963649
dsum = [1 for _ in range(max_number+1)]
for i in range(2, max_number + 1):
    for j in range(i + i, max_number + 1, i):
        if (dsum[j] == j): get_div_cnt(j)
        dsum[j] += i
done = perf_counter() - start
print(f"Done in: {done:.3f}s")
Output:
      24 equals the sum of its first 6 divisors
    2016 equals the sum of its first 31 divisors
    8190 equals the sum of its first 43 divisors
   42336 equals the sum of its first 66 divisors
   45864 equals the sum of its first 66 divisors
  714240 equals the sum of its first 113 divisors
  392448 equals the sum of its first 68 divisors
 1571328 equals the sum of its first 115 divisors
61900800 equals the sum of its first 280 divisors
91963648 equals the sum of its first 142 divisors
Done in: 10.478s

Raku

use Prime::Factor;

sub is-Erdős-Nicolas ($n) {
    my @divisors = $n.&proper-divisors: :s;
    ((@divisors.sum > $n) && (my $key = ([\+] @divisors).first: $n, :k)) ?? 1 + $key !! False
}

my $count;

(1..*).hyper(:2000batch).map( * × 2 ).map: {
    if my $key = .&is-Erdős-Nicolas {
        printf "%8d == sum of its first %3d divisors\n", $_, $key;
        exit if ++$count >= 8;
    }
}
Output:
      24 == sum of its first   6 divisors
    2016 == sum of its first  31 divisors
    8190 == sum of its first  43 divisors
   42336 == sum of its first  66 divisors
   45864 == sum of its first  66 divisors
  392448 == sum of its first  68 divisors
  714240 == sum of its first 113 divisors
 1571328 == sum of its first 115 divisors

Ring

see "works..." + nl
erdos = []
limit = 1600000
for n = 1 to limit
    num = 0
    sum = 0
    erdos = []
    for m = 1 to n/2
        if n%m = 0
           add(erdos,m)
        ok
    next
    lenErdos = len(erdos)
    for p = 1 to lenErdos 
        sum = sum + erdos[p]
        if sum = n and p < lenErdos 
           num++
           see "" + n + " equals the sum of its first " + p + " divisors" + nl
           exit
        ok
    next
    if num = 8
       exit
    ok
next
see "done..." + nl
Output:
        24 equals the sum of its first 6 divisors
      2016 equals the sum of its first 31 divisors
      8190 equals the sum of its first 43 divisors
     42336 equals the sum of its first 66 divisors
     45864 equals the sum of its first 66 divisors
    714240 equals the sum of its first 113 divisors
    392448 equals the sum of its first 68 divisors
   1571328 equals the sum of its first 115 divisors

Sidef

func is_Erdős_Nicolas(n) {

    n.is_abundant || return false

    var sum   = 0
    var count = 0

    n.divisors.each {|d|
        ++count;         sum += d
        return count if (sum == n)
        return false if (sum >  n)
    }
}

var count = 8   # how many terms to compute

^Inf -> by(2).each {|n|
    if (is_Erdős_Nicolas(n)) { |v|
        say "#{'%8s'%n} is the sum of its first #{'%3s'%v} divisors"
        --count || break
    }
}
Output:
      24 is the sum of its first   6 divisors
    2016 is the sum of its first  31 divisors
    8190 is the sum of its first  43 divisors
   42336 is the sum of its first  66 divisors
   45864 is the sum of its first  66 divisors
  392448 is the sum of its first  68 divisors
  714240 is the sum of its first 113 divisors
 1571328 is the sum of its first 115 divisors

Wren

Version 1

Library: Wren-math
import "./math" for Int

var erdosNicolas = Fn.new { |n|
    var divisors = Int.properDivisors(n) // excludes n itself
    var dc = divisors.count
    if (dc < 3) return 0
    var sum = divisors[0] + divisors[1]
    for (i in 2...dc-1) {
        sum = sum + divisors[i]
        if (sum == n) return i + 1
        if (sum > n)  break
    }
    return 0
}

var limit = 8
var n = 2
var count = 0
while (true) {
    var k = erdosNicolas.call(n)
    if (k > 0) {
        System.print("%(n) from %(k)")
        count = count + 1
        if (count == limit) return
    }
    n = n + 2
}
Output:
24 from 6
2016 from 31
8190 from 43
42336 from 66
45864 from 66
392448 from 68
714240 from 113
1571328 from 115

Version 2

Translation of: C++
Library: Wren-fmt

This takes 7 minutes to run (about 10 times slower than C++) but finds 2 more numbers, albeit not in strictly increasing order.

If `maxNum` is set to 2 million, then it finds the first 8 in about 5.2 seconds which is more than 10 times faster than Version 1's 58 seconds.

import "./fmt" for Fmt

var maxNum = 1e8
var dsum = List.filled(maxNum+1,1)
var dcount = List.filled(maxNum+1, 1)
for (i in 2..maxNum) {
    var j = i + i
    while (j <= maxNum) {
        if (dsum[j] == j) {
            Fmt.print("$8d equals the sum of its first $d divisors", j, dcount[j])
        }
        dsum[j] = dsum[j] + i
        dcount[j] = dcount[j] + 1
        j = j + i
    }
}
Output:
      24 equals the sum of its first 6 divisors
    2016 equals the sum of its first 31 divisors
    8190 equals the sum of its first 43 divisors
   42336 equals the sum of its first 66 divisors
   45864 equals the sum of its first 66 divisors
  714240 equals the sum of its first 113 divisors
  392448 equals the sum of its first 68 divisors
 1571328 equals the sum of its first 115 divisors
61900800 equals the sum of its first 280 divisors
91963648 equals the sum of its first 142 divisors

XPL0

Translation of Algol 68 and C++.

def Max = 2_000_000;
int DSum(Max+1), DCount(Max+1);
int I, J;
[for I:= 0 to Max do
    [DSum(I):= 1;
    DCount(I):= 1;
    ];
Format(8, 0);
for I:= 2 to Max do
    [J:= I+I;
    while J <= Max do
        [if DSum(J) = J then
            [RlOut(0, float(J));
            Text(0, " equals the sum of its first ");
            IntOut(0, DCount(J));
            Text(0, " divisors^j^m");
            ];
        DSum(J):= DSum(J)+I;
        DCount(J):= DCount(J)+1;
        J:= J+I;
        ]
    ]
]
Output:
      24 equals the sum of its first 6 divisors
    2016 equals the sum of its first 31 divisors
    8190 equals the sum of its first 43 divisors
   42336 equals the sum of its first 66 divisors
   45864 equals the sum of its first 66 divisors
  714240 equals the sum of its first 113 divisors
  392448 equals the sum of its first 68 divisors
 1571328 equals the sum of its first 115 divisors