function [Fn,Fnp1] = fibmod(n,p)
% fibmod - compute the nth fibonacci number, as a remainder mod p
% usage: Fn = fibmod(n,p)
%
% An O(log2(n)) algorithm is employed, in the sense that the number of
% operations will be O(log2(n)). If those computations are done using sym
% or vpij objects, the overhead of computations using those large integer
% forms is ignored.
% 
% arguments: (input)
%  n - scalar, non-negative integer
%  p - scalar integer, p >= 2
%
% arguments: (output)
%  Fn - The desired Fibonacci number, modulo p
%        The class of Fn will be the same as p.
%
%  The basic algorithm uses the following identities
%    F(2k) = F(k)*(2F(k+1) - F(k))
%    F(2k+1) = F(k+1)^2 + F(k)^2
%  then essentially doubling the index of the Fibonacci number at each
%  iteration.
%
% Examples:
% Compute the remainder: mod(fibonacci(9931),23)
% fibmod(sym(9931),23)
% ans =
%      5
%  Test (i.e., computing the Fibonacci number directly (which has 2076 decimal digits), then computing the remainder)
% mod(fibonacci(sym(9931)),23)
% ans =
% 5


% was p supplied?
if (nargin < 2) || isempty(p) || (p < 2) || (n <= 0)
  error('both n and p must be supplied, p >= 2, n >= 0, is required')
end

% insure the result will be of the same class as is p
pclass = class(p);
Fn = zeros(1,1,pclass);
Fnp1 = ones(1,1,pclass);

% can we stop early?
if n == 0
  Fn = zeros(1,1,pclass);
  Fnp1 = ones(1,1,pclass);
  return
elseif n == 1
  Fn = ones(1,1,pclass);
  Fnp1 = Fn;
  return
end
% n must be at least 2 if we drop through to here.

% we need n in binary form
nbin = dec2bin(n) - '0';
nbits = numel(nbin);

% The highest order bit of n will ALWAYS be 1. 
% what actually matters in the loop will be the next lower order bit.
% (From a different perspective, the highest order bit we have not yet used.)
% That is, the two highest order bits will be either '10' or '11'. at any
% step of this iteration, we need to worry about the next lower order
% bit of n.
% If the next lower order bit is a 0, then we will map 
%     f(k) and f(k+1) into f(2k), f(2k+1)
% If the next lower bit is a 1, then we wil map
%     f(k) and f(k+1) into f(2k+1), f(2k+2)

% loop over the bits of n. Since n must be at least 2, then nbits will
% always be at least 2.
for i = 1:nbits
  if nbin(i)
    % the incoming bit of n is a 1. We use these rules to step up:
    %    F(2k+1) = F(k+1)^2 + F(k)^2
    %    F(2k+2) = F(2k) + F(2k+1) = 2*F(k)*F(k+1) + F(k+1)^2
    [Fn,Fnp1] = deal(mod(Fnp1^2 + Fn^2 , p) , mod(Fnp1*(Fnp1 + 2*Fn) , p));

  else
    % The incoming bit of n is a 0, Just step up using 
    %    F(2k) = F(k)*(2F(k+1) - F(k))
    %    F(2k+1) = F(k+1)^2 + F(k)^2
    % putting them into Fn and Fnp1, repectively
    [Fn,Fnp1] = deal(mod(Fn*(2*Fnp1-Fn),p) , mod(Fnp1^2 + Fn^2 , p));

  end
end

end