Exploring Computational Number Theory (Part 1)

• Download project files - 21 MB

Introduction

This article explores computational number theory and the relationships between the various classical number theorists’ theories by using Euler Pseudoprimes to bridge these theories in a base 2 computational environment while providing the user the information in a human readable form to enhance understanding in an exploratory environment.

To accomplish this in a VS2005 .NET environment involved the use and modification of the following:

1. "BigInteger Library" by Mihnea Rădulescu, 23 Sep 2011 at http://www.codeproject.com/Articles/60108/BigInteger-Library. This work was modified to support arbitrary size big integers and includes additional algorithms to support this project.
2. "C# BigInteger Class" by Chew Keong TAN, 28 Sep 2002 at http://www.codeproject.com/Articles/2728/C-BigInteger-Class. This work supplied most of the additional algorithms for supplementing the "BigInteger Library".
3. "Using C# .NET User Defined Functions (UDF) in Excel" by Adam Tibi, 13 Jun 2013 at http://www.codeproject.com/Articles/606446/UsingplusC-plus-NETplusUserplusDefinedplusFuncti was used to link the computational environment to the Excel exploratory environment.
4. Excel-DNA Version 0.32 can be obtained from here http://exceldna.codeplex.com/releases/view/119190 .
5. "Performance Tests: Precise Run Time Measurements with System.Diagnostics.Stopwatch" by Thomas Maierhofer (Tom), 28 Feb 2010 at http://www.codeproject.com/Articles/61964/Performance-Tests-Precise-Run-Time-Measurements-wi were also used to provide the user with an analytical capability within this environment.

Background

The following web posting in 2001, "All Fermat and Mersenne numbers are squarefree" by Jacques BASALDÚA have deeply influenced this article as well as the well known classical number theorists:

	Marin Mersenne 1588-1648 Definition of Mersenne Number: A Mersenne Number is of the form, Mn = 2n – 1 where n is a positive integer. Definition of Mersenne Prime: Mersenne Primes are of the form, Mp = 2p - 1 where p and 2p - 1 is prime.
	Pierre de Fermat 1601-1665 Definition of Fermat Number: A Fermat number is of the form Fn = 22^n+1 where n is a positive integer. Fermat's Little Theorem. Let p be prime and suppose p | b. Then bp-1 º 1 (mod p)
	Leonhard Euler 1707-1783 Definition of Euler Pseudoprime: An Euler pseudoprime to base 2 is of the form En = 1 mod n and is the largest number of the form En = c*n+1 where c is a positive integer and n is the least positive exponent base 2. All factors are of this form. (see Germain) Euler-Fermat Theorem. If gcd(b, n) = 1 then bf(n) º 1(mod n).
	Sophie Germain 1776-1831 Definition of Germain Pseudoprime: A Germain number Gn = c*n+1 where c is a positive integer and n is the least positive exponent base 2. If n is odd or n = 2x, then Gn = 2k*n+1 or c = 2k. If n is even, excepting n = 2x, then Gn = c * n + 1. All factors are of this form.
	Francois Lucas 1842-1891 Definition of Fermat-Lucas Number: A Fermat-Lucas number is of the form FLn = 2n+1 where n is a positive integer. Definition of Order (ord): Let gcd(b, n) = 1. The least positive integer m which bm º 1(mod n) is the order of b modulo n, denoted ordnb.

Many of us from our high school and college days are familiar with the classical number theorists and their formulas. This article will codify the relationship of their formulas and their Euler psuedoprime attributes. For instance, M_p = E_p (see Exhibit 1 below) which has one prime base, but has different prime bases which creates a binary sequence of (p-1) ... 0. Similarly, F_n = 2^{2^n}+1 is equivalent to the class of Euler psuedoprimes as J = 2*2ⁿ or 2ⁿ⁺¹ then E_j = F_n has one prime base of 2 and creates a binary sequence of bit n + bit 0 where bit n is a multiplicative function of 2ⁿ⁺¹. This also applies to other prime bases. For example, E₉ = 73 has one prime base of 3 which has 2 bits set {1, 0} in a computational environment. However, in a mathematical environment the two bits are {2, 1}. When multiplied by 3 because 9 = 3², the sequence becomes bits {6, 3} + bit 0. This approach works with multiple prime bases as well. For example, E₁₅ = 151 has two prime bases, three and five, and the prime base sequence is intertwined in the following sequence pattern. The sequence initiation is 15 = 3 * 5 and 3 + 5 = 8 so the first bit set is 15 – 8 = 7 then 7 – 3 = 4 then 4 – 3 = 1 then 7 – 5 = 2. The sequence becomes bits {7, 4, 1, 2} + bit 0. Note: E_2n = FL_n

