Exponentiation by squaring #math #recursion

exp.cpp

#include "exp.h"
#include <stack>

// Calculate the exponentiation by squaring:
//   k ^ n = (k ^ (n / 2)) ^ 2           if n is even
//        or k * (k ^ ((n - 1) / 2))^ 2  if n is odd
//   or
//   k ^ n = (k ^ 2) ^ (n / 2)           if n is even
//        or k * (k ^ 2) ^ ((n - 1) / 2) if n is odd
//
// https://en.wikipedia.org/wiki/Exponentiation_by_squaring
uint64_t pow1(unsigned int k, unsigned int n)
{
  if (!n) {
    return 1;
  }

  if (n % 2) {
    return k * pow1(k, (n - 1) / 2) * pow1(k, (n - 1) / 2);
  } else {
    return pow1(k, n / 2) * pow1(k, n / 2);
  }
}

uint64_t pow2(unsigned int k, unsigned int n)
{
  if (!n) {
    return 1;
  }

  if (n % 2) {
    return k * pow2(k * k, (n - 1) / 2);
  } else {
    return pow2(k * k, n / 2);
  }
}

uint64_t pow3(unsigned int k, unsigned int n)
{
  if (!n) {
    return 1;
  }

  // uint64_t x = pow3(k, n / 2);
  // x *= x;
  uint64_t x = pow3(k * k, n / 2);
  // so x = (k^2)^(n/2),      if n is even
  //     or (k^2)^((n-1)/2),  if n is odd

  if (n % 2) {
    x = x * k; // so x = k * (k^2)^((n-1)/2)
  }
  return x;
}

uint64_t pow4(unsigned int k, unsigned int n)
{
  if (!n) {
    return 1;
  }

  uint64_t x = pow4(k * k, /*n / 2*/n >> 1);
  return (/*n % 2*/ n & 1) ? x * k : x;
}

uint64_t pow5(unsigned int k, unsigned int n)
{
  if (!n) {
    return 1;
  }

  uint64_t r = 1; // The remaining part for the squaring.
  while (n > 1) {
    if (/*n % 2*/ n & 1) {
      r *= k;
      k *= k;
      n = (n - 1) / 2;
    } else {
      k *= k;
      n = n / 2;
    }
  }

  return r * k;
}

uint64_t pow6(unsigned int k, unsigned int n)
{
  // The `r` should be the remaining part for the squaring(in pow5).
  // However, we notice that the `r * k` is duplicated in pow5. We will get
  // the answer by `r * k` when n = 1. If we keep looping when n = 1,
  // `r` is our answer. Nevertheless, we will waste time to do `k *= k`
  // when n = 1.
  uint64_t r = 1;

  while (n) {
    if (/*n % 2*/ n & 1) {
      r *= k;
      k *= k;
      n = (n - 1) / 2;
    } else {
      k *= k;
      n = n / 2;
    }
  }

  return r;
}

uint64_t pow7(unsigned int k, unsigned int n)
{
  uint64_t r = 1;
  while (n) {
    if (/*n % 2*/ n & 1) {
      r *= k;
    }
    k *= k;
    /*n /= 2*/n >>= 1;
  }

  return r;
}

// uint64_t pow8(unsigned int k, unsigned int n)
// {
//   std::stack<unsigned int> s;
//   while(n) {
//     s.push(n);
//     n >>= 1;
//   }

//   uint64_t a = 1; // a = k^0 = 1
//   while (!s.empty()) {
//     unsigned int x = s.top(); s.pop();
//     // Let y = floor(x/2), y = x/2 if x is even, y = (x-1)/2 if x is odd.
//     // then a = k^y now.
//     if (x % 2) {      // x is odd:
//       a = k * a * a;  //   a = k^x = k^(2y+1) = k * (k^y)^2
//     } else {          // x is even:
//       a = a * a;      //   a = k^x = k^(2y) = (k^y)^2
//     }
//   }

//   return a;
// }

uint64_t pow8(unsigned int k, unsigned int n)
{
  std::stack<unsigned int> s;
  while(n) {
    s.push(n);
    n >>= 1;
  }

  uint64_t a = 1; // a = k^0 = 1
  while (!s.empty()) {
    unsigned int x = s.top(); s.pop();
    // Let y = floor(x/2), y = x/2 if x is even, y = (x-1)/2 if x is odd.
    // then a = k^y now.
    a *= a; // a = (k^y)^2 = k^(2y)
                  // x is even:
                  //   a = k^x = k^(2y)
    if (x % 2) {  // x is odd:
      a *= k;     //   a = k^x = k^(2y+1) = k * k^(2y)
    }
  }

  return a;
}

