\\ PARI code for A343867

\\ The following program is only suitable for Latin squares of prime order.
\\ Since only odd orders that are not divisible by 3 have solutions, the first difficult case is for order 25 (n=12).

\\ Enumerates all solutions of the n-queens problem on a n X n toroidal board with a queen in the (0,0) cell.
\\ The number of solutions is A071607.
\\ Each solution is passed to the user defined function f and the sum of the ouputs from that function returned.
\\ This method is not particularly efficient.
Queens(n, f(v)=1)={
  local(v=vector(n));
  my(recurse(k, S, d1, d2)=if(k==n, f(v), 
     sum(j=1, n-1, if(!bittest(S,j) && !bittest(d1, (j-k)%n) && !bittest(d2, (j+k)%n),
                    v[1+k]=j; self()(1+k, S+(1<<j), d1+(1<<((j-k)%n)), d2+(1<<((j+k)%n))) )))
  );
  recurse(1, 1, 1, 1);
}

\\ Tests if a n-queens solution corresponds with a cyclic Latin square (not counted by this sequence)
IsCyclic(v)={my(n=#v, d=(v[1]-v[n])%n); for(i=2, n, if((v[i]-v[i-1])%n <> d, return(0))); 1}

\\ Test for a semi-cyclic square with step (1,d).
\\ With squares which are not of prime order there may also be steps that are not of this form.
IsSemi(v, d) = {my(n=#v, b=0); for(i=0, n-1, my(j=(v[i+1]-d*i)%n); if(bittest(b, j), return(0)); b=bitor(b, 1<<j)); 1}

\\ Returns number of semicylic squares corresponding to an n-queens solution.
\\ Valid only for prime orders.
Semis(v) = {my(n=#v); if(!IsCyclic(v), 2 + sum(d=1, n-1, IsSemi(v, d)))}

\\ Function for A343867.
\\ Too slow to produce more than the first 3 nonzero terms.
\\ a(6) takes about 250ms, a(8) takes approx 2 mins and a(9) many hours.
a343867(n)={
  my(m=2*n+1);
  if (m%3==0||m==1, 0, if(isprime(m), Queens(m,Semis), error("order ", m, " is not prime")))
}

\\ End