Exhibit 2 below displays the result of the process described above with one and two prime bases and with an exponent of one for each base. It also indicates the relationship of Mersenne numbers when they are reduced to an Euler pseudoprime. Exhibit 3 & 4 below displays similar data as the first chart with exponents greater than or equal to one and indicates Fermat numbers to other bases and reduced to an Euler pseudoprime. Exhibit 5 below displays FermatLucas numbers, a larger class than Fermat numbers, and how they can be reduced to an Euler pseudoprime. In all these exhibits the black set characters are the binary sequence and the colored characters are the integer results.

The code below generates the information in the exhibits above for one and two base sequences to any exponent within the resources of the computer.

private static BigInteger NonOrderedSeq(long p1, long p2, long multiplier)
{
    long h, n;
    BigInteger result = new BigInteger(0);
    BigInteger temp1 = new BigInteger(0);
    BigInteger temp2 = new BigInteger(0);
 
    if (p1 == 1)
        h = p2;
    else
        h = (p1 - 1) * p2;
 
    while (h > 0)
    {
        n = h;
        while (n >= p1)
        {
            n = n - p1;
            if (p1 == 1)
                result = result + BigInteger.setBit(n * multiplier);
            else
                temp1 = temp1 + BigInteger.setBit(n * multiplier);
        }
 
        h = h - p2;
    }
 
    if (result == 0)
        if (multiplier == 1)
            result = temp1 + BigInteger.setBit(0);
        else
        {
            temp2 = temp1 << Convert.ToInt32(multiplier);
            result = temp2 - temp1;
            result = result + BigInteger.setBit(0);
        }
    return result;
}

In Exhibit 5 above, you may have noticed that in some two base cases the sequence has not been fully reduced to an Euler pseudoprime. For example, E₆ = 1 but is indicating 3 and E₂₀ = 41 but is indicating 205. There is an order factor for all bases greater than 1. Exhibit 6 below is a small set of these factors. The factor is determined as follows:

1. Order_factor22 = Number/Factor2^e then divide the sequence by Factor2
2. Order_factor32 = Number/Factor3^e then divide the sequence by Factor3
3. Order_factor42 = Number/Factor4^e then divide the sequence by Factor4
4. etcetera

The following Exhibit 7 provides the general formula for determining the Euler pseudoprime by the number of prime bases.

The following Exhibit 8 & 9 show how the general formula is applied to specific cases.

Notice in the binary sequences in Exhibit 9 above that the multiplicative nature of the results before the division will grow very large. This will diminish the range of exploration in a finite computational environment and slow the results as n of E increases. Exhibit 10 below shows how this can be done through division. Exhibit 11 below shows how the division groupings break into a Pascal Triangle. Exhibit 12 below shows the Pascal Triangle where division moves from the bases of the triangle to the apex which is pairing the sequences for division within the triangle/pyramid structure and is equivalent to multiplication and division though the matrix. However, in the case of A₇ in Exhibit 9 above, the calculation is bounded with a maximum size of bit 24023.

The following code creates the Pascal Triangle for the divisions.

