\section{Introduction}

\looseness-1
Let~$X$
be a quasi-projective variety over a finite
field~${\mathbbm F}_{\!q}$
of
characteristic~$p > 0$.
One~of the most natural questions concerning the arithmetic
of~$X$
is certainly to pinpoint the number of points that
$X$~has
over its base
field~${\mathbbm F}_{\!q}$.
For~example, let
$X := V(f) \subset {\text{\bf P}}^N$
be given
for~$f \in {\mathbbm F}_{\!q}[T_0,\ldots,T_N]$
a homogeneous polynomial. Then~this means just to count the number of solutions in
${\mathbbm F}_{\!q}^{N+1}$
of the equation
$f(T_0,\ldots,T_N) = 0$,
ignoring about
$(0,\ldots,0)$
and up to~scaling.

More generally, one asks for the sequence
$(\#X({\mathbbm F}_{\!q^i}))_{i\ge1}$
of integers associated
with~$X$
that is composed of the numbers of points
on~$X$
that are defined over the extension
fields~${\mathbbm F}_{\!q^i}$.
This~sequence is of practical as well as theoretical interest. For~example, it is known since the days when E.\ Artin wrote his Ph.D.~thesis~\cite[\S22]{Ar} that it is wise to form the generating~function
\begin{equation}
\label{zeta}
Z_X(t) := \exp\left(\sum_{i=1}^\infty \#X({\mathbbm F}_{\!q^i}) \frac{t^i}i\right) ,
\end{equation}
which is called the zeta function of the
variety~$X$.
In~fact,
$Z_X$
is always a rational~function. According~to \cite[Exp.~XV, \S3, n$^\circ$2 and~3]{SGA5}, it is the alternating product of the Weil
poly\-no\-mi\-als~$\chi_i$
associated
with~$X$,
the characteristic polynomials of the operation of the geometric Frobenius
$\mathop{\text{\rm Frob}}\nolimits$
on the
\mbox{$l$-adic}
cohomology vector spaces
$H_{c,{{\text{\rm {\'e}t}}}}^i(X_{\overline{\mathbbm F}_{\!q}}, {\mathbbm Q}_l)$
with compact~support,
\begin{equation}
\label{ratzeta}
Z_X(t) = \!\!\prod_{j=0}^{2\dim X}\!\!\!\! \det \big(1-t\mathop{\text{\rm Frob}}\nolimits | H_{c,{{\text{\rm {\'e}t}}}}^j(X_{\overline{\mathbbm F}_{\!q}}, {\mathbbm Q}_l)\big)^{(-1)^{j+1}} =: \frac{\chi_1(t)\chi_3(t)\cdots\chi_{2\dim X-1}(t)}{\chi_0(t)\chi_2(t)\chi_4(t)\cdots\chi_{2\dim X}(t)} \, .
\end{equation}
Moreover,~at least when
$X$
is proper and smooth, the Weil polynomials
$\chi_j$
have remarkable properties. For~example, as was conjectured by A.~Weil and proven by P.~Deligne in \cite[Th\'eo\-r\`eme~(1.6)]{DeW1} and \cite[Corollaire~(3.3.9)]{DeW2}, every complex root
of~$\chi_j$
is of absolute value
$q^{j/2}$.
In~particular, no cancellations occur in formula~(\ref{ratzeta}).

Assume~now that the geometry of the variety
$X_{\overline{\mathbbm F}_{\!q}}$
is well understood. By~which we mean that at least all its
\mbox{$l$-adic}
Betti numbers are~known. This~is usually the case in~practice, as the~varieties considered are curves of known genus, abelian varieties of known dimension, non-singular complete intersections of known multidegree,
$K3$~surfaces,~etc.
From~the algorithmic point of view, formula~(\ref{ratzeta}) then shows that it is easy to calculate
$\#X({\mathbbm F}_{\!q^i})$
for some very
large~$i$,
provided that
$\#X({\mathbbm F}_{\!q}), \#X({\mathbbm F}_{\!q^2}), \ldots, \#X({\mathbbm F}_{\!q^k})$
have been found. Thereby,~the
bound~$k$
only depends on the geometry
of~$X_{\overline{\mathbbm F}_{\!q}}$.
Indeed,~once
$k$~is
properly chosen, (\ref{zeta}) and~(\ref{ratzeta}) allow to compute all coefficients of the
polynomials~$\chi_j$
from the numbers of points known. Then~to pinpoint
$\#X({\mathbbm F}_{\!q^i})$
is just the calculation of the
\mbox{$i$-th}
power series coefficient of
$\sum_j (-1)^{j+1} \log\chi_j$.

Unfortunately,~for most nontrivial types of varieties, the
bound~$k$
turns out to be too large, such that it is not feasible to count
$\#X({\mathbbm F}_{\!q}), \#X({\mathbbm F}_{\!q^2}), \ldots, \#X({\mathbbm F}_{\!q^k})$
in a naive way, even for medium-sized values
of~$q$.
For~instance, when
$X$~is
a
$K3$~surface,
things are relatively~simple. In~fact,
$\chi_1(t) = \chi_3(t) = 1$,
$\chi_0(t) = 1-t$,
$\chi_4(t) = 1-q^2t$,
and only
$\chi_2$
varies,
which is of
degree~$22$.
Nevertheless,~one may only hope that
$k=10$
is~sufficient. It~is known that there are cases where
$k$~has
to be chosen larger~\cite[Table~6]{EJants9}.\smallskip

{\em Our motivation. Real and complex multiplication for
$K3$~surfaces.}
Our motivation comes from projective varieties
$X$
that are defined over a number
field~$K$.
As~was first noticed by R.~van Luijk, from point counts on the reductions
$X_{{\mathfrak{p}}_1}$,
$X_{{\mathfrak{p}}_2}$
modulo two primes~\mbox{\cite{vL05,vL07,EJants8,EJ12}} of good reduction or sometimes only one~\cite{EJ11}, one may determine the geometric N\'eron-Severi rank
of~$X$.
This~applies well even to varieties of general type~\cite{EK}, but a nontrivial case is provided already by
$K3$~surfaces.

The very basic idea behind van Luijk's method is that, in the case of a
$K3$~surface,
$\smash{\mathop{\text{\rm rk}}\nolimits\mathop{\text{\rm Pic}}\nolimits X_{\overline{\mathbbm F}_{\!{\mathfrak{p}}}}}$
is always even and one has
$\smash{\mathop{\text{\rm rk}}\nolimits\mathop{\text{\rm Pic}}\nolimits X_{\overline{\mathbbm F}_{\!{\mathfrak{p}}}} \geq \mathop{\text{\rm rk}}\nolimits\mathop{\text{\rm Pic}}\nolimits X_{\overline{K}}}$.
Moreover,~one hopes to find a prime such that
$\smash{\mathop{\text{\rm rk}}\nolimits\mathop{\text{\rm Pic}}\nolimits X_{\overline{\mathbbm F}_{\!{\mathfrak{p}}}}}$
is actually equal to
$\mathop{\text{\rm rk}}\nolimits\mathop{\text{\rm Pic}}\nolimits X_{\overline{K}}$
or~$\mathop{\text{\rm rk}}\nolimits\mathop{\text{\rm Pic}}\nolimits X_{\overline{K}} + 1$.
As~was observed by F.~Charles~\cite{Ch}, the existence of such primes is related to whether
$X$
has real multiplication or~not. More~precisely, existence is provided unless
$X$
has real multiplication and
$(22 - \mathop{\text{\rm rk}}\nolimits\mathop{\text{\rm Pic}}\nolimits X_{\overline{K}})/[E:{\mathbbm Q}]$
is odd, for
$E$
the endomorphism~field.

We~say that a
$K3$~surface~$X$
has real or complex multiplication when the endomorphism algebra of the transcendental part
$T \subset H^2(X({\mathbbm C}), {\mathbbm Q})$,
considered as a pure
\mbox{${\mathbbm Q}$-Hodge}
structure, is strictly larger
than~${\mathbbm Q}$.
It~is then either a totally real field or a CM~field \cite[Theorem~1.6.a) and Theorem~1.5.1]{Za83}. Starting~from this definition, it seems, however, hard to conclude anything for concrete~examples. One~may deduce~\cite{vGe} at least that having RM or CM is a property of positive codimension in the analytic moduli space of
$K3$~surfaces
of fixed Picard rank (unless the Picard rank is~20, in which case every
$K3$~surface
has~CM). Thus,~having RM or CM should be thought of as being exceptional.

In~particular, it seems that one does not find examples of RM or CM~surfaces
just by~accident. Our~approach is therefore to perform our searches through the {\em arithmetic consequences\/} of RM and CM. These~include that, for all
primes~${\mathfrak{p}}$
that are (at least partially) inert
in~$E$,
the reduction
$X_{\mathfrak{p}}$
is non-ordinary. I.e.,~that one has
$\#X_{\mathfrak{p}}({\mathbbm F}_{\!{\mathfrak{p}}}) \equiv 1 \pmod p$
for
$p$
the prime number
below~${\mathfrak{p}}$.
Moreover,~the transcendental factors of the corresponding Weil polynomials are all partially split over the same field, the endomorphism
field~$E$.
Cf.~the top of Section~\ref{RMCM} for a more precise~statement.
In~order to detect the very few surfaces showing such a behaviour within a large family, efficient methods for point counting are asked~for.\smallskip

{\em Harvey's
\mbox{$p$-adic}
method for counting points.}
The most efficient methods for counting points on a variety over a finite field that are known today are
\mbox{$p$-adic}
in nature. One~of them was originally developed by K.~Kedlaya in~\cite{Ke}. Kedlaya's point of view was cohomological. This~means, he computed the characteristic polynomial of the Frobenius operation on the Monsky-Washnitzer cohomology, a well-behaved cohomology theory with
\mbox{$p$-adic}
coefficients for affine varieties in
characteristic~$p$.
This~approach seems to be very natural. It~is, however, still limited to curves, and in fact to particular types of~them. T.~Satoh's
\mbox{$p$-adic}
algorithm~\cite{Sa} for ordinary elliptic curves could perhaps be seen as a predecessor of this~method.

On~the other hand, as was shown by D.~Harvey~\cite{Ha}, an entirely elementary approach is possible, as well as a generalisation to arbitrary schemes of finite type over a finite~field. Harvey's method is, at least partially, based on earlier~ideas. We~are most likely unaware of some of its predecessors, but the work of D.~Wan~\cite{Wa99,Wa08}, as well as that of A.~Lauder and D.~Wan~\cite{LW}, certainly should be mentioned. Probably,~B.\ Dwork~\cite{Dw60a,Dw60b} is the very first, to whom the ideas behind~\cite{Ha} may be traced~back.

We~report on a variation of Harvey's
\mbox{$p$-adic}
method for double covers of the projective space, which we implemented as far as required for the needs of our search for
$K3$~surfaces
of degree two having RM or~CM. Although~it was never mentioned in print, this variation was certainly known, or at least obvious, to D.~Harvey~before. Generally~speaking, it seems that the algorithm described in~\cite{Ha}, although general in theory, may practically be implemented only in a form specialised to particular types of~varieties. Our code will be available within {\tt magma}~\cite{BCP}, version~2.22.\medskip

{\em Acknowledgement.}
We wish to thank David Harvey for valuable one-to-one talks, as well as the two anonymous referees for their suggestions on how to improve this~note.

\section{Counting points on varieties over finite fields--Elementary methods}

Let us assume that
$\smash{X \subset {\text{\bf A}}^N_{{\mathbbm F}_{\!q}}}$
is an affine hypersurface. I.e., that
$X = V(f)$
for an arbitrary polynomial
$f \in {\mathbbm F}_{\!q}[T_1,\ldots,T_N]$.
At~least in theory, this is not a serious restriction as every variety over a base field is birationally equivalent to a hypersurface~\cite[Chap.~I, Prop.~4.9]{Hart}, cf.\ \cite[\S1.3]{Ha}. In~practice, of course, too complicated transformations are undesirable, as they obviously reduce the speed of the~method. Even~if that concerns only a constant~factor.

\begin{ttt}[Determination of the points]
The most naive approach to point counting is certainly to {\em determine\/} all the points. If,~for whatever application, the actual points are asked for then there is no alternative to such a~method. Moreover,~to have some code available that realises a naive approach is useful for testing more advanced~algorithms.

In~order to determine the
${\mathbbm F}_{\!q^i}$-rational
points
on~$X$,
one has to run an iterated loop over all
$(N-1)$-tuples
$\smash{(x_1,\ldots,x_{N-1}) \in {\mathbbm F}_{\!q^i}^{N-1}}$.
Each~time, the roots in
${\mathbbm F}_{\!q^i}$
of a univariate
polynomial~$g$
have to be found, which just means to compute
$\smash{\gcd(g, T^{q^i}-T)}$.
Thus,~the complexity to determine the
${\mathbbm F}_{\!q^i}$-rational
points is essentially
$O(q^{(N-1)i})$.
\end{ttt}

\begin{rem}
When~$X$
is defined over
${\mathbbm F}_{\!q}$,
but points over
${\mathbbm F}_{\!q^i}$
for
$i>1$
are sought for, one may gain a factor
of~$i$
by dealing with Frobenius orbits instead of~points.
\end{rem}

\begin{ttt}[FFT point counting]
Assume that
$X$
is given by a {\em decoupled\/}
polynomial~$f$.
I.e., that
$f = f_1(T_1,\ldots,T_M) + f_2(T_{M+1},\ldots,T_N)$
for some
$1 < M < N$.
Then~$\#X({\mathbbm F}_{\!q^i}) = (c_1 * c_2)(0)$
in terms of the two counting functions
$c_1$
and~$c_2$,
given by
\begin{align*}
c_1(c) & := \#\{(x_1,\ldots,x_M) \in {\mathbbm F}_{\!q^i}^M \mid f_1(x_1,\ldots,x_M) = c \} \tag*{and\hspace{1.5cm}}\\
c_2(c) & := \#\{(x_{M+1},\ldots,x_N) \in {\mathbbm F}_{\!q^i}^{N-M} \mid f_2(x_{M+1},\ldots,x_N) = c \} \, .
\end{align*}
In~order to compute
$\#X({\mathbbm F}_{\!q^i})$,
one first compiles look-up tables for
$c_1$
and~$c_2$
and then calculates
$(c_1 * c_2)(0)$,
according to the~definition. The complexity of this method to count the
${\mathbbm F}_{\!q^i}$-rational
points is
$O(q^{\max\{M,N-M\}i})$.

Similarly,~for the double cover of
${\text{\bf A}}^N$
given by
$W^2 = f$,
one~has
$$\#X({\mathbbm F}_{\!q^i}) = q^{iN} + \sum_{\smash{c \in {\mathbbm F}_{\!q^i}}} \!\chi(c) \!\cdot\! (c_1 * c_2)(c) \, ,$$
where~$\chi\colon {\mathbbm F}_{\!q^i} \to \{-1,0,1\}$
denotes the quadratic~character. The~convolution
$c_1 * c_2$
then has to be computed using the FFT method~\cite[S\"atze~20.2 and~20.3]{Fo} or \cite[Theorem~32.8]{CLR}. Here,~again, the complexity is essentially
$O(q^{\max\{M,N-M\}i})$,
with an additional
\mbox{$\log$-factor}
coming from the Fourier~transforms.
Cf.~\cite[Algorithm~17]{EJants8} for more details on this~method.
\end{ttt}

\begin{ex}[cf.~{\cite[Example~3.13]{EJ15}}]
Consider the
$K3$~surface
$X$
of
degree~$2$
over~${\mathbbm F}_{\!7}$,
given~by
\begin{align*}
W^2 = 6 T_0^6 + 6T_0^5T_1 + 2T_0^5T_2 + 6T_0^4T_1^2 + 5T_0^4T_2^2 + 5T_0^3T_1^3 + T_0^2T_1^4 + 6T_0T_1^5 +{}& 5T_0T_2^5 \\[-1mm]
&{} + 3T_1^6 + 5T_2^6 \,.
\end{align*}
The polynomial on the right hand side does not contain any monomial involving both
$T_1$
and~$T_2$.
Thus,~over the affine plane
${\text{\bf P}}^2_{\overline{\mathbbm F}_{\!7}} \!\setminus\! V(T_0)$,
the double cover
$X$~is
given by a decoupled~polynomial.
The~numbers of points over
${\mathbbm F}_{7},\ldots,{\mathbbm F}_{7^{10}}$
are
\begin{align*}
60,~
2\,488,~
118\,587,~
5\,765\,828,~
282\,498\,600,&~
13\,841\,656\,159,~
678\,225\,676\,496,\\
33\,232\,936\,342\,644,~
&1\,628\,413\,665\,268\,026,~
79\,792\,266\,679\,604\,918.
\end{align*}
Using~the FFT method, it took about 2 hours of CPU time and required approximately 5\,GB of memory to compute these numbers on an AMD 248 Opteron processor running at 2.2\,GHz.
\end{ex}

\begin{ttt}[Making \'etale cohomology explicit]
For a general variety, it is probably a tough ask to make its
\mbox{$l$-adic}
cohomology vector spaces explicit, including the Frobenius operation on~them. There~are, however, certain special types of varieties, for which this is possible. For~example, for an elliptic curve, one has
$\smash{H_{{\text{\rm {\'e}t}}}^1(X_{\overline{\mathbbm F}_{\!q}}, {\mathbbm Q}_l) \cong \mathop{\underleftarrow{\rm lim}}\limits_n E(\overline{\mathbbm F}_{\!q})[l^n]} \!\otimes_{{\mathbbm Z}_l}\! {\mathbbm Q}_l$,
which is the starting point of Schoof's algorithm~\cite{Sch}.

A~different kind are those varieties, for which all the cycle maps~\cite[Cycle, \S2.2.10]{SGA4h}
$$\mathop{\text{\rm cl}}\nolimits\colon \mathop{\text{\rm CH}}\nolimits^j(X_{\overline{\mathbbm F}_{\!q}}) \!\otimes_{\mathbbm Z}\! {\mathbbm Q}_l \longrightarrow H_{{\text{\rm {\'e}t}}}^{2j}(X_{\overline{\mathbbm F}_{\!q}}, {\mathbbm Q}_l(j))$$
are~surjective and, moreover,
$H_{{\text{\rm {\'e}t}}}^j(X_{\overline{\mathbbm F}_{\!q}}, {\mathbbm Q}_l) = 0$
for every odd
integer~$j$.
These~assumptions are fulfilled, for instance, for surfaces that are geometrically~rational.

A~case for which this approach is implemented in {\tt magma} is that of a cubic~surface.
Here,~$\mathop{\text{\rm CH}}\nolimits^1(X_{\overline{\mathbbm F}_{\!q}}) = \mathop{\text{\rm Pic}}\nolimits(X_{\overline{\mathbbm F}_{\!q}}) \cong {\mathbbm Z}^7$
is generated by the 27~lines. One~has to determine the~lines by a Gr\"obner base calculation. From~the operation
of~$\mathop{\text{\rm Frob}}\nolimits$
on the lines, one deduces the operation
on~$H_{{\text{\rm {\'e}t}}}^2(X_{\overline{\mathbbm F}_{\!q}}, {\mathbbm Q}_l(1))$.
\end{ttt}

\begin{ex}
\leavevmode
\vspace{-.8mm}
{\footnotesize
\begin{verbatim}
> p := NextPrime(31^31);
> p;
17069174130723235958610643029059314756044734489
> rr<x,y,z,w> := PolynomialRing(GF(p),4);
> gl := x^3 + 2*y^3 + 3*z^3 + 5*w^3 - 7*(x+y+z+w)^3;
> time NumberOfPointsOnCubicSurface(gl);
291356705504951337929728147720395643095420441743546148074332970578622743757660955298550825611
Time: 1.180
\end{verbatim}}
\end{ex}

\begin{rem}
The~list of methods given in the section is by no means meant to be exhaustive. There~is at least one further idea, which we can mention only in~passing, namely the use of a~fibration. For~more information, we advise the reader to consult the articles \cite{La} of A.~Lauder and~\cite{PT} of S.~Pancratz and J.~Tuitman.
\end{rem}

\section{Harvey's
\mbox{$p$-adic}
method}

As~above, we assume that
${\mathbbm F}_{\!q}$
is a finite field of
characteristic~$p>0$.
The~method is
\mbox{$p$-adic}
in~nature. It~finds an approximation
of~$\#X({\mathbbm F}_{\!q^i})$
with respect to the
\mbox{$p$-adic}
valuation. Thus,~in~order to pinpoint
$\#X({\mathbbm F}_{\!q^i})$
exactly, one needs, in addition, an estimate for the number of~points.\medskip\smallskip

{\em Characteristic functions.}

\begin{nota}
\label{quasi}

\begin{iii}
\item
\looseness-1
We~continue to assume that
$X$
is an affine hypersurface. However,~we now suppose that
$\smash{X := V(f) \cap {\mathbbm G}_{m,{\mathbbm F}_{\!q}}^N \subset {\mathbbm G}_{m,{\mathbbm F}_{\!q}}^N}$
is a subscheme of an affine~torus. Again,~this restriction is not serious as
${\text{\bf A}}^N_{{\mathbbm F}_{\!q}}$
is a finite union of affine~tori. For~varieties more general than hypersurfaces, one might, in principle, simplify the situation by iteratively projecting away from points, cf.~\cite[\S1.3]{Ha}.
\item
For~every natural
number~$i$,
let us denote by
$\mu_{q^i-1}$
the group of all
\mbox{$(q^i-1)$-th}
roots~of unity and write
${\mathbbm Z}_{q^i}$
for the integer ring of the local field
${\mathbbm Q}_p(\mu_{q^i-1})$,
which is an unramified extension
of~${\mathbbm Q}_p$.
${\mathbbm Z}_{q^i}$~is
thus a complete discrete valuation ring with residue
field~${\mathbbm F}_{\!q^i}$.
Furthermore,~$\mu_{q^i-1} \subset {\mathbbm Z}_{q^i}$
is bijectively mapped onto
$\smash{{\mathbbm F}_{\!q^i}^*}$
under the residue map
$$\pi = \pi_i\colon {\mathbbm Z}_{q^i} \twoheadrightarrow {\mathbbm F}_{\!q^i} \, .$$
This~is, of course, an instance of taking Teichm\"uller representatives.
\item
Moreover,~one has
${\mathbbm Z}_{q^i} \subseteq {\mathbbm Z}_{q^{i'}}$
if and only if
${\mathbbm F}_{\!q^i} \subseteq {\mathbbm F}_{\!q^{i'}}$.
In~this case, we introduce the {\em quasi norm\/} map
\begin{eqnarray*}
N_{{\mathbbm Z}_{q^{i'}}/{\mathbbm Z}_{q^i}}\colon {\mathbbm Z}_{q^i}[T_1,\ldots,T_N] &\longrightarrow& {\mathbbm Z}_{q^i}[T_1,\ldots,T_N] \, , \\[-1mm]
g &\mapsto& g \!\cdot\! g^{(q^i)} \!\cdot\! g^{(q^{2i})} \cdots g^{(q^{i'-i})} \, ,
\end{eqnarray*}
for
$g^{(n)}$
the {\em quasi power,} given by
$g^{(n)}(T_1,\ldots,T_N) := g(T_1^n,\ldots,T_N^n)$.

The quasi norm is compatible with the usual norm map
$\smash{\mathop{\text{\rm N}}\nolimits_{{\mathbbm F}_{\!q^{i'}}/{\mathbbm F}_{\!q^i}}\colon {\mathbbm F}_{\!q^{i'}} \to {\mathbbm F}_{\!q^i}}$
between the finite residue fields in the sense that, for arbitrary
$g \in {\mathbbm Z}_{q^i}[T_1,\ldots,T_N]$
and
$z_1,\ldots,z_N \in {\mathbbm Z}_{q^{i'}}$,
one~has
$N_{{\mathbbm Z}_{q^{i'}}/{\mathbbm Z}_{q^i}}g(z_1,\ldots,z_N) \in {\mathbbm Z}_{q^i}$
and
$$\pi_i(N_{{\mathbbm Z}_{q^{i'}}/{\mathbbm Z}_{q^i}}g(z_1,\ldots,z_N)) = \mathop{\text{\rm N}}\nolimits_{{\mathbbm F}_{\!q^{i'}}/{\mathbbm F}_{\!q^i}}(\overline{g}(\pi_{i'}(z_1),\ldots,\pi_{i'}(z_N))) \, .$$
Here,~$\overline{g} \in {\mathbbm F}_{\!q^i}[T_1,\ldots,T_N]$
denotes the reduction
of~$g$
modulo~$(p)$.
\item
We~choose a lift
$F \in {\mathbbm Z}_q[T_1,\ldots,T_N]$
of~$f$.
\end{iii}
\end{nota}

\begin{defi}
We~call a function
$\Phi\colon \mu_{q^i-1}^N \to {\mathbbm Z}_q$
a {\em characteristic function\/}
of~$X({\mathbbm F}_{\!q^i})$
{\em modulo\/}~$p^l$
if, for every
$(z_1,\ldots,z_N) \in \mu_{q^i-1}^N$,
\begin{align*}
&\Phi(z_1,\ldots,z_N) \equiv 0 \pmod {p^l} \;\;\text{ if } (\pi(z_1),\ldots,\pi(z_N)) \not\in X({\mathbbm F}_{\!q^i}) \quad\quad\text{ and} \\
&\Phi(z_1,\ldots,z_N) \equiv 1 \pmod {p^l} \;\;\text{ if } (\pi(z_1),\ldots,\pi(z_N)) \in X({\mathbbm F}_{\!q^i}) \, .
\end{align*}
\end{defi}

\begin{prop}
\label{charfkt}
Let\/~$a$
be a positive integer such that\/
$a(q-1) \geq l$
and\/
$i \in {\mathbbm N}$
be~arbitrary. Then
\begin{eqnarray*}
\Phi & := & (1-N_{{\mathbbm Z}_{q^i}/{\mathbbm Z}_q}F^{a(q-1)})^l \\
 & = & 1 + \sum\limits_{k=1}^l (-1)^k \binom{l}k N_{{\mathbbm Z}_{q^i}/{\mathbbm Z}_q}F^{ka(q-1)}
\end{eqnarray*}
is a characteristic function
of\/~$X({\mathbbm F}_{\!q^i})$
modulo\/~$p^l$.\medskip

\noindent
{\bf Proof.}
{\em
By~definition,
$(\pi(z_1),\ldots,\pi(z_N)) \in X({\mathbbm F}_{\!q^i})$
if and only if
$f(\pi(z_1),\ldots,\pi(z_N)) = 0 \in {\mathbbm F}_{\!q^i}$,
which is equivalent to
$\mathop{\text{\rm N}}\nolimits_{{\mathbbm F}_{\!q^i}/{\mathbbm F}_{\!q}}(f(\pi(z_1),\ldots,\pi(z_N))) = 0 \in {\mathbbm F}_{\!q}$
and
$$N_{{\mathbbm Z}_{q^i}/{\mathbbm Z}_q}F(z_1,\ldots,z_N) \equiv 0 \pmod p \, .$$
If~this is true then
$\smash{N_{{\mathbbm Z}_{q^i}/{\mathbbm Z}_q}F^{a(q-1)}(z_1,\ldots,z_N) \equiv 0 \pmod {p^{a(q-1)}}}$,
hence the same modulo
$p^l$,
which implies
$\smash{1-N_{{\mathbbm Z}_{q^i}/{\mathbbm Z}_q}F^{a(q-1)}(z_1,\ldots,z_N) \equiv 1 \pmod {p^l}}$,
and
$$(1-N_{{\mathbbm Z}_{q^i}/{\mathbbm Z}_q}F^{a(q-1)}(z_1,\ldots,z_N))^l \equiv 1 \pmod {p^l} \, .$$
On~the other hand, in the opposite case, one has
$\smash{N_{{\mathbbm Z}_{q^i}/{\mathbbm Z}_q}F^{q-1}(z_1,\ldots,z_N) \equiv 1 \pmod p}$
and, consequently,
$1-N_{{\mathbbm Z}_{q^i}/{\mathbbm Z}_q}F^{a(q-1)}(z_1,\ldots,z_N) \equiv 0 \pmod p$,
which shows
$$(1-N_{{\mathbbm Z}_{q^i}/{\mathbbm Z}_q}F^{a(q-1)}(z_1,\ldots,z_N))^l \equiv 0 \pmod {p^l} \, , $$
as required.
\ifmmode\rule[-2pt]{0pt}{1pt}\ifinner\tag*{$\square$}\else\eqno{\square}\fi\else\vadjust{}\hfill$\square$\fi
}
\end{prop}

\begin{rems}

\begin{iii}
\item
According to the definition, if
$\Phi$
is a characteristic function
of~$X({\mathbbm F}_{\!q^i})$
modulo~$p^l$
then
\begin{equation}
\label{summ}
\#X({\mathbbm F}_{\!q^i}) \equiv \!\!\!\!\!\!\!\!\!\sum_{(z_1,\ldots,z_N) \in \mu_{q^i-1}^N}\!\!\!\!\!\!\!\!\!\!\! \Phi(z_1,\ldots,z_N) \pmod {p^l} \, .
\end{equation}
\item
One~has
$$\sum_{z \in \mu_{q^i-1}} \!\!\!\!z^e = 0$$
unless
$e$
is a multiple
of~$q^i-1$.
Therefore,~the right hand side of~(\ref{summ}) is
$(q^i-1)^N$~times
the sum of the coefficients
of~$\Phi$
at the monomials of the form
$\smash{(T_1^{e_1} \cdots T_N^{e_N})^{q^i-1}}$,
for~$e_1, \ldots, e_N \in {\mathbbm Z}_{\geq0}$.
\item
For~$l=1$,
these facts are known for a long time, cf.~\cite[Sec.~1.1, Theorem~4]{BS}.
\item
One might want to take the step from
$l=1$
to a larger value in a more naive manner as~follows.
For~$(\pi(z_1),\ldots,\pi(z_N)) \not\in X({\mathbbm F}_{\!q^i})$,
one has
$N_{{\mathbbm Z}_{q^i}/{\mathbbm Z}_q}F^{a(q-1)}(z_1,\ldots,z_N) = 1 + sp$,
for some
unknown~$s \in {\mathbbm Z}_{q^i}$.
This~yields
\begin{align*}
N_{{\mathbbm Z}_{q^i}/{\mathbbm Z}_q}F^{2a(q-1)}(z_1,\ldots,z_N) & = 1 + 2sp +  s^2p^2 \, , \\
N_{{\mathbbm Z}_{q^i}/{\mathbbm Z}_q}F^{3a(q-1)}(z_1,\ldots,z_N) & = 1 + 3sp + 3s^2p^2 + s^3p^3 \, , \\[-2mm]
& \hspace{1.8mm}\vdots \\
N_{{\mathbbm Z}_{q^i}/{\mathbbm Z}_q}F^{la(q-1)}(z_1,\ldots,z_N) & = \textstyle 1 + lsp + \binom{l}2 s^2p^2 + \binom{l}3 s^3p^3 + \cdots + \binom{l}{l-1} s^{l-1}p^{l-1} + \binom{l}l s^lp^l \, .
\end{align*}
There~is a linear combination of these expressions that eliminates all terms containing
$p, p^2, \ldots, p^{l-1}$.
Namely,~one needs to take as the coefficient vector the first row of the Pascal-like matrix
$$
\left(
\begin{array}{ccccc}
1 & 1 & 0 & \cdots & 0 \\
1 & 2 & 1 & \cdots & 0 \\[-1mm]
\vdots & \vdots & \vdots & \ddots & \vdots \\
1 & l\!-\!1 & \binom{l-1}2 & \cdots & \binom{l-1}{l-1} \\
1 & l & \binom{l}2 & \cdots & \binom{l}{l-1}
\end{array}
\right)^{-1} .
$$
However,~well-known relations between binomial coefficients reveal that the first row of this matrix is exactly
$(l, -\binom{l}2, \binom{l}3, \ldots, (-1)^{l+1}\binom{l}l)$,
in agreement with Proposition~\ref{charfkt}.
\end{iii}
\end{rems}

{\em Linear operators.}

\begin{nota}

\begin{iii}
\item
We~consider the polynomial ring
${\mathbbm Z}_q[T_1,\ldots,T_N]$
as a free
\mbox{${\mathbbm Z}_q$-module}
and equip it with the symmetric bilinear form
$\langle.,.\rangle$,
given by
$$\left\langle \sum_{i_1,\ldots,i_N}\!\!\! a_{i_1,\ldots,i_N} T_1^{i_1}\cdots T_N^{i_N}, \sum_{i_1,\ldots,i_N}\!\!\! b_{i_1,\ldots,i_N} T_1^{i_1}\cdots T_N^{i_N} \right\rangle := \sum_{i_1,\ldots,i_N}\!\!\! a_{i_1,\ldots,i_N}b_{i_1,\ldots,i_N} \, .$$
Then~the monomials form an orthonormal~basis. Furthermore,~the bilinear form
$\langle .,. \rangle$
yields a
norm~$\|.\|$
on~${\mathbbm Z}_q[T_1,\ldots,T_N]$
such that
$\|g\| := [\|\langle g, g\rangle \|_p]^{1/2}$,
for
$\|.\|_p$
the usual normalised
\mbox{$p$-adic}~valuation.
The~completion
${\mathbbm Z}_q\{\{T_1,\ldots,T_N\}\}$
is a
\mbox{${\mathbbm Z}_q$-Hilbert}
space.
\item
There is the bounded linear operator
$\varphi\colon {\mathbbm Z}_q\{\{T_1,\ldots,T_N\}\} \to {\mathbbm Z}_q\{\{T_1,\ldots,T_N\}\}$,
provided by taking the quasi power
$h \mapsto h^{(q)}$.
Its~adjoint
$$\kappa\colon {\mathbbm Z}_q\{\{T_1,\ldots,T_N\}\} \longrightarrow {\mathbbm Z}_q\{\{T_1,\ldots,T_N\}\}$$
maps the monomial
$(T_1^{e_1}\cdots T_N^{e_N})^q$
to~$T_1^{e_1}\cdots T_N^{e_N}$
and all monomials not being of this particular type to~zero. In~particular,
$\kappa$
is a retraction for~$\varphi$.
I.e., one has
$\kappa \!\circ\! \varphi = {{\text{\rm id}}}$.

On~the other hand, for
$i \in {\mathbbm N}$,
the composition
$\varphi^i \!\circ\! \kappa^i \colon {\mathbbm Z}_q\{\{T_1,\ldots,T_N\}\} \to {\mathbbm Z}_q\{\{T_1,\ldots,T_N\}\}$
is the orthogonal projection to the
sub-${\mathbbm Z}_q$-Hilbert
space spanned by all monomials of the type
$\smash{(T_1^{e_1}\cdots T_N^{e_N})^{q^i}}$.
\item
For~every polynomial
$g \in {\mathbbm Z}_q[T_1,\ldots,T_N]$,
the multiplication map
$$m_g\colon {\mathbbm Z}_q\{\{T_1,\ldots,T_N\}\} \longrightarrow {\mathbbm Z}_q\{\{T_1,\ldots,T_N\}\} \, , \quad h \mapsto gh \, ,$$
is a bounded linear~operator.\vspace{-2mm}
\end{iii}
\end{nota}

\begin{fac}
\label{cmp}
Let\/~$g \in {\mathbbm Z}_q[T_1,\ldots,T_N]$
and\/~$i \in {\mathbbm N}$
be arbitrary. Then
$$\kappa^i \!\circ\! m_{g^{(q^i)}} \!\circ\! \varphi^i = m_g \, . \ifmmode\rule[-2pt]{0pt}{1pt}\ifinner\tag*{$\square$}\else\eqno{\square}\fi\else\vadjust{}\hfill$\square$\fi\vspace{-5mm}$$
\end{fac}

\begin{defi}[Linear operator associated with a polynomial]
\label{tr_op}
Given~a positive
integer~$i$
and a polynomial
$g \in {\mathbbm Z}_q[T_1,\ldots,T_N]$,
we consider the linear operator
$$\mathop{\text{\rm M}}\nolimits_{g,i} := \kappa^i \!\circ\! m_g \,.$$
\end{defi}

This operator is of interest because of the two results below.\vspace{-3mm}

\begin{lem}
For~every positive integer\/
$i$
and every
polynomial\/~$g \in {\mathbbm Z}_q[T_0,\ldots,T_N]$,
one~has
\begin{equation}
\label{mat_prod}
\mathop{\text{\rm M}}\nolimits_{\mathop{\text{\rm N}}\nolimits_{{\mathbbm Z}_{q^i}/{\mathbbm Z}_q} \!g,\,i} = (\mathop{\text{\rm M}}\nolimits_{g,1})^i \, .
\end{equation}
{\bf Proof.}
{\em
We~will show this inductively. The~assertion is clearly true for
$i=1$.
For~the inductive step, let us assume that equation~(\ref{mat_prod}) is true for a certain value
of~$i$.
Then,~by \ref{quasi}.iii), we have
$\mathop{\text{\rm N}}\nolimits_{{\mathbbm Z}_{q^{i+1}}/{\mathbbm Z}_q} g = g^{(q^i)} \!\cdot\! \mathop{\text{\rm N}}\nolimits_{{\mathbbm Z}_{q^i}/{\mathbbm Z}_q} g$.
Hence,~Definition~\ref{tr_op} shows~that
\begin{eqnarray*}
\mathop{\text{\rm M}}\nolimits_{\mathop{\text{\rm N}}\nolimits_{{\mathbbm Z}_{q^{i+1}}/{\mathbbm Z}_q} \!g,\,i+1} & = & \kappa^{i+1} \!\circ\! m_{\mathop{\text{\rm N}}\nolimits_{{\mathbbm Z}_{q^{i+1}}/{\mathbbm Z}_q} \!g} \\
 & = & \kappa^{i+1} \!\circ\! m_{g^{(q^i)}} \!\circ\! m_{\mathop{\text{\rm N}}\nolimits_{{\mathbbm Z}_{q^i}/{\mathbbm Z}_q} \!g} \\
 & = & \kappa^{i+1} \!\circ\! m_{g^{(q^i)}} \!\circ\! \varphi^i \!\circ\! \kappa^i \!\circ\! m_{\mathop{\text{\rm N}}\nolimits_{{\mathbbm Z}_{q^i}/{\mathbbm Z}_q} \!g} \, .
\end{eqnarray*}
Observe~that the projection
$\varphi^i \!\circ\! \kappa^i$
has no effect here. Indeed,~it is followed by the multiplication with a polynomial, all of whose monomials have exponents only divisible
by~$q^i$,
and the
re\-trac\-tion~$\kappa^i$.
Thus,~according to Fact~\ref{cmp}, Definition~\ref{tr_op}, and the induction hypothesis,
\begin{eqnarray*}
\mathop{\text{\rm M}}\nolimits_{\mathop{\text{\rm N}}\nolimits_{{\mathbbm Z}_{q^{i+1}}/{\mathbbm Z}_q} \!g,\,i+1} & = & \kappa \!\circ\! m_g \!\circ\! \mathop{\text{\rm M}}\nolimits_{\mathop{\text{\rm N}}\nolimits_{{\mathbbm Z}_{q^i}/{\mathbbm Z}_q} \!g,\,i} \\
 & = & \mathop{\text{\rm M}}\nolimits_{g,1} \!\circ\, (\mathop{\text{\rm M}}\nolimits_{g,1})^i \\[1.4mm]
 & = & (\mathop{\text{\rm M}}\nolimits_{g,1})^{i+1} \, ,
\end{eqnarray*}
as required.
\ifmmode\rule[-2pt]{0pt}{1pt}\ifinner\tag*{$\square$}\else\eqno{\square}\fi\else\vadjust{}\hfill$\square$\fi
}
\end{lem}

\begin{prop}
Let\/~$i$
be a positive integer and\/
$g \in {\mathbbm Z}_q[T_0,\ldots,T_N]$
be any~polynomial.

\begin{abc}
\item
Then the linear operator\/
$\mathop{\text{\rm M}}\nolimits_{g,i}$
is trace~class.
\item
Its trace\/
$\mathop{\text{\rm tr}}\nolimits \mathop{\text{\rm M}}\nolimits_{g,i}$
is equal to the sum of all coefficients
of\/~$g$
at the monomials of the form\/
$(T_1^{e_1} \cdots T_N^{e_N})^{q^i-1}$.
\end{abc}\smallskip

\noindent
{\bf Proof.}
{\em
We~will only prove b), as this immediately implies~a). By~definition,
$$\mathop{\text{\rm tr}}\nolimits \mathop{\text{\rm M}}\nolimits_{g,i} = \sum_{k\in K} \langle \mathop{\text{\rm M}}\nolimits_{g,i}(v_k), v_k \rangle \, ,$$
for
$(v_k)_{k\in K}$
any complete orthonormal~system. In~particular, for the system of all~monomials.

For~the monomial
$T_1^{e_1} \cdots T_N^{e_N}$,
the pairing
$\langle \mathop{\text{\rm M}}\nolimits_{g,i}(T_1^{e_1} \cdots T_N^{e_N}), T_1^{e_1} \cdots T_N^{e_N} \rangle$
yields the coefficient of
$\mathop{\text{\rm M}}\nolimits_{g,i}(T_1^{e_1} \cdots T_N^{e_N}) = \kappa^i (gT_1^{e_1} \cdots T_N^{e_N})$
at the monomial
$T_1^{e_1} \cdots T_N^{e_N}$.
This~is nothing but the coefficient of
$gT_1^{e_1} \cdots T_N^{e_N}$
at
$\smash{(T_1^{e_1} \cdots T_N^{e_N})^{q^i}}$
and, hence, the coefficient
of~$g$
at
$\smash{(T_1^{e_1} \cdots T_N^{e_N})^{q^i-1}}$.
The assertion~follows.
}
\ifmmode\rule[-2pt]{0pt}{1pt}\ifinner\tag*{$\square$}\else\eqno{\square}\fi\else\vadjust{}\hfill$\square$\fi
\end{prop}

\begin{rem}
For~general information on
$p$-adic
Hilbert spaces, linear operators, and the concept of the trace, the interested reader might compare Chapters 2 and~3 of the textbook~\cite{Di} of T.~Diagana.
\end{rem}

{\em An algorithm.}

\begin{algo}
\label{Harvey_alg}
Given a projective hypersurface
$X = V(f) \subset {\text{\bf P}}^N_{{\mathbbm F}_{\!q}}$
and a positive
integer~$n$,
this algorithm computes the numbers of points
$\#X({\mathbbm F}_{\!q})$,
$\#X({\mathbbm F}_{\!q^2})$,
\ldots,~$\#X({\mathbbm F}_{\!q^n})$.

\begin{iii}
\item
Choose a lift
$\smash{\widetilde{F} \in {\mathbbm Z}_q[T_0,\ldots,T_N]}$
of~$f$.
Decompose~$X$
into~$X^c := [V(T_0) \cup\ldots\cup V(T_N)] \cap X$
and~$X^0 := X \!\setminus\! X^c$.
For~the lower dimensional closed subscheme
$X^c$,
either use naive point counting methods or apply the method recursively.
For~$X^0$,
do the following.
\item
Use elementary arguments or the Weil conjectures or both to find estimates for
$\#X^0({\mathbbm F}_{\!q})$,
$\#X^0({\mathbbm F}_{\!q^2})$,
\ldots,~$\#X^0({\mathbbm F}_{\!q^n})$.
From~these, deduce a
\mbox{$p$-adic}
precision~$l$
that suffices to determine these numbers~exactly. Put,~finally,
$a := \lceil l/(q-1) \rceil$.
\item
Define
$F \in {\mathbbm Z}_q[T_1,\ldots,T_N]$
by putting
$\smash{F(T_1,\ldots,T_N) := \widetilde{F}(1,T_1,\ldots,T_N)}$.
\item
Find matrix representations of the operators
$\mathop{\text{\rm M}}\nolimits_{F^{a(q-1)},1}$,
$\mathop{\text{\rm M}}\nolimits_{F^{2a(q-1)},1}$,
\ldots,
$\mathop{\text{\rm M}}\nolimits_{F^{la(q-1)},1}$.
Work~on finite dimensional subspaces that are large enough to get the traces~right. For~a polynomial of total
degree~$d$,
as a basis one may use all monomials of total
degree~$\leq\!\! \lceil d/(q-1) \rceil$.
\item
For~$i := 1$
to~$n$,
do the following.

Calculate the matrix powers
$\mathop{\text{\rm M}}\nolimits_{F^{a(q-1)},1}^i$,
$\mathop{\text{\rm M}}\nolimits_{F^{2a(q-1)},1}^i$,
\ldots,
$\mathop{\text{\rm M}}\nolimits_{F^{la(q-1)},1}^i$
and eventually their traces
$\smash{\mathop{\text{\rm tr}}\nolimits \mathop{\text{\rm M}}\nolimits_{F^{ka(q-1)},1}^i}$.
Then
\begin{equation}
\label{point_count}
\#X^0({\mathbbm F}_{\!q^i}) \equiv (q^i-1)^N \left[ 1 + \sum\limits_{k=1}^l (-1)^k \binom{l}k \mathop{\text{\rm tr}}\nolimits \mathop{\text{\rm M}}\nolimits_{F^{ka(q-1)},1}^i \right] \pmod {p^l} \, .
\end{equation}
Determine the exact value of
$\#X^0({\mathbbm F}_{\!q^i})$
from this congruence and the estimate established in step~ii).
Add~$\#X^c({\mathbbm F}_{\!q^i})$
and output the integer~found.
\end{iii}
\end{algo}

\begin{rems}

\begin{iii}
\item
Formula~(\ref{point_count}) is essentially the same as the ``trace formula'' in \cite[Theorem~3.1]{Ha}.
\item[ii) (Cohomological interpretation) ]
The matrix
$\mathop{\text{\rm M}}\nolimits_{F^{q-1},1}$
appears for the first time in \cite[\S2]{Dw60a} and later in \cite[Corol\-lar\-ies~1 and~2]{Mi}, where it is shown to be the natural generalisation of the classical Hasse-Witt matrix~\cite{HW}. In~other words, it represents the Frobenius homomorphism on the middle coherent cohomology
$H^{N-1}(X,{\mathscr{O}}_X)$.
The~characteristic polynomial
of~$\mathop{\text{\rm M}}\nolimits_{F^{q-1},1}$
is hence congruent
modulo~$p$
to the Weil polynomial
$\chi_{N-1}$
\cite[Exp.\,XXII, Th\'eor\`eme~3.1]{SGA7}. Unfortunately,~the algorithm above does not provide a canonical lift representing the crystalline Frobenius, but only successive
\mbox{$p$-adic}
approximations of the proper eigenvalues, so that (\ref{point_count}) finally holds due to many~cancellations.
\end{iii}
\end{rems}

\section{A variation adapted to double covers}

{\em The variation.}

\begin{nota}

\begin{iii}
\item
We~assume in this section that
${\mathbbm F}_{\!q}$
is a finite field of
characteristic~$p \neq 2$.
\item
Furthermore,~let
$X$
be a double cover of an affine torus
$\smash{{\mathbbm G}_{m,{\mathbbm F}_{\!q}}^N}$,
given by an equation of the form
$w^2 = f$
for a polynomial
$f \in {\mathbbm F}_{\!q}[T_1,\ldots,T_N]$.
Our~interest stems from double covers of projective space
of~$\smash{{\text{\bf P}}_{{\mathbbm F}_{\!q}}^N}$,
but, as before, going over from projective space to an affine torus is not a serious~modification. One~may at the very end cover
$\smash{{\text{\bf P}}_{{\mathbbm F}_{\!q}}^N}$
by affine~tori.
\item
We~choose a lift
$F \in {\mathbbm Z}_q[T_1,\ldots,T_N]$
of~$f$.
\item
For a positive
integer~$l$
and
$k = 1,\ldots,l$,
let us put
\begin{equation}
\label{ali}
A^{(l)}_k := \frac{(2l)!}{2^{2l-1}l!(l-1)!} \!\cdot\! \frac{(-1)^{k+1}}{2k-1} \binom{l-1}{k-1} \, .
\end{equation}
\end{iii}
\end{nota}

\begin{lem}
\label{Pascal}
Let\/~$l > 0$
and\/~$b$
be~integers.

\begin{abc}
\item
Then~the Pascal-like~matrix
\begin{equation}
\label{Pasc_mat}
P_b^{(l)} :=
\left(
\begin{array}{lccccc}
1 & b & \binom{b}2 & \binom{b}3 &\cdots & \binom{b}{l-1} \\
1 & b+2 & \binom{b+2}2 & \binom{b+2}3 &\cdots & \binom{b+2}{l-1} \\
\vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\
1 & b\!+\!2l\!-\!4 & \binom{b\!+\!2l\!-\!4}2 & \binom{b\!+\!2l\!-\!4}3 &\cdots & \binom{b\!+\!2l\!-\!4}{l-1} \\
1 & b\!+\!2l\!-\!2 & \binom{b\!+\!2l\!-\!2}2 & \binom{b\!+\!2l\!-\!2}3 &\cdots & \binom{b\!+\!2l\!-\!2}{l-1}
\end{array}
\right) \in M_{l\times l}({\mathbbm Q})
\end{equation}
is~invertible.
\item
All coefficients of its inverse\/
$\smash{\big( P_b^{(l)} \big)^{-1}}$
are\/
\mbox{$p$-adic}
integers.
\item
The first row of\/
$\smash{\big( P_1^{(l)} \big)^{-1}}$
is\/
$(A^{(l)}_1, A^{(l)}_2, \ldots, A^{(l)}_l)$.
\end{abc}\smallskip

\noindent
{\bf Proof.}
{\em
a) and~b)
Column transformations convert
$\smash{P_b^{(l)}}$
essentially into a Vandermonde matrix. This~shows
$\smash{\det P_b^{(l)} = \frac{(2l-2)!!(2l-4)!!\ldots 2!!}{1!2!\ldots (l-1)!} = 2^{(l-1)+(l-2)+\cdots+1} = 2^{\binom{l}2} \neq 0}$,
which immediately implies~a). Assertion~b) follows from this together with Cramer's~rule.\smallskip

\noindent
c)
The assertion is that the linear combination with coefficients
$(A^{(l)}_1, A^{(l)}_2, \ldots, A^{(l)}_l)$
of the rows of the matrix
$\smash{P_1^{(l)}}$
yields~$(1, 0, \ldots, 0)$.
Here,~the occurrence of the zeroes is a consequence of the standard relation
$\sum_{k=1}^l (-1)^k \binom{l-1}{k-1} k^e = 0$
among binomial coefficients,
for~$e = 0, \ldots, l-2$.
It~is a trivial, but tedious work to adjust the constant factor.
}
\ifmmode\rule[-2pt]{0pt}{1pt}\ifinner\tag*{$\square$}\else\eqno{\square}\fi\else\vadjust{}\hfill$\square$\fi
\end{lem}

\begin{defi}
We~call a function
$\Phi\colon \mu_{q^i-1}^N \to {\mathbbm Z}_q$
an {\em Eulerian function\/}
for~$X({\mathbbm F}_{\!q^i})$
{\em modulo\/}~$p^l$
if, for every
$(z_1,\ldots,z_N) \in \mu_{q^i-1}^N$,
\begin{align*}
&\Phi(z_1,\ldots,z_N) \equiv \phantom{-}1 \pmod {p^l} \quad\text{ if } f(\pi(z_1),\ldots,\pi(z_N)) \text{ is a square in } {\mathbbm F}_{\!q^i}^* \, , \\
&\Phi(z_1,\ldots,z_N) \equiv -1 \pmod {p^l} \quad\text{ if } f(\pi(z_1),\ldots,\pi(z_N)) \text{ is a non-square in } {\mathbbm F}_{\!q^i}^* \, , \text{ and}\\
&\Phi(z_1,\ldots,z_N) \equiv \phantom{-}0 \pmod {p^l} \quad\text{ if } f(\pi(z_1),\ldots,\pi(z_N)) = 0 \, .
\end{align*}
\end{defi}

\begin{theo}
Suppose that\/
$p$
is an odd prime and that\/
$b$
is odd such that\/
$b(q-1)/2 \geq l$.
Put\/~$A^{(l)}_{b,k} := \big(\big( P_b^{(l)} \big)^{-1}\big)_{1,k}$.
Then
$$\Phi := \sum\limits_{k=1}^l A^{(l)}_{b,k} \!\cdot\! N_{{\mathbbm Z}_{q^i}/{\mathbbm Z}_q}F^{(b+2k-2)(q-1)/2}$$
is an Eulerian function
for\/~$X({\mathbbm F}_{\!q^i})$
modulo\/~$p^l$.\medskip

\noindent
{\bf Proof.}
{\em
Recall~that each of the coefficients
$\smash{A^{(l)}_{b,1} := \big(\big( P_b^{(l)} \big)^{-1}\big)_{1,1}}$,
\ldots,
$\smash{A^{(l)}_{b,l} := \big(\big( P_b^{(l)} \big)^{-1}\big)_{1,l}}$
is a
\mbox{$p$-adic}
integer. Furthermore, one has that
$f(\pi(z_1),\ldots,\pi(z_N)) = 0 \in {\mathbbm F}_{\!q^i}$
if and only if
$\smash{\mathop{\text{\rm N}}\nolimits_{{\mathbbm F}_{\!q^i}/{\mathbbm F}_{\!q}}(f(\pi(z_1),\ldots,\pi(z_N))) = 0 \in {\mathbbm F}_{\!q}}$,
which, in turn, is equivalent to
$\smash{N_{{\mathbbm Z}_{q^i}/{\mathbbm Z}_q}F(z_1, \ldots, z_N)}$
being a non~unit
in~${\mathbbm Z}_q$.
If~this is true then
$\smash{N_{{\mathbbm Z}_{q^i}/{\mathbbm Z}_q}F^{(b+2k-2)(q-1)/2}(z_1,\ldots,z_N)}$
is of
\mbox{$p$-adic}
valuation at least
$(b+2k-2)(q-1)/2 \geq b(q-1)/2 \geq l$.
Consequently,
$$\sum_{k=1}^l A^{(l)}_{b,k} \!\cdot\! N_{{\mathbbm Z}_{q^i}/{\mathbbm Z}_q}F^{(b+2k-2)(q-1)/2} \equiv 0 \pmod {p^l} \, .$$

Otherwise,
$\smash{f(\pi(z_1),\ldots,\pi(z_N)) \in ({\mathbbm F}_{\!q^i}^*)^2}$
precisely when
$\smash{\mathop{\text{\rm N}}\nolimits_{{\mathbbm F}_{\!q^i}/{\mathbbm F}_{\!q}}(f(\pi(z_1),\ldots,\pi(z_N))) \in {\mathbbm F}_{\!q}^*}$
is a square. Then
$$\mathop{\text{\rm N}}\nolimits_{{\mathbbm F}_{\!q^i}/{\mathbbm F}_{\!q}}(f^{(q-1)/2}(\pi(z_1),\ldots,\pi(z_N))) = 1 \in {\mathbbm F}_{\!q} \, ,$$
while the same expression is equal
to~$(-1)$
in the case of a non-square. By~construction, this~means
$\smash{N_{{\mathbbm Z}_{q^i}/{\mathbbm Z}_q}F^{(q-1)/2}(z_1, \ldots, z_N) \equiv \pm1 \pmod p}$,
respectively.

Write
$N_{{\mathbbm Z}_{q^i}/{\mathbbm Z}_q}F^{(q-1)/2}(z_1, \ldots, z_N) = \pm(1 + sp)$
for some
unknown~$s \in {\mathbbm Z}_{q^i}$.
This~yields
\begin{eqnarray*}
N_{{\mathbbm Z}_{q^i}/{\mathbbm Z}_q}F^{3(q-1)/2}(z_1,\ldots,z_N) & = & \pm\big(1 + 3sp + 3s^2p^2 + s^3p^3\big) \, ,\\
N_{{\mathbbm Z}_{q^i}/{\mathbbm Z}_q}F^{5(q-1)/2}(z_1,\ldots,z_N) & = & \pm\big(1 + 5sp + 10s^2p^2 + 10s^3p^3 + 5s^4p^4 + s^5p^5\big) \, , \\[-2mm]
&\vdots& \\
N_{{\mathbbm Z}_{q^i}/{\mathbbm Z}_q}F^{(b+2l-2)(q-1)/2}(z_1,\ldots,z_N) & = & \textstyle \pm\big(1 + (b\!+\!2l\!-\!2)sp + \binom{b+2l-2}2 s^2p^2 + \\[-1mm]
 & & \textstyle \hspace{1.45cm} \binom{b+2l-2}3 s^3p^3 + \cdots + \binom{b+2l-2}{l-1} s^{l-1}p^{l-1} + \ldots \big) \, .
\end{eqnarray*}
Lemma~\ref{Pascal}.a) and~b) now show
$\smash{\sum\limits_{k=1}^l A^{(l)}_{b,k} \!\cdot\! N_{{\mathbbm Z}_{q^i}/{\mathbbm Z}_q}F^{(b+2k-2)(q-1)/2}(z_1,\ldots,z_N) \equiv \pm1 \pmod {p^l}}$,
as required.%
}%
\ifmmode\rule[-2pt]{0pt}{1pt}\ifinner\tag*{$\square$}\else\eqno{\square}\fi\else\vadjust{}\hfill$\square$\fi
\end{theo}

\begin{rem}
According to the definition, if
$\Phi$
is an Eulerian function
for~$X({\mathbbm F}_{\!q^i})$
modulo~$p^l$
then
\begin{equation}
\#X({\mathbbm F}_{\!q^i}) \equiv (q^i-1)^N + \!\!\!\!\!\!\!\!\!\!\!\sum_{(z_1,\ldots,z_N) \in \mu_{q^i-1}^N}\!\!\!\!\!\!\!\!\!\!\! \Phi(z_1,\ldots,z_N) \pmod {p^l} \, .
\end{equation}
This leads to the algorithm~below.
\end{rem}

\begin{algo}
Given a polynomial
$f \in {\mathbbm F}_{\!q}[T_0,\ldots,T_N]$
and a positive
integer~$n$,
this algorithm computes the numbers of points
$\#X({\mathbbm F}_{\!q})$,
$\#X({\mathbbm F}_{\!q^2})$,
\ldots,~$\#X({\mathbbm F}_{\!q^n})$
on the double cover
of~$\smash{{\text{\bf P}}^N_{{\mathbbm F}_{\!q}}}$,
given
by~$w^2 = f$.

\begin{iii}
\item
Choose a lift
$\smash{\widetilde{F} \in {\mathbbm Z}_q[T_0,\ldots,T_N]}$
of~$f$.
Decompose~$\smash{{\text{\bf P}}^N_{{\mathbbm F}_{\!q}}}$
into~$P^c := V(T_0) \cup\ldots\cup V(T_N)$
and~$\smash{P^0 := {\text{\bf P}}^N_{{\mathbbm F}_{\!q}} \!\setminus\! P^c}$.
Then~$\smash{P^0 \cong {\mathbbm G}_{m,{\mathbbm F}_{\!q}}^N}$.
For~the
restriction~$X^c$
of the double cover to the lower dimensional closed
subscheme~$P^c$,
either use naive point counting methods or apply the method recursively.
For~$X^0$,
the restriction
to~$P^0$,
do the~following.
\item
Use elementary arguments or the Weil conjectures or both to find estimates for
$\#X^0({\mathbbm F}_{\!q})$,
$\#X^0({\mathbbm F}_{\!q^2})$,
\ldots,~$\#X^0({\mathbbm F}_{\!q^n})$.
From~these, deduce a
\mbox{$p$-adic}
precision~$l$
that suffices to determine these numbers~exactly. Put, finally,
$b := \lfloor 2l/(q-1)\rfloor$.
\item
If~$b > 1$
then determine the coefficient vector
$\smash{(A^{(l)}_{b,1}, \ldots, A^{(l)}_{b,l})}$
by inverting the matrix
$P_b^{(l)}$
from~(\ref{Pasc_mat}). Otherwise,~calculate first
$$C_l := \frac{(2l)!}{2^{2l-1}l!(l-1)!}$$
\mbox{and then the coefficients
$\smash{\!A^{(l)}_k \!:=\! C_l \!\cdot\! \frac{(-1)^{k+1}}{2k-1} \binom{l-1}{k-1}}$.
Finally,
set~$\smash{\!(A^{(l)}_{1,1}, \ldots, A^{(l)}_{1,l}) \!:=\! (A^{(l)}_1, \ldots, A^{(l)}_l)}$.}
\item
Define
$F \in {\mathbbm Z}_q[T_1,\ldots,T_N]$
by putting
$\smash{F(T_1,\ldots,T_N) := \widetilde{F}(1,T_1,\ldots,T_N)}$.
\item
\!Find \!matrix \!representations \!of \!the \!operators
$\mathop{\text{\rm M}}\nolimits_{F^{(q-1)/2},1}$,
$\mathop{\text{\rm M}}\nolimits_{F^{3(q-1/2)},1}$,
\ldots,
$\mathop{\text{\rm M}}\nolimits_{F^{(b+2l-2)(q-1/2)},1}$.
Work~on finite dimensional subspaces that are large enough to get the traces~right. For~a polynomial of total
degree~$d$,
as a basis one may use all monomials of total
degree~$\leq\!\! \lceil d/(q-1) \rceil$.
\item
For~$i := 1$
to~$n$,
do the following.

Calculate the matrix powers
$\mathop{\text{\rm M}}\nolimits_{F^{(q-1)/2},1}^i$,
$\mathop{\text{\rm M}}\nolimits_{F^{3(q-1)/2},1}^i$,
\ldots,
$\mathop{\text{\rm M}}\nolimits_{F^{(b+2l-2)(q-1)/2},1}^i$
and eventually their traces
$\smash{\mathop{\text{\rm tr}}\nolimits \mathop{\text{\rm M}}\nolimits_{F^{(b+2k-2)(q-1)/2},1}^i}$.
Then
$$\#X^0({\mathbbm F}_{\!q^i}) \equiv (q^i-1)^N \left[ 1 + \sum\limits_{k=1}^l A^{(l)}_{b,k} \!\cdot\! \mathop{\text{\rm tr}}\nolimits \mathop{\text{\rm M}}\nolimits_{F^{(b+2k-2)(q-1)/2},1}^i \right] \pmod {p^l} \, .$$
Determine the exact value of
$\#X^0({\mathbbm F}_{\!q^i})$
from this congruence and the estimate established in step~ii).
Add~$\#X^c({\mathbbm F}_{\!q^i})$
and output the integer~found.
\end{iii}
\end{algo}

\begin{rems}

\begin{iii}
\item (Complexity.)
The algorithm contains two time consuming operations. The~first is the computation of the powers
of~$F$
and the second is the multiplication of~matrices. We~tested our implementation in {\tt magma} on one core of an Intel i5-4690 processor running at~3.5\,GHz.

For~a
$K3$~surface
of degree two, in order to compute the whole Weil polynomial, the algorithm has to work with a
\mbox{$p$-adic}
precision of at least 11~digits. For~this, it has to compute
$F^{k(p-1)/2}$
for all odd integers
$1 \leq k \leq 21$.
This~results in matrices of size up to
$2080 \times 2080$
having entries
in~${\mathbbm Z}/p^{11}{\mathbbm Z}$.

The table below lists the time used for a randomly chosen surface in a test modulo four representative primes of good~reduction. The computation required about 13\,GB of memory.
\begin{table}[H]
\small
\centerline{
\begin{tabular}{c|c|c|c}
$p$ & powers of $F$ & matrix build & matrix operations \\
\hline
 31  &  12.65 & 21.05 & 55.36 \\
 61  &  76.91 & 21.63 & 71.00 \\
 97  & 236.92 & 22.30 & 73.53 \\
127  & 489.92 & 22.36 & 73.97
\end{tabular}}
\caption{Time (in s) for the steps of the algorithm}
\end{table}\vspace{-3mm}

The~figures show that the time for the computation of the powers
of~$F$
is  approximately cubic
in~$p$.
This~is explained by the fact that the number of terms of the final result is quadratic
in~$p$,
while the number of multiplications necessary to approach it is linear
in~$p$.
The time for the matrix operations is almost constant as the size of the
matrices is constant. The variation of the time is explained completely by the
change of the ring the matrix entries are taken~from.

In~theory, there exists an algorithm that is quadratic
in~$p$.
This~is, however, not yet implemented in {\tt magma} and, reportedly, it is not easy to~implement.
\item
We~incorporated the variable
$b$
only for completeness. In~practice, the algorithm is of interest particularly when
$q$
is not very~small. Then
$q-1 \geq 2l$
and~$b = 1$.
The~same applies to the variable
$a$
in Algorithm~\ref{Harvey_alg}.
\item
One might as well determine the coefficient vector
$\smash{(A^{(l)}_{b,1}, \ldots, A^{(l)}_{b,l})}$
by a matrix inversion, also when
$b=1$.
This~step is in fact not time-critical. Nevertheless, we find it interesting that it is possible to give an explicit~formula.
\end{iii}
\end{rems}

\begin{rem}[The moving simplex idea]
The~idea behind most of the operator calculus above is that only very few of the coefficients of the powers
of~$F$
are actually~needed. An~extreme case is the~following. Assume~that a
$K3$~surface
$X$
is given as a double cover
of~${\text{\bf P}}^2_{\mathbbm Q}$,
by the equation
$W^2 = F(T_0,T_1,T_2)$
for~$F$
a homogeneous polynomial of
degree~$6$.
Moreover,~one only wants to decide for which primes
$p$
up to a certain
bound~$B$
the reduction
$X_p$
is ordinary. This~means that only
$(\#X_p({\mathbbm F}_{\!p}) \bmod p)$
is asked~for.
Thus,~for each prime, just a single coefficient is to be computed, that of
$\smash{F^{\frac{p-1}2}}$
at
$T_0^{p-1} T_1^{p-1} T_2^{p-1}$.

\looseness-1
This~might be done as~follows. Consider~just a small triangle of coefficients
of~$F^e$,
those at
$T_0^{a_0} T_1^{a_1} T_2^{a_2}$
such that
$\smash{|a_i-\frac{a_0+a_1+a_2}3| \leq c_i}$
for
$i = 0, 1, 2$
and some constants
$c_0$,
$c_1$,
and~$c_2$.
These triangles may be calculated inductively for
$e = 1, \ldots, (B-1)/2$.
For~$F$
homogeneous of
degree~$d$,
any choice such that
$c_0+c_1+c_2 := 2(d-1)$
should~work.

\looseness-1
Indeed,~$F$
is given and therefore
$\frac{dF}{dT_0}$,
$\frac{dF}{dT_1}$,
and~$\frac{dF}{dT_2}$
are~known. When~the triangle for
$F^e$
has been established, computing
$F^{e+1}$
naively one would lose some of the~coefficients. But~the coefficients of
$F^{e+1}$
have massive linear~relations. A~single coefficient of
$F^{e+1}$
appears as a linear combination of the coefficients
of~$F^e$,
when evaluating the formula
$F^{e+1} = F \!\cdot\! F^e$.
It~appears as well in
$\smash{\frac{dF^{e+1}}{dT_0} \!=\! (e+1) F^e \frac{dF}{dT_0}}$
and the analogous formulae involving
$T_1$
and~$T_2$.
Hence,~each hypothetical coefficient yields three linear~relations. One~may thus reconstruct some neighbouring coefficients using the linear relations and therefore keep the size of the triangle constant under~iteration. This~approach stems from \cite[\S4.1]{Ha} and is implemented in the function {\tt NonOrdinaryPrimes}.

\looseness-1
On~the other hand, {\tt WeilPolynomialOfDegree2K3Surface} does not use the moving simplex idea, but the approach described above in this~section. The~point is that, for example for
$p = 127$,
one has to compute
$F^{63 \cdot 21}$,
which results in a sum of 31\,517\,830~terms. Thereby,~from the coefficients, about 14\% are used as matrix~entries. Thus, computing only the used coefficients of the powers
of~$F$
would speed up the algorithm only for considerably larger values
of~$p$.\vspace{-3mm}
\end{rem}

\begin{rem}
There are improvements to the moving simplex idea, which we did not implement. The point is that instead of keeping track of the coefficients
of~$F^e$
on the entire simplex, one may sometimes, and probably often, get by with working only on a relatively small~subset. For~example, the idea has been worked out for smooth space quartics, in which case 64 coefficients suffice out of the 220, one would naively~have. This~is called the controlled reduction Abbott--Kedlaya--Roe algorithm \cite[Algorithm~3.4.10]{AKR}. Cf.~\cite{Co}.
\end{rem}

\section{Families of
$K3$~surfaces
that are highly likely to have real or complex multiplication}
\label{RMCM}

\begin{ttt}
\label{arith_cons}
Based~on the Mumford-Tate conjecture, which was proven for
$K3$~surfaces
by S.\,G.~Tankeev~\cite{Ta90,Ta95}, we established in~\cite{EJ14} certain arithmetic consequences of real or complex multiplication for
$K3$~surfaces
defined
over~${\mathbbm Q}$.
Assume~for simplicity that the endomorphism field
$E$
is a quadratic number~field. Then~these consequences include the~following.

\begin{iii}
\item
If~$p$
is inert
in~$E$
then
$\#X_p({\mathbbm F}_{\!p}) \equiv 1 \pmod p$.
I.e.~$X_p$
is non-ordinary.
\item
For every prime
$p$
of good reduction, the transcendental factor
$\chi_p^{{\text{\rm tr}}}$
of the Weil polynomial
of~$X_p$
either splits over
$E$
into two factors conjugate to each other or becomes a square under raising all its roots to the
\mbox{$f$-th}
power, for
some~$f > 0$.
\end{iii}
During the last few years, we systematically searched for
$K3$~surfaces,
the reductions of which show such an unusually regular behaviour. So~far, we have the following conjectural list of suspicious~surfaces.
\end{ttt}

\begin{cons}
\label{list_ex}
Consider~the following\/
$K3$~surfaces
and families of~such,
\begin{align*}
V_a^{(2)} \colon
    & W^2 = \textstyle
 [(\frac18 a^2 \!-\! \frac12 a \!+\! \frac14)T_1^2 + (a^2 \!-\! 2a \!+\! 2)T_1T_2 + (a^2 \!-\! 4a \!+\! 2)T_2^2] \\[-.5mm]
    &  \textstyle \hspace{1.7mm}
 [(\frac18 a^2 \!+\! \frac12 a \!+\! \frac14)T_0^2 + (a^2 \!+\! 2a \!+\! 2)T_0T_2 + (a^2 \!+\! 4a \!+\! 2)T_2^2]
 [2T_0^2 + (a^2 \!+\! 2)T_0T_1 + a^2T_1^2] , \\
V_a^{(5)} \colon
    & W^2 = \textstyle
 [T_1^2 + aT_1T_2 + (\frac5{16} a^2 \!+\! \frac54 a \!+\! \frac54)T_2^2]
 [T_0^2 + T_0T_2 + (\frac1{320} a^2 \!+\! \frac1{16} a \!+\! \frac5{16}) T_2^2] \\[-.5mm]
    &   \textstyle\hspace{97.2mm}
 [T_0^2 + T_0T_1 + \frac1{20} T_1^2] , \\
V^{(13)} \colon
    & W^2 = (25T_1^2 + 26T_1T_2 + 13T_2^2)(T_0^2 + 2T_0T_2 + 13T_2^2)(9T_0^2 + 26T_0T_1 + 13T_1^2) \, , \\
V_{a,b}^{(2)} \colon & W^2 = T_0 T_1 T_2 f_{a,b} \quad\text{\rm for} \\[-.5mm]
    &   \hspace{1cm} f_{a,b} \!:= a (-T_0^2 T_1 \!+\! T_0^2 T_2 \!+\! 2 T_0 T_1^2 \!-\! 3 T_0 T_1 T_2 \!+\! T_0 T_2^2 \!+\! T_1^3 \!-\! 4 T_1^2 T_2 \!+\! 5 T_1 T_2^2 \!-\! 2 T_2^3) \\[-.5mm]
    &   \hspace{2cm} {}+ b (T_0^3 - 2 T_0^2 T_1 - T_0 T_1^2 + 3 T_0 T_1 T_2 - T_0 T_2^2 + T_1^2 T_2 - T_1 T_2^2) \, , \\
V_{a,b}^{(3)} \colon & W^2 = T_0 T_1 T_2 f_{a,b} \quad\text{\rm for} \\[-.5mm]
    &   \hspace{1cm} f_{a,b} := (T_0^3 - 2 T_0^2 T_2 - T_0 T_1^2 - T_0 T_2^2 - 2 T_1^2 T_2 + 2 T_2^3) a^2  \\[-.5mm]
    &   \hspace{2cm} {}+ (6 T_0^2 T_1 + 6 T_0^2 T_2 + 6 T_0 T_1^2 + 6 T_1^2 T_2 - 6 T_2^3) a b \\[-.5mm]
    &   \hspace{2cm} {}+ (-3 T_0^2 T_1 - 6 T_0^2 T_2 + 3 T_1^3 - 6 T_1^2 T_2 - 3 T_1 T_2^2 + 6 T_2^3) b^2 , \\
V_{{\bf a},{\bf b}}^{(-1)} \colon
    & W^2 = T_0 T_1 T_2 (T_0+T_1+T_2)(a_1 T_0 + a_2 T_1 + a_3 T_2)(b_1 T_0 + b_2 T_1 + b_3 T_2) \\[-.8mm]
    &   \hspace{1.3cm} \text{\rm for } {\bf a}, {\bf b} \in \smash{{\mathbbm C}^3} \text{\rm ~such that } a_1 b_3 \!+\! a_2 b_1 \!-\! 2 a_3 b_1 = 0 \text{\rm ~and } a_1 b_2 \!+\! a_2 b_3 \!-\! 2 a_3 b_2 = 0 \, . \\
V^{(-1,\mu_7)} \colon
    & W^2 = T_0T_1T_2(7 T_0^3 - 7 T_0^2 T_1 + 49 T_0^2 T_2 - 21 T_0 T_1 T_2 + 98 T_0 T_2^2 + T_1^3 - 7 T_1^2 T_2 \!+\! 49 T_2^3) \, , \\
V^{(-1,\mu_9)} \colon
    & W^2 = T_0 T_1 T_2 (T_0^3 - 3 T_0^2 T_2 - 3 T_0 T_1^2 - 3 T_0 T_1 T_2 + T_1^3 + 9 T_1^2 T_2 + 6 T_1 T_2^2 + T_2^3) \, , \\
V^{\!(-1,\mu_{19})} \colon
    & W^2 = T_0 T_1 T_2 (49T_0^3 - 304T_0^2T_1 + 570T_0^2T_2 + 361T_0T_1^2 - 2793T_0T_1T_2 + 2033T_0T_2^2 \\[-.5mm]
    & \hspace{57mm} {} + 361T_1^3 + 2888T_1^2T_2 - 5415T_1T_2^2 + 2299T_2^3) \, .
\end{align*}%
\begin{abc}
\item
Then the generic fibre of each family, as well as each of the four individual surfaces has geometric Picard
rank\/~$16$.
\item
Moreover,~the generic fibre
of\/
$\smash{V_a^{(2)}}$
has real multiplication
by\/
${\mathbbm Q}(\sqrt{2})$,
that
of\/
$\smash{V_a^{(5)}}$
RM
by\/
${\mathbbm Q}(\sqrt{5})$,
that
of\/
$\smash{V_{a,b}^{(2)}}$
RM
by\/
${\mathbbm Q}(\sqrt{2})$,
that
of\/
$\smash{V_{a,b}^{(3)}}$
RM
by\/
${\mathbbm Q}(\sqrt{3})$,
and that
of\/
$\smash{V_{{\bf a},{\bf b}}^{(-1)}}$
complex multiplication
by\/~${\mathbbm Q}(\sqrt{-1})$.

Finally,~$\smash{V^{(13)}}$
has RM
by\/~${\mathbbm Q}(\sqrt{13})$,
$\smash{V^{(-1, \mu_7)}}$
CM
by\/~${\mathbbm Q}(\zeta_{28} + \zeta_{28}^{13}) = {\mathbbm Q}(i, \zeta_7 + \zeta_7^{-1})$,
$\smash{V^{(-1, \mu_9)}}$
CM
by\/~${\mathbbm Q}(\zeta_{36} + \zeta_{36}^{17}) = {\mathbbm Q}(i, \zeta_9 + \zeta_9^{-1})$,
and\/
$\smash{V^{(-1, \mu_{19})}}$
CM
by\/
$L(i)$,
for
$L \subset {\mathbbm Q}(\mu_{19})$
the unique cubic~subfield.
\end{abc}
\end{cons}

\begin{rems}

\begin{iii}
\item
The~equations given actually describe singular models of the
$K3$~surfaces
to be~considered.
\item
(Evidence.)
Thanks~to Harvey's
\mbox{$p$-adic}
point counting method, we can now give a lot more numerical evidence for our conjectures than before~\cite{EJ14}.
\item
(Proven cases.)
The first three examples were already published in~\cite{EJ14}. For the
family~$\smash{V_a^{(2)}}$,
both a) and~b) were proven in \cite[Theorem~6.6]{EJ14}. The~same method applies to
$\smash{V_{{\bf a},{\bf b}}^{(-1)}}$
and provides a proof also for this~family. As~far as only~a) is concerned, using van Luijk's method one shows in each case that the geometric Picard rank is
$16$
or~$17$,
and~$16$
as soon as~b) is~true.
\item
In~particular, for each family or surface, we explicitly know 16 divisors that are linearly independent in the Picard~group. In~neither case, all of them are defined
over~${\mathbbm Q}$.
For~instance,
for~$\smash{V^{(-1, \mu_7)}}$,
$\smash{V^{(-1, \mu_9)}}$,
and
$\smash{V^{(-1, \mu_{19})}}$,
the fields of definition are
$\smash{{\mathbbm Q}(\zeta_7 + \zeta_7^{-1})}$,
$\smash{{\mathbbm Q}(\zeta_9 + \zeta_9^{-1})}$,
and the unique cubic subfield of
$\smash{{\mathbbm Q}(\mu_{19})}$,
respectively.

On~the other hand, for
$\smash{V_a^{(2)}}$,
the field of definition of the known divisors
is~${\mathbbm Q}(\sqrt{2})$,
while it
is~${\mathbbm Q}(\sqrt{5}, \sqrt{a^2-20a-20})$
for
$\smash{V_a^{(5)}}$,
${\mathbbm Q}(\sqrt{13}, \sqrt{-3})$
for
$V^{(13)}$,
$\smash{{\mathbbm Q}(\sqrt{2}, \sqrt{a}, \sqrt{b}, \sqrt{a^2-6ab+b^2})}$
for~$\smash{V_{a,b}^{(2)}}$
and
$\smash{{\mathbbm Q}(\sqrt{3}, \sqrt{a^4 - 12a^3b + 30a^2b^2 - 36ab^3 + 9b^4}, \sqrt{2a^2 - 6ab + 6b^2})}$
for
$\smash{V_{a,b}^{(3)}}$.
We~do not know whether the following is more than a~coincidence.
\end{iii}
\end{rems}

\begin{eo}
All the examples listed in Conjectures~\ref{list_ex} that are supposed to have real multiplication have the property that the endomorphism field is contained in the field of definition of the Picard~group.
\end{eo}

{\em Numerical evidence. Data for the real multiplication examples.}\medskip

For each prime
$p \in \{19, \ldots, 499\}$,
we inspected all the possible specialisations of the families after reduction
modulo~$p$.
Whenever~this resulted in a non-singular surface, we computed the Weil polynomial and derived the geometric Picard~rank. This~confirmed properties \ref{arith_cons}.i) and~ii) in every~case. In~particular, only the Picard~ranks 18 and 22 are~occurring. There~is the following~statistics.

\begin{table}[H]
\small
\centerline{
\begin{tabular}{p{1.3cm}|c|c|p{1.4cm}|c|c|c}
Family & \multicolumn{6}{c}{\# relative frequency (in\,\%) of rank 22 per prime} \\[-.7mm]
       & \multicolumn{3}{c}{inert primes} & \multicolumn{3}{c}{split primes} \\[-.7mm]
       & $\min$ & average & $\max$ & $\min$ & average & $\max$ \\[-.7mm]\hline
$V_a^{(2)}$     & 0.00 &  7.42 & 25.00 & 0.00 & 6.14 & 25.00 \\
$V_a^{(5)}$     & 2.33 &  9.32 & 24.24 & 0.00 & 5.84 & 16.00 \\
$V_{a,b}^{(2)}$ & 4.75 & 12.66 & 38.89 & 1.50 & 7.39 & 22.73 \\[.85mm]
$\smash{V_{a,b}^{(3)}}$ & 0.00 &  0.00 & 0.000 & 2.33 & 7.40 & 20.59 \\
\end{tabular}}
\caption{Frequency of reduction to geometric Picard rank 22}
\end{table}

One~might want to look more closely at the primes that result in reduction to geometric Picard
rank~$22$.
In~the case of the
family~$\smash{V_a^{(2)}}$,
77 of 88 primes are~occurring. The~missing ones are
$23$,
$29$,
$31$,
$47$,
$97$,
$127$,
$193$,
$241$,
$401$,
$433$,
and~$449$.
The reductions of these primes
modulo~$8$
are
$1^6, 5, 7^4$.
The~only one being inert in the endomorphism field
${\mathbbm Q}(\sqrt{2})$
is~$29$.
For~the family
$\smash{V_{a,b}^{(2)}}$,
all 88 inspected primes do occur, while for
$\smash{V_a^{(5)}}$
87 of the 88 primes occur and only
$29$
is~missing.
However,~in the example
$\smash{V_{a,b}^{(3)}}$,
only 44 of the 88 primes result in reductions to rank 22 and in fact {\em exactly those splitting\/} in the endomorphism field
${\mathbbm Q}(\sqrt{3})$,
a phenomenon we have no explanation~for.\medskip

{\em Primes of reduction to rank 18.}
In~the case of reduction to Picard
rank~$18$,
the transcendental factor split off the Weil polynomial is of degree~four. In~our families, it was never a perfect power, cf.~\cite[Theorem~1.1]{Za93}, but always turned out to be~irreducible. It~was either the norm of a quadratic polynomial
$g \in E[t]$
over the endomorphism field or turned into a square under the operation of squaring all its~roots. The~first alternative occurred precisely at the primes split
in~$E$,
while the second option came up at the inert~primes.\bigskip

{\em Numerical evidence. Data for the complex multiplication examples.}\medskip

Let~us report only on the numerical evidence for the three isolated~examples.
These~are contained in the family
$\smash{V_{{\bf a},{\bf b}}^{(-1)}}$.
Thus,~at least CM
by~${\mathbbm Q}(\sqrt{-1})$
is~proven.
As the surfaces have geometric Picard
rank~$16$,
the transcendental part of the cohomology is of relative dimension~one. Thus,~one should expect phenomena that are very close to those known to occur for elliptic curves with complex~multiplication.

To~verify this, we run the point counting for all the primes
below~$1000$.
It~turns out that all the reductions have geometric Picard rank either
$16$
or~$22$,
where
rank~$16$
appears exactly for the primes
$1 \bmod 4$.
A~prime is ordinary if and only if it completely splits in the endomorphism~field.
Furthermore, the Frobenius eigenvalues on the transcendental part of
$\smash{H^2_{{\text{\rm {\'e}t}}}(V_{\overline{\mathbbm F}_{\!p}}, {\mathbbm Q}_l)}$
are closely related to the endomorphism
field~$E$.
More~precisely,

\begin{abc}
\item
In the case that
$p$
completely splits
in~$E$,
the eigenvalues
of~$\mathop{\text{\rm Frob}}\nolimits$
are contained
in~$E$.
They~are of the~form
$\pm p \frac{\overline{\pi}}{\pi}$
for~$(\pi)$
a prime
above~$p$.
Note that the endomorphism fields have class
number~$1$,
such that the ideals
above~$p$
are~principal.
\item
In the case that
$p$
splits into three primes, the characteristic polynomial
of~$\mathop{\text{\rm Frob}}\nolimits$
on the transcendental lattice is
$t^6 - 3 p^2 t^4 + 3 p^4 t^2 - p^6 = (t^2 - p^2)^3$.
\item
In the case that
$p$
splits into two primes, the characteristic polynomial
of~$\mathop{\text{\rm Frob}}\nolimits$
on the transcendental lattice is of the form
$t^6 + a p^2 t^3 + p^6$,
for~$a \in {\mathbbm Z}$
having the property that the discriminant of
$t^2 + a t + p^2$
is minus a~square. Thus, the Weil polynomial
of~$\smash{X_{{\mathbbm F}_{\!p^3}}}$
is a cube of a quadratic polynomial having its roots in
${\mathbbm Q}(i)$.
These~are of the form
$p^2 \varrho^2$,
for
$\varrho \in {\mathbbm Q}(i)$
a prime
above~$p$.
Let~us note here once again that
${\mathbbm Q}(i) \subset E$.
\item
In the case that
$p$
is totally inert, the degree six transcendental factor of the Weil polynomial is equal
to~$t^6 - p^6$.
\end{abc}
The total computation took about 110 days of CPU-time on one
Intel Xeon E5-4650 processor running at 2.7\,GHz.

\frenchspacing


\section{Introduction}

The Sum-Product Network (SPN) \citep{Poon} is a tractable and interpretable deep model. An advantage of SPNs over other graphical models such as Bayesian Networks is that they allow efficient exact inference in linear time with network size.
An SPN represents a multivariate probability distribution with a directed acyclic graph consisting of sum nodes (clusters over instances), product nodes (partitions over features), and leaf nodes (univariate distributions over features), as shown in Figure \ref{fig:spn}.

The standard algorithms for learning SPN structure assume discrete data with no missingness, and mostly test on the same set of benchmark data sets that satisfy these criteria \citep{Rooshenas}. This is not a reasonable assumption when dealing with messy data sets in real applications. The Google Knowledge Graph (KG) is a semantic network of facts, based on Freebase \citep{Freebase}, used to generate Knowledge Panels in Google Search. KG data is quite heterogeneous, with a lot of it missing, since much more is known about some entities in the graph than others. High missingness rates can also worsen the impact of discretizing continuous variables before doing structure learning, which results in losing more of the already scarce covariance information.

Applications like the KG are common, and there is a need for an SPN learning algorithm that can handle this kind of data. In this paper, we present MiniSPN, a simplification of a state-of-the-art SPN learning algorithm that improves its applicability, performance and speed. We demonstrate the performance of MiniSPN on KG data and on standard benchmark data sets.

\begin{figure}[t]
\centering
\begin{minipage}{.5\textwidth}
\includegraphics[scale=0.6]{spn}
\caption{Example of an SPN structure \\(figure from \cite{Zhao})}\label{fig:spn}
\end{minipage
\begin{minipage}{.5\textwidth}
\includegraphics[scale=0.17]{learnspn}
\caption{Recursive partitioning process in the LearnSPN algorithm \\(figure from \cite{Gens})}\label{fig:learnspn}
\end{minipage}
\end{figure}

\section{Variation on the LearnSPN algorithm}

LearnSPN \citep{Gens} is a greedy algorithm that performs co-clustering by recursively partitioning variables into approximately independent sets and partitioning the training data into clusters of similar instances, as shown in Figure \ref{fig:learnspn}. The variable and instance partitioning steps are applied to data slices (subsets of instances and variables) produced by previous steps.
The variable partition step uses pairwise independence tests on the variables, and the approximately independent sets are the connected components in the resulting dependency graph.
The instance clustering step uses a naive Bayes mixture model for the clusters, where the variables in each cluster are assumed independent. The clusters are learned using hard EM with restarts, avoiding overfitting using an exponential prior on the number of clusters.
The splitting process continues until the data slice has too few instances to test for independence, at which point all the variables in that slice are considered independent. The end result is a tree-structured SPN.

The standard LearnSPN algorithm assumes that all the variables are discrete and there is no missing data. Hyperparameter values for the cluster penalty and the independence test critical value are determined using grid search. The clustering step seems unnecessarily complex, involving a penalty prior, EM restarts, and hyperparameter tuning. It is by far the most complicated part of the algorithm in a way that seems difficult to justify, and likely the most time-consuming due to the restarts and hyperparameter tuning.
We propose a variation on LearnSPN called MiniSPN that handles missing data, performs lazy discretization of continuous data in variable partition step, simplifies the model in the instance clustering step, and does not require hyperparameter search.

\begin{table}[t]
\centering
\small
\caption{Average log likelihood and runtime comparison on KG data sets (best performing methods are shown in bold).}
\label{kg-results}
\begin{tabular}{l|lll|lll}\hline
& \multicolumn{3}{c}{Test set log likelihood} & \multicolumn{3}{|c}{Runtime (seconds)} \\\hline
Data set         & Pareto  & Hybrid          & MiniSPN         & Pareto   & Hybrid   & MiniSPN         \\\hline
Professions-10K  & -10.2   & -6.2            & \textbf{-6.09}  & 5.3      & 3.7      & \textbf{0.4}    \\
Professions-100K & -6.61   & -6.53           & \textbf{-6.44}  & 72       & 131      & \textbf{7.2}    \\
Dates-10K        & -8.66   & \textbf{-8.53}  & -8.68           & 1.7      & 2.4      & \textbf{0.26}   \\
Dates-100K       & -17.1   & -16.7           & \textbf{-16.5}  & 29       & 566      & \textbf{5.4}  \\\hline
\end{tabular}
\vspace{1em}
\centering
\caption{Average log likelihood and runtime comparison on literature data sets (best performing methods are shown in bold).}
\label{lit-results}
\small
\begin{tabular}{l|llll|llll}\hline
& \multicolumn{4}{c}{Test set log likelihood} & \multicolumn{4}{|c}{Runtime (seconds)} \\\hline
Data set               & Pareto & Hybrid         & MiniSPN        & LearnSPN        & Pareto       & Hybrid & MiniSPN      & LearnSPN \\\hline
NLTCS        & -6.33  & \textbf{-6.03} & -6.12          & -6.1            & \textbf{4.8} & 35     & \textbf{1.4} & 60       \\
MSNBC        & -6.54  & -6.4           & -6.61          & \textbf{-6.11}  & 61           & 212    & \textbf{5.6} & 2400     \\
KDDCup  & -2.17  & \textbf{-2.13} & \textbf{-2.14} & -2.21           & 152          & 2080   & \textbf{23}  & 400      \\
Plants       & -17.3  & \textbf{-13.1} & \textbf{-13.2} & \textbf{-13}    & 28           & 780    & \textbf{11}  & 160      \\
Audio        & -41.9  & \textbf{-39.9} & \textbf{-40}   & -40.5           & 28           & 556    & \textbf{12}  & 955      \\
Jester       & -54.6  & \textbf{-52.9} & \textbf{-53}   & -53.4           & 13           & 193    & \textbf{6.7} & 1190     \\
Netflix      & -59.5  & \textbf{-56.7} & \textbf{-56.8} & -57.3           & 27           & 766    & \textbf{14}  & 1230     \\
Accidents    & -40.4  & -32.5          & -32.6          & \textbf{-30.3}  & 31           & 1140   & \textbf{18}  & 330      \\
Retail       & -11.1  & \textbf{-11}   & -11.1          & \textbf{-11.09} & 25           & 63     & \textbf{7.3} & 100      \\
Pumsb-star   & -40.8  & -28.4          & -28.3          & \textbf{-25}    & 47           & 1100   & \textbf{22}  & 350      \\
DNA          & -98.1  & -91.5          & -93.9          & \textbf{-89}    & \textbf{6.3} & 45     & \textbf{3}   & 300      \\
Kosarek      & -11.2  & \textbf{-10.8} & \textbf{-10.9} & \textbf{-11}    & 90           & 537    & \textbf{22}  & 200      \\
MSWeb        & -10.7  & \textbf{-9.94} & \textbf{-10.1} & -10.26          & 75           & 572    & \textbf{34}  & 260      \\
Book         & -35.1  & \textbf{-34.7} & \textbf{-34.7} & -36.4           & 83           & 181    & \textbf{32}  & 350      \\
EachMovie    & -55    & \textbf{-52.3} & \textbf{-52.2} & -52.5           & 62           & 218    & \textbf{22}  & 220      \\
WebKB        & -161   & \textbf{-155}  & \textbf{-155}  & -162            & 37           & 169    & \textbf{38}  & 900      \\
Reuters-52   & -92    & \textbf{-85.2} & \textbf{-84.7} & -86.5           & \textbf{76}  & 656    & 95           & 2900     \\
Newsgroup & -156   & \textbf{-152}  & \textbf{-152}  & -160.5          & 181          & 1190   & \textbf{139} & 28000    \\
BBC          & -258   & \textbf{-250}  & \textbf{-249}  & \textbf{-250}   & \textbf{33}  & 123    & \textbf{42}  & 900      \\
Ad           & -52.3  & -49.5          & -49.2          & \textbf{-22}    & \textbf{58}  & 92     & \textbf{50}  & 300        \\\hline
\end{tabular}
\end{table}

We simplify the naive Bayes mixture model in the instance clustering step by attempting a split into two clusters at any given point, and eliminating the cluster penalty prior, which results in a more greedy approach than in LearnSPN that does not require restarts or hyperparameter tuning. This seems like a natural choice of simplification - an extension of the greedy approach used at the top level of the LearnSPN algorithm.
We compare a partition into univariate leaves to a mixture of two partitions into univariate leaves (generated using hard EM), and the split succeeds if the two-cluster version has higher validation set likelihood.
If the split succeeds, we apply it to each of the two resulting data slices, and only move on to a variable partition step after the clustering step fails. The greedy approach is similar to the one used in the SPN-B method \citep{Vergari}, which however alternates between variable and instance splits by default, thus building even deeper SPNs.

In the variable partition step, we perform an independence test using the subset of rows where both variables are not missing, and conclude independence if the number of such rows is below threshold.
We apply binary binning to each continuous variable, using its median in the data slice as a cutoff.

We compare to the ``Pareto" algorithm, previously used for learning SPN models in KG, inspired by the work of \cite{Grosse}.
It produces a Pareto-optimal set of models, trading off between degrees of freedom and validation set log likelihood score. At each iteration, production rules are randomly applied to add partition and mixture splits to the models in the current model set, and the new models are added to the model set. If a model in the model set has both lower degrees of freedom and higher log likelihood score than another model, the inferior model is removed from the set. The algorithm returns the model from the Pareto model set with the highest validation log likelihood. We also compare to a hybrid method, with the Pareto algorithm initialized by MiniSPN.

\section{Summary of experiments}
We use two data sets from the Knowledge Graph People collection. In the KG Professions data set, most of the variables are boolean indicators of whether each person belongs to a particular profession. There are 83 boolean variables and 4 continuous variables. In the KG Dates data set, there are 14 continuous variables representing dates of life events for each person and their spouse(s), with around 95\% of the data missing. We use subsets of 10000 and 100000 instances from each of these data sets, and randomly split the data sets into a training and test set.

On the KG data sets, we compare MiniSPN, Pareto and Hybrid algorithms. We were not able to apply the standard LearnSPN algorithm on these data sets, since they contain missing data. Table \ref{kg-results} shows log likelihood performance on the test set and runtime performance. MiniSPN does better than Pareto, both in terms of log likelihood and runtime. Hybrid performs comparably to MiniSPN, but is usually the slowest of the three.

We use 20 benchmark data sets from the literature (exactly the same ones used in the
LearnSPN paper \citep{Gens}) to compare the performance of MiniSPN with the standard LearnSPN algorithm. We are particularly interested in the effect of MiniSPN's simple two-cluster instance split relative to the more complex instance split with the exponential prior and EM restarts used in the standard LearnSPN. Table \ref{lit-results} shows log likelihood performance on the test set and runtime performance. Like on the KG data, we find that MiniSPN uniformly outperforms Pareto, and performs similarly to Hybrid and LearnSPN but runs much faster (on the most time-intensive data set, Newsgroup, MiniSPN takes 2 minutes while LearnSPN takes 8 hours).

\section{Conclusion}

Sum-product networks have been receiving increasing attention from researchers due to their expressiveness, efficient inference and interpretability, and many learning algorithms have been developed in the past few years. While recent developments have mostly focused on improving performance on benchmark data sets, our variation on a classical learning algorithm is simple yet has a large impact on usability, by improving speed and making it possible to apply to messy real data sets.


\section{Introduction}
\label{sec:intro}
 We have recently undertaken a spectroscopic survey of all apparently slowly rotating bright early A stars (A0-A1V) and late B stars (B8-B9V) observable from the northern hemisphere. The incentive is
to search for rapid rotators seen pole-on or new chemically peculiar B and A stars which have thus far remained unnoticed.
This project addresses fundamental questions of the physics of late-B and early-A stars: i) can we find new instances of rapid rotators seen pole-on (other than Vega) and study their physical properties (gradient of temperature across the disk, limb and gravity darkening), ii) is our census of Chemically Peculiar stars complete up to the magnitude limits we adopted ? If not, what are the  physical properties of the newly found CP stars?.
   The abundance results for the A0-A1V sample have been published in \cite{Royer}.
   The selection criteria were: a declination higher than $-15$\degr, spectral class A0 or A1 and luminosity class V and IV and magnitudes $V$ brighter than 6.65 and a \ensuremath{v_{{\mathrm e}}\sin i} less than 65 km~s$^{-1}$. The B8-9 sample employs the same criteria, except for the V magnitude brighter than 7.85 as these B stars are intrinsically brighter in the V band where SOPHIE reaches its maximum efficiency. Most of the stars of that B8-9 sample (40 stars) have been observed in December 2014 with SOPHIE, the \'echelle high-resolution spectrograph
   at Observatoire de Haute Provence yielding spectra coving the 3900 \AA-6800 \AA\ spectral ranger over 39 orders at a resolving power R = 75000.
   A careful abundance analysis of the high resolution high signal-to-noise ratio spectra of the A stars sample has allowed to sort out the sample of 47 A stars into 17 chemically normal stars (ie. whose abundances do not depart by more  than $\pm 0.20$ dex from solar values), 12 spectroscopic binaries and 13 chemically peculiar stars (CPs) among which five are new   CP stars. The status of these new CP stars still needs to be fully specified by spectropolarimetric observations to address their magnetic nature or by exploring new spectral ranges which we had not explored in this first study. Indeed, the abundance analysis of the A stars sample in \cite{Royer} relied only on four spectral regions: 4150--4300\,\AA, 4400--4790\,\AA, 4920--5850\,\AA, and 6000--6275\,\AA, avoiding Balmer lines and atmospheric telluric lines.
   \\
   We have now started to examine the B9-B8V sample using the full wavelength coverage provided by SOPHIE to search for new Chemically Peculiar stars. We have already reported on the discovery of 4 new HgMn stars \citep{Monier} whose spectra display strong
 \ion{Hg}{2} lines at 3984\,\AA\ and strong \ion{Mn}{2} lines. These new HgMn stars are  HD 18104, HD 30085, HD 32867, HD 53588. In the process of our analysis of the B9-B8V sample, we have just found that HD 67044, currently classified B8 in SIMBAD, is actually another new CP star, most likely another new HgMn star and that HD 42035, classified B9V,  actually is a new spectroscopic binary containing a slow rotator which we will refer to as HD 42035 S (S for "sharp" lines) and a fast rotator, HD 42035 B (B for "broad" lines).
 A bibliographic query of the CDS for HD 67044 and HD 42035 actually reveals only 4 and 26 publications respectively.  Neither HD 67044 nor HD 42035 appear in \citeauthor{Cowley1972}'s classification (\citeyear{Cowley1972}) of the bright B8 stars.
The purpose of this paper is to report on the detection of strong   \ion{Ti}{2},  \ion{Cr}{2}, \ion{Mn}{2}, \ion{Sr}{2}, \ion{Y}{2},
\ion{Zr}{2}, \ion{Ba}{2} lines  and a weak \ion{Hg}{2} line in the spectrum of HD 67044 and perform a line identification of all the very sharp lines of HD 42035 S. We also have determined the abundances of 17 chemical elements for HD 67044 using spectrum synthesis to quantify the enhancements and depletions of these elements. As several of the lines of HD 42035 S are blended with those of HD 42035 B, we have only been able to derive provisional abundances for HD 42035 S using lines little affected by the fast rotator.

\section{Observations and reduction}
\label{sec:obs-reduc}

HD 67044 and HD 42035 have been observed at Observatoire de Haute Provence using the high resolution ($R = 75000$) mode of the SOPHIE \'echelle spectrograph \citep{Perruchot} in December 2014.
  The  $\frac{S}{N}$ ratio of the spectra are about 130 and 174 at 5500 \AA\ respectively. The observations log is displayed in Table\,\ref{table:1}.
 The data are automatically reduced to produce 1D extracted and wavelength calibrated \'echelle orders. Each reduced order was normalised separately using a Chebychev polynomial fit with sigma clipping, rejecting points above or below 1 $\sigma$ of the local continuum. Normalized orders were merged together, corrected  by the blaze function and resampled into a constant wavelength step of about 0.02\,\AA\ \cite[see][for more details]{Royer}.

The radial velocities of HD 67044 and HD 42035 were derived from cross-correlation techniques, avoiding the Balmer lines and the atmospheric telluric lines. The normalized spectrum was cross-correlated with a synthetic template extracted from the POLLUX database\footnote{\url{http://pollux.graal.univ-montp2.fr}}  \citep{Palacios} corresponding to the parameters $T_\mathrm{eff}=11000$\,K, $\log g=4$ and solar abundances. A parabolic fit of the upper part of the resulting cross-correlation function yields the Doppler shift, which is then used to shift spectra to rest wavelengths. The radial velocities of HD 67044 and HD 42035 are collected in Table 2.

\section{The line spectrum of HD 67044}
\label{sec:linespec}

 Several spectral regions have been used to establish the chemical peculiarity of HD 67044. First, the red wing of H${\epsilon}$, which lies in order 3, harbors the \ion{Hg}{2} $\lambda$ 3983.93\,\AA\ line and several \ion{Zr}{2} and \ion{Y}{2} lines likely to be strengthened in CP stars. After proper correction for the stellar radial velocity, we found that HD 67044
does  show the \ion{Hg}{2} 3983.93\,\AA\ line as a feature absorbing about 2\% , and also strong lines of \ion{Y}{2} at 3982.59 \AA\ and of \ion{Zr}{2} at 3991.13 \AA\ and 3998.97 \AA.
The \ion{Y}{2},  \ion{Zr}{2} and \ion{Hg}{2} lines in the spectral range 3980 to 4000 \AA\ are displayed in figure \,\ref{HD67044ZrY} together with the synthetic spectrum fiting best this spectral interval.
Several other lines of \ion{Y}{2} and \ion{Zr}{2} are strong absorbers in the spectrum of HD 67044 and have been used to derive the abundances in these elements.
Second, we examined the region from 4125\,\AA\ to 4145\,\AA\  (order 6) for the \ion{Si}{2} lines at 4128.054 \AA\ and 4130.894 \AA\
and  the  \ion{Mn}{2} line at 4136.92\,\AA. In an HgMn star, the \ion{Mn}{2} line at 4136.92\,\AA\ should be strong whereas it should be absent in any comparison normal late B-type star.
Furthermore, the  lines of \ion{Mn}{2} at 4206.37\,\AA\ and 4252.96\,\AA\ should also be enhanced in the spectra of HgMn stars \citep{Gray}. These three \ion{Mn}{2} lines correspond to moderately strong  features in HD 67044, they usually are broad lines absorbing  5\% of the local continuum. Several other strong  \ion{Mn}{2} lines could be found and have been used to derive the abundance of manganese.

We find that the strongest expected lines of  \ion{Ti}{2}, \ion{Cr}{2}, \ion{Mn}{2} \ion{Sr}{2}, \ion{Y}{2}, \ion{Zr}{2}, \ion{Ba}{2} are all indeed strong features in the spectrum of HD 67044. In Table 3, we give a list of the strongest unblended lines of these species together with abundance determinations.

The presence of the \ion{Hg}{2} 3983.93\,\AA\ line and of several strong \ion{Mn}{2} and \ion{Sr}{2}, \ion{Y}{2} and \ion{Zr}{2}
lines led us to conclude that HD 67044 is actually another new HgMn star
 and should be reclassified as such. The abundance determinations presented in next paragraph do confirm this proposal.

 \begin{table}
\tiny
\caption{Observation log}
\label{table:1}
\centering
\begin{tabular}{cccccc}
\hline\hline
Star ID & Spectral & V &Observation & Exposure & S/N \\
& type&  & Date & time (s ) &\\
\hline
HD 67044 & B 8    & 7.48 & 2014-12-16 & 2400 & 130 \\
HD 42035  & B9V  & 6.55 & 2014-12-18 & 750 &  174 \\
\hline
\end{tabular}
\end{table}

 \begin{figure}[h!]
\vskip 0.5cm
   \centering
      \includegraphics[scale=0.33]{YZrHg_HD67044.eps}
   \caption{Comparison of the observed \ion{Y}{2} line at 3991.13 \,\AA\ , \ion{Hg}{2}  at 3983.93 \AA\ and \ion{Zr}{2} line at 3998.96 \,\AA\ of HD 67044 to a synthetic profile computed for overabundances of 1125 , 200 and 5000 times solar for yttrium, zirconium and mercury respectively.}
   \label{HD67044ZrY}
 \end{figure}

\section{The line spectrum of HD 42035}

The high resolution spectrum of HD 42035 is characterised by very sharp lines originating in HD 42035 S. A first estimate of the projected equatorial rotational velocity
of HD 42035 S
has been obtained by Fourier transform analysis (Fig.\,\ref{FT}), using the two \ion{Fe}{2} lines at 4508.29 and 4515.34\,\AA. From the position of the first zero of the Fourier transform, we derived a value of $\ensuremath{v_{{\mathrm e}}\sin i} = 3.7\pm 0.2$ km\,s$^{-1}$.
A close inspection of the spectral region around the \ion{Mg}{2} triplet at 4481 \AA\ reveals that the spectrum of HD 42035 is actually a composite spectrum
containing shallow and broad lines coming from HD 42035 B onto which are overimposed redshifted very sharp lines coming from HD 42035 S. In the SOPHIE spectrum of HD 42035 obtained in December 2014, the sharp lines of the
\ion{Mg}{2} triplet are displaced, after correction for the barycentric velocity of the Earth, by about +53.0
km\,s$^{-1}$ with respect to their laboratory positions. A comparison in figure \, \ref{CompHD42035}
of our SOPHIE spectrum of HD 42035 with an archival ELODIE spectrum (R=42000) obtained in December 2000 (also corrected for barycentric velocity of the Earth)
reveals large radial velocity variations of the sharp lines of the \ion{Mg}{2} triplet and possibly a change of the asymmetry of the broad  and shallow \ion{Mg}{2} line formed in HD 42035 B. Indeed the shallow line  had an extended blue wing in December 2014 whereas it had an extended red wing in December 2000. Recent spectroscopy obtained at DAO by E. Griffin confirms large radial velocity variations of the sharp lines and support the binary nature of this star.
The composite nature of the spectrum is seen in many features, we show in figure \, \ref{OIHD42035} the  region of Multiplet 7 of \ion{O}{1}.

   \begin{figure}[!ht]

 \centering
   \includegraphics[scale=0.83]{ff.eps}
      \caption{Fourier transforms of the \ion{Fe}{2} lines at 4508.29 and 4515.34\,\AA (solid lines) and of a synthetic profile with a $\ensuremath{v_{{\mathrm e}}\sin i}=3.7$ km\,s$^{-1}$ at a spectral resolution of 75000. On the velocity displayed on the x-axis, the position of the first zero yields the projected equatorial rotational velocity.}
         \label{FT}
   \end{figure}

   \begin{figure}[th!]
\vskip 1cm
 \centering
   \includegraphics[scale=0.35]{MgII_HD42035.eps}
      \caption{Comparison of  the composite \ion{Mg}{2} triplet profile  of HD 42035 taken 14 years apart: SOPHIE spectrum (solid line) and archival ELODIE spectrum (dashed line)
      .The vertical lines depict the laboratory wavelengths of the \ion{Mg}{2} triplet.}
         \label{CompHD42035}
   \end{figure}

   \begin{figure}[th!]
\vskip 1cm
 \centering
   \includegraphics[scale=0.35]{OI_HD42035.eps}
      \caption{Composite profile of Multiplet 7 of \ion{O}{1} for HD 42035
      The vertical lines depict the laboratory wavelengths of the \ion{O}{1} lines.}
         \label{OIHD42035}
   \end{figure}

We have measured the centroids of all sharp lines of HD 42035 S absorbing more than 2 \% of the normalized flux by adjusting gaussians. This yields wavelengths accurate to
$\pm 0.02$ \AA.
Seven hundred and thirteen line centers were thus measured.
In order to identify these lines, a model atmosphere was computed using ATLAS9 \citep{Kurucz} with 72 layers for the effective temperature and surface gravity of HD 42035 (see par. 5.1) and a solar composition. A synthetic spectrum was then computed using SYNSPEC48 (Hubeny \& Lanz, 1992) over the entire range 3900 \AA\ to 6800 \AA\ for solar abundances
and a microturbulent velocity 0.0  km\,s$^{-1}$. It was further convolved with the ROTIN3 routine (provided with SYNSPEC48) for the FWHM of SOPHIE and the $\ensuremath{v_{{\mathrm e}}\sin i}$ of HD 42035 S.
We found that the following species are definitely present in HD 42035 S: \ion{He}{1}, \ion{C}{1}, \ion{C}{2}, \ion{Mg}{2}, \ion{N}{1}, \ion{O}{1},
\ion{Na}{1}, \ion{Mg}{1}, \ion{Mg}{2}, \ion{Al}{1}, \ion{Al}{2}, \ion{Si}{1}, \ion{Si}{2}, \ion{S}{2}, \ion{Ca}{1}, \ion{Ca}{2},
\ion{Sc}{1}, \ion{Sc}{2}, \ion{Ti}{2}, \ion{V}{2}, \ion{Cr}{1}, \ion{Cr}{2}, \ion{Mn}{2}, \ion{Fe}{1}, \ion{Fe}{2},
\ion{Co}{2}, \ion{Sr}{2} and \ion{Ba}{2} (see Table 5).
In contrast, the following elements have very few lines identified and therefore their presence is not confirmed:   \ion{P}{2} (only one line)
and  \ion{Nd}{2} (only one line).
We note the absence of  \ion{Y}{2} lines and that of lines of  \ion{Hg}{2} and  \ion{Pt}{2} characteristic of HgMn stars.
Also the lines of \ion{Zr}{2} are weak and we could only use one unblended line, $\lambda$ 4443.008 \AA\, to derive the zirconium abundance.
A few sharp lines have triangular profiles, they all can be identified as \ion{Fe}{2} lines\\
In table 5, we have indicated as ``broad'' and sometimes ``asymmetric'' the lines which show a broad shallow component originating in HD 42035 B which can sometimes show an assymetry.
These broad and shallow lines are due to \ion{N}{1}, \ion{O}{1}, \ion{Mg}{2},  \ion{Si}{2},  \ion{Ti}{2} and  \ion{Fe}{2}.
From these identifications, we can infer that HD 42035 B probably has an effective temperature which cannot differ much from that of HD 42035 S.A rough estimate of its projected rotational equatorial velocity has been obtained by adjusting the broad component of the composite \ion{Mg}{2} triplet profile at 4481 \AA\ with a model of
effective temperature 10740 K, $\log g$ =3.80, solar metallicity and yields $\ensuremath{v_{{\mathrm e}}\sin i}$ about 150 $\pm$ 25 km\,s$^{-1}$.

 \section{Abundance determinations}

 \subsection{Fundamental parameters determinations}

For HD 67044, we have adopted the effective temperature $T_{\rm{eff}}$ = 10519 K and surface gravity $\log g$ = 3.72 derived by \cite{Huang} from fitting the H${\gamma}$ profiles to prediction of model atmospheres. Indeed HD 67044 does not have Str\"{o}mgren's photometry which precludes applying Napiwotzky's (\citeyear{Napiwotzki}) UVBYBETA procedure to derive its fundamental parameters.
 A spectrum synthesis of the H${\gamma}$ profile was run to confirm these parameters. For HD 42035, the effective temperature, $T_{\rm{eff}}$ = 10740 K, and surface gravity, $\log g$ = 3.80,  were also taken from \cite{Huang}. Applying Napiwotzky's procedure to the Str\"{o}mgren's photometry of HD 42035 yields very similar fundamental parameters. The effective temperature obtained in this manner is probably a complex mean of the individual effective temperatures of
 HD 42035 S and HD 42035 B as the Str\"{o}mgren's  photometry most likely measured the combination of the lights of both stars.
 The adopted effective temperatures, surface gravities, projected equatorial velocities and radial velocities of HD 67044 and HD 42035 S are collected in Table\,\ref{table:2}.

 \begin{table}
\caption{Fundamental parameters}
\label{table:2}
\centering
\begin{tabular}{c c c c c}
\hline\hline
Star ID & $T_\mathrm{eff}$ & $\log g$  & $v\sin i$ & $V_{rad}$\\
            &               &               &  (km\,s$^{-1}$) & (km\,s$^{-1}$)\\
\hline
   HD 67044 & 10519 & 3.72 & 45.0 & 14.50 \\
   HD 42035 S & 10740 & 3.80 & 3.7  & 53.00 \\
\hline
\end{tabular}
\end{table}

 \subsection{Model atmospheres and spectrum synthesis calculations}

 Plane parallel model atmospheres assuming radiative equilibrium and hydrostatic equilibrium were computed using the ATLAS9 code \citep{Kurucz}. The linelist was built from \cite{Kurucz} gfhyperall.dat and its revision gfall18sep15.dat  \footnote{\url{http://kurucz.harvard.edu/linelists/gfnew/gfall28sep15.dat}} which includes hyperfine splitting levels. A grid of synthetic spectra was computed with SYNSPEC48 \citep{Hubeny} to model the
 lines of \ion{He}{1}, \ion{C}{2}, \ion{O}{1}, \ion{Mg}{2}, \ion{Al}{1} and \ion{Al}{2},
  \ion{Si}{2}, \ion{Ca}{2}, \ion{Sc}{2}, \ion{Ti}{2}, \ion{Cr}{2}, \ion{Mn}{2}, \ion{Fe}{2}, \ion{Sr}{2},  \ion{Y}{2}, \ion{Zr}{2}, \ion{Ba}{2}  and \ion{Hg}{2} lines. Computations were iterated varying the unknown abundance until minimization of the chi-square between the observed and synthetic spectrum was achieved. The microturbulent velocity was first assumed to be  1.5 km\,s$^{-1}$ (in agreement with the run of $v_{micr}$ with $T_\mathrm{eff}$ we established in \cite{2014psce.conf..193G}). The synthesis of the \ion{Fe}{2} lines of various strengths however imposed a lower $v_{micr} = 0$ km\,s$^{-1}$ for HD 67044  typical of an HgMn star. A null microturbulent velocity was also finally adopted for HD 42035 S after unsuccessful attempts to model the \ion{Fe}{2} lines with higher velocities.

\subsection{The derived abundances of HD 67044}

We have used only unblended lines to derive the abundances. There are actually few unblended lines as the projected equatorial velocity of HD 67044 broadens significantly the lines.
For a given element, the final abundance is a weighted mean of the abundances derived for each transition. These final abundances and their estimated uncertainties for HD 67044 are collected in Table 3 \cite[the determination of the uncertainties is discussed in][]{Royer}.
Table 3 contains for each analysed species the adopted laboratory wavelength, logarithm of
oscillator strength, its source, the logarithm of the absolute abundance normalised to that of hydrogen (on a scale where $\log(N_{H}) = 12$) for each transition, and the final abundance and estimated uncertainty.
In this work, we adopted \cite{asplund2009} abundances for the Sun as a reference.
\\
The iron abundance, which is found to be about solar in HD 67044,  has been derived mostly by using several \ion{Fe}{2}  lines of multiplets 37, 38 and 186  in the range 4500--4600\,\AA\ whose atomic parameters are critically assessed in NIST\footnote{\url{http://www.nist.gov}} (these are C+ and D quality lines). These lines are widely spaced and the continuum is fairly easy to trace in this spectral region. Their synthesis always yields consistent iron abundances from the various transitions with very little dispersion. The iron abundance is probably the most accurately determined of the abundances derived here. The NLTE abundance correction for the $\lambda$ 4471.48 \AA\  \ion{He}{1} line is very small as shown by \cite{lemke} for early A-type stars and we have not corrected for it.
We find that the following species are underabundant: He, C, O, Si, Ca, Sc while Mg, V, Fe and Ni  show mild enhancements (less than 5 times solar) and Ti, Cr, Mn, Sr, Y, Zr, Ba and Hg  show pronounced overabundances (larger than 5 times solar)., the largest overabundance being for Hg. The Sr-Y-Zr triad is inverted, yttrium beeing more abundant than strontium and zirconium. The general trend is that the heaviest elements are the most overabundant.

\subsection{The provisional abundances for HD 42035 S}

Assuming an effective temperature of 10740 K, $\log g$ = 3.80 and a solar composition, for HD 42035 S, we have derived
provisional  abundances for HD 42035 S which are collected in Table 4 The effective temperature  of HD 42035 S is probably different from this mean value so that the derived abundances are rough estimates only.
Furthermore, placing a continuum level in the composite spectrum of HD 42035 turned out to be difficult because of the presence of the rotationally broadened lines and of the continuum of HD 42035 B. The provisional abundances we derive here are therefore affected by systematic errors due to (at least) the possibility of an improper placement of the continuum
(too low  compared to what it actually is). As a consequence, the residual fluxes $\frac{F}{F_{c}}$ are likely to be too large, placing the continuum higher would decrease
these residual fluxes and therefore lead to higher abundances. We therefore consider that the provisional abundances derived here are lower limits. The actual
abundances of HD 42035 S must be larger than these lower limits.
To minimize this effect  we chose to synthesize sharp lines which were as little as possible blended with the broad lines of HD 42035 B.
Realistic abundances will only be derived if we manage to disentangle the spectrum of HD 42035 S from that of HD 42035 B.
\\
We find that the light elements up to Z=21, C, O, Mg, Al, Si, S, Ca and Sc are all underabundant. The scandium deficiency is large (about 4 \% of the solar scandium abundance).
The elements heavier than Z=21 are overabundant except for Ti, Fe and Sr which are underabundant. The zirconium overabundance has been derived from only one unblended line and thus should be taken with caution.
The absence of detection of the two strongest \ion{Y}{2} lines at 3982.60 \AA\ and 5662.92 \AA\ implies that yttrium must be solar or underabundant.

\section{Conclusions}
\label{sec:conclusion}

We find that HD 67044, hitherto classified a normal B8 star,
displays underabundances of He, C, O, Si, Ca and Sc, mild overabundances of  Mg, V, Fe and Ni and pronounced enhancements of Ti, Cr, Mn, Sr, Y, Zr, Ba and Hg..
We therefore propose that HD 67044 actually is a new HgMn star and should be reclassified as such.
Current monitoring of HD 67044 suggests that the \ion{Hg}{2} line at 3983.93 \AA\ is probably variable. One possible interpretation is that HD 67044 has one or several spots of overabundant Hg over its surface. More observations are planed to address this issue.
\\
The high resolution spectrum of HD 42035 appears to be a composite spectrum where redshifted very sharp lines originating from HD 42035 S are overimposed onto the broad shallow lines originating from HD 42035 B.
The presence of a very slow rotator and a fairly fast rotator inside the same system suggests that the semi-axis of the orbit must be large.
A FT analysis of the sharp lines of HD 42035 S yields
a very low projected equatorial velocity ($\ensuremath{v_{{\mathrm e}}\sin i} = 3.7$ km\,s$^{-1}$) which makes this star one of the very few late-B type slow rotators.
The identification of 713 very sharp lines reveals that most of the lines are due to elements up to Nickel. Strontium and Barium are also present.
Assuming a similar effective temperature and surface gravity for HD 42035 S and HD 42035 B, we have derived provisional abundances for  HD 42035 S. Light elements up to Z=22 appear to be underabundant in HD 42035 S whereas heavier elements are overabundant with the exception of Fe and Sr.
These preliminary abundances are lower limits of the actual abundances of HD 42035 S because we are probably placing the continuum too low.
A complete follow-up of this new spectroscopic binary will hopefully enable a proper disentangling of the spectra of HD 42035 S and HD 42035 B and help characterize fully the nature of each of these components, in particular the nature and abundance pattern of HD 42035 S.

\begin{acknowledgements}

We thank the referee, Dr. Fiorella Castelli, for her insightful comments which resulted into many improvements.
We also thank the OHP night assistants for their helpful support during the three observing runs.
We acknowledge the use of the ELODIE archive at OHP.
\end{acknowledgements}