Modular Exponentiation

Question

I am trying to use this method in order to break down bases with large exponents because data types in the standard C++ library do not store numbers that large.

The problem is in the last loop where I use the fmod() function to mod my large numbers. The answer is supposed to be 1 but I am getting 16. Does someone see a problem?

#include <iostream>
#include <vector>
#include <math.h>

using namespace std;

typedef vector<int> ivec;
ivec binStorage, expStorage;

void exponents()
{
    for (int j=binStorage.size(); j>=0; j--)
        if(binStorage[binStorage.size()-j-1]!=0)
            expStorage.push_back(pow(2, j));
}

void binary(int number)
{
    int remainder;

    if(number <= 1)
    {
        cout << number;
        return;
    }

    remainder = number%2;
    binary(number >> 1);
    cout << remainder;
    binStorage.push_back(remainder);
}

int main()
{
    int num = 117;
    int message = 5;
    int mod = 19;
    int prod = 1;
    binary(num);
    cout << endl;
    exponents();

    cout << "\nExponents: " << endl;
    for (int i=0; i<expStorage.size(); i++)
        cout << expStorage[i] << " " ;

    cout << endl;
    cout << "\nMessage" << "-" << "Exponent" << endl;
    for (int i=0; i<expStorage.size(); i++)
    {
        cout << message << "-" << expStorage[i] << endl;
        prod *= fmod(pow(message, expStorage[i]), mod);
    }
    cout << "\nAnswer: " << fmod(prod, mod) << endl;
    return 0;
}

Here are my results:

1110101

Exponents:
64 32 16 4 1

Message-Exponent
5-64
5-32
5-16
5-4
5-1

Answer: 16

Process returned 0 (0x0)   execution time : 0.085 s
Press any key to continue.

Edit: Here is the problem loop.

for (int i=0; i<expStorage.size(); i++)
    {
        cout << message << "-" << expStorage[i] << endl;
        prod *= fmod(pow(message, expStorage[i]), mod);
    }

You couldn't narrow this down to one pair of input and output? I'm sure we don't need this whole program and all those irrelevant calculations. — Lightness Races in Orbit
Does that do 5 to the power of 64 at some point? ie Overflow, which misses the point of the substitution in the link — doctorlove
I didn't mean for you to dump a random loop on us. I meant for you to present a testcase, a program containing only includes, a small main, one or two constant input values, the problematic function call, and perhaps some output to demonstrate that the output isn't what it should be. This should have been your first debugging step. — Lightness Races in Orbit

101010 101010 · Accepted Answer · 2014-05-29T15:21:58

The algorithm you posted is the modular exponentiation algorithm. Following the steps in the link you posted the algorithm reduces down to the following piece of code:

#include <iostream>
#include <cmath>

// B : Base
// E : Exponent
// M : Modulo
constexpr int powermod(int const B, int const E, int const M) {
  return ((E > 1) ? (powermod(B, E / 2, M) * powermod(B, E / 2, M) * powermod(B, E % 2, M)) % M
                  : (E == 0) ? 1 : B % M);
}

int main() {
  int const e = 117;
  int const b = 5;
  int const m = 19;

  std::cout << "Answer: " << powermod(b, e, m) << std::endl;

  return 0;
}

Note, that I used constexpr. If your compiler doesn't support it, you can remove it. Using constexpr and provided that the input arguments are constant expressions, like the example above, the computation of the power exponential will take place at compile time.

Now regarding the code you posted:

It seems that fmod doesn't work well with big numbers like (5^32 and 5^64) and gives false results.
Also your code suffers from compile errors and runtime errors so I corrected it.
I coded an algorithm that computes the modulo based on recursion. Basically is a variation of the algorithm I posted above with a safe guard power of 4. (see function safemod() below):

#include <iostream>
#include <cmath>
#include <vector>

using namespace std;

typedef vector<int> ivec;

// B : Base     (e.g., 5)
// E : Exponent (e.g., 32)
// M : Modulo   (e.g., 19)
double safemod(double B, double E, double M) {
  return ((E > 4) ? fmod(safemod(B, E / 2, M) * safemod(B, E / 2, M), M)
    :
    fmod(pow(B, E), M));
}

void exponents(ivec const &binStorage, ivec &expStorage) {
  int j(pow(2.0, binStorage.size() - 1));
  for (vector<int>::const_iterator it(binStorage.begin()), ite(binStorage.end()); it != ite; ++it) {
    if (*it != 0) expStorage.push_back(j);
    j /= 2;
  }
}

void binary(int const number, ivec &binStorage) {
  if (number > 0) {
    int remainder = number % 2;
    binary(number / 2, binStorage);
    binStorage.push_back(remainder);
  }
}

int main() {
  int num     = 117;
  int message = 5;
  int mod     = 19;
  int prod    = 1;
  ivec binStorage, expStorage;

  binary(num, binStorage);
  for (size_t i(0); i < binStorage.size(); ++i) cout << binStorage[i];
  cout << endl;

  exponents(binStorage, expStorage);
  cout << "\nExponents: " << endl;
  for (size_t i(0); i<expStorage.size(); ++i) cout << expStorage[i] << " ";
  cout << endl;

  cout << "\nMessage" << "-" << "Exponent" << endl;
  for (size_t i(0); i<expStorage.size(); ++i) {
    cout << message << "-" << expStorage[i] << endl;
    prod *= safemod(message, expStorage[i], mod);
  }

  cout << "\nAnswer: " << fmod(prod, mod) << endl;

  return 0;
}

Output:

1110101

Exponents: 64 32 16 4 1

Message-Exponent

5 - 64

5 - 32

5 - 16

5 - 1

Answer: 1

Modular Exponentiation

1 Answers

Output: