Antimatroid, The

thoughts on computer science, electronics, mathematics

Archive for the ‘Number Theory’ Category

Integer Factorization by Dynamic Programming with Number Theoretic Applications

with 7 comments

Having been a participant of a number of mathematical programming competitions over the years, I’ve had to find a number of efficient ways of implementing many common Number Theoretic Functions. In this write up, I’m going to go over a method I’ve found useful for easily factoring numbers using a sieving method, go over some the implementation of a few number theory functions along with time complexity analysis of each. The cornerstone to many of these implementations relies on the ability to quickly factor integers and find primes.

One of the most popular methods of for finding primes is the Sieve of Eratosthenes. The algorithm starts by populating a table with every positive integer from 2 to a ceiling value. Then find the first integer not yet crossed off, in the case 2, and eliminate every multiple of 2 from the table then return to 2 and find the next positive integer not yet crossed off and repeat the procedure until the end of the table is reached. The method is fine and all, but a lot of really great information is lost in that computation. Here is a sample implementation:

bool[] isPrime = new bool[400];
for (uint n = 2; n < isPrime.Length; n++)
    isPrime[n] = true;
for (uint n = 2; n < isPrime.Length; n++)
    if (isPrime[n])
        for (uint m = 2, c = 0; (c = m * n) < isPrime.Length; m++)
            isPrime[c] = false;

On the other hand, say we approach sieve a little differently. Create an empty table as large as the ceiling value. Start at 2 and for every multiple of 2, create a record that has two parts: 2 and half of the multiple. Return to 2 and find the next integer in the table that has yet to be recorded, in this case 3. For every multiple of 3, create a record that has two parts: 3 and third of the multiple (only if the multiple was not previously recorded. E.g., 6 because 2 previously recorded the record). Return to 3 and find the next integer in the table that has yet to be recorded so on and so forth until every integer in the table has been recorded.

The following graphic demonstrates this process for a ceiling values of 25. If we wish to factor 16, we go to 16’s record (2, 8), follow to 8’s record (2,4), again follow 4’s record (2, 2) and finally 2’s record (2, λ). Thus the prime factorization of 16 is 2, 2, 2, 2.

factor_table_polar

It should be apparent that this algorithm is a simple dynamic programming solution that yields two major results:

  1. We have factored every positive integer up to a ceiling value.
  2. We have found every positive prime integer up to a ceiling value.

And one major draw back

  1. Uses a lot of memory as a trade off for speed.

Let’s get into the C# implementation. To start off, we need a record class that’ll store the information about the first prime that divides an entry and the composite to jump to if the record corresponds to a composite.

public class Record {
    public uint Prime {get; set;}
    public uint? JumpTo {get; set;}
}

We’ll have a class called NumberTheory and assume that it has the following structure. If you want, you could make this a Singleton class but I felt it was unnecessary for the scope of this write up.

public class NumberTheory {
    private Record[] table;

    ...
}

It makes sense to put the core algorithm in the constructor and then have member methods for each of the functions we’d like to have. It should be assumed that for the lifetime of the class that the largest value ever called on the methods will be N otherwise an exception should be thrown by the methods (omitted here for brevity).

public NumberTheory(uint N){
    uint c = 0;
    table = new Record[N+1];
    for(uint n = 2; n < table.Length; n++) {
        if(table[n] != null)
            continue;
        table[n] = new Record() { Prime = n };
        for(uint m = 2; (c = n * m) < table.Length; m++)
            if(table[c] == null)
                table[c] = new Record() { JumpTo = m, Prime = n };
    }
}

The time complexity of the implementation can be derived using some analysis and by having some knowledge of certain identities. Starting from 2 there are N – 2 numbers to check of which 1/2 will be visited by the interior loop, starting from 3 there are N – 3 numbers to check of which 1/3 will be visited visited by the interior loop, so on and so forth leading to the following summation:
\displaystyle T(n) = \sum_{p \le n} \left( \frac{1}{p}(n - p) \right)
\displaystyle T(n) = \sum_{p \le n} \left( \frac{n}{p} - 1 \right)

If we separate the summation into the harmonic series of primes (HSP) and prime counting function (aggregate 1 for every prime less than n) we get:
\displaystyle T(n) = n \sum_{p \le n} \left( \frac{1}{p} \right) - \pi(n)

Asymptotically, the HSP tends towards \ln \ln n and \pi(n) towards \frac{n}{\ln n}.
\displaystyle T(n) = n \left( \ln \ln n \right) - \pi(n)
\displaystyle T(n) = n \left( \ln \ln n \right) - \frac{n}{\ln n}
\displaystyle T(n) = n \frac{ \ln n \cdot \ln \ln(n) - 1}{\ln n}

Giving us our final asymptotic time complexity of O(n \ln \ln n).

