3

I'd like to verify whether

pow(a, b) % b == a

is true in C, with 2 ≤ b ≤ 32768 (215) and 2 ≤ a ≤ b with a and b being integers.

However, directly computing pow(a, b) % b with b being a large number, this will quickly cause C to overflow. What would be a trick/efficient way of verifying whether this condition holds?

This question is based on finding a witness for Fermat's little theorem, which states that if this condition is false, b is not prime.

Also, I am also limited in the time it may take, it can't be too slow (near or over 2 seconds). The biggest Carmichael number, a number b that's not prime but also doesn't satisfy pow(a, b)% b == a with 2 <= a <= b (with b <= 32768) is 29341. Thus the method for checking pow(a, b) % b == a with 2 <= a <= 29341 shouldn't be too slow.

Isaiah
  • 1,852
  • 4
  • 23
  • 48
  • 3
    By `a ^ b` you mean `a` raised to the power `b` or `bitwise XOR` of `a` and `b`? – Ajay Brahmakshatriya Sep 20 '17 at 08:14
  • @AjayBrahmakshatriya i thinks he means power because else it is much harder to get an overflow. However `^` in C is normally XOR, so pretty hard to say for sure. – izlin Sep 20 '17 at 08:18
  • @izlin I figured he meant power but I was asking just to clarify. – Ajay Brahmakshatriya Sep 20 '17 at 08:19
  • Look at the iterative square-and-multiply algorithm for a general way of computing the exponentiation modulo n. (https://en.wikipedia.org/wiki/Exponentiation_by_squaring). This can be even more efficient by using the k-ary- the improved k-ary or the sliding window method. Multiplications can be computed more efficiently when using the Karatsuba-Multiplication. Karatsuba + Sliding Window yields the fastest algorithm. – Maximilian Gerhardt Sep 20 '17 at 08:31
  • @AjayBrahmakshatriya, I did, thanks! Fixed it now! – Isaiah Sep 20 '17 at 08:37
  • @MaximilianGerhardt sorry I keep messing it up :] – Isaiah Sep 20 '17 at 08:40
  • @MaximilianGerhardt you don't need to do the exponentiation as that would need a lot of time and memory for large powers – phuclv Sep 20 '17 at 08:43
  • Possible duplicate of [Calculate (a^b)%c where 0<=a,b,c<=10^18](https://stackoverflow.com/questions/32485750/calculate-abc-where-0-a-b-c-1018) – phuclv Sep 20 '17 at 08:43
  • 2
    Post should clearly state `a,b` are integers. – chux - Reinstate Monica Sep 20 '17 at 11:00
  • the `pow()` function expects `double` for parameter types. Not all integer values can be exactly represented via a `double` so some integer values will not result in the same value after being run through this code. – user3629249 Sep 22 '17 at 01:11
  • Possible duplicate of [Calculating (a^b)%MOD](https://stackoverflow.com/questions/11272437/calculating-abmod) – phuclv May 10 '19 at 14:51
  • [Raising large number to large power and mod it by a large number?](https://stackoverflow.com/q/27153665/995714), [Calculate (a^b)%c where 0<=a,b,c<=10^18](https://stackoverflow.com/q/32485750/995714), [Calculating (a^b)%MOD](https://stackoverflow.com/q/11272437/995714), [a to power b modulus k](https://stackoverflow.com/q/30138020/995714)... – phuclv May 10 '19 at 14:52

2 Answers2

5

You can use the Exponentiation by squaring method.

The idea is the following:

  • Decompose b in binary form and decompose the product
  • Notice that we always use %b which is below 32768, so the result will always fit in a 32 bit number.

So the C code is:

/*
 * this function computes (num ** pow) % mod
 */
int pow_mod(int num, int pow, int mod)
{
    int res = 1

    while (pow>0)
    {
        if (pow & 1)
        {
            res = (res*num) % mod;
        }
        pow /= 2;
        num = (num*num)%mod;
    }

    return res;
}
Mathieu
  • 8,840
  • 7
  • 32
  • 45
  • This works indeed! I had looked at this already but had difficulty implementing it, this seems to work! (but you forgot ; after res = 1 :). Thank you! – Isaiah Sep 20 '17 at 08:55
  • 1
    Apart from data types (`short` or `unsigned short`, for input parameters, return type and `res` datatype to `unsigned long int` or `long int`), this looks good. Maybe also braces around `res = (res*num) % mod`. +1. – Maximilian Gerhardt Sep 20 '17 at 08:58
  • Corner case failure: This `pow_mod(num, 0, 1)` --> 1 and should be 0. Suggest change: `int res = 1%mod;`. But then OP did say: 2 ≤ a ≤ b. – chux - Reinstate Monica Sep 20 '17 at 11:06
3

You are doing modular arithmetic in Z/bZ.

Note that, in a quotient ring, the n-th power of the class of an element is the class of the n-th power of the element, so we have the following result:

(a^b) mod b = ((((a mod b) * a) mod b) * a) mod b [...] (b times)

So, you do not need a big integer library.

You can simply write a C program using the following algorithm (pseudo-code):

  • declare your variables a and b as integers.
  • use a temporary variable temp that is initialized with a.
  • do a loop with b steps, and compute (temp * a) mod b at each step, to get the new temp value.
  • compare the result with a.

With this formula, you can see that the highest value for temp is 32768, so you can choose an integer to store temp.

Alexandre Fenyo
  • 4,526
  • 1
  • 17
  • 24
  • 1
    The number of iterations for this algorithm would be `b` multiplications, which is very bad. The square-and-multiply algorithm achieves a complexity of `log2(b)` squarings and `hammingWeight(b)` (also at max `log2(b)`) multiplications. When `b=32768` your algorithm needs `32768` iterations, whereas the simple SQ+MUL needs 15 iterations. – Maximilian Gerhardt Sep 20 '17 at 08:28
  • Yes, you are totally right. I wanted to show the **simplest** algorithm that achieve the main goal in the question: not having an overflow in C. Of course, this is not the fastest algorithm, as you said. – Alexandre Fenyo Sep 20 '17 at 08:34
  • @AlexandreFenyo thank you, this indeed works! However, I am also limited in the time it may take; the biggest Carmichael number (a number b that's not prime but also doesn't satify pow(a,b)%b == a with 2<=a<=b) within the limit for b is 29341, and this method would then have to be repeated 29341 as well, being considerably slow. – Isaiah Sep 20 '17 at 08:40
  • "highest value for temp is 32767 ^ 2" --> I'd say the highest value for temp is `32768`. The `(temp * a) mod b` needs to be done with about 2x wider math. – chux - Reinstate Monica Sep 20 '17 at 11:12
  • Yes your are right, I've patched my answer. Thanks. – Alexandre Fenyo Sep 20 '17 at 11:18