1

I have written a code to go around the big numbers: The problem is, there is a slight problem I cant seem to catch it. It is accurate till exponent or b is 2, then 3-4, it is slightly off, then 5-6 it just starts to deviate from the true answer.

#include <iostream>
#include <conio.h>
using namespace std;

unsigned long mul_mod(unsigned long b, unsigned long n, unsigned long m);

unsigned long exponentV2(unsigned long a, unsigned long b, unsigned long m);

int main()
{
    cout << exponentV2(16807, 3, 2147483647);
    getch();
}

unsigned long exponentV2(unsigned long a, unsigned long b, unsigned long m)
{
   unsigned long result = (unsigned long)1;
   unsigned long mask = b;    // masking

   a = a % m;
   b = b % m;

   unsigned long pow = a;

   // bits to exponent
   while(mask)
   {
     //Serial.print("Binary: ");  // If you want to see the binary representation,     uncomment this and the one below
     //Serial.println(maskResult, BIN);
      //delay(1000);  // If you want to see this in slow motion 
     if(mask & 1)
     {
            result = (mul_mod(result%m, a%m, m))%m;

           //Serial.println(result);  // to see the step by step answer, uncomment this
     }
     a = (mul_mod((a%m), (a%m), m))%m;
     //Serial.print("a is ");
     //Serial.println(a);
     mask = mask >> 1;          // shift 1 bit to the left

   }
   return result;
}

unsigned long add_mod(unsigned long a, unsigned long b, unsigned long m)
{
    a = a%m;
    b = b%m;
    return (a+b)%m;
}
High Performance Mark
  • 77,191
  • 7
  • 105
  • 161
  • What answer do you get, and what do you think it should be? – Beta Oct 06 '12 at 14:16
  • i've been comparing it with wolfram alpha, and trying to track down the very sly errors. I cant give numbers right now cause I'm in the middle of that, stepping through everything. – Joey Arnold Andres Oct 06 '12 at 14:26
  • Are you sure that `unsigned long` is 64 bits everywhere the code should run? (That may not be necessary if `mul_mod` is analogous to `exponentV2`, but it would if it is basically `return ((a%m)*(b%m))%m;`.) – Daniel Fischer Oct 06 '12 at 22:35

1 Answers1

0

i just looked at your code and noticed few things:

in function:

unsigned long exponentV2(unsigned long a, unsigned long b, unsigned long m);
  • i understand that this function return a^b mod m
  • in the initialization you destroy exponent (b=b%m) this invalidates result !!! remove that line

in function:

unsigned long add_mod(unsigned long a, unsigned long b, unsigned long m);
  • you do not handle overflow (what if a+b is bigger than long ?)
  • in that case (a+b)%m is wrong,...
  • you should in case of overflow substract m*x from result and than do the modulo.
  • x must be as big so m*x is almost the size of long to eliminate overflow ...
  • so (a+b)-(m*x) also fits to the long variable

for more insight look at this: Modular arithmetics and NTT (finite field DFT) optimizations

hope it helps

Community
  • 1
  • 1
Spektre
  • 49,595
  • 11
  • 110
  • 380