The first member method we’ll implement is a trivial check to see if a given number is prime by checking the table’s Record’s JumpTo property for null.

public bool IsPrime(uint N){
    return !table[N].JumpTo.HasValue;
}

The time complexity here is a simple O(1).

From IsPrime, we can easily implement a function that will get every single prime up to a given value by iterating over the table.

public void PrimesLessThan(uint value, Action<uint> actOnPrime) {
    for(int n = 2; n < value; n++)
	if(IsPrime(n))
		actOnPrime(n);
}

We can get the prime counting function \pi(n) but utilizing the PrimeLessThan function.

public uint CountPrimes(uint n) {
	uint count = 0;
	PrimesLessThan(n, (p) => {count++;});
	return count;
}

Here the time complexity is the same as PrimesLessThan: O(n).

We can get the prime factorization of a composite easily. To do so, we simply iterate over the records until we reach a Record with no JumpTo value.

public void PrimeFactorsOf(uint composite, Action actOnFactor) {
    Record temp = table[composite];
    while(temp != null) {
        actOnFactor(temp.Prime);
        if(temp.JumpTo.HasValue)
            temp = table[temp.JumpTo.Value];
        else
            temp = null;
    }
}

The time complexity of this implementation relies on the the prime omega function \Omega(n) which is the number of prime factors (not necessarily distinct) of n. The function tends to O( \ln \ln n).

From PrimeFactorsOf, we can also easily implement Euler’s Totient Function \phi(n)– the function tells us how many positive integers less than n are coprime to n. It is defined as:
\displaystyle \phi(n) = n \prod_{p|n} \left( 1 - \frac{1}{p} \right)
Which essentially states that if you multiply all of the repeated prime factors of n together by all of the non-repeat prime factors – 1 of n together, you will have the result of \phi(n).

public uint EulerTotient(uint n){
    uint phi = 1, last = 0;
    PrimeFactorsOf(n, (p) => {
        if(p != last) {
            phi *= p - 1;
            last = p;
        } else {
            phi *= p;
        }
    });
    return phi;
}

A similar function known as the Dedekind \psi(n) Function defined as
\displaystyle \psi(n) = n \prod_{p|n} \left ( 1 + \frac{1}{p} \right )
can be implemented in a similar way as \phi(n):

public uint DedekindPsi(uint n){
    uint phi = 1, last = 0;
    PrimeFactorsOf(n, (p) => {
        if(p != last) {
            phi *= p + 1;
            last = p;
        } else {
            phi *= p;
        }
    });
    return phi;
}

The Von Mangoldt Function \Lambda(n) is another interesting function, unfortunately I haven’t had a chance to use it, but it is trivial to implement so I will include it here for completeness. It is defined as
\displaystyle \Lambda (n) = \begin{cases} \ln{p} & \text {if n = prime to some positive integer power} \\ 0 & \text{otherwise} \end{cases}

public double VonMangoldt(uint n){
    uint P = 0;
    PrimeFactorsOf(n, (p) => {
        if(P == 0) {
            P = p;
        } else if (P != p) {
            P = 1;
        }
    });
    return Math.Log(P);
}

The Möbius Function \mu(n) is a handy little function for determining if a number is square free or not (among many more interesting things). It is defined as:
\displaystyle \mu(n) = \begin{cases} 0 & \text{if n has one or more replicated factors} \\ 1 & \text{if n = 1} \\ (-1)^{k} & \text{if n is a product of k distinct primes} \end{cases}

public int MoebiusFunction(uint N){
    if(N == 1)
        return 1;
    bool distinct = true;
    uint last = 0, k = 0;
    PrimeFactorsOf(N, (p) => {
        if(p == last) {
            distinct = false;
        } else {
            k++;
            last = p;
        }
    });

    if(distinct) 
        return ((k & 1) == 0) ? 1 : -1; 
    return 0;
}

Since EulerTotient, DedekindPsi, VonMangoldt and Moebius each use PrimeFactorsOf without any additional lifting, their time complexities are the same as PrimeFactorsOf – O(\ln \ln n).

The last function I’ll implement is the Mertens Function M(n) which is simply defined as \displaystyle M(n) = \sum_{k = 1}^{n} \mu(k):

public int Mertens(uint n){
    int m = 0;
    for(uint k = 1; k <= n; k++)
        m += MoebiusFunction(k);
    return m;
}

The implementation’s time complexity is simply O(n \ln \ln n).

Advertisements

Written by lewellen

2009-04-01 at 12:00 am

The ring of Gaussian Integers

leave a comment »

