function TF = isProbablePrimeFLT(x,w)
  % isProbablePrimeFLT: Applies a Fermat primality test to x, with witness(s) w
  % usage: TF = isProbablePrimeFLT(x,w)
  % 
  % The Fermat test compares
  %    mod(w^x-1,x) == 1
  % for a positive indication of primality. If this fails, then x is not
  % prime.
  %
  % Such failures are not terribly common for the Fermat test,
  % although they do exist, and for some values of x, ALL witnesses are
  % Fermat liars.
  % 
  % A pretest is done to insure that gcd(w,x) ==1, as otherwise x cannot
  % be prime.
  %
  % When isProbablePrimeFLT does predict x is composite, it is ALWAYS
  % correct. There are no false negatives, only false positives for being
  % prime.
  % 
  % Arguments:
  %  x - any positive integer, x >= 3, or array thereof to be tested for primality
  %  
  %  w - (optional) integer, generally scalar, with 2 <= w <= x-2
  %      If multiple witnesses are provided, then the test is performed
  %      sequentially for each witness. Note that if w does not lie in the
  %      indicated interval, then it will be reduced modulo x.
  %
  %       Default value: w = 210
  %
  % TF - boolean array of the same size and shape as x. 1 indicates x may
  %       be prime, unless w is a Fermat liar.
  %
  %  Examples:
  % 37 is prime, and 13 is coprime with 37
  %   isProbablePrimeFLT(37,13)
  %  ans =
  %    logical
  %    1
  %
  % 143 is composite, and it is coprime with 14
  % isProbablePrimeFLT(143,14)
  % ans =
  %   logical
  %    0
  %
  % 561 is composite, and it is coprime with 7, HOWEVER, 561 is a
  % Carmichael number. It will fail the Fermat test, producing a false positive.
  % Luckily these Fermat false positives are somewhat rare.
  % isProbablePrimeFLT(143,14)
  % ans =
  %   logical
  %    1
  %
  % 100 is clearly not prime, but if x and the witness w are not co-prime,
  % (relatively prime) the test does not apply.
  %   isProbablePrimeFLT(100,210)
  %     Error using isProbablePrimeFLT (line 62)
  %     Fermat test requires x to be coprime with the given witness
  %
  % See also: isprime
  %
  % John D'Errico, woodchips@rochester.rr.com, 1/15/2025

  % the default witness is 210
  if (nargin < 2) || isempty(w)
    w = 210;
  end
   
  TF = false(size(x));
  for i = 1:numel(x)
    xi = x(i);
    if xi == 2
      TF(i) = true;
    elseif mod(xi,2) ~= 0
      if any(w == xi)
        error("A witness was identical to x.")
      end
      if any(gcd(w,xi) ~= 1)
        % w and x(i) are not coprime, which insures that x(i) is composite
        TF(i) = false;
      else
        % apply the Fermat test for each witness
        for j = 1:numel(w)
          wj = w(j);
          % insure all witnesses are reduced modulo xi
          wj = mod(wj,xi);
          if powermod(wj,xi-1,xi) == 1
            % A positive indication  was found, but in case other witnesses
            % were provided, then we need to check all witnesses.
            TF(i) = true;
          else
            % if we ever get a composite signal, then we need not search
            % further for this x. This is true even if other witnesses
            % indicated primality.
            TF(i) = false;
            break
          end % if
        end % for j
      end % if 
    end % if
  end % for i
end % function