Modular exponentiation using only addition and subtraction

Question

The challenge

In the least number of source code characters, in a language of your choosing, devise a function which raises base a to power b, modulo m, using purely integer arithmetic, where a, b, and m are all unsigned 64-bit values or larger. You may use signed arithmetic if your language does not have unsigned arithemetic.

This is the standard number-theory operation a^b (mod m), and not a^(b mod m). Conceptually, and perhaps literally, depending on your implementation, all of your intermediate results should conform to the modulus m. That is to say, you won't actually be raising a to the b power and then computing the remainder after dividing by m — that would require as much as two trillion gigabytes of RAM to store a number that big. Instead, you'll be raising a to the b power while staying within the constraints of modulus m, or straying from it only a bit (no pun intended).

This is actually a lot easier than it sounds!—Raising numbers to astronomically large powers can actually be a surprisingly fast operation. If you do it right, you can compute the result in a handful of microseconds, or perhaps few milliseconds if you're using a slow language. However, you are free to use any algorithm you like, even if it uses an increment-counter in a hard loop running for quadrillions of years. The only requirement with regard to execution time is that your program provably halt in a finite amount of time, providing the correct answer.

Input assumptions

You can assume that a ≥ 0, b ≥ 0, m > 0, and a < m, but you should not assume that b < m. That is, you should assume that both b and m can be any value as large as your integer type supports (this will be 2⁶³–1 for signed 64-bit integers and 2⁶⁴–1 for unsigned 64-bit integers).

If you are given a = b = 0, you can return whatever value you like, since 0^0 is undefined. Similarly, if you are given m = 0, you can return whatever you like, since 0 is not a valid modulus.

Test values

Your program should work correctly for all inputs, but here are a few samples to help you with validation:

              Base             Exponent               Modulus                Result
                 2                    8                   255                     1
                 2                    8                   256                     0
                 2                    8                   257                   256
                 3                    7                 10000                  2187
                 2                 2046                  2047                     1
               123                  456                   789                   699
                 3                 1000  18446744073709551615  12311760789144243126
             86400             22157322        48519018006822        40149207423504
           8675309         100018327824   8621993634251008000   3858055581225668161
325284989554104320  1508436685178379520   8582294829391072256   6354230931838838784

Counting source characters

Spaces and tabs each count as 1 character. Newlines do not count, unless the newline character adds syntactic whitespace in order to permit parsing. For example, the newline at the end of the first line here counts as 1 character because it is required to permit proper parsing:

#define foo 7
int f(){...}

But the newline character at the end of the first line here does not count:

int f(){...}
int g(){...}

because parsing would still be possible if the two lines were adjoined:

int f(){...}int g(){...}

This allows you to present your source code with line-breaks at natural locations without penalty, so that you don't have to try to cram everything onto a single line. Finally, if you're using C/C++ and you put a \ at the end of a line, then the newline also does not count because you're using the \ only for readability of presentation.

Restrictions & Allowances:

1. Arithmetic: You may use only integer addition, subtraction, logical operations, and conditionals. Comparisons count as a form of subtraction. Multiplication and division are not allowed. (However, the specific case of left- and right-shifting may be written using *2, *=2, /2, /=2 if these are syntactically shorter than the bit-shifting operators in your language, e.g., <<1, <<=1, >>1, >>=1, etc.) Use of a built-in exponentiation operator (e.g., **) is not allowed.
2. Floating-point arithmetic is not allowed. However, you may use floating-point variables provided you only use them to hold integers. (Floating-point operations aren't really going to help you on this challenge anyway, so don't be discouraged by this restriction.)
3. Your function is not required to have a name (it can be an anonymous or lamba function if your language supports that), but it does have to be a bona fide function that you can call an arbitrary number of times.
4. No global variables. (Exception: If there is no way in your language to provide a solution without global variables, then you may use global variables. But only if there is no other way.)
5. You may ignore any and all compiler, interpreter, or runtime warnings, so long as your program computes correct results.
6. If you are using a language which requires importation of type declarations in order to guarantee a 64-bit integer size (e.g, <stdint.h> in C), you may assume that this is already included. Note that uint64_t is considerably fewer characters than unsigned long long.

