\\ A400659 -- self-contained verification program, PARI/GP (2.13 or later).
\\
\\   a(n) = least k >= 1 such that A002110(n-1)*prime(n)^k - 1 is prime.
\\
\\ Run:   gp -q A400659.gp            (checks n = 1..100, the b-file range)
\\ For n = 1..30 only: delete the last two lines (verify(); quit;), then in gp
\\ type  \r A400659.gp  followed by  verify(30).
\\
\\ For every n the program finds a(n) itself and proves it, using nothing but
\\ PARI's integer arithmetic:
\\
\\  * every k < a(n): N = A002110(n-1)*prime(n)^k - 1 is shown composite, either
\\    by a prime factor p < 10^5 (with p < N), or by a Fermat witness b, i.e.
\\    b^(N-1) != 1 (mod N), which no prime N can satisfy.
\\    (N = 1, at n = 1, k = 1, is not prime by definition.)
\\
\\  * k = a(n): N is proved prime by the Lucas N+1 test, which applies because
\\    N + 1 = A002110(n-1)*prime(n)^k is completely factored.
\\
\\ The Lucas N+1 test used (Lucas; Lehmer; Brillhart-Lehmer-Selfridge 1975).
\\ Let N >= 3 be odd and let F = N + 1 be completely factored. Suppose there is
\\ one integer D with gcd(D, N) = 1 and Jacobi(D, N) = -1, and for each prime
\\ q dividing F an integer P_q, P_q == D (mod 2), Q_q = (P_q^2 - D)/4, with
\\     gcd(2*Q_q, N) = 1,
\\     U_{N+1}(P_q, Q_q) == 0 (mod N),
\\     gcd(U_{(N+1)/q}(P_q, Q_q), N) = 1,
\\ where U is the Lucas sequence U_0 = 0, U_1 = 1, U_{m+1} = P*U_m - Q*U_{m-1}.
\\ Then every prime r | N satisfies (N+1) | r - Jacobi(D, r), so r >= N and N is
\\ prime. The discriminant D must be the SAME for every q (only P may change);
\\ that is what makes the conditions for the different q combine.
\\
\\ U_m is read off x^m in (Z/NZ)[x]/(x^2 - P*x + Q), where x^m = U_m*x - Q*U_{m-1}.
\\ The powers x^((N+1)/q) for all q are computed together by a product tree, at
\\ the cost of about log2(n) full-size powers instead of n.
\\
\\ For each q the program uses the first P in 1, 3, 5, ... that satisfies the
\\ conditions, and it tries D in the order 5, -7, 9, -11, ... This is the same
\\ search order as the Python program A400659.py, which was written separately,
\\ so the two print the same lines apart from the timings.
\\
\\ The terms found are compared with the b-file, which is embedded below.

if(default(parisizemax) < 2^30, default(parisizemax, 2^30));

A400659 = [2, 1, 1, 3, 1, 1, 2, 5, 5, 4, 16, 4, 1, 12, 9, 2, 5, 3, 6, 2, 34, 3, 2, 1, 6, 65, 3, 5, 6, 3, 229, 11, 17, 5, 7, 16, 19, 15, 46, 11, 14, 20, 10, 207, 12, 12, 32, 63, 6, 93, 4, 16, 96, 5, 5, 329, 172, 34, 68, 4, 395, 150, 20, 2, 20, 1, 54, 1, 106, 5, 103, 11, 20, 34, 22, 63, 139, 48, 71, 25, 55, 45, 27, 22, 119, 9, 124, 470, 20, 51, 6, 71, 74, 95, 83, 589, 171, 30, 118, 1279];

\\ The prime factors of N + 1 = A002110(n-1)*prime(n)^k: the first n primes.
yfactors(n) = primes(n);

yN(n, k) = prod(j = 1, n - 1, prime(j)) * prime(n)^k - 1;

