fibs2(n)

function [Fn,Ln] = fibs2(n)
% fibs2: vpi tool to efficiently compute the n'th Fibonacci number and the n'th Lucas number
% usage: [Fn,Ln] = fibs2(n)
%
% the "2*n" identities are employed here:
%  F(2*m) = F(m)*L(m)
%  L(2*m) = (5*F(m)^2 + L(m)^2)/2
%
% coupled with the addition identities
%
%  F(m+1) = (F(m) + L(m))/2
%  L(m+1) = (5*F(m) + L(m))/2


if (nargin~=1) || (numel(n) ~= 1) || (n~=round(n)) || (abs(n) > 2^53)
  error('n must be scalar and an integer, <= 2^53 in absolute value')
end

% much faster if n is not a vpi but a double.
% also, we don't need to worry about n being
% larger than 2^53, as this would have a vast
% number of digits.
n = double(n);

% ensure that n is positive, handle the negative
% cases too here.
if n < 0
  n = abs(n);
  [Fn,Ln] = fibs2(n);
  if iseven(n)
    Fn = -Fn;
  else
    Ln = -Ln;
  end
  return
end

% a few special cases to end the recursion,
% depending where we came in.
if n <= 3
  switch n
    case 2
      Fn = vpi(1);
      Ln = vpi(3);
    case 3
      Fn = vpi(2);
      Ln = vpi(4);
    case 1
      Fn = vpi(1);
      Ln = vpi(1);
    case 0
      Fn = vpi(0);
      Ln = vpi(2);
  end
  return
end
% n is now known to be at least 4.
%
% Next, we will use several identities relating
% F(2*n) and L(2*n) to F(n) and L(n), and between
% the Fibonacci and Lucas numbers.

% Is n a multiple of 2?
if iseven(n)
  % even, so use 
  [Fnover2,Lnover2] = fibs2(n./2);
  
  % Use the double argument formulas
  Fn = Fnover2.*Lnover2;
  Ln = (5 .*Fnover2.*Fnover2 + Lnover2.*Lnover2)./2;
  
else
  % n must be odd. So if we subtract 1, then
  % n-1 will be even. I could have added 1
  % instead, but I want to make n as small
  % as possible, as quickly as I can.
  [Fnover2,Lnover2] = fibs2((n-1)./2);
  
  % Use the double argument identities
  Fnminus1 = Fnover2.*Lnover2;
  Lnminus1 = (5 .*Fnover2.*Fnover2 + Lnover2.*Lnover2)./2;
  
  % this has only gotten us to F(n-1) and
  % L(n-1). use the addition identities to
  % move up one more step.
  Fn = (Lnminus1 + Fnminus1)./2;
  Ln = (5 .*Fnminus1 + Lnminus1)./2;
end
  
function result = iseven(n)
% tests if a scalar value is an even integer
if isnumeric(n)
  result = (mod(n,2) == 0);
elseif isa(n,'vpi')
  % must have been a vpi
  result = (mod(trailingdigit(n,1),2) == 0);
else
  error('n must be either numeric or vpi')
end