Cheating:

1. You may optionally provide a cheat solution, but only if you are also providing a legal solution and your cheat solution computes correct results as well (please post the cheat version below the legal version and do not count the cheat version as your score).
2. If you do provide a cheat version that uses a built-in exponentiation operator or function call (e.g., ** in AWK or pow() in C), beware of rounding and truncation issues. Just because you're providing a cheat version doesn't mean you're allowed to produce incorrect results. :-) For example, writing uint64_t modpow(uint64_t a,uint64_t b,uint64_t m){return pow(a,b)%m;} in C hasn't even the remotest chance of being correct for anything but the smallest values, so don't do that. (Actually, it might even be a challenge to use built-in non-modular exponentiation. I'd be interested to see if anyone actually manages to pull that off. But even though I'm posing that as a sub-challenge, it's still considered a cheat for the purposes of the larger challenge.)

If you need help getting started:

Here are a couple resources with simple algorithms for computing modular exponentiation:

1. Modular exponentiation [Wikipedia.org]
2. Fastest modular exponentiation [CodeGolf.StackExchange.com]

Lastly:

I will be posting a solution in C99 that I worked out earlier today as a byproduct of writing a Miller–Rabin primality test for another purpose. My answer should not be treated as an eligible entry to the challenge.

I am really looking forward to seeing what you come up with in Golfscript, Befunge, Haskell, Clojure, and more!

Answer 1

Python 2 (70)

def f(a,b,m,n=1):
 exec("n=0"+"+n"*a+";")*b
 while n>=m:n-=m
 return n

Very, very slow. The multiplications are string operations, not arithmetic ones, so I hope they're allowed.

The line with exec generates and executes code like (with a=4, b=3):

n=0+n+n+n+n;n=0+n+n+n+n;n=0+n+n+n+n;

Edit: Here's two alternative solutions that use list/string multiplication in a cheaper way by representing the current number as a list of ones. They're both exactly the same length as the original (70 chars).

def f(a,b,m,s=[1]):
 exec"s*=a;"*b;n=len(s)
 while n>=m:n-=m
 return n

def f(a,b,m,s=[1]):
 exec"s*=a;"*b
 while s[m:]:s=s[m:]
 return len(s)

More edit: A 66-char solution that uses string splitting to do the modulo m

def f(a,b,m,s="1"):
 exec"s*=a;"*b
 return len(s.split("1"*m)[-1])

and a 65-char recursive one-liner

f=lambda a,b,m:b and len((a*f(a,b-1,m)*"1").split("1"*m)[-1])or 1

Answer 2

Haskell (78)(67)(61)

p a b m|b==0=1|0<1=until(<m)(+(-m))$sum$replicate a$p a(b-1)m

The function f is a simple modulus function implemented using repeated subtraction.

The modular exponentiation function p first tries to match b==0. If this is true, then it returns 1. Otherwise, it calls p recursively, "multiplies" it with a, and reduces it mod m.

It work with all the inputs until p 86400 22157322 48519018006822, where the stack overflows.

EDIT: Thank's Matt! Using Matt's idea, I was able to shorten the mod function and have it in place. Still overflows though.

EDIT: Realized foldl(+)0 is just a sum. Doh!

Answer 3

Haskell: 87 bytes

Performs modular exponentiation x^y mod n in a number of operations which is logarithmic in the exponent y. The code

(x!y)m|odd y=z q|y>0=q|0<1=1where z 0=0;z n=until(<m)(+(-m))$x+z(n-1);q=(z x!div y 2)m

defines an operator (!) to do modular exponentiation, so you can use it like