// uint64_t pow9(unsigned int k, unsigned int n)
// {
//   // The position of the highest bit of n.
//   // So we need to loop `h` times to get the answer.
//   // Example: n = (Dec)50 = (Bin)00110010, then h = 6.
//   //                               ^ 6th bit from right side
//   unsigned int h = 0;
//   for (unsigned int i = n ; i ; ++h, i >>= 1);

//   uint64_t a = 1; // a = k^0 = 1
//   // There is only one `1` in the bits of `mask`. The `1`'s position is same as
//   // the highest bit of n(mask = 2^(h-1) at first), and it will be shifted right
//   // iteratively to do `AND` operation with `n` to check `n_j` is odd or even,
//   // where n_j is defined below.
//   for (unsigned int i = 1, mask = 1 << (h-1) ; i <= h ; ++i, mask >>= 1) {
//     // Let j = h-i (looping from i = 1 to i = h), n_j = floor(n / 2^j) = n >> j
//     // (n_j = n when j = 0), x = floor(n_j / 2), then a = k^x now.
//     if (n & mask) {   // n_j is odd: x = (n_j - 1) / 2 => n_j = 2x + 1
//       a = k * a * a;  //   a = k^(n_j) = k^(2x+1) = k * (k^x)^2
//     } else {          // n_j is even: x = n_j / 2 => n_j = 2x
//       a = a * a;      //   a = k^(n_j) = k^(2x) = (k^x)^2
//     }
//   }

//   return a;
// }

// uint64_t pow9(unsigned int k, unsigned int n)
// {
//   // The position of the highest bit of n.
//   // So we need to loop `h` times to get the answer.
//   // Example: n = (Dec)50 = (Bin)00110010, then h = 6.
//   //                               ^ 6th bit from right side
//   unsigned int h = 0;
//   for (unsigned int i = n ; i ; ++h, i >>= 1);

//   uint64_t a = 1; // a = k^0 = 1
//   // There is only one `1` in the bits of `mask`. The `1`'s position is same as
//   // the highest bit of n(mask = 2^(h-1) at first), and it will be shifted right
//   // iteratively to do `AND` operation with `n` to check `n_j` is odd or even,
//   // where n_j is defined below.
//   for (unsigned int i = 1, mask = 1 << (h-1) ; i <= h ; ++i, mask >>= 1) {
//     // Let j = h-i (looping from i = 1 to i = h), n_j = floor(n / 2^j) = n >> j
//     // (n_j = n when j = 0), x = floor(n_j / 2), then a = k^x now.
//     a *= a; // a = (k^x)^2 = k^(2x)
//                     // n_j is even: x = n_j / 2 => n_j = 2x
//                     //   a = k^(n_j) = k^(2x)
//     if (n & mask) { // n_j is odd: x = (n_j - 1) / 2 => n_j = 2x + 1
//       a *= k;       //   a = k^(n_j) = k^(2x+1) = k * k^(2x)
//     }
//   }

//   return a;
// }

uint64_t pow9(unsigned int k, unsigned int n)
{
  // The position of the highest bit of n.
  // So we need to loop `h` times to get the answer.
  // Example: n = (Dec)50 = (Bin)00110010, then h = 6.
  //                               ^ 6th bit from right side
  unsigned int h = 0;
  for (unsigned int i = n ; i ; ++h, i >>= 1);

  uint64_t a = 1; // a = k^0 = 1
  // There is only one `1` in the bits of `mask`. The `1`'s position is same as
  // the highest bit of n(mask = 2^(h-1) at first), and it will be shifted right
  // iteratively to do `AND` operation with `n` to check `n_j` is odd or even,
  // where n_j is defined below.
  for (unsigned int mask = 1 << (h - 1) ; mask ; mask >>= 1) { // Run h times!
    // Let j = h-i (looping from i = 1 to i = h), n_j = floor(n / 2^j) = n >> j
    // (n_j = n when j = 0), x = floor(n_j / 2), then a = k^x now.
    a *= a; // a = (k^x)^2 = k^(2x)
    // n_j is even: x = n_j / 2 => n_j = 2x
    //   a = k^(n_j) = k^(2x)
    if (n & mask) { // n_j is odd: x = (n_j - 1) / 2 => n_j = 2x + 1
      a *= k;       //   a = k^(n_j) = k^(2x+1) = k * k^(2x)
    }
  }

  return a;
}

exp.h

// Exponentiation by squaring
#ifndef EXP_BY_SQUARING
#define EXP_BY_SQUARING

#include <cstdint>  // uint64_t

