Пусть у нас есть некоторое число a, которое требуется возвести в степень k(возможно, по модулю n). Можно просто умножить a само на себя k раз, но при больших размерах чисел это довольно сложная и медленная операция. Рассмотрим алгоритм, вычисляющий a^k за O(log k) операций. В квадратных скобках указаны изменения, необдимые для вычисления степени по модулю.

1. Представим d в двоичной системе счисления d = d₀2^r + ... + d_r-12 + d_r, где d_i, цифры в двоичном представлении, равны 0 или 1, d₀=1.

2. Положим a₀ =a и затем для i=1, ... ,r вычислим
a_i = a_i-1² * a^di [(mod n)].

3. a_r есть искомый вычет a^d [(mod n)].

Исходник на Си

/* Author:  Pate Williams (c) 1997 */

#include <stdio.h>

#define BITS_PER_LONG 32l
#define DEBUG

// Перевод числа в двоичную систему
long long_to_binary(long K, long *k)
{
  int found = 0;
  long a = K, i, l = 0, length;

  while (!found && l < BITS_PER_LONG) {
    found = ((a & 0x80000000l) >> 31) == 1;
    if (!found) a <<= 1, l++;
  }

  length = BITS_PER_LONG - l;
  for (i = 0; i < length; i++) k[i] = K & 1, K >>= 1;
  return length;
}

// Возведение a в степень K по модулю n
// Если модуль не нужен - уберите все вхождения n
long powmod(long a, long K, long n)
{
  long A = a, b = 1, i, k[32];
  long t = long_to_binary(K, k);

  if (K == 0) return b;
  if (k[0] == 1) b = a;

  #ifdef DEBUG
  printf("-------------\n");
  printf("i k   A    B \n");
  printf("-------------\n");
  printf("%ld %ld %4ld %4ld\n", i = 0, k[i], A, b);
  #endif

  for (i = 1; i < t; i++) {
    A = (A * A) % n;
    if (k[i]) b = (A * b) % n;
    #ifdef DEBUG
    printf("%ld %ld %4ld %4ld\n", i, k[i], A, b);
    #endif
  }

  #ifdef DEBUG
  printf("-------------\n");
  #endif

  return b;
}

int main(void)
{

  long a = 5, K = 596, n = 1234;

  printf ("%ld\n",powmod(a, K, n));

  return 0;
}

А вот - другой вид функции возведения в степень, неявно использующий двоичное представление(вычисляющий его на лету) и ненуждающийся в long2binary():

long powmod(long a, long k, long n)
{
  long b=1;

  while (k) {
    if (k%2==0) {
      k /= 2;
      a *= a; // [ a = (a*a)%n; ]
      }
    else {
      k--;
      b *= a; // [ b = (b*a)%n; ]
      }
  }
  return b;
}