Numeric Algorithms

Basics

A numeric algorithm does some computation given one or more numeric values. For the purposes of measuring complexity, the size of a number is the number of bits (or digits) in the numbers, not the value of the numbers themselves!

Note the base of the numerals does not matter when computing asymptotic complexity. There is always a linear relationship between the number of digits required to represent a number in base $a$ and base $b$.

Wikipedia has a nice collection of links to articles on numeric algorithms.

Addition

The addition algorithm you already know works in any base $\geq 2$, because the sum of any three digits is at most a two-digit number (in any base!).

The algorithm goes from right to left, propagating carries. In binary:

 101011010
 100010111
1001110001

In decimal:

 939972
 654108
1594080

The complexity is $\Theta(n)$. We can’t do better because we have to touch each bit anyway.

Subtraction is similar; you can do it in linear time.

Multiplication and Division

You could multiply by repeated addition...

$$ \begin{eqnarray} 39053 \times 46987 & = & \underbrace{39053 + 39053 + 39053 + \ldots + 39053}_{46987\:terms} \\ \end{eqnarray} $$ ...but uggggh! That was 46986 additions of 5-digit integers!

You probably learned this algorithm before (shown here in base 2):

    1011
  × 1101
    1011
   0000
  1011
 1011   
10001111

The complexity of that algorithm is $\Theta(n^2)$, because multiplying an $n$-bit number by a one-bit number can be done in linear time, and we basically do that $n$ times. Then we do $n+1$ additions of $2n$-bit numbers: also quadratic.

Did you notice that all we did was just shifting and adding? And that shifting is really multiplying by 2? Look at the above example: it is saying $11 \times 13 = 11 + 44 + 88 = 143$. Some people actually carry out this binary-ish algorithm in decimal numbers like this:

 11    13
 22     6
 44     3
 88     1
143

Let’s generalize this to an actual algorithm. Here it is in Python, assuming $y \geq 0$:

def times(x, y):
    result = 0
    while y > 0:
        if y & 1:            # accumulate only when y is odd
            result += x
        y >>= 1              # halve y
        x <<= 1              # double x
    return result

There is also a quadratic time division algorithm. As long as $y \geq 1$, this Python function returns the quotient and the remainder:

def divide(x, y):
    if x == 0:
        return (0, 0)
    q, r = divide(x >> 1, y)
    q, r = (q << 1, r << 1)
    if x & 1:
        r += 1
    if r >= y:
        q, r = (q + 1, r - y)
    return (q, r)

Exponentiation

You can do exponentiation by repeated multiplication:

$$ \begin{eqnarray} 31^{400} & = & \underbrace{31 \times 31 \times 31 \times \ldots \times 31}_{400\:terms} \end{eqnarray} $$

Do we really have to do 399 multiplications? No! Let’s divide an conquer. First notice that, for example:

$$ x^{400} = x^{200}x^{200} $$

So we just have to compute $x^{200}$ once! Also note:

$$ x^{125} = x^{62}x^{62}x $$

In general:

$$ x^n = \begin{cases} 1 & n = 0\\ x^{\lfloor n/2\rfloor} \cdot x^{\lfloor n/2\rfloor} & n\;\mathrm{is\;even} \\ (x^{\lfloor n/2\rfloor} \cdot x^{\lfloor n/2\rfloor}) \cdot x & n\;\mathrm{is\;odd} \end{cases} $$

So as you recurse, halving the exponent, you do one multiplication when you get an even value, and two with an odd:

