Variable Precision Integer Arithmetic: powermod(a,d,n) - File Exchange

Code covered by the BSD License

Highlights from
Variable Precision Integer Arithmetic

demo_vpi
base2vpi(B,base) bin2vpi: converts an integer in an arbitrary base into vpi (decimal) form
bin2vpi(B) bin2vpi: converts a binary representation of an integer into vpi (decimal) form
binomfactors(n,k) binomfactors: list all factors of the binomial coefficient nchoosek(n,k)
catdigits(N,M) catdigits: concatenates the digits of N and M into an aggregate number
createPrimesList createPrimesList - For users of older matlab releases, this function will generate a compatible _primeslist_ file
factorialfactors(n) factorialfactors: efficient computation of the prime factors of factorial(n)
fibonacci(n) fibonacci: vpi tool to efficiently compute the n'th Fibonacci number and the n'th Lucas number
getprimeslist loads the primeslist file, and decompresses it, returning the list of primes up to 2^26
ispalindrome(N) ispalindrome: test if the number N (vpi or numeric, or a digit string as a vector) is a palindrome
iszero(INT) vpi/iszero: test to see if a numeric object is zero
legendresymbol(a,p) legendresymbol: computes the legendre symbol (a/p) for prime p
lineardiophantine(A,B,C) lineardiophantine: solve the linear Diophantine equation, A*x + B*y = C
mersenne(p) mersenne: identify whether 2^p-1 is a Mersenne prime, using the Lucas-Lehmer test
minv(a,p)
modfibonacci(n,modulus) fibonacci: compute the n'th Fibonacci number and the n'th Lucas number, all modulo a given value
modrank(A,p) modrank: compute the rank of an integer array, modulo p
modroot(a,p)
modsolve(A,rhs,p)
nextprime(N,direction,kprimes) nextprime: finds the next larger prime number directly above (or below) N
numberOfPartitions(N) numberOfPartitions: compute the number of partitions of the positive integer n
powermod(a,d,n) vpi/powermod: Compute mod(a^d,n)
quadraticresidues(N) quadraticresidues: returns a list of the possible quadratic residues of the integer N
quotient(numerator,denominato... quotient: divides two integers, computing a quotient and remainder
subfactorial(N) subfactorial: The subfactorial of an integer (or integers) N, known as !N
totient(N) vpi/totient: the number of positive integers less than N that are coprime to N
vpi(N) vpi: Creator function for a variable precision integer
View all files

from Variable Precision Integer Arithmetic by John D'Errico
Arithmetic with integers of fully arbitrary size. Arrays and vectors of vpi numbers are supported.

powermod(a,d,n)

function R = powermod(a,d,n)
% vpi/powermod: Compute mod(a^d,n)
% usage: R = powermod(a,d,n)
% 
% powermod is MUCH faster than direct exponentiation
% with mod for large numbers. powermod does NOT
% suppoort array or vector inputs, only scalar inputs.
%
% arguments: (input)
%  a,d,n - vpi SCALAR integers, or numeric values
%
% arguments: (output)
%  R - a vpi scalar integer, representing mod(a^d,n)
%
% Example:
%  Compare exponentiation plus mod to
%  the direct application of powermod:
%
%  tic,M = powermod(vpi(123),200,497);toc
%  Elapsed time is 0.044618 seconds.
%
%  tic,M = mod(vpi(123)^200,497);toc
%  Elapsed time is 0.971667 seconds.
%
%
%  See also: power, mod, rem, quotient
%  
% 
%  Author: John D'Errico
%  e-mail: woodchips@rochester.rr.com
%  Release: 1.0
%  Release date: 1/19/09


% convert d to binary, either from a vpi
% or a double
if isnumeric(d)
  db = dec2bin(d);
else
  db = vpi2bin(d);
end
db = fliplr(db == '1');

% if a is too large, the repeated squarings will
% cause flint overflow as a double
if (a > 2^26) || (n > 2^26)
  a = vpi(a);
end

if isnumeric(a) && isnumeric(d) && isnumeric(n)
  % pure numeric
  % use the binary expansion of d to form the
  % desired power as efficiently as possible,
  % repeatedly squaring a on each pass.
  if db(1)
    R = mod(a,n);
  else
    R = 1;
  end
  for i = 2:length(db)
    if i > 2
      a2 = mod(a2*a2,n);
    else
      a2 = mod(a*a,n);
    end
    
    % do we need to multiply this power
    % of a into the result?
    if db(i)
      % take the mod on each pass through
      R = mod(R*a2,n);
    end
  end
  
else
  % use the binary expansion of d to form the
  % desired power as efficiently as possible,
  % repeatedly squaring a on each pass.
  if db(1)
    R = mod(vpi(a),n);
  else
    R = vpi(1);
  end
  for i = 2:length(db)
    if i > 2
      a2 = mod(a2*a2,n);
    else
      a2 = mod(a*a,n);
    end
    
    % do we need to multiply this power
    % of a into the result?
    if db(i)
      % take the mod on each pass through
      R = mod(R*a2,n);
    end
  end
end

Highlights from Variable Precision Integer Arithmetic

Highlights from
Variable Precision Integer Arithmetic