(2!3) 5 == 3

The program works by using the two relations x^{2 n} = (x²)ⁿ and x^{2 n + 1} = x (x²)ⁿ to accumulate the result, halving the exponent at each step. Intermediate results are taken mod m by repeated subtraction. Here's an ungolfed version:

(x ^% y) m
    | y == 0     =  1
    | odd y      =  x `times` (squaredX ^% halfY) m
    | otherwise  =  (squaredX ^% halfY) m
    where squaredX = x `times` x
          halfY = y `div` 2
          a `times` 0 = 0
          a `times` b = reduce $ a + a `times` (b - 1)
          reduce n = if n < m then n else reduce (n - m)

Note: div is used to carry out a right shift operation only.

Answer 4

Javascript (92) 89

function f(a,b,m){
    c=1;
    while(b--){
        o=0;
        for(i=0;i<a;i++){
            o+=c-m;
            if(o<0)
                o+=m;
        }
        c=o;
    }
    return c;
}

Thanks to Todd for the while loop trick. I didn't know Javascript could do that.

Each iteration through the outer loop "multiplies" the result by a. Essentially it adds the previous iteration's result to itself a times (i.e. multiplying itself by a). It also subtracts m each time, and if it's less than 0, it adds it back.

This function is not very efficient, especially for large inputs, but it should eventually return the correct response.

Answer 5

C, 175 characters

Answering my own question. This answer should not be considered an eligible entry to the challenge.

#define N uint64_t
#define X(F,G,d)\
N F(N m,N a,N b){N c=d;for(;b;b/=2){if(b&1)c=G(m,c,a);a=G(m,a,a);}return c;}
N A(N m,N a,N b){return a+=b,a-(a<m&&a>=b?0:m);}X(M,A,0)X(E,M,1)

How it works: The macro invocation X(E,M,1) generates a function E(m,a,b), and the invocation X(M,A,0) generates a function M(m,a,b) using the same macro template with slightly different parameters. E(m,a,b) is the modular Exponentiation function. It calls the modular Multiplication function M(m,a,b), which in turn calls the modular Addition function A(m,a,b). The modular addition function has subtle logic to detect and handle both 64-bit overflow and modular overflow correctly, which is transparently leveraged by the higher-level functions.

Performance: Worst-case performance of this is approximately (log₂m)² steps, or about 4096 steps for large 64-bit m.

Runtime on a 3.4 GHz Intel Core i7 to compute 1 million different values of

a ^ 18446744073709551615 (mod 18446744073709551615)

for a ∈ [1,10⁶] is 51.8 seconds, or 51.8 microseconds per exponentiation invocation, where 18446744073709551615 = 2⁶⁴ – 1, which is the most taxing 64-bit value for this method.

At the other end of the spectrum, using much smaller values, the runtime to compute 1 million different values of

a ^ a (mod 2a + 1)

for a ∈ [1,10⁶] is only 3.41 seconds, or 3.41 microseconds per exponentiation operation.

---

C, 110 characters

Here is formulation using nested loops:

#define N uint64_t
N E(N M,N A,N B){N C=1;while(B--){N c=0,a=A;while(a--)c+=C,c>=M||c<C?c-=M:0;C=c;}return C;}

How it works: The capital variables A, B, and M are the original values passed in; the lowercase variables c and a are temporary work variable copies of their uppercase counterparts. The inner loop performs repeated addition to produce multiplication. The outer loop performs repeated multiplication (using the inner loop) to produce exponentiation.

Performance: This formulation is reasonably fast for very small input values, but performs extremely poorly for even moderately large values. The runtime is proportional to the product of a and b and is independent of m; in other words, this is an O(ab) algorithm. Worst-case performance on 64-bit input is 2¹²⁸ steps, which on my system would require 4.2 × 10²⁹ seconds or about 13 thousand million million million years. For worst-case input, this is 10³⁴ times slower than the more efficient version earlier.