I recently started taking a foray into Complex Analysis as a means of filling in the gaps of my undergraduate mathematics knowledge. After reading about holomorphic functions, the Cauchy-Riemann equations and how to model ideal fluid flows I decided to take a break to digest it all. During this period I was reflecting on what I had done with Number Theory the previous summer and the question popped in my head: what is the complex equivalent of w \equiv v \pmod{z}? Or for that matter, are integer primes also primes in the complex domain (and vice versa)? And, are all complex z = \displaystyle\prod_{k = 0}^{\infty} p_{k}^{e_{k}} factorizations unique?

After looking around the internet for a while I came across a rather well written paper [pdf] by Assistant Professor Keith Conrad of the University of Connecticut that answered all of the questions brewing in my head along with the ones that weren’t. I’m going to surmise the basic ideas behind the solutions Conrad wrote about, if you have the free time it’s worthwhile to read the original text in its entirety.

If we consider the set of complex numbers of the form n + mi and restrict n, m \in \mathbb{Z} we have the set of the Gaussian Integers \mathbb{Z}[i]. If we take two values z, w \in \mathbb{Z}[i] we can say that z|w only if w = zu where u \in \mathbb{Z}[i]. Thus, \displaystyle\bar{z}w = \bar{z}zu \Rightarrow u = \frac{\bar{z}w}{\bar{z}z}. With this last step we’ve reduced the problem to satisfying two conditions: \displaystyle\bar{z}z|\Re{\bar{z}w} \wedge \bar{z}z|\Im{\bar{z}w}. For example let z = 1 + 2i, w = -5 + 10i \Rightarrow \bar{z}z = 5, \bar{z}w = 15 + 20i which satisfies the conditions 5|15 \wedge 5|20 \therefore z|w. For the congruence w \equiv v \pmod{z} with w,v,z \in \mathbb{Z}[i] to hold true, it must be the case that w = zu + v \Rightarrow z|(w - v).

Under \mathbb{Z} any integer m that is divisible by any integer other than a unit \pm 1 or unit multiple \pm m is said to be composite, otherwise m is said to be prime. The same definition holds in \mathbb{Z}[i] but our units are now \pm 1, \pm i and our unit multiples are now \pm z, \pm zi. As one might have guessed all primes in \mathbb{Z} are primes in \mathbb{Z}[i] however, the converse does not hold. For example 5 is prime but can be written as (1 + 2i)(1 - 2i) and is thus a composite number under \mathbb{Z}[i]. It is also useful to mention that if \bar{z}z \in \mathbb{P} then z is prime. By way of the Fundamental Theorem of Arithmetic, every composite integer can be written as the product of prime integers, this is also the case under \mathbb{Z}[i] and each factorization is unique.

The mechanics of finding the prime factorization of m \in \mathbb{Z}^{+} can be na├»vely done by trial division from n = 2 to \sqrt{m}. Once n|m reset n = 2 and m_{k} = m_{k - 1} / n and continue this procedure until m_{k} = 1. That is all well and fine for \mathbb{Z}^{+} but it isn’t immediately useful for doing the same under \mathbb{Z}[i].

For example, let’s consider z = -1395 - 12410i, since we are uncertain about how to factor z, lets think about how to factor \bar{z}z as we already know how to factor under \mathbb{Z} and see where that takes us. \bar{z}z = 155954125 = 5^{3} \cdot 61 \cdot 113 \cdot 181 We went from \mathbb{Z}[i] \to \mathbb{Z} so we’d hope that we could go the other way around by finding p = (n + mi)(n - mi) for each of the prime factors of \bar{z}z. In other words, we want to know p = n^{2} + m^{2} which will produce four Gaussian primes n \pm mi, m \pm ni. Thus, we find 5 = 1 + 2^{2}, 61 = 5^{2} + 6^{2}, 113 = 7^{2} + 8^{2}, 181 = 9^{2} + 10^{2}. Now, according to Conrad we should be able to pick any one of the four possible Guassian primes for each of the prime factors and multiply each Gaussian prime across, in doing so we get z = (2 + i)^{2} (2 - i) (6-5i) (8 - 7i) (10 - 9i). Now that’s pretty damn cool.

The general procedure for factoring a Gaussian Integer z into it’s prime components is to find the integer factorization of \bar{z}z and for each prime factor find p = n^{2} + m^{2}. Next, explore the Cartesian space formed by each factor’s four Gaussian primes until you come across a product coordinate that equals z and output said coordinate. As one can image this is a woefully poor algorithm for find the prime factorization of z. If anyone has ideas or knows of more efficient methods post them in the comments.

There is clearly much, much more to be said about Gaussian Integers, but this feels like a good stopping point. If you want to find out more about how some of the traditional Number Theory constructs are defined under \mathbb{Z}[i] you’ll want to read the entirety of Conrad’s paper or jump on Google and see what’s out there.

Written by lewellen

2008-08-10 at 8:05 pm

Posted in Number Theory

Tagged with