private static object[,] PyramidPairing(long numOfBases)
{
    long j, k;
    object[,] elements = new object[numOfBases, numOfBases];
    for (j = 0; j <= numOfBases-1; j++)
    {
        if (j == 0)
            elements[j, j] = 1;
        else
            for (k = 0; k <= j; k++)
            {
                if (k == 0)
                    elements[j,k] = elements[j-1,k];
                if (k > 0 && k < j)
                    elements[j,k] = Convert.ToInt64(elements[j-1,k-1]) +
                                    Convert.ToInt64(elements[j-1,k]);
                if (k == j)
                    elements[j,k] = elements[j-1,k-1];
            }
    }
    return elements;
}

Euler pseudoprime code snippets

Setup

public static BigInteger Euler(object[,] ArrIn)
{
    #region Euler setup

    object[,] Arr1;
    object[,] Arr2;
    object[,] element;
    BigInteger[] division;
    BigInteger number = new BigInteger(1);
    long i, j, k, nGroup, numDiv, iCount;
    long p1, p2;

    Arr1 = ArrIn;

    if (Arr1.GetUpperBound(0) < 1)
    {
        p1 = 1;
        p2 = Convert.ToInt32(Arr1[0, 0]);
    }
    else
    {
        p1 = Convert.ToInt32(Arr1[0, 0]);
        p2 = Convert.ToInt32(Arr1[1, 0]);
    }

    #endregion Euler setup

BaseA Sequence

#region BaseA Sequence

///
/// Creates the matrixed sequences
/// If one or two base sequence then
/// division[i] with i = 0 contains the result.
///
Arr2 = AExp(Arr1);
division = new BigInteger[Arr2.GetUpperBound(0) + 1];
for (i = 0; i <= Arr2.GetUpperBound(0); i++)
    division[i] = new BigInteger(Number_Theory.NumberTheory.baseA2(p1.ToString(),
                  p2.ToString(), Convert.ToInt64(Arr2[Arr2.GetUpperBound(0)-i, 0])));

#endregion BaseA Sequence

Prime base greater than 2

#region Prime Base Greater Than 2

//  Prime bases greater than 2

if (Arr1.GetUpperBound(0) > 1)
{
    element = PyramidPairing(Arr1.GetUpperBound(0) - 1);
    numDiv = Convert.ToInt64(BigInteger.Power(new BigInteger(2),
             Arr1.GetUpperBound(0) - 1).ToString());
    Boolean[] bDivides = new Boolean[numDiv/2];
    //division = new BigInteger[numDiv];
    numDiv = numDiv - 1;

    ///
    /// Divides up the pyramid from the base to the apex
    /// with Division[i] where i = 0 is the apex.
    ///
    for (i = element.GetUpperBound(0); i >= 0; i--)
    {
        iCount = 0;
        numDiv = 0;
        bDivides = new Boolean[division.GetLength(0) / 2];
        BigInteger[] hold = new BigInteger[division.GetLength(0) / 2];
        for (j = 0; j <= i; j++)
        {
            nGroup = Convert.ToInt64(element[i, j]);
            for (k = 0; k <= nGroup - 1; k++)
            {
                bDivides[iCount] = false;
                if (division[numDiv] % division[numDiv + nGroup] == 0)
                    bDivides[iCount] = true;

                hold[iCount] = division[numDiv] / division[numDiv + nGroup];
                numDiv = numDiv + 1;
                iCount = iCount + 1;
            }
            numDiv = numDiv + nGroup;
        }
        division = new BigInteger[hold.GetLength(0)];
        for (j = 0; j <= hold.GetUpperBound(0); j++)
            division[j] = hold[j];
    }
}

#endregion Prime Base Greater Than 2

Order base factor?

    #region Order Base Factor?
    ///
    /// Checks for the order base factor
    /// for bases greater than one and
    /// reduces the result if factor exists.
    ///
    if (Arr1.GetUpperBound(0) > 0)
    {
        for (j = 0; j <= Arr1.GetUpperBound(0); j++)
            number = number * BigInteger.Power(new BigInteger(Convert.ToInt64(Arr1[j, 0])),
                     Convert.ToInt64(Arr1[j, 1]));

        if (division[0] % number != division[0].One())
            division[0] =
               division[0] / new BigInteger(Convert.ToInt64(Arr1[Arr1.GetUpperBound(0), 0]));
    }
    #endregion Order Base Factor?

    return division[0];
}