uint64_t pow1(unsigned int k, unsigned int n);
uint64_t pow2(unsigned int k, unsigned int n);
uint64_t pow3(unsigned int k, unsigned int n);
uint64_t pow4(unsigned int k, unsigned int n);
uint64_t pow5(unsigned int k, unsigned int n);
uint64_t pow6(unsigned int k, unsigned int n);
uint64_t pow7(unsigned int k, unsigned int n);
uint64_t pow8(unsigned int k, unsigned int n);
uint64_t pow9(unsigned int k, unsigned int n);

#endif // EXP_BY_SQUARING

makefile

CC=g++
CFLAGS=-Wall --std=c++11

SOURCES=test.cpp\
        exp.cpp
OBJECTS=$(SOURCES:.cpp=.o)

# Name of the executable program
EXECUTABLE=run_test

all: $(EXECUTABLE)

# $@ is same as $(EXECUTABLE)
$(EXECUTABLE): $(OBJECTS)
	$(CC) $(CFLAGS) $(OBJECTS) -o $@

# ".cpp.o" is a suffix rule telling make how to turn file.cpp into file.o
# for an arbitrary file.
#
# $< is an automatic variable referencing the source file,
# file.cpp in the case of the suffix rule.
#
# $@ is an automatic variable referencing the target file, file.o.
# file.o in the case of the suffix rule.
#
# Use -c to generatethe .o file
.cpp.o:
	$(CC) -c $(CFLAGS) $< -o $@

clean:
	rm *.o $(EXECUTABLE)

test.cpp

#include "exp.h"
#include <cassert>  // assert
#include <ctime>    // std::clock
#include <cstdint>  // uint64_t
#include <iostream> // std::cout, std::cin, std::endl

typedef struct Time
{
  double cpu;
  double wall;
} Time;

typedef struct Result
{
  uint64_t value;
  Time time;
} Result;

template <typename Function, typename... Args>
Result calculate(Function aFunction, Args&&... aArgs)
{
  std::clock_t cpuStart = std::clock();
  auto wallStart = std::chrono::high_resolution_clock::now();

  uint64_t r = aFunction(std::forward<Args>(aArgs)...);

  std::clock_t cpuEnd = std::clock();
  auto wallEnd = std::chrono::high_resolution_clock::now();
  return Result { r, Time { 1000.0 * (cpuEnd - cpuStart) / CLOCKS_PER_SEC,
                            std::chrono::duration<double, std::milli>
                              (wallEnd - wallStart).count()
                          }
                };
}

int main()
{
  unsigned int k = 0, n = 0;
  std::cout << "Enter k, n to get k^n: ";
  std::cin >> k >> n;
  Result p1 = calculate(pow1, k, n);
  Result p2 = calculate(pow2, k, n);
  Result p3 = calculate(pow3, k, n);
  Result p4 = calculate(pow4, k, n);
  Result p5 = calculate(pow5, k, n);
  Result p6 = calculate(pow6, k, n);
  Result p7 = calculate(pow7, k, n);
  Result p8 = calculate(pow8, k, n);
  Result p9 = calculate(pow9, k, n);
  assert(p1.value == p2.value &&
         p2.value == p3.value &&
         p3.value == p4.value &&
         p4.value == p5.value &&
         p5.value == p6.value &&
         p6.value == p7.value &&
         p7.value == p8.value &&
         p8.value == p9.value);
  std::cout << k << "^" << n << " = " << p1.value << std::endl;

  std::cout << "\nTime\n-----" << std::endl;
  std::cout << "1 \tcpu: " << p1.time.cpu << "ms, wall: " << p1.time.wall << " ms" << std::endl;
  std::cout << "2 \tcpu: " << p2.time.cpu << "ms, wall: " << p2.time.wall << " ms" << std::endl;
  std::cout << "3 \tcpu: " << p3.time.cpu << "ms, wall: " << p3.time.wall << " ms" << std::endl;
  std::cout << "4 \tcpu: " << p4.time.cpu << "ms, wall: " << p4.time.wall << " ms" << std::endl;
  std::cout << "5 \tcpu: " << p5.time.cpu << "ms, wall: " << p5.time.wall << " ms" << std::endl;
  std::cout << "6 \tcpu: " << p6.time.cpu << "ms, wall: " << p6.time.wall << " ms" << std::endl;
  std::cout << "7 \tcpu: " << p7.time.cpu << "ms, wall: " << p7.time.wall << " ms" << std::endl;
  std::cout << "8 \tcpu: " << p8.time.cpu << "ms, wall: " << p8.time.wall << " ms" << std::endl;
  std::cout << "9 \tcpu: " << p9.time.cpu << "ms, wall: " << p9.time.wall << " ms" << std::endl;

  return 0;
}