The Holy-Grail Formula For Primes Generation

Students at the University of Buenos Aires have published a very-short paper this month (November, 2020) providing an astonishing formula to generate all primes - in sequence.

This has been the Holy Grail in Number Theory: is there a way to generate all the prime numbers one by one with a deterministic formula? The students have proved that there is. Their paper is here: 2010.15882.pdf (arxiv.org)

The formula is actually very simple. Here it is:

And F1 should be a constant which "magically" is

F1 = 2.9200509773161350318170795639

This constant is being called now the "Buenos Aires Prime Constant". The only problem is how to find this constant without knowing a priori all the primes. This is the current open problem. If someone can extend this constant to a very-very-very large number of decimal points, the problem of primes generation will be solved. With the current one, though, it is good enough to generate few primes, in order. The code below (using Miller-Rabin Primality Test) makes a codification of the paper, and the initial results are indeed surprising. Check out the results, and the code. Cheers, ACC. 

Buenos-Aires Prime Constant: 2.9200509773161350318170795639

Prime: 2

Prime: 3

Prime: 5

Prime: 7

Prime: 11

Prime: 13

Prime: 17

Prime: 19

Prime: 23

Prime: 29

Prime: 31

Prime: 37

Prime: 41

Prime: 43

And so on...

using System;
using System.Numerics;

namespace PrimeGenerator
{
    class Program
    {
        static void Main(string[] args)
        {
            Decimal f = FindingBuenosAiresConstant();
            Console.WriteLine("Buenos-Aires Prime Constant: {0}", f);
            BigInteger number = 0;
            while(Check(f, out number))
            {
                Console.WriteLine("Prime: {0}", number);
                f = Math.Floor(f) * (f - Math.Floor(f) + 1);
            }
            Console.WriteLine("And so on...");
        }

        public static Decimal FindingBuenosAiresConstant()
        {
            int MAX = 2000;
            BigInteger[] primes = new BigInteger[MAX];
            int index = 0;
            BigInteger candidate = 2;
            while (index < primes.Length)
            {
                if (IsPrimeMillerRabin(candidate))
                {
                    primes[index++] = candidate;
                }
                candidate++;
            }

            //Calculate
            Decimal bac = (Decimal)primes[0] - 1;
            BigInteger den = 1;
            for (int i = 1; i < MAX; i++)
            {
                Decimal partial = (Decimal)primes[i] - 1;
                den *= primes[i - 1];
                partial = (Decimal)((Double)partial / (Double)den);
                bac += partial;
            }

            return bac;
        }

        public static bool Check(Decimal d, out BigInteger number)
        {
            number = BigInteger.Parse(d.ToString().Substring(0, d.ToString().IndexOf(".")));
            return IsPrimeMillerRabin(number);
        }

        public static bool IsPrimeMillerRabin(BigInteger n)
        {
            //It does not work well for smaller numbers, hence this check
            int SMALL_NUMBER = 1000;

            if (n <= SMALL_NUMBER)
            {
                return IsPrime(n);
            }

            int MAX_WITNESS = 100;
            for (long i = 2; i <= MAX_WITNESS; i++)
            {
                if (IsPrime(i) && Witness(i, n) == 1)
                {
                    return false;
                }
            }

            return true;
        }

        public static BigInteger SqRtN(BigInteger N)
        {
            /*++
             *  Using Newton Raphson method we calculate the
             *  square root (N/g + g)/2
             */
            BigInteger rootN = N;
            int count = 0;
            int bitLength = 1;
            while (rootN / 2 != 0)
            {
                rootN /= 2;
                bitLength++;
            }
            bitLength = (bitLength + 1) / 2;
            rootN = N >> bitLength;

            BigInteger lastRoot = BigInteger.Zero;
            do
            {
                if (lastRoot > rootN)
                {
                    if (count++ > 1000)                   // Work around for the bug where it gets into an infinite loop
                    {
                        return rootN;
                    }
                }
                lastRoot = rootN;
                rootN = (BigInteger.Divide(N, rootN) + rootN) >> 1;
            }
            while (!((rootN ^ lastRoot).ToString() == "0"));
            return rootN;
        }

        public static bool IsPrime(BigInteger n)
        {
            if (n <= 1)
            {
                return false;
            }

            if (n == 2)
            {
                return true;
            }

            if (n % 2 == 0)
            {
                return false;
            }

            for (int i = 3; i <= SqRtN(n) + 1; i += 2)
            {
                if (n % i == 0)
                {
                    return false;
                }
            }
            return true;
        }

        private static int Witness(long a, BigInteger n)
        {
            BigInteger t, u;
            BigInteger prev, curr = 0;
            BigInteger i;
            BigInteger lln = n;

            u = n / 2;
            t = 1;
            while (u % 2 == 0)
            {
                u /= 2;
                t++;
            }

            prev = BigInteger.ModPow(a, u, n);
            for (i = 1; i <= t; i++)
            {
                curr = BigInteger.ModPow(prev, 2, lln);
                if ((curr == 1) && (prev != 1) && (prev != lln - 1)) return 1;
                prev = curr;
            }
            if (curr != 1) return 1;
            return 0;
        }
    }
}