Exhibit 13 shows the time for E_n to 10,000 in tenths of a second by the number of bases.

Exhibit 14 shows the time for E_n to 3,000 decimal digits in tenths of a second by number of bases.

Exhibit 16 below is an example of the impact of register size on the calculation of Euler pseudo-primes.

Factors of Euler pseudoprimes are Germain primes/pseudoprimes and factors of Germain pseudoprimes are Germain primes/pseudoprimes. Exhibit 17 below shows this relationship. However, finding Germain pseudoprimes from large Euler numbers is still very time consuming and unlike Euler numbers, at this time there seems to be no way to calculate Germain numbers directly.

Problem: When linking to the Access database, one may get the following error: "The 'Microsoft.ACE.OLEDB.12.0' provider is not registered on the local machine." To correct this problem and/or for a limited explanation of the problem go to the following link: https://social.msdn.microsoft.com/Forums/en-US/1d5c04c7-157f-4955-a14b-41d912d50a64/how-to-fix-error-the-microsoftaceoledb120-provider-is-not-registered-on-the-local-machine?forum=vstsdb

Time32 Update Data

The "Time32 Update Data" routine updates the Access database from Excel to time the generation of Euler pseudoprimes.

Using the Code

This project was built as a concept prototype and as such does not have user friendly error checking of a more finished product. The C# code was built for a 64 bit system. Microsoft Office 2007 is a 32 bit environment requiring Excel-DNA to be compiled for 32 bits. If you want this project to run on a Microsoft Office version that is 64 bit, you will want to recompile Excel-DNA for 64 bits. The BigInterger Library has been modified for arbitrary size integers and routines from C# BigInteger Class have been used to expand the algorithms required to support this project. I have also modified the BigInteger Library with a shift divide for greater accuracy but with some loss of speed of calculation. There have been several other support routines in BigInteger Library that I have modified to expand their accuracy or to support this project. The code should be ready to go for experimenting in C# projects, Excel VBA and Excel spreadsheets. To use in Excel open BigInteger-packed.xll in a spreadsheet and there will be two function libraries, "BigInteger" and "Number Theory". A reference to BigInteger-packed.xll in the VBA environment is required for some of the macros to build the spreadsheets used in this project or to be used in other VBA applications. Microsoft Access and Excel require references to the following libraries:

1. Visual Basic For Applications
2. Microsoft Access 12.0 Object Library
3. OLE Automation
4. Microsoft Office 12.0 Access database engine Object
5. Microsoft ActiveX Data Objects 6.1 Library
6. BigIntExcel

If you are using .NET 4.0 or later you may want to alter this project code to use Microsoft’s Numeric library. From my research, the division algorithm is the key to accuracy and speed.

Points of Interest

This code is not thread safe because of the "RegisterSize" function which changes the values of static fields. By removing this function the process should be fairly straight forward in making this project thread safe. At present I have not explored the use of parallel processing in a multi-CPU environment. The maximum size of a BigInteger in Microsoft’s Excel 2007 is 32,767 decimal digits and will cause errors in calculations if this threshold is exceeded. See https://support.office.com/en-us/article/Excel-specifications-and-limits-16c69c74-3d6a-4aaf-ba35-e6eb276e8eaa.

There are two glaring mathematical omissions this article does not cover. The algorithm for finding the Order of a BigInteger is slow for finding an integers order with a large e. Second, as I stated earlier in this article, there is no direct way known (to me) to find Sophie Germain factors. However, there are techniques for finding Sophie Germain factors that are beyond the scope of this article and these techniques are relatively slow for finding large Sophie Germain factors. Note: Links to Jacques BASALDÚA paper are no longer working. I have not re-posted them. He will have to do that. Also, D:\BaseASeq\ProgramsC#\ExcelDna-0.32\Distribution\ExcelDnaPack.noexe should be changed to .exe before compiling the C# code. I have deleted all .exes in ExcelDNA so it will need to be recompiled or unzipped from the original zip file if you need them.

History

This is the first release of this software.