def pow(x, y):
    if y == 0:
        return 1
    h = pow(x, y // 2)
    result = h * h
    if y & 1:
        result *= x
    return result

This lets us do $31^{400}$ in 11 multiplcations, not 399. Divide and Conquer FTW:

$$ \begin{eqnarray} 31^{400} & = & 31^{200} \cdot 31^{200} \\ 31^{200} & = & 31^{100} \cdot 31^{100} \\ 31^{100} & = & 31^{50} \cdot 31^{50} \\ 31^{50} & = & 31^{25} \cdot 31^{25} \\ 31^{25} & = & 31^{12} \cdot 31^{12} \cdot 31 \\ 31^{12} & = & 31^6 \cdot 31^6 \\ 31^6 & = & 31^3 \cdot 31^3 \\ 31^3 & = & 31^1 \cdot 31^1 \cdot 31 \\ 31^1 & = & 31^0 \cdot 31^0 \cdot 31 \\ \end{eqnarray} $$

BUT WAIT! THIS IS NOT TAIL RECURSIVE! Even though you only recurse 7 levels deep, you are still stacking up partial results. This extra space can be significant! We can do better by using an iterative approach, based on the following observation:

$$ \begin{eqnarray} 31^{400} = 31^{16} + 31^{128} + 31^{256} \end{eqnarray} $$

I got that from the binary representation of 400 ($110010000_2$). Now walk through it, right to left; accumulate the result into $a$:

a = 1      #              a == 31^0
x = 31     # x == 31^1
x *= x     # x == 31^2
x *= x     # x == 31^4
x *= x     # x == 31^8
x *= x     # x == 31^16
a *= x     #              a == 31^16
x *= x     # x == 31^32
x *= x     # x == 31^64
x *= x     # x == 31^128
a *= x     #              a == 31^144
x *= x     # x == 31^256
a *= x     #              a == 31^400

Here it is in Python. Note the halving-and-doubling:

def pow(x, y):
    result = 1
    while y > 0:
        if y & 1:
            result *= x
        y >>= 1
        x *= x
    return result

Analysis Time! The number of multiplications is linear in the number of bits of the exponent (it can’t be over $2n$, right?) but what’s the complexity of the whole thing? You still have to compute (that is, at least write out), the whole answer, and how many bits are in that? First let’s find out what the answer is. Python says:

>>> 31**400
35049153521902777114708986131739468304466934398225961899541896904817914996246
23354936652173201312485109671587265639812785215391138217944131481131092898719
74685894570874283007085730451124921973944672872304632607933438457819117259406
11939966676143051741472986147114119144281388547152293994389062523413380947354
01786426523926432321776970921620521452184049782817872783269059536200968913247
54765133653979081932865951491386227358704384254668823602537295120143399193347
26030963677780786780733041458754894712578986705032469879343454522360163745820
4353222027948981902971920726384470581625600419245578432001

Okay now lets estimate the number of bits:

31400 ≅ (24.95)400 ≅ 21980

Practically a 2000-bit number (more or less). Hey, that’s pretty big, considering 31 is a 5-bit number, and 400 is a 9-bit number!

So writing out the result of an exponentiation is just too darn big! We’re usually a lot more interested in modular exponentation! This is coming up soon.

Integer Square Root

For a positive integer $n$, $\mathrm{isqrt}(n) =_{def}\lfloor \sqrt{n}\rfloor$.

The naïve approach is to “binary search it”. Let’s try to find the integer square root of 149.

• $149 \cdot 149 = 22201$, too high
• $74 \cdot 74 = 5476$, too high
• $37 \cdot 37 = 1369$, too high
• $18 \cdot 18 = 324$, too high
• $9 \cdot 9 = 81$, too low
• $13 \cdot 13 = 169$, too high
• $11 \cdot 11 = 121$, too low
• $12 \cdot 12 = 144$, low, but 13 was too high, so answer is 12.

We can do a little better by:

• Choosing a better starting value.
• Use a better way to choose the next guess, using the algorithm at Wikipedia.

Factorial

This was covered on the first day of class. Review if you need to.

Fibonacci

This was covered on the first day of class. Review if you need to.

Greatest Common Divisor

Euclid’s algorithm is, for $y \geq 0$, shown here in Python:

def gcd(x, y):
    return x if y == 0 else gcd(y, x % y)

Hey, did ya notice? It’s TAIL RECURISVE WOOO!

Examples

    gcd(253, 113)
    = gcd(113, 27)
    = gcd(27, 5)
    = gcd(5, 2)
    = gcd(2, 1)
    = gcd(1, 0)
    = 1

    gcd(2535, 1130)
    = gcd(1130, 275)
    = gcd(275, 30)
    = gcd(30, 5)
    = gcd(5, 0)
    = 5

    gcd(33031, 182845)
    = gcd(182845, 33031)
    = gcd(33031, 17690)
    = gcd(17690, 15341)
    = gcd(15341, 2349)
    = gcd(2349, 1247)
    = gcd(1247, 1102)
    = gcd(1102, 145)
    = gcd(145, 87)
    = gcd(87, 58)
    = gcd(58, 29)
    = gcd(29, 0)
    = 29

Least Common Multiple

It turns out that, for positive integers:

$$ \mathrm{lcm}(x,y) = \frac{xy}{\mathrm{gcd}(x,y)} $$

Since multiplying $x$ and $y$ might lead to a super large number, it is much more efficient to compute

$$ \mathrm{lcm}(x,y) = \left(\frac{x}{\mathrm{gcd}(x,y)}\right) y $$

In other words:

def gcd(x, y):
    return x if y == 0 else gcd(y, x % y)

def lcm(x, y):
    return (x / gcd(x, y)) * y

Modular Arithmetic

When numbers get too big, start counting over.

Let’s look at numbers mod 6:

-12	-11	-10	-9	-8	-7	-6	-5	-4	-3	-2	-1	-0	1	2	3	4	5	6	7	8	9	10	11	12	13	14
0	1	2	3	4	5	0	1	2	3	4	5	0	1	2	3	4	5	0	1	2	3	4	5	0	1	2

In the “mod 6 world”, both 2 and 10 map to 2. So 2 and 10 are in some sense equivalent in that world. We now have the notation “x ≡ y (mod N)”. More examples:

• 2 ≡ −10 (mod 6)      ...because 2 mod 6 = −10 mod 6
• 14 ≡ 8 (mod 6)
• 16 ≡ −2 (mod 6)

What do addition, multiplication, and exponentiation look like in a modular context?

Modular Addition

If the values to be added are in $0 ... N-1$, you don’t have to divide to get the remainder:

If adding a lot of numbers, do the modulos on the intermediate results:

  (30 + 26 + 8 + 29 + 28 + 5) mod 31
= (((((30 + 26) mod 31 + 8) mod 31 + 29) mod 31 + 28) mod 31 + 5) mod 31
= ((((25 + 8) mod 31 + 29) mod 31 + 28) mod 31 + 5) mod 31
= (((2 + 29) mod 31 + 28) mod 31 + 5) mod 31
= ((0 + 28) mod 31 + 5) mod 31
= (28 + 5) mod 31
= 2

Modular Multiplication

To compute $xy\:\mathrm{mod}\:N$, multiply (that’s quadratic) then divide to get the remainder (also quadratic). Again, do the modulos on intermediate results:

  (30 × 26 × 8 × 29 × 28 × 5) mod 31
  (((((30 × 26) mod 31 × 8) mod 31 × 29) mod 31 × 28) mod 31 × 5) mod 31
= ((((5 × 8) mod 31 × 29) mod 31 × 28) mod 31 × 5) mod 31
= (((9 × 29) mod 31 × 28) mod 31 × 5) mod 31
= ((13 × 28) mod 31 × 5) mod 31
= (23 × 5) mod 31
= 22

Modular Exponentiation

To compute $x^y\:\mathrm{mod}\:N$, we can use the exponentiation algorithm above, doing the modulos for all intermediate results.

Let’s try 31⁴⁰⁰ mod 100:

a = 1              #              a == 1
p = 31             # p == 31^1 == 31
(p *= p) % 100     # p == 31^2 == 61
(p *= p) % 100     # p == 31^4 == 21
(p *= p) % 100     # p == 31^8 == 41
(p *= p) % 100     # p == 31^16 == 81
(a *= p) % 100     #              a == 31^16 == 81
(p *= p) % 100     # p == 31^32 == 61
(p *= p) % 100     # p == 31^64 == 21
(p *= p) % 100     # p == 31^128 == 41
(a *= p) % 100     #              a == 31^144 == 21
(p *= p) % 100     # p == 31^256 == 81
(a *= p) % 100     #              a == 31^400 == 1

Also see the Wikipedia article; it includes a nice code fragment. Here it is in Python:

def modpow(x, y, n):
    result = 1
    while y > 0:
        if y & 1:
            result = (result * x) % n
        y >>= 1
        x = (x * x) % n
    return result

That algorithm is pretty cool, and you should know how to do it, but in practice there will be libraries.

Modular Exponentiation is built into many languages

Let’s compute $5438394857757488^{3424255654452323}\:\mathrm{mod}\:234235256666421$ in multiple languages!

Python

>>> pow(5438394857757488,3424255654452323,234235256666421)
189325407946433

Julia

julia> powermod(5438394857757488,3424255654452323,234235256666421)
189325407946433

Scala

scala> BigInt("5438394857757488").modPow(BigInt("3424255654452323"), BigInt("234235256666421"))
res0: scala.math.BigInt = 189325407946433

Ruby

>> require 'openssl'
=> true
>> 5438394857757488.to_bn.mod_exp(3424255654452323,234235256666421).to_i
=> 189325407946433

Clojure

user=> (.modPow (biginteger 5438394857757488) (biginteger 3424255654452323) (biginteger 234235256666421))
189325407946433

JavaScript

Not built-in; you have to install a big integer library. Fortunately that’s not too hard. You can get big-integer:

> var bigInt = require("big-integer");
> bigInt(5438394857757488).modPow(bigInt(3424255654452323), bigInt(234235256666421)).toString()
'189325407946433'

Java

Oh, Java doesn’t have a REPL? Booooooo.

Doing modular exponentiation in your head

The fact that you can do the modulos at any intermediate results lets you do some cool computations without a calculator:

2345 mod 31
= (25)69 mod 31
= 3269 mod 31
= (32 mod 31)69 mod 31
= 169 mod 31
= (1 mod 31)69 mod 31
= 1

Once you get good at this stuff you just "pull out" the mod 31 and write only this:

2345 ≡ (25)69 ≡ 3269 ≡ 169 ≡ 1 (mod 31)

Fermat’s Little Theorem for more fun!

Fermat’s Little Theorem is:

• If $p$ is prime, then for any integer $a$, $a\,^p\:\mathrm{mod}\:p = a$

A corollary is:

• If $p$ is prime and integer $a$ is not divisible by $p$, then $a\,^{p-1}\:\mathrm{mod}\:p = 1$

This means we can do even more modular exponentiation in our head.

  6^36 mod 31
  = (6^(30) * 6^6) mod 31
  = ((6^30 mod 31) * (6^6 mod 31)) mod 31
  = (1 * (6^6 mod 31)) mod 31
  = (6^2 * 6^2 * 6^2) mod 31
  = ((6^2 mod 31) * (6^2 mod 31) * (6^2 mod 31)) mod 31
  = (5 * 5 * 5) mod 31
  = 125 mod 31
  = 1

Discrete Logaarithm

Remember, in real numbers, when $x^y=z$ we say that $y$ is the logaritm of $z$ (base $x$)? Well in the modular integer world, whenever $b^k \equiv g\;(mod\:N)$, then we say $k$ is the discrete logarithm of $g$ in base $b$ modulo $N$.

Modular Exponentiation is probably One-Way

IMPORTANT We’ve just seen a polynomial time algorithm for computing $z = x^y\:\mathrm{mod}\:N$, but there currently exists no known efficient algorithm for finding $y$ given $x$, $z$, and $N$. Finding the discrete logarithm is in general computationally intractable. So modular exponentiation is a kind of one-way function.

Modular Multiplicative Inverse

The (regular) multiplicative inverse of a non-zero real number $x$ is the number $y$ such that $xy = 1$. We usually write it as $x^{-1}$ or $\frac{1}{x}$.

In the modular world the same idea holds.

But wait.... You don’t always have a modular inverse. There is no number $y$ such that $12y\:\mathrm{mod}\:34 = 1$. There is a theorem that says $x$ has an inverse mod $N$ if and only if $x$ is relatively prime to $N$. ($x$ and $y$ are relatively prime means that gcd($x$, $y$) = 1.)

Code to find the modular inverse is built into many languages, or at least easy to find in a library.

Julia

julia> invmod(234535235234,2345665654331)
146170270779

Scala

scala> BigInt("234535235234").modInverse(BigInt("2345665654331"))
res0: scala.math.BigInt = 146170270779

Clojure

user=> (.modInverse (biginteger 234535235234) (biginteger 2345665654331))
146170270779

To compute the modular inverse, yourself, from scratch (for educational purposes), you can use the Extended Euclidean algorithm. Details are in any good book on number theory, algorithms, or cryptography.

Primality Testing

There are lots of algorithms out there to determine if a number is prime. Two main classes of algorithms:

• Deterministic: The algorithm always produces a correct answer. It answers yes if and only if the candidate is prime and answers no if and only if the candidate is composite.
• Probabilistic: The algorithm answers yes for all primes, but also for some composites. If it answers no, you definitely have a composite. If it answers yes for an input $n$, then $n$ is called a probable prime. If $n$ is a composite for which the algorithm answers yes, then $n$ is called a pseudoprime.

Even the fastest deterministic algorithms are really, really, slow compared to the probabilistic ones. Besides, you can get probabilistic algorithms to have a probability of error less than that of hardware failures or getting hit by lightning and winning the lottery on the same day.

Some Algorithms

• Deterministic
  • Trial division: exponential bit complexity, but fine for small numbers!
  • APR (Adelman-Pomerance-Rumley, 1979): Runs in $O(n^{c \log \log n})$ time for some $c>0$.
  • AKS (Agrawal-Kayal-Saxena, 2002): The first published general, polynomial-time, deterministic, and unconditional test. Runs in bitwise $O(n^{12+\epsilon})$ time.
  • Lenstra-Pomerance (2003): Improvement on AKS, sped up to $O(n^{6+\epsilon})$ time.
  • ECPP (a.k.a Elliptic Curve Primality Proving): A family of algorithms for determining primality. Fast and widely used. Las Vegas algorithm, so there are some nasty inputs that take too long.
• Probabilistic
  • Fermat: Simple, but fails spectacularly on Carmichael numbers, so do not use
  • Solovay-Strassen: Fairly simple, but replaced by Miller-Rabin.
  • Lehmann: Also simple, but Miller-Rabin is more popular.
  • Miller-Rabin: Seems to be most popular of the probabilisitic methods.
  • Baillie-PSW (a.k.a. Baille-Pomerance-Selfridge-Wagstaff): Based on the idea that Fermat pseudoprimes and Lucas pseudoprimes don’t overlap. For this algorithm, there are no pseudoprimes below $2^{64}$.

Actually we care about a few other attributes of these algorithms:

Algorithm	General?	Deterministic?	Polynomial?	Unconditional?
Trial Division	Yes	Yes	No, not at all ☠	Yes
Miller	I think so	Maybe	Yes	No
Miller-Rabin	Yes	No	Yes, $O(n^3)$, or $O(n^{2+\epsilon})$ with FFT	Yes
Baillie-PSW	Yes	No	Yes	Yes
ECPP	Yes	Yes	Expected, $O(n^{5+\epsilon})$ for most inputs	TODO
Fast ECPP	Yes	Yes	Expected, $O(n^{4+\epsilon})$ for most inputs	TODO
AKS	Yes	Yes	Yes, $O(n^{12+\epsilon})$	Yes
Lenstra-Pomerance	Yes	Yes	Yes, $O(n^{6+\epsilon})$	Yes

Solovay-Strassen Test

It turns out that, for an odd positive $n$, if $n$ is prime then

a(n-1)/2 ≡ Jacobi(a,n) (mod n)

is true for all $a \in \{1 ... n-1\}$, and if $n$ is composite then that congruence holds true for at most half of the values in $\{1 ... n-1\}$.

So if you pick a random $a$ less than $n$, and the congruence fails, you definitely know $n$ is composite. But if it holds, $n$ has a 50% chance (1 out of 2) of being composite. So pick another value for $a$. If it holds again, you have a 1/4 chance of being composite. If another one passes, your odds of being composite go down to 1 in 8. Pass 30 straight tests and you have only a 1 in a billion chance.

Lehmann Test

This one doesn’t require computing the Jacobi symbol. You simply compute $a^{(p-1)/2}\:\mathrm{mod}\:p$ for $k$ randomly selected $a$’s (all less than $p$). If you ever get a value other than $1$ or $p-1$, then $p$ is definitely composite. Otherwise, if you always get $1$ and $p-1$ and there is at least one $p-1$ then the probability that $p$ is prime is $1-2^{-k}$.

Miller-Rabin Test

Read the Wikipedia article.

Also some cool implementations on Rosetta Code.

Integer Factorization

No one as yet knows any efficient general-purpose algorithm for determining the prime factors of an integer, but a number of special-purpose algorithms (besides the simple-minded "direct search" approach) have been invented.

MathWorld has good articles on Prime Factorization and Prime Factorization Algorithms. Start your self-study tour there.

Julia has a factor function built-in:

julia> factor(227275645766597)
Dict{Int64,Int64} with 3 entries:
  23  => 3
  211 => 1
  97  => 4

julia> factor(2398472839472983423499238211)
Dict{Int128,Int64} with 5 entries:
  7                 => 1
  911063            => 1
  83                => 1
  113               => 1
  40098840940603649 => 1

julia> factor(2398472839472983423499238212398472304972342300923841)
Dict{BigInt,Int64} with 5 entries:
  125573593                             => 1
  4327                                  => 1
  3                                     => 2
  367                                   => 1
  1336413995720810139514854155187173177 => 1

Why care? The fact we can generate primes efficiently but not (yet???) factor integers efficiently gives us one way to implement public-key cryptography.

Oh yeah, you might want to read about Shor’s Algorithm.

Floating-Point Algorithms

These are covered in a different class (CMSI 284).