\\ Elements of (Z/NZ)[x]/(x^2 - P*x + Q); the coefficient of x in x^m is U_m(P, Q).
lucasRingX(P, Q, N) = Mod(Mod(1, N) * 'x, 'x^2 - P * 'x + Q);
lucasU(z) = lift(polcoef(lift(z), 1));

\\ 1 if N has a prime factor p < 10^5 with p < N, else 0. One gcd with the
\\ product of all primes below 10^5. When that gcd is N itself, N is a product
\\ of such primes and has a smaller one unless N is prime.
SIEVE = prod(i = 1, primepi(10^5 - 1), prime(i));
hasSmallFactor(N) = my(g = gcd(N, SIEVE)); g > 1 && (g < N || !isprime(N));

\\ A Fermat witness b <= 97 for composite N, or 0 if none of these bases is one.
fermatWitness(N) =
{
  forprime(b = 2, 97, if(N % b && Mod(b, N)^(N - 1) != 1, return(b)));
  0;
}

\\ [B^(F/q_1), ..., B^(F/q_m)] for S = [[q_1, e_1], ...], F = prod q_i^e_i, by a
\\ product tree: about log2(m) full-size powers instead of m.
cofactorPowers(B, S) =
{
  my(m = #S, mid, F1 = 1, F2 = 1);
  if(m == 1, return([B^(S[1][1]^(S[1][2] - 1))]));
  mid = m \ 2;
  for(i = 1, mid, F1 *= S[i][1]^S[i][2]);
  for(i = mid + 1, m, F2 *= S[i][1]^S[i][2]);
  concat(cofactorPowers(B^F2, S[1..mid]), cofactorPowers(B^F1, S[mid+1..m]));
}

\\ Lucas N+1 proof for N with N + 1 = A002110(n-1)*prime(n)^k.
\\ Returns [D, P, overrides] (overrides = [q, P_q] where P_q differs from P),
\\ or 0 if no witness was found within the search limits (never a pass).
\\ For each q the first P = 1, 3, 5, ... that works is used; all still-open q
\\ are handled together for each P.
lucasProof(N, n) =
{
  my(F = N + 1, qs = yfactors(n), qe, D, d = 5, Ps, P0, Q, X, sub, fsub, zs, pending, still, over);
  if(N < 3 || N % 2 == 0, return(0));
  qe = vector(#qs, i, [qs[i], valuation(F, qs[i])]);
  if(F != prod(i = 1, #qe, qe[i][1]^qe[i][2]) || vecmin(vector(#qe, i, qe[i][2])) == 0,
    error("N+1 is not A002110(n-1)*prime(n)^k"));
  \\ Selfridge order for D: 5, -7, 9, -11, 13, ...
  for(t = 1, 200,
    D = if(t % 2, d, -d); d += 2;
    if(gcd(D, N) != 1, if(gcd(D, N) < N, return(0)); next);
    if(kronecker(D, N) == -1, break);
    if(t == 200, return(0)));
  Ps = vector(#qe);
  pending = vector(#qe, i, i);
  forstep(P = 1, 2001, 2,                          \\ D is odd, so P must be odd
    if(#pending == 0, break);
    Q = (P^2 - D) / 4;                             \\ an integer: every D tried is 1 mod 4
    if(type(Q) != "t_INT", error("Q is not an integer"));
    if(Q == 0 || gcd(2 * Q, N) != 1, next);
    X = lucasRingX(P, Q, N);
    sub = vector(#pending, j, qe[pending[j]]);
    fsub = prod(j = 1, #sub, sub[j][1]^sub[j][2]);
    zs = cofactorPowers(X^(F / fsub), sub);        \\ x^((N+1)/q) for each open q
    if(lucasU(zs[1]^sub[1][1]) != 0, return(0));   \\ U_{N+1} != 0: N is composite
    still = List();
    for(j = 1, #pending,
      if(gcd(lucasU(zs[j]), N) == 1, Ps[pending[j]] = P,
        listput(still, pending[j])));              \\ this P does not separate q
    pending = Vec(still));
  if(#pending, return(0));
  \\ Report the most common P (the smallest, on a tie) and the exceptions.
  my(c = matreduce(Ps), best = 0);
  for(j = 1, #c~, if(c[j, 2] > best, best = c[j, 2]; P0 = c[j, 1]));
  over = List();
  for(i = 1, #qe, if(Ps[i] != P0, listput(over, [qe[i][1], Ps[i]])));
  [D, P0, Vec(over)];
}

\\ Find and prove a(n). Prints one line; returns a(n), or -1 on failure.
proveTerm(n) =
{
  my(k = 0, N, b, w, bases = List(), one = 0, byf = 0, byp = 0, parts, how, below, t0 = getabstime());
  while(1,
    k++;
    N = yN(n, k);
    if(N == 1, one = 1; next);                     \\ n = 1, k = 1: 1 is not prime
    if(hasSmallFactor(N), byf++; next);
    b = fermatWitness(N);
    if(b, listput(bases, b); next);
    \\ No witness of compositeness: N must now be proved prime, or the run stops.
    w = lucasProof(N, n);
    if(w == 0,
      if(isprime(N), printf("n=%d k=%d: prime by isprime, no Lucas witness found\n", n, k); return(-1));
      printf("n=%d k=%d: composite by isprime (no other witness found)\n", n, k);
      byp++; next);
    parts = List();
    if(one, listput(parts, "k=1 gives N=1 (not prime)"));
    if(byf || #bases || byp,
      how = List();
      if(byf, listput(how, Str(byf, " with a prime factor below 10^5")));
      if(#bases, listput(how, Str(k - 1 - one - byf - byp, " by a Fermat test, bases ", Set(Vec(bases)))));
      if(byp, listput(how, Str(byp, " by PARI's isprime")));
      listput(parts, Str("k<", k, " composite (", strjoin(Vec(how), ", "), ")")));
    below = if(#parts, strjoin(Vec(parts), "; "), "no smaller k");
    printf("n=%d  a(n)=%d  digits=%d  %s  |  k=%d prime: Lucas N+1, D=%d, P=%d%s  [%d ms]\n",
      n, k, #Str(N), below, k, w[1], w[2],
      if(#w[3], Str(", except (q,P)=", w[3]), ""), getabstime() - t0);
    return(k));
}

verify(nmax = #A400659) =
{
  my(r, bad = 0, t0 = getabstime());
  if(nmax < 1 || nmax > #A400659, error("the b-file covers n = 1..", #A400659));
  for(n = 1, nmax,
    r = proveTerm(n);
    if(r != A400659[n], bad++; printf("MISMATCH at n=%d: found %d, b-file has %d\n", n, r, A400659[n])));
  if(bad, printf("FAILED: %d of %d terms not confirmed.\n", bad, nmax),
    printf("OK: a(1)..a(%d) proved and equal to the b-file. Total %d s.\n", nmax, (getabstime() - t0) \ 1000));
  !bad;
}

\\ Optional cross-check with PARI's own primality prover, which trusts nothing
\\ above: it re-proves that N is prime at k = a(n) (not that smaller k give
\\ composites). Slow for the largest terms: crossCheck(1, 60)
crossCheck(n1, n2) =
{
  for(n = n1, n2,
    my(t0 = getabstime(), ok = isprime(yN(n, A400659[n])));
    printf("n=%d a(n)=%d isprime=%d [%d ms]\n", n, A400659[n], ok, getabstime() - t0));
}

\\ When the file is run directly (gp -q A400659.gp), check the full range and quit.
\\ To load the functions without running, comment out the next two lines.
verify();
quit;