function [asol,pasol,pafsol] = emirpsearchpar(nd,amax)
% emirpsearchpar - finds reversible prime pairs of the form 10^nd + a
%   This code employs the parallel processing toolbox, in the form of a
%   parfor loop.
% 
% If amax is a scalar, then we assume a will be searched over the domain
%   [1,amax]
%
% if amax is a vector of length 2, then it indicates the search domain as:
%   [amax(1),amax(2)]

tic

if isscalar(amax)
  amax = [1,amax];
end

% set of primes to test for roughness
% sqrt(flintmax('double')) is 9.49e7
% sqrt(intmax('int64')) is 3e9
roughset = int64(primes(2e9))';
np = numel(roughset);

% compute powers of 10 at and just below 10^nd, all modulo roughset.
expbin = dec2bin(nd-12) - '0';
T = ones(np,1,'int64');
pow10 = mod(int64(10),roughset);

for i = 1:numel(expbin)
  if expbin(end - i + 1)
    T = mod(T.*pow10,roughset);
  end
  pow10 = mod(pow10.*pow10,roughset);
end

Tenmods(:,1) = T;
% at this point, Tenmods(:,1) = mod(10^(nd-12),roughset)
% now we can work upwards through the columns of Tenmods,
% so that Tenmods(:,12) = mod(10^nd,roughset)

for i = 1:12
  Tenmods(:,i+1) = mod(Tenmods(:,i)*int64(10),roughset);
end

% get the entire list of rough increments beyond 10^nd.
TenNDmod = Tenmods(:,end);
ar = nextRough(amax,roughset,mod(TenNDmod + amax(1) - 1,roughset));
na = numel(ar);

disp("Number of rough increments found for a: " + na)

% next, which of those rough increments are also rough after reversal?
% get the digits, then flip.
ahatDigits = flip(dec2base(ar,10),2);
ahat = int64(base2dec(ahatDigits,10))';
nDigits = size(ahatDigits,2);

Tenmodsf = Tenmods(:,12 - nDigits + 2);

% process them in blocks, but keeping a RAM limit used for the array as
RAMlim = 12; % GB

pblock = min(np,RAMlim*2^30/(8*na));
pind = 1:pblock;
b0 = pind(end);

Hwb = waitbar(0);
while b0 <= np
  waitbar(b0/np,Hwb,"Rough search, # remaining rough candidates: " + na)

  ahatrem = mod(Tenmodsf(pind).*ahat + 1,roughset(pind));
  roughind = all(ahatrem,1);
  ahat = ahat(roughind);
  ar = ar(roughind);
  na = numel(ar);
  
  pblock = min(np,RAMlim*2^30/(8*na));
  pind = (b0 + 1):(b0+pblock);
  b0 = pind(end);

end

toc

disp(" Begin searching for rough pairs at " + (nd+1) + " digits")

% next, find pairs of rough reversible numbers, where both
% elements of the number and its reversed cousin in the pair
% are rough
roughpaircount = numel(ar);

disp("Total rough pairs found: " + roughpaircount)

% do one Fermat test, to get a prediction of the time to be done.
t1 = datetime;
pa = sym(10)^nd + ar(1);
res = isProbablePrimeFLT(pa,210);
t2 = datetime;
dt = t2 - t1;

disp("Time (in seconds) for one Fermat test: " + seconds(dt))

asol = [];pasol = [];pafsol = [];
t0 = datetime;
close(Hwb)

poolobj = parpool; % maximum number of workers possible
poolsize = poolobj.NumWorkers;

fintime = datestr(t0 + roughpaircount*dt/poolsize);
disp("Predicted time this run will finish: " + fintime)

parfor i = 1:roughpaircount
 
  % at this point, we have rough pairs in the vectors ar and ahat
  % disp("ai = " + ai + " and its cousin were rough")
  tnow = datetime;

  pa = sym(10)^nd + ar(i);

  % test using little fermat, with 210 as a witness
  if isProbablePrimeFLT(pa,210)

    disp("10^n+a passed Fermat, with a = " + ar(i))
    emirpdisplay(pa)

    paf = sym(10)^(nd - nDigits +1)*ahat(i) + 1;

    if isProbablePrimeFLT(paf,210)
      disp("ahat*10^(N-d)+1 also passed Fermat, with ahat = " + ahat(i))
      emirpdisplay(paf)

      if isProbablePrimeJ(pa) && isProbablePrimeJ(paf)
        disp("$$$$$$$$$$$$$$$$$$$$$$$$$$ An emirp prime pair has been verified $$$$$$$$$$$$$$$$$$$$$$$$$$$$$")

        asol = [asol;ar(i)];
        pasol = [pasol;pa];
        pafsol = [pafsol;paf];
        
      else
        disp("False alarm, they were Fermat pseudo-primes")
      end
  
    end % if isprobableprime

  end % for

end

delete(poolobj)

end

function roughInc = nextRough(Ndelta,roughSet,N0rem)
% nextRough - finds the next rough numbers in sequence after N0, using a simple Euclidean sieve.
% 
% arguments: 

% Set up a boolean vector to sieve 
Nbool = true(1,diff(Ndelta)+1);

nr = numel(roughSet);
for i = 1:nr
  P_i = roughSet(i);
  
  % Start point for the sieve
  Nbool(P_i - N0rem(i):P_i:end) = false;
end

roughInc = int64(find(Nbool)) + int64(Ndelta(1) - 1);

end

function TF = isProbablePrimeJ(N)
% isProbablePrimeJ - tests if N is prime, using a Miller-Rabin test
% 
% isProbablePrimeJ converts the number to a Java.Math.BigInteger form,
% then uses isprime for that class, which performs a Miller-Rabin test.
%
% N is assumed to be a symbolic toolbox integer, or an array of
% such integers.

% convert N into a character form. no test is done to insure that N is
% indeed purely numeric.
if ischar(N)
  CN = N;
else
  CN = char(N);
end

JNum = java.math.BigInteger(CN);

% now just pass the call to isProbablePrime
TF = JNum.isProbablePrime(10);

end

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


function emirpdisplay(X)
% display an emirp number, dropping out the zero digits in the middle
% 

% Get the digits as a character string
DX = char(X);
nDX = numel(DX);

% Find the middle string of zeros. It will be by far the longest string of
% zeros. So BX is true for non-zeros.
BX = DX ~= '0';

startloc = strfind(char(BX+'0'),'10');
endloc = strfind(char(BX+'0'),'01');

% what is the longest string of zeros?
zlen = endloc - startloc;
[zlen,ind] = max(zlen);
startloc = startloc(ind);
endloc = endloc(ind);

% If the zero string is at least 50 zeros in length, then zap away most of
% them. If not, then just return the entire number.
if zlen > 50
  nrem = zlen - 10;
  cDX = [DX(1:(startloc + 5)) ,'... (',num2str(nrem),' zeros redacted) ...', DX((endloc-4):end)];
else
  cDX = DX;
end

% get the name of the variable as it was passed in
Xname = inputname(1);

disp([Xname,' ='])
disp(cDX)
disp(' ')

end