\section{Introduction}
\label{sec:intro}

JPEG lossy compression is one of the most common coding techniques to store images. It uses a block based coding scheme in frequency domain, e.g. discrete cosine transform (DCT), for compression. Since $B \times B$ blocks are coded independent of each other, blocking artifacts are visible in JPEG compressed images specially under low bit rate compression. Several no-reference image quality assessment models (NR-IQAs) have been proposed to objectively assess the quality of the JPEG compressed images \cite{impair1997, freq2000, harmonic2000, Wang2000, DCT2001, Wang2002, orientation2004, Perra2005, Jmath2007, NBM2008, eurasip2009, PCM2010, new2012, NJQA, Learning2014, Tchebichef2014, GridSAR}. NR-IQAs do not need any information of the reference image. NR-IQAs are of high interest because in most present and emerging practical real-world applications, the reference signals are not available \cite{SPM2011}. In the following, we will have an overview on NR-IQAs for JPEG compressed images.

In \cite{impair1997} for each block, horizontal and vertical difference at block boundaries are used to measure horizontal and vertical blockiness, respectively. The authors in \cite{freq2000} proposed a blockiness metric via analysis of harmonics. They used both the amplitude and the phase information of harmonics to compute a quality score. Harmonic analysis was also used to model another blockiness metric in \cite{harmonic2000}.

Wang et. al. \cite{Wang2000} modeled the blocky image as a non-blocky image interfered with a pure blocky signal. Energy of the blocky signal is then used to calculate a quality score. In DCT domain, a metric was proposed in \cite{DCT2001} that models the blocking artifacts by a 2-D step function. The quality score is calculated following the human vision measurement of block impairments. The metric proposed in \cite{Jmath2007} measures blockiness artifact in both the pixel and the DCT domains. In \cite{NJQA}, zero values DCT coefficients within each block are counted and a relevance map is estimated that distinguishes between naturally uniform blocks and compressed uniform blocks. For this end, an analysis in both DFT and DCT domains is conducted.

Wang et. al. \cite{Wang2002} proposed an efficient metric that measures blockiness via horizontally and vertically computed features. These features are average differences across block boundaries, average absolute difference between in-block image samples, and zero crossing rate. Using a set of subjective scores, five parameters of this model are estimated via nonlinear regression analysis. In \cite{orientation2004}, the edge orientation changes of blocks were used to measure severity of blockiness artifacts. Perra et. al. \cite{Perra2005} analyzed the horizontal, vertical and intra-block sets of $8\times8$ blocks after applying the Sobel operator to the JPEG compressed images.

The difference of block boundaries plus luminance adaptation and texture masking were used in \cite{NBM2008} to form a noticeable blockiness map (NBM). From which, the quality score is calculated by a Minkowski summation pooling. In \cite{eurasip2009}, 1-D signal profile of gradient image is used to extract block sizes and then priodic peaks in DCT domain are analyzed to calculate a quality score. Chen et. al. \cite{PCM2010} proposed a very similar metric.

In \cite{Learning2014}, three features including the corners, block boundaries (horizontal, vertical and intra-block), and color changes, together with the subjective scores are used to train a support vector regression (SVR) model. Li et. al. \cite{Tchebichef2014} measured the blocking artifacts through weighting a set of blockiness scores calculated by Tchebichef moments of different orders.

Lee and Park \cite{new2012} proposed a blockiness metric that first identifies candidates of having blockiness artifacts. The degree of blockiness of these candidates is then used to compute a quality score. Recently a blockiness metric is proposed that performs in three steps \cite{GridSAR}. Block grids are extracted in the spatial domain and their strength and regularity is measured. Afterwards, a masking function is used that gives different weights to the smooth and textured regions.

The aforementioned indices have at least one of the following drawbacks. They might not be robust to block size and block misalignment (examples are \cite{Wang2002, orientation2004, Perra2005, NJQA, Tchebichef2014, NBM2008}). They are complex (examples are \cite{DCT2001, NJQA, GridSAR, Tchebichef2014, Learning2014, eurasip2009}), or have many parameters to set (\cite{Wang2002, NJQA, GridSAR, eurasip2009, Tchebichef2014, Learning2014}). Indices like NJQA \cite{NJQA} and GridSAR \cite{GridSAR} are too much slow. Some indices need training (\cite{Wang2002, Learning2014}). Also, the range of quality scores provided by some of the indices like \cite{Wang2002} is not well defined, or they show other numerical issues \cite{GridSAR}.

In this letter, we propose a quality assessment model for JPEG compressed images that overcomes all aforementioned drawbacks. The proposed index is very simple and efficient, it is parameterless, and robust to block size and misalignment. The proposed metric called MUG is based on two simple facts about blockiness artifact. As a result of more JPEG compression, the number of unique gradient magnitude values decreases, and the median value of unique gradient magnitude values increases. The proposed blockiness metric MUG uses these two simple facts to provide accurate quality predictions for JPEG compressed images. Unlike other metrics that presume position of blocks beforehand or localize the position of blocks, MUG is not a local model and hence does not need any information on the position of blocks.

\section{Proposed Metric (MUG)}
\label{method}

The proposed index called MUG predicts the quality of JPEG compressed images as follows. Given the JPEG distorted image $\mathcal{D}$, the Scharr gradient operator is used to approximate horizontal $G_x$ and vertical $G_y$ gradients of $\mathcal{D}$: $G_x = h_x \ast \mathcal{D}$ and $G_y = h_y \ast \mathcal{D}$, where $h_x$ and $h_y$ are horizontal and vertical gradient operators, and $\ast$ denotes the convolution. From which, the gradient magnitude is computed as $G(\textbf{x})=\sqrt{G^2_x(\textbf{x})+G^2_y(\textbf{x})}$. It is worth to mention that within the context of the proposed metrics, the Scharr operator performs better than the Sobel and Prewitt operators. The proposed metric works directly on the gradient magnitude instead of directional gradients. Let's denote $uG$ as the unique numerical values of $G(\textbf{x})$. We show in the following that two properties of $uG$ can be used to predict quality of JPEG compressed images: i) number of values in $uG$, and ii) median of $uG$ values.

\subsection{Number of unique gradients (NUG)}
\label{NUG}

The number of unique gradients (NUG), e.g. the number of values in vector $uG$, indicates how many distinct edge strengths exist in JPEG compressed image $\mathcal{D}$. It is very likely that a JPEG compressed image with blocking artifacts has smaller values of NUG than its uncompressed version. To verify this statement, JPEG compressed images of TID2013 dataset were chosen. For each of the 25 distortion-free images in TID2013, there are five JPEG compressed images of different distortion levels. The values of NUG for each of the 25 sets are found inversely proportional to the amount of distortion:

\begin{equation}
  \ \text{Compression rate} \propto \frac{1}{\text{NUG}}
  \label{equ:NUG}
\end{equation}

In other words, the Spearman Rank-order Correlation coefficient (SRCC) between NUG values and mean opinion score (MOS) values is equal to 1 for each of the 25 sets. This experiment shows that aforementioned statement holds true. Fig. \ref{fig1} shows scatter plot of NUG scores against the subjective MOS on the LIVE dataset \cite{LIVE} (see experimental results section to see how this plot is drawn). This plot shows that there is noticeable correlation between NUG scores and MOS on this dataset. Unfortunately, NUG does not take into account the content of original images. An image may originally have less edge strengths variation than another. Therefore, there are cases that NUG can not fairly judge images having different contents. This issue is solved through including median of unique gradients ($MUG$) into the proposed model.

\begin{figure}[htb]
\scriptsize
\begin{minipage}[b]{0.99\linewidth}
  \centering
  \centerline{\includegraphics[height=3.5cm]{NUGplotLIVEcrop.png}}
\end{minipage}
\caption{Scatter plot of NUG scores against the subjective MOS on the LIVE dataset. The Pearson linear Correlation Coefficient (PLCC) is equal to 0.9105.}
\label{fig1}
\end{figure}

\subsection{Median of unique gradients (MUG)}
\label{MUG}

As mentioned above, the image content is a factor that needs to be taken into account. Let's repeat the same experiment on JPEG compressed images of the TID2013 dataset, but this time for the median of unique gradients ($MUG$). The experiments show that the same statement holds true, e.g. the values of $MUG$ for each of the 25 sets are proportional\footnote{Except for one case where SRCC is equal to 0.6, not 1.} to the amount of distortion:

\begin{equation}
  \ \text{Compression rate} \propto {MUG}
  \label{equ:MUG}
\end{equation}

In fact, $MUG$ determines how strong is the middle value of unique gradients which helps in taking into account the content of images. However, the values of $MUG$ are not always reliable because image quality is not only related to the edge strengths. The distribution of the unique gradients $uG$ is another factor that can not be considered by direct median value. Therefore, a simple standard deviation normalization was applied on the $uG$ values before median value being computed:

\begin{equation}
  \ uG' = \frac{uG}{\sqrt{ \sigma(uG) }}
  \label{equ:uG}
\end{equation}

Unique gradients vector $uG$ has different behavior for images having naturally uniform regions and block uniform regions. For images with mostly naturally uniform regions, the standard deviation in general decreases by more compression. In contrast, standard deviation value in general increases by more compression for images having less naturally uniform regions. Therefore, median of $uG'$ takes into account the \textit{content} of images. The effect of standard deviation normalization is visually shown in the scatter plots of Fig. \ref{fig2}.

\begin{figure}[htb]
\scriptsize
\begin{minipage}[b]{0.49\linewidth}
  \centering
  \centerline{\includegraphics[height=3.4cm]{MUGplotLIVEcrop.png}}
\end{minipage}
\begin{minipage}[b]{.49\linewidth}
  \centering
  \centerline{\includegraphics[height=3.4cm]{MUGnormplotLIVEcrop.png}}
\end{minipage}
\caption{Scatter plots of $MUG$ scores against the subjective MOS on the LIVE dataset. Left: $MUG$ without normalization (PLCC = 0.8422), and right: $MUG$ with standard deviation normalization (PLCC = 0.8768).}
\label{fig2}
\end{figure}

The proposed quality assessment model for JPEG compressed images (MUG) can be written by combining relations (\ref{equ:NUG}) and (\ref{equ:MUG}):

 %
\begin{equation}
  \ \text{MUG} = \frac{MUG}{\text{NUG}}
  \label{equ:Q}
\end{equation}

where, $MUG$ (in italic) is the median value of $uG'$. It can be seen that the proposed metric is parameterless. To the best of our knowledge, MUG is the only parameterless metric in the literature. MUG is therefore completely independent to the misalignment. This advantage is shown in the experimental results. Since the proposed metric is parameterless, it should be invariant to the block size as well. However, no dataset is available to experimentally verify this statement. It is worth to mention that when the input image is in color, MUG converts it to a luminance channel: L = 0.06R + 0.63G + 0.27B. According to \cite{invariance2001}, this conversion may be imperfect, but it is likely to offer accurate estimates of differential measurements. Therefore, image gradient computation from L should yield more accurate results. Since MUG only uses the median value of unique gradients, it might not be very accurate for some cases. In the following, MUG is modified by adding a few more unique gradient values.

\subsection{Stable MUG (MUG$^+$)}
\label{stability}

The distributions of unique gradient values can be very different for images having diverse edge information. This distribution might be skewed (usually right-skewed), bimodal, etc. Median value can suddenly be shifted to the left or right if some adjacent values of median are not considered. Therefore, the MUG index can become more stable by adding a few more values to the median. These values must be smaller than the median value because larger values than median have much more variations and might be unreliable. Suppose that $uG'$ values are sorted from smallest to largest. In this case NUG/2 is the index of median value in $uG'$. One easy way to add a few values that mentioned above is to use corresponding values of these indices: NUG$/i$, $i\in \{2, 3, ..., M+1\}$, where $i=2$ is the index of median and $M$ is the total number of values used ($M=19$ in this paper). In fact by adding these extra values, the proposed metric becomes numerically more stable. Moreover, there are cases that there are not $M$ unique values in the vector $uG'$. This property often happens when the majority of the input image or the whole image is naturally uniform or textured. Suppose that there are $1 \leq N \leq M$ of such values available. The stable MUG (called MUG$^+$) takes into account this behavior by the following formulation:

\begin{equation}
  \ \text{MUG}^+ = \frac{\text{MUG}}{M-N+1}
  \label{equ:Q}
\end{equation}
where MUG$^+ = $ MUG for $N=M$.

Apart from the block misalignment problem, several JPEG quality assessment models like \cite{Wang2002, NJQA} provide quite wrong predictions in special cases that image has large amount of naturally uniform regions and/or it is textured. Fig. \ref{chess} shows a high quality image of chessboard. This image has a very bad quality according to the \cite{Wang2002} (Q = -245.89). NJQA \cite{NJQA} likewise assessed this image as being of bad quality (Q = 0.3414). GridSAR \cite{GridSAR} was not able to provide a numerical value. MUG is equal to 0.8060 (very bad quality) which also provides wrong assessment. In contrast, MUG$^+$ = 0.0448 which truly means that chessboard image has a very good quality. This is another advantage of the proposed index MUG$^+$. Note that the datasets used in this paper do not have any image sample with this behavior.

\begin{figure}[htb]
\scriptsize
\begin{minipage}[b]{0.99\linewidth}
  \centering
  \centerline{\includegraphics[height=2.9cm]{chess.jpg}}
\end{minipage}
\caption{A high quality image of chessboard with naturally uniform and textured regions. The image size is 1024$\times$1024 and block sizes are all 128$\times$128.}
\label{chess}
\end{figure}

\section{Experimental results}
\label{results}

In the experiments, six standard datasets of natural images and a benchmark dataset of synthetic images are used. The TID2013 \cite{TID2013} dataset contains 125 JPEG compressed images in total. The CSIQ dataset \cite{MAD} has 150, LIVE dataset \cite{LIVE} has 175, VCL dataset \cite{VCL} has 138, and the MICT dataset \cite{MICT} has 84 JPEG compressed images. ESPL dataset \cite{ESPL} is a synthetic dataset which contains 100 JPEG compressed images. The TID2008 dataset \cite{TID2008} is another dataset with 100 JPEG compressed images which is in fact a subset of TID2013.

For objective evaluation, two evaluation metrics were used in the experiments: the Spearman Rank-order Correlation coefficient (SRCC), and the Pearson linear Correlation Coefficient (PLCC). The SRCC and PLCC metrics measure prediction monotonicity and prediction linearity, respectively.

To get a visual observation, the scatter plots of the proposed NR-IQA models MUG and MUG$^+$ on the LIVE dataset are shown in Fig. \ref{scatter}. The logistic function suggested in \cite{LIVE} was used to fit a curve on each plot:

\begin{equation}
  \ f(x) = \beta_1\Big(\frac{1}{2}-\frac{1}{1-e^{\beta_2(x-\beta_3)}}\Big)+\beta_4x+\beta_5
  \label{equ:REG}
\end{equation}

where $\beta_1$, $\beta_2$, $\beta_3$, $\beta_4$ and $\beta_5$ are fitting parameters computed by minimizing the mean square error between quality predictions $x$ and subjective scores MOS.

\begin{figure}[htb]
\scriptsize
\begin{minipage}[b]{0.49\linewidth}
  \centering
  \centerline{\includegraphics[height=3.4cm]{MUGLIVE.png}}
\end{minipage}
\begin{minipage}[b]{.49\linewidth}
  \centering
  \centerline{\includegraphics[height=3.4cm]{MUGplusLIVE.png}}
\end{minipage}
\caption{Scatter plots of MUG and MUG$^+$ scores against the subjective MOS on the LIVE dataset. Left: MUG (PLCC = 0.9649), and right: MUG$^+$ (PLCC = 0.9730).}
\label{scatter}
\end{figure}

\begin{table}[htb]
\centering
\scriptsize
\caption{PERFORMANCE COMPARISON OF THE IQA MODELS ON JPEG COMPRESSION DISTORTION TYPE OF SEVEN DATASETS IN TERMS OF SRCC AND PLCC}
\begin{tabular}{|l|c|c|ccccc|}
\hline
\multicolumn{2}{|c|}{Index}                                               & SSIM            & \cite{Wang2002}            & NJQA            & \cite{GridSAR}         & MUG    & MUG$^+$            \\
\multicolumn{2}{|c|}{Type}                                                & FR              & NR              & NR              & NR              & NR     & NR              \\ \hline
\multirow{2}{*}{\begin{tabular}[c]{@{}l@{}}TID\\ 2008\end{tabular}} & PLCC & \textbf{0.9540} & 0.9518          & 0.9442          & 0.9511          & 0.9408 & \textbf{0.9529} \\
                                                                    & SRCC & \textbf{0.9252} & 0.9129          & 0.8993          & 0.9166          & 0.9169 & \textbf{0.9239} \\ \hline
\multirow{2}{*}{\begin{tabular}[c]{@{}l@{}}TID\\ 2013\end{tabular}} & PLCC & 0.9544          & 0.9530          & 0.9477          & \textbf{0.9545} & 0.9419 & \textbf{0.9546} \\
                                                                    & SRCC & 0.9200          & \textbf{0.9267} & 0.8860          & \textbf{0.9309} & 0.9077 & 0.9185          \\ \hline
\multirow{2}{*}{CSIQ}                                               & PLCC & \textbf{0.9786} & 0.9751          & 0.9539          & \textbf{0.9788} & 0.9674 & 0.9717          \\
                                                                    & SRCC & 0.9546          & \textbf{0.9551} & 0.9249          & \textbf{0.9565} & 0.9304 & 0.9372          \\ \hline
\multirow{2}{*}{LIVE}                                               & PLCC & \textbf{0.9790} & \textbf{0.9787} & 0.9562          & 0.9756          & 0.9649 & 0.9730          \\
                                                                    & SRCC & \textbf{0.9764} & \textbf{0.9735} & 0.9562          & 0.9726          & 0.9596 & 0.9677          \\ \hline
\multirow{2}{*}{VCL}                                                & PLCC & 0.9257          & \textbf{0.9433} & 0.8611          & \textbf{0.9304} & 0.8683 & 0.8868          \\
                                                                    & SRCC & 0.9236          & \textbf{0.9403} & 0.8445          & \textbf{0.9313} & 0.8659 & 0.8850          \\ \hline
\multirow{2}{*}{MICT}                                               & PLCC & 0.8664          & \textbf{0.8876} & \textbf{0.8746} & 0.8305          & 0.8341 & 0.8503          \\
                                                                    & SRCC & 0.8590          & \textbf{0.8829} & \textbf{0.8728} & 0.8333          & 0.8263 & 0.8513          \\ \hline
\multirow{2}{*}{ESPL}                                               & PLCC & 0.9431          & \textbf{0.9599} & 0.8089          & \textbf{0.9623} & 0.9398 & 0.9370          \\
                                                                    & SRCC & 0.9042          & \textbf{0.9327} & 0.7388          & \textbf{0.9331} & 0.9284 & 0.9265          \\ \hline
\end{tabular}
\label{results1}
\end{table}

SSIM \cite{SSIM} as an FR-IQA, as well as five NR-IQAs including \cite{Wang2002}, NJQA \cite{NJQA}, GridSAR \cite{GridSAR}, and the proposed indices MUG and MUG$^+$ were used in the experiments. \cite{Wang2002} was chosen because it shows outstanding performance, and NJQA because it follows a different approach with promising performance. GridSAR is recently introduced blockiness metric which is also able to handle block misalignment. Table \ref{results1} provides a performance comparison between the six aforementioned FR/NR-IQAs in terms of SRCC and PLCC. The same experiment is repeated on JPEG compressed images with misaligned blocks. JPEG compressed images with misaligned blocks are generated by cutting one pixel from the borders of the images. Since only one pixel width is cropped from image borders, the MOS values should remain unchanged. When block positions are known beforehand, the NR-IQA of \cite{Wang2002} shows the best overall performance for the seven datasets. The proposed indices show consistent prediction accuracy over different datasets and comparable to the GridSAR and SSIM. The proposed indices in general outperform NJQA \cite{NJQA}. When block positions are not known, it can be seen from the Table \ref{results2} that the proposed indices, e.g. MUG and MUG$^+$, and GridSAR show almost the same prediction accuracy as in Table \ref{results1}. This means that they are robust to the block misalignment. In contrast, \cite{Wang2002} provides predictions with low accuracy.

While GridSAR performs better than MUG$^+$ on more considered datasets, it should be noted that GridSAR is a complex metric with several parameters to set. It is also computationally inefficient and numerically unstable.

\begin{table}[htb]
\centering
\scriptsize
\caption{PERFORMANCE COMPARISON OF THE IQA MODELS ON JPEG COMPRESSION DISTORTION TYPE OF SEVEN DATASETS WITH BLOCK MISALIGNMENT IN TERMS OF SRCC AND PLCC}
\begin{tabular}{|l|c|c|ccccc|}
\hline
\multicolumn{2}{|c|}{Index}                                               & SSIM            & \cite{Wang2002}   & NJQA   & \cite{GridSAR}         & MUG             & MUG$^+$            \\
\multicolumn{2}{|c|}{Type}                                                & FR              & NR     & NR     & NR              & NR              & NR              \\ \hline
\multirow{2}{*}{\begin{tabular}[c]{@{}l@{}}TID\\ 2008\end{tabular}} & PLCC & 0.9247          & 0.3742 & 0.8499 & \textbf{0.9540} & 0.9407          & \textbf{0.9528} \\
                                                                    & SRCC & 0.8989          & 0.3146 & 0.8128 & \textbf{0.9197} & 0.9171          & \textbf{0.9242} \\ \hline
\multirow{2}{*}{\begin{tabular}[c]{@{}l@{}}TID\\ 2013\end{tabular}} & PLCC & 0.9328          & 0.5087 & 0.8540 & \textbf{0.9566} & 0.9418          & \textbf{0.9545} \\
                                                                    & SRCC & 0.9096          & 0.2372 & 0.8107 & \textbf{0.9317} & 0.9075          & \textbf{0.9177} \\ \hline
\multirow{2}{*}{CSIQ}                                               & PLCC & \textbf{0.9750} & 0.6350 & 0.8899 & \textbf{0.9790} & 0.9676          & 0.9718          \\
                                                                    & SRCC & \textbf{0.9504} & 0.5642 & 0.8694 & \textbf{0.9560} & 0.9303          & 0.9370          \\ \hline
\multirow{2}{*}{LIVE}                                               & PLCC & \textbf{0.9761} & 0.5667 & 0.9214 & \textbf{0.9762} & 0.9646          & 0.9728          \\
                                                                    & SRCC & \textbf{0.9722} & 0.4088 & 0.9131 & \textbf{0.9727} & 0.9593          & 0.9673          \\ \hline
\multirow{2}{*}{VCL}                                                & PLCC & \textbf{0.9043} & 0.2949 & 0.6816 & \textbf{0.9265} & 0.8683          & 0.8867          \\
                                                                    & SRCC & \textbf{0.9017} & 0.1923 & 0.6498 & \textbf{0.9268} & 0.8652          & 0.8845          \\ \hline
\multirow{2}{*}{MICT}                                               & PLCC & 0.7967          & 0.4646 & 0.7647 & 0.8189          & \textbf{0.8316} & \textbf{0.8475} \\
                                                                    & SRCC & 0.7865          & 0.4443 & 0.7450 & 0.8217          & \textbf{0.8248} & \textbf{0.8474} \\ \hline
\multirow{2}{*}{ESPL}                                               & PLCC & \textbf{0.9510} & 0.6458 & 0.9414 & \textbf{0.9626} & 0.9398          & 0.9370          \\
                                                                    & SRCC & 0.9144          & 0.6412 & 0.9154 & \textbf{0.9333} & \textbf{0.9285} & 0.9265          \\ \hline
\end{tabular}
\label{results2}
\end{table}

\subsection{Complexity}

To show the efficiency of the proposed indices, a run-time comparison between six IQAs is performed and shown in Table \ref{time}. The experiments were performed on a Corei7 3.40 GHz CPU with 16 GB of RAM. The IQA model was implemented in MATLAB 2013b running on Windows 7. It can be seen that MUG and MUG$^+$ have satisfactory run-times. Compared to the competing metric GridSAR, the proposed metric is about 250 times faster.

\begin{table}[htb]
\centering
\caption{Run time comparison of six IQA models when applied on an image of 1080$\times$1920 size.}
\label{time}
\begin{tabular}{l|c
>{\columncolor[HTML]{656565}}c l|c}
Index & Time (ms) &  & Index    & Time (ms) \\ \hline
JPEGind \cite{Wang2002}   & 140.21      &  & MUG$^+$ & 225.52    \\
SSIM \cite{SSIM}  & 187.85      &  & GridSAR \cite{GridSAR}  & 56810.53    \\
MUG  & 222.06      &  & NJQA \cite{NJQA}     & 79983.76    \\
\hline
\end{tabular}
\end{table}

\section{Conclusion}
\label{conclusion}
In this letter, two novel image quality assessment models for JPEG compressed images were proposed. The proposed indices are very simple and do not need training. They are based on the two simple facts of gradient magnitude of JPEG compressed images. As a result of more
JPEG compression, the number of unique gradient magnitude
values decreases and the median value of unique gradient
magnitude values increases. The extensive experimental results shown that the proposed indices are robust to block misalignment and have consistent performance on seven benchmark datasets.

\section*{Acknowledgments}
The authors thank the NSERC of Canada for their financial support under Grants RGPDD 451272-13 and RGPIN 138344-14.

\bibliographystyle{IEEEtran}


\section{Introduction}

The point of this note is to prove that a language is in the complexity class
$\mathsf{PP}$ if and only if the strings of the language encode valid inferences
in a Bayesian network defined using function-free first-order logic with equality.
Before this statement can be made precise, a number of definitions are needed.
Section \ref{section:Background} summarizes the necessary background and
Section \ref{section:Characters} defines first-order Bayesian network specifications
and the complexity class $\mathsf{PP}$.
Section \ref{section:CapturePP} states and proves the former captures the latter.

\section{Background}\label{section:Background}

We collect a number of definitions here \cite{Gradel2007,Libkin2012},
so as to fix our terminology and notation.

We consider input strings in the alphabet $\{0,1\}$; that is, a {\em string} is a sequence
of $0$s and $1$s. A {\em language} is a set of strings; a {\em complexity class} is
a set of languages.
A language is {\em decided} by a Turing machine if the machine accepts each string
in the language, and rejects each string not in the language. The complexity class
$\mathsf{NP}$ contains each languages that can be decided by a nondeterministic
Turing machine with a polynomial time bound.

We focus on function-free first-order logic with equality (denoted by
$\mathsf{FFFO}$). That is, all formulas we contemplate are well-formed
formulas of first-order logic with equality but without functions, containing
predicates, negation ($\neg$), conjunction ($\wedge$), disjunction ($\vee$),
implication ($\Rightarrow$), equivalence ($\Leftrightarrow$),
 existential quatification ($\exists$)  and universal
quantification ($\forall$).
The set of predicates is the {\em vocabulary}.

A formula $\phi$ in existential function-free second-order logic
(denoted by $\mathsf{ESO}$) is a formula
of the form $\exists r_1 \dots \exists r_m \phi'$, where $\phi'$ is a sentence
of $\mathsf{FFFO}$ containing predicates $r_1,\dots,r_m$. Such a sentence allows
existential quantification over the predicates themselves. Note that again we
have equality in the language (that is, the built-in predicate $=$ is always
available).

For a given vocabulary, a structure $\mathfrak{A}$ is a pair consisting of a {\em domain} and
an {\em interpretation}. A domain is simply a set. An {\em interpretation} is
a truth assignment for every grounding of every predicate that is not
existentially quantified. As an example, consider the following formula
of $\mathsf{ESO}$ as discussed by Gr\"adel \cite{Gradel2007}:
\[
\exists\mathsf{partition}: \forall \logvar{x}: \forall \logvar{y}:
\left( \mathsf{edge}(\logvar{x},\logvar{y}) \Rightarrow
(\mathsf{partition}(\logvar{x}) \Leftrightarrow \neg \mathsf{partition}(\logvar{y}))
\right).
\]
A domain is then a set that can be
taken as the set of nodes of an input graph. An interpretation is a truth
assignment for the $\mathsf{edge}$ predicate and can be taken as the set
of edges of the input graph. The formula is satisfied if and only if it is possible
to partition the vertices into two subsets such that if a node is in one subset,
it is not in the other. That is, the formula is satisfied if and only if the input
graph is bipartite.

We only consider finite vocabularies and finite domains in this note.
If a formula $\phi(\hat{\logvar{x}})$ has free logical variables $\hat{\logvar{x}}$, then
denote by $\mathfrak{A} \models \phi(\hat{a})$ the fact that formula
$\phi(\hat{a})$ is true in structure $\mathfrak{A}$ when the logical
variables $\hat{\logvar{x}}$ are replaced by elements of the domain $\hat{a}$.
In this  case say that $\mathfrak{A}$ is a {\em model} of $\phi(\hat{a})$.

Note that if $\phi(\hat{\logvar{x}})$ is a formula in $\mathsf{ESO}$
 as in the previous paragraphs, then its interpretations runs over the groundings
of the non-quantified predicates; that is, if $\phi$ contains predicates
$r_1,\dots,r_m$ and $s_1,\dots,s_M$, but $r_1,\dots,r_m$ are all
existentially quantified, then a model for $\phi$ contains an intepretation
for $s_1,\dots,s_M$.

There is an {\em isomorphism} between structures $\mathfrak{A}_1$ and $\mathfrak{A}_2$
when there is a bijective mapping $g$ between the domains such that if
$r(a_1,\dots,a_k)$ is true in $\mathfrak{A}_1$, then $r(g(a_1),\dots,g(a_k))$ is true
in $\mathfrak{A}_2$, and moreover if $r(a_1,\dots,a_k)$ is true in $\mathfrak{A}_2$,
then $r(g^{1}(a_1),\dots,g^{-1}(a_k))$ is true in $\mathfrak{A}_1$ (where $g^{-1}$
denotes the inverse of $g$). A set of structures is {\em isomorphism-closed} if whenever
a structure is in the set, all structures that are isomorphic to it are also in the set.

We assume that every structure is given as a string, encoded as follows for a fixed
vocabulary \cite[Section 6.1]{Libkin2012}.
First, if the domain contains elements $a_1,\dots,a_n$, then the string begins with
$n$ symbols $0$ followed by $1$. The vocabulary is fixed, so we take some order for
the predicates, $r_1,\dots,r_m$. We then append, in this order, the encoding of
the interpretation of each predicate. Focus on predicate $r_i$ of arity $k$.
To encode it with respect to a domain, we need to order the elements of the domain,
say $a_1 < a_2 < \dots < a_n$. This total ordering is assumed for now to be always
available; it will be important later to check that the ordering itself can be defined. In any
case, with a total ordering we can enumerate lexicographically all $k$-tuples over the
domain. Now suppose $\hat{a}_j$ is the $j$th tuple in this enumeration; then the $j$the
bit of the encoding of $r_i$ is $1$ if $r(\hat{a}_j)$ is true in the given interpretation, and
$0$ otherwise. Thus the encoding is a string containing
$n+1+\sum_{i=1}^m n^{\mathrm{arity}(r_i)}$ symbols (either $0$ or $1$).

We can now state Fagin's theorem:

\begin{Theorem}
Let $\mathcal{S}$ be an isomorphism-closed set of finite structures of some
non-empty finite vocabulary. Then $\mathcal{S}$ is in $\mathsf{NP}$ if and only if
$\mathcal{S}$ is the class of finite models of a sentence in existential function-free
second-order logic.
\end{Theorem}

Denote by $\mathsf{CHECK}$ the problem of deciding whether an input structure is
a model of a fixed existential function-free second-order sentence.
Fagin theorem means first that $\mathsf{CHECK}$ is in $\mathsf{NP}$ (this is the easy
part of the theorem). Second, the theorem means that every  language that can be decided
by a polynomial-time nonderministic Turing machine can be exactly encoded as the set
of models for a sentence in existential second-order logic (this is the surprising part of
the theorem).
This implies that $\mathsf{CHECK}$ is $\mathsf{NP}$-hard, but the theorem is much
more elegant (because it says that there is no need for some polynomial processing outside
of the specification provided by existential second-order logic).

The significance of Fagin's theorem is that it offers a definition of $\mathsf{NP}$
that is not tied to any computational model; rather, it is tied to the expressivity
of the language that is used to specify problems. Any language
that can be decided by a polynomial nondeterministic Turing machine can
equivalently be
be decided using first-order logic with some added quantification over predicates.

\section{First-order Bayesian network specifications and the complexity class $\mathsf{PP}$}
\label{section:Characters}

We start by defining our two main characters: on one side we have Bayesian networks
that are specified using $\mathsf{FFFO}$; on the other side we have the complexity
class $\mathsf{PP}$.

It will now be convenient to view each grounded predicate $r(\hat{a})$ as a random
variable once we have a fixed vocabulary and domain. So, given a domain $\mathcal{D}$,
we understand $r(\hat{a})$ as a function over all possible interpretations of the
vocabulary, so that $r(\hat{a})(\mathbb{I})$ yields $1$ if $r(\hat{a})$ is true in interpretation
$\mathbb{I}$, and $0$ otherwise.

\subsection{First-order Bayesian network specifications}

A {\em first-order Bayesian network specification} is a directed graph where each
node is a predicate, and where each root node $r$ is associated with a probabilistic
assessment
\[
\pr{r(\hat{\logvar{x}})=1}=\alpha,
\]
while each non-root node $s(\hat{\logvar{x}})$ is associated with a formula (called the
{\em definition} of $s$)
\[
s(\hat{x}) \Leftrightarrow \phi(\hat{\logvar{x}}),
\]
where $\phi(\hat{\logvar{x}})$ is a formula in $\mathsf{FFFO}$ with free variables $\hat{\logvar{x}}$.

Given a domain, a first-order Bayesian network specification can be grounded into
a unique Bayesian network. This is done:
\begin{enumerate}
\item by producing every grounding of the predicates,
\item by associating with each grounding $r(\hat{a})$ of a root predicate the
grounded assessment $\pr{r(\hat{a})=1}=\alpha$;
\item by associating with each grounding $s(\hat{a})$ of a non-root predicate the
grounded definition $s(\hat{a}) \Leftrightarrow \phi(\hat{a})$;
\item finally, by drawing a graph where each node is a grounded predicate and where
there is an edge into each grounded non-root predicate $s(\hat{a})$ from each
grounding of a predicate that appears in the grounded definition of $s(\hat{a})$.
\end{enumerate}

Consider, as an example, the following model of asymmetric friendship, where
an individual is always a friend of herself, and where two
individuals are friends if they are both fans (of some writer, say) or if there is some
``other'' reason for it:
\begin{eqnarray}
\label{equation:Friends}
\pr{\mathsf{fan}(\logvar{x})} & = & 0.2, \nonumber \\
\pr{\mathsf{friends}(\logvar{x},\logvar{y})} & \Leftrightarrow & (\logvar{x} = \logvar{y}) \vee \nonumber \\
& & (\mathsf{fan}(\logvar{x}) \wedge \mathsf{fan}(\logvar{y})) \vee \\
& & \mathsf{other}(\logvar{x},\logvar{y}), \nonumber \\
\pr{\mathsf{other}(\logvar{x},\logvar{y})} & = & 0.1. \nonumber
\end{eqnarray}
Suppose we have domain $\mathcal{D}=\{a,b,c\}$. Figure \ref{figure:Friends} depicts
 the Bayesian
network generated by $\mathcal{D}$ and Expression (\ref{equation:Friends}).

\begin{figure}
\begin{center}
\begin{tikzpicture}
\node[rectangle,rounded corners,draw,fill=yellow] (fa) at (4,2.5) {$\mathsf{fan}(a)$};
\node[rectangle,rounded corners,draw,fill=yellow] (fb) at (6,2.5) {$\mathsf{fan}(b)$};
\node[rectangle,rounded corners,draw,fill=yellow] (fc) at (8,2.5) {$\mathsf{fan}(c)$};
\node[rectangle,rounded corners,draw,fill=yellow] (frab) at (1,1) {$\mathsf{friends}(a,b)$};
\node[rectangle,rounded corners,draw,fill=yellow] (frac) at (3,1) {$\mathsf{friends}(a,c)$};
\node[rectangle,rounded corners,draw,fill=yellow] (frba) at (5,1) {$\mathsf{friends}(b,a)$};
\node[rectangle,rounded corners,draw,fill=yellow] (frbc) at (7,1) {$\mathsf{friends}(b,c)$};
\node[rectangle,rounded corners,draw,fill=yellow] (frca) at (9,1) {$\mathsf{friends}(c,a)$};
\node[rectangle,rounded corners,draw,fill=yellow] (frcb) at (11,1) {$\mathsf{friends}(c,b)$};
\node[rectangle,rounded corners,draw,fill=yellow] (oab) at (1,0) {$\mathsf{other}(a,b)$};
\node[rectangle,rounded corners,draw,fill=yellow] (oac) at (3,0) {$\mathsf{other}(a,c)$};
\node[rectangle,rounded corners,draw,fill=yellow] (oba) at (5,0) {$\mathsf{other}(b,a)$};
\node[rectangle,rounded corners,draw,fill=yellow] (obc) at (7,0) {$\mathsf{other}(b,c)$};
\node[rectangle,rounded corners,draw,fill=yellow] (oca) at (9,0) {$\mathsf{other}(c,a)$};
\node[rectangle,rounded corners,draw,fill=yellow] (ocb) at (11,0) {$\mathsf{other}(c,b)$};
\draw[->,>=latex,thick] (fa)--(frab);
\draw[->,>=latex,thick] (fa)--(frac);
\draw[->,>=latex,thick] (fa)--(frba);
\draw[->,>=latex,thick] (fa)--(frca);
\draw[->,>=latex,thick] (fb)--(frba);
\draw[->,>=latex,thick] (fb)--(frbc);
\draw[->,>=latex,thick] (fb)--(frab);
\draw[->,>=latex,thick] (fb)--(frcb);
\draw[->,>=latex,thick] (fc)--(frca);
\draw[->,>=latex,thick] (fc)--(frcb);
\draw[->,>=latex,thick] (fc)--(frac);
\draw[->,>=latex,thick] (fc)--(frbc);
\draw[->,>=latex,thick] (oab)--(frab);
\draw[->,>=latex,thick] (oac)--(frac);
\draw[->,>=latex,thick] (oba)--(frba);
\draw[->,>=latex,thick] (obc)--(frbc);
\draw[->,>=latex,thick] (oca)--(frca);
\draw[->,>=latex,thick] (ocb)--(frcb);
\node[rectangle,rounded corners,draw,fill=yellow] (fraa) at (4,3.5) {$\mathsf{friends}(a,a)$};
\node[rectangle,rounded corners,draw,fill=yellow] (oaa) at (4,4.5) {$\mathsf{other}(a,a)$};
\draw[->,>=latex,thick] (fa)--(fraa);
\draw[->,>=latex,thick] (oaa)--(fraa);
\node[rectangle,rounded corners,draw,fill=yellow] (frbb) at (6,3.5) {$\mathsf{friends}(b,b)$};
\node[rectangle,rounded corners,draw,fill=yellow] (obb) at (6,4.5) {$\mathsf{other}(b,b)$};
\draw[->,>=latex,thick] (fb)--(frbb);
\draw[->,>=latex,thick] (obb)--(frbb);
\node[rectangle,rounded corners,draw,fill=yellow] (frcc) at (8,3.5) {$\mathsf{friends}(c,c)$};
\node[rectangle,rounded corners,draw,fill=yellow] (occ) at (8,4.5) {$\mathsf{other}(c,c)$};
\draw[->,>=latex,thick] (fc)--(frcc);
\draw[->,>=latex,thick] (occ)--(frcc);
\end{tikzpicture}
\end{center}
\caption{The Bayesian network generated by Expression (\ref{equation:Friends}) and domain
$\mathcal{D}=\{a,b,c\}$.}
\label{figure:Friends}
\end{figure}

For a given Bayesian network specification $\tau$ and a domain $\mathcal{D}$, denote
by $\mathbb{B}(\tau,\mathcal{D})$ the Bayesian network obtained by grouning $\tau$
with respect to $\mathcal{D}$.
The set of all first-order Bayesian network specifications is denoted by $\mathcal{B}(\mathsf{FFFO})$.

\subsection{Probabilistic Turing machines and the complexity class $\mathsf{PP}$}

If a Turing machine is such that, whenever its transition function maps to a
non-singleton set, the transition is selected with uniform probability
within that set, then the Turing machine is a {\em probabilistic Turing machine}.
The complexity class $\mathsf{PP}$ is the set of languages that are decided
by a probabilistic Turing machine in polynomial time, with an error probability
strictly less than $1/2$ for all input strings.

Intuitively, $\mathsf{PP}$ represents the complexity of computing probabilities
for a phenomenon that can be simulated by a polynomial probabilistic Turing machine.

This complexity class can be equivalently defined as follows:
a language is in $\mathsf{PP}$ if and only if there is a polynomial nondeterministic
Turing machine such that a string is  in the language if and only if  more than
half of the computation paths of the machine end in the accepting state when the
string is the input.
We can imagine that there is a special class of nondeterministic Turing
machines that, given an input, not only accept it or not, but actually write in some
special tape whether that input is accepted in the majority of computation paths.
Such a special machine could then be used directly to decide a language in
$\mathsf{PP}$.

\section{$\mathcal{B}(\mathsf{FFFO})$ captures $\mathsf{PP}$}
\label{section:CapturePP}

Given a first-order Bayesian network specification and a domain, an {\em evidence piece}
$\mathbf{E}$ is a partial interpretation; that is, an evidence piece assigns a truth value for
some groundings of predicates.

We encode a pair domain/evidence $(\mathcal{D},\mathbf{E})$ using the same strategy
used before to  encode a structure; however, we must take into account the fact that
a particular grounding of a predicate can be either assigned true or false or be left
without assignment. So we use a pair of symbols in $\{0,1\}$ to encode each
grounding; we assume that $00$
means ``false'' and $11$ means ``true'', while say $01$ means lack of assignment.

Say there is an {\em isomorphism} between pairs $(\mathcal{D}_1,\mathbf{E}_1)$
and $(\mathcal{D}_2,\mathbf{E}_2)$
when there is a bijective mapping $g$ between the domains such that if
$r(a_1,\dots,a_k)$ is true in $\mathbf{E}_1$, then $r(g(a_1),\dots,g(a_k))$ is true
in $\mathbf{E}_2$, and moreover if $r(a_1,\dots,a_k)$ is true in $\mathbf{E}_2$,
then $r(g^{1}(a_1),\dots,g^{-1}(a_k))$ is true in $\mathbf{E}_1$ (where again $g^{-1}$
denotes the inverse of $g$). A set of pairs domain/evidence is {\em isomorphism-closed}
if whenever a pair is in the set, all pairs that are isomorphic to it are also in the set.

Suppose a set of pairs domain/evidence
is given with respect to a fixed vocabulary $\sigma$. Once encoded, these pairs
form a language $\mathcal{L}$ that can for instance belong to $\mathsf{NP}$ or to $\mathsf{PP}$.
One can imagine building a Bayesian network specification $\tau$ on an extended
vocabulary consisting of $\sigma$ plus some additional predicates, so as to
decide this language $\mathcal{L}$ of domain/evidence pairs. For a given input pair
$(\mathcal{D},\mathbf{E})$, the Bayesian network specification and the domain
lead to a Bayesian network $\mathbb{B}(\tau,\mathcal{D})$; this
network can be used  to compute the probability of some groundings, and that
probabiility in turn can be
used to accept/reject the input. This is the sort of strategy we pursue.

The point is that we must determine some prescription by which, given a Bayesian
network and an evidence piece, one can generate an actual decision so as to accept/reject
the input pair domain/evidence. Suppose we take the following strategy. Assume that
in the extended vocabulary of $\tau$ there are two sets of distinguished auxiliary
predicates $A_1,\dots,A_{m'}$ and $B_1,\dots,B_{m''}$ that are not in $\sigma$.
We can use the Bayesian network $\mathbb{B}(\tau,\mathcal{D})$ to
compute the probability $\pr{\mathbf{A}|\mathbf{B},\mathbf{E}}$  where
$\mathbf{A}$ and $\mathbf{B}$ are interpretations of $A_1,\dots,A_{m'}$ and
$B_1,\dots,B_{m''}$ respectively. And then we might accept/reject the input on
the basis of $\pr{\mathbf{A}|\mathbf{B},\mathbf{E}}$. However, we cannot
specify particular intepretations $\mathbf{A}$ and $\mathbf{B}$ as the related
predicates are not in the vocabulary $\sigma$. Thus the sensible strategy is
to fix attention to some selected pair of intepretations;
we simply take the interpretations that assign true to every grounding.

In short: use the Bayesian network $\mathbb{B}(\tau,\mathcal{D})$
to determine whether or not
$\pr{\mathbf{A}|\mathbf{B},\mathbf{E}}>1/2$, where
$\mathbf{A}$ assigns true to every grounding of $A_1,\dots,A_{m'}$,
and
$\mathbf{B}$ assigns true to every grounding of $B_1,\dots,B_{m''}$.
If this inequality is satisfied, the input pair is {\em accepted}; if not,
the input pair is {\em rejected}.

We refer to $A_1,\dots,A_{m'}$ as the {\em conditioned predicates}, and
to $B_1,\dots,B_{m''}$ as the {\em conditioning predicates}.

Here is the main result:

\begin{Theorem}\label{theorem:CapturePP}
Let $\mathcal{S}$ be an isomorphism-closed set of pairs domain/evidence of some
non-empty finite vocabulary, where all domains are finite.
Then $\mathcal{S}$ is in $\mathsf{PP}$ if and only if $\mathcal{S}$ is the class
of domain/evidence pairs that are accepted by a fixed first-order Bayesian network
specification with fixed conditioned and conditioning   predicates.
\end{Theorem}
\begin{proof}
First, if $\mathcal{S}$ is a class of domain/query pairs that are accepted by a
fixed first-order Bayesian network specification, they can be decided by
a polynomial time probability Turing machine. To see that, note that we
can build a nondeterministic Turing machine that guesses the truth value of
all groundings that do not appear in the query (that is, not in
$\mathbf{A}\cup\mathbf{B}\cup\mathbf{E}$), and then verify whether the resulting
complete interpretation is a model of the first-order Bayesian network specification
(as model checking of a fixed first-order sentence is in $\mathsf{P}$ \cite{Libkin2012}).

To prove the other direction,
we must adapt the proof of Fagin's theorem as described by Gr\"adel \cite{Gradel2007},
along the same lines as the proof of Theorem 1 by Saluja et al.\ \cite{Saluja95}.
So, suppose that $\mathcal{L}$ is a language decided by some probabilistic
Turing machine. So equivalently there is a nondeterministic Turing machine
that determines whether the majority of its computation paths accept an input,
and accepts/rejects the input accordingly. By the mentioned proof of Fagin's
theorem, there is a first-order sentence $\phi'$ with vocabulary consisting of the vocabulary
of the input plus additional auxiliary predicates, such that each interpretation of
this joint vocabulary is a model of the sentence if it is encodes a computation
path of the Turing machine, as long as there is an available additional predicate
that is guaranteed to be a linear order on the domain. Denote by $A$ the
zero arity predicate with associated definition $A \Leftrightarrow \phi'$.
Suppose a linear order
is indeed available; then by creating a first-order Bayesian network specification
where all groundings are associated with probability $1/2$, and where a
non-root node is associated with the sentence in the proof of Fagin's theorem,
we have that the probability of the query is larger than $1/2$ iff the
majority of computation paths accept. The challenge is to encode a linear
order. To do so, introduce a new predicate $<$ and  the first-order sentence $\phi''$
that forces $<$ to be a total order, and a zero arity predicate $B$ that is
associated with definition  $B \Leftrightarrow \phi''$. Now an input domain/pair
$(\mathcal{D},\mathbf{E})$ is accepted by
the majority of computation paths in the Turing machine if and only if
we have $\pr{A|B,\mathbf{E}}>1/2$. Note that there are actually $n!$
linear orders that satisfy $B$, but for each one of these linear orders we
have the same assignments for all other predicates, hence the ratio between
accepting computations and all computations is as desired.
\end{proof}

We might picture this as follows. There is always a Turing machine $\mathbb{TM}$
and a corresponding triple $(\tau,A,B)$ such that for any
pair $(\mathcal{D},\mathbf{E})$, we have
\[
(\mathcal{D},\mathbf{E}) \mbox{ as input to } \mathbb{TM} \mbox{ with output given by }
\pr{\mathbb{TM} \mbox{ accepts } (\mathcal{D},\mathbf{E})} > 1/2,
\]
if and only if
\[
(\mathcal{D},\mathbf{E}) \mbox{ as ``input'' to }  (\tau,A,B) \mbox{ with ``output'' given by }
\mathbb{P}_{\tau,\mathcal{D}}(A|B,\mathbf{E})>1/2,
\]
where $\mathbb{P}_{\tau,\mathcal{D}}(A|B,\mathbf{E})$ denotes   probability
with respect to $\mathbb{B}(\tau,\mathcal{D})$. (Of course, there is no need
to use only zero-arity predicates $A$ and $B$, as Theorem \ref{theorem:CapturePP}
allows for sets of predicates.)

Note that the same result could be proved if every evidence piece was taken
to be a complete interpretation for the vocabulary $\sigma$. In that case we could
directly speak of structures as inputs, and then the result would more closely
mirror Fagin's theorem. However it is very appropriate, and entirely in line
with practical use, to take the inputs to a Bayesian network as the
groundings of a partially observed interpretation. Hence we have preferred
to present our main result as stated in Theorem \ref{theorem:CapturePP}.


\section{Introduction}

For a prime $p$, denote by $\mathbb{F}_p$ the field of $p$ elements and always assume that it is represented by the first $p$-many non-negative integers $\{0,1,\dots, p-1\}$.

Let $\mathbf{E}$ be an elliptic curve defined over $\mathbb{F}_p$ given by an \emph{affine Weierstrass equation}, which for $\gcd(p,6)=1$ takes the form
\begin{equation}\label{eq:weierstrass}
 y^2=x^3+Ax+B,
\end{equation}
for some $A,B\in \mathbb{F}_p$ with non-zero discriminant $4A^3+27B^2\neq 0$.

The $\mathbb{F}_p$-rational points $\E(\fp)$ of $\mathbf{E}$ form an Abelian group (with respect to the usual addition  denoted by $\oplus$) with the point at infinity $\mathcal{O}$ as the neutral element. Let $x(\cdot)$ and $y(\cdot)$ be the coordinate functions, then for a point $P\in\E(\fp)$, $P\neq \mathcal{O}$ with affine components $P=(x_P,y_P)$ we have $x(P)=x_P$ and $y(P)=y_P$.

For a given point $G\in\E(\fp)$ and \emph{initial point} $W_0\in\E(\fp)$ the \emph{elliptic curve congruential generator}  is the sequence $(W_n)$ of points of $\E(\fp)$ satisfying the recurrence relation
\begin{equation*
 W_n=G \oplus W_{n-1}=nG\oplus W_0, \quad n=1,2,\dots
\end{equation*}
We also define the \emph{elliptic curve congruential generator with respect to the first coordinate} as the sequence $(x_n)$ in $\mathbb{F}_p$ as
\begin{equation}\label{eq:def-x}
 x_n=x(W_n)=x(nG\oplus W_0), \quad n=1,2,\dots
\end{equation}

The elliptic curve congruential generator has been widely studied, many positive results have been proven about its randomness, see \cite{BD,Clegendre,CGP,CLX,ElMahassniShparlinski,GBS,HessShparlinski,HHF,Liu,LWZ,merai-1,merai-2,TopWinterhof} and see also the survey paper \cite{shparlinski-EC-survey}.
In particular, El Mahassni and Shparlinski~\cite{ElMahassniShparlinski} showed that  $(W_n)$, and so the sequence $(x_n)$,  is well-distributed.  Hess and Shparlinski~\cite{HessShparlinski}, and Topuzo\u glu and Winterhof~\cite{TopWinterhof} provided lower bounds to the linear complexity profile of the sequence $(x_n)$.

However, it is clear, that when the curve $\mathbf{E}$ is given, the sequence $(W_n)$ is predictable from two consecutive points $W_n$, $W_{n+1}$. In \cite{GutierrezIbeas}, Gutierrez and Ibeas showed that when the prime $p$ and $G$ are known, then the sequence $(W_n)$ is predictable even if just an approximation of $W_n$, $W_{n+1}$ are revealed (even if the curve $\mathbf{E}$ is private). These show that the security of the point sequence $(W_n)$ is not well-established. However these attacks use the assumption
that the prime and the curve (resp. the point $G$) are given which assumption is quite optimistic (in the viewpoint of the attacker).

In our cryptographic settings, all the parameters, as the initial point $W_0=(x_0,y_0)$, the generator $G=(x_G,y_G)$, the parameters of the curve $A,B$ and the prime $p$ are assumed to be secret and just the output of the generator $x_1,x_2,\dots$ represented as non-negative integers, are used. The main contribution of this paper is that if some consecutive elements  are revealed, then one can compute in polynomial time (polynomial in $\log p$) an elliptic curve congruential generator (where the curve possibly defined over the rational or over a residue ring) such that the generated sequence is identical to $(x_n)$ in the revealed segment. It turns out that in practice, all the secret parameters, and thus the whole sequence $(x_n)$, can be computed from eight consecutive elements.

The result suggests that for cryptographic applications the elliptic curve congruential generator should be used with great care.

In Section~\ref{sec:background} we summarize some basic facts about elliptic curves. In Section~\ref{sec:main} we present the algorithm, and in Section \ref{sec:num} we discuss results of numerical tests.

\section{Background}\label{sec:background}

In this section we summarize some basic facts about elliptic curves. First we recall the definition of the group operation $\oplus$ of $\mathbf{E}$ defined over arbitrary field $k$. Then we extend the notion of elliptic curve to the case where it is defined over a ring.

\subsection{The group law on the curve}
Let $k$ be a field of characteristic different from $2,3$. Let $\mathbf{E}$ be an elliptic curve defined over $k$ given by an affine Weierstrass equation \eqref{eq:weierstrass} with $A,B\in k$, $4A^3+27B^2\neq 0$. The group operation $\oplus$ in $\mathbf{E}$ is defined in the following way.

\begin{definition}\label{def:1}
The operation $\oplus$ is defined over $\mathbf{E}$ as follows. If $P=(x_P,y_P)$ and $Q=(x_Q,y_Q)$ are in $\mathbf{E}(k)$, then
\[
 P\oplus Q=R=(x_R,y_R),
\]
where
\begin{enumerate}[(i)]
 \item \label{def:-i} if $x_P\neq x_Q$, then
 \[
  x_R=s^2-x_P-x_Q, \quad y_R=s(x_P-x_R)-y_P, \quad \text{where } s=\frac{y_Q-y_P}{x_Q-x_P};
 \]
 \item \label{def:-ii} if $x_P=  x_Q$ but $y_P\neq y_Q$, then $P\oplus Q=\mathcal{O}$;
 \item  if $P=Q$ and $y_P\neq 0$, then
 \[
  x_R=s^2-2x_P, \quad y_R=s(x_P-x_R)-y_P, \quad \text{where } s=\frac{3x_P^2+A}{2y_P};
 \]
 \item  if $P=Q$ and $y_P= 0$, then $P\oplus Q=\mathcal{O}$.
\end{enumerate}

\end{definition}
\subsection{Elliptic curves over $\mathbb{Z}_m$}

If $m$ is a composite integer with $\gcd(m,6)=1$, elliptic curve $\mathbf{E}$ can be also defined over $\mathbb{Z}_m$ via the \emph{projective Weierstrass equation}
\begin{equation*}
 y^2z=x^3+Axz^2+Bz^3
\end{equation*}
with $A,B\in\mathbb{Z}_m$, $\gcd(4A^3+27B^2,m)=1$. The $\mathbb{Z}_m$-rational points $\mathbf{E}(\mathbb{Z}_m)$ of $\mathbf{E}$ with projective coordinates can be represented as a triple $(x:y:z)$ such that $\gcd(m,x,y,z)=1$ (but not necessarily $z=0$ or $1$).

As in the field case, group operation can be defined on $\mathbf{E}$ (see \cite{lenstra,washington}) whose formulas correspond to Definition \ref{def:1}  if the divisor in (\ref{def:-i}) or (\ref{def:-ii})  is co-prime to $m$.

We remark that for integers $m_1,m_2$, $\gcd(m_1,m_2)=1$ we have
\[
 \mathbf{E}(\mathbb{Z}_{m_1m_2})\cong \mathbf{E}(\mathbb{Z}_{m_1}) \otimes \mathbf{E}(\mathbb{Z}_{m_2})
\]
as groups. Moreover, if $\mathbf{E}$ is an elliptic curve over $\mathbb{Q}$ defined by \eqref{eq:weierstrass} with integers $A,B$, then for  $m$ with $\gcd(4A^3+27B^2,m)=1$, the map $ \mathbf{E}(\mathbb{Q})  \rightarrow \mathbf{E}(\mathbb{Z}_m)$ defined by
\[
 (x:y:z) \mapsto  (x \bmod m:y\bmod m:z\bmod m)
\]
is a group homomorphism (where the representation $(x:y:z)$ is chosen as $x,y,z\in\mathbb{Z}$ and $\gcd(x,y,z)=1$).

Finally, for arbitrary integers $m_1,m_2$ (not necessarily co-primes), the the map $\mathbf{E}(\mathbb{Z}_{m_1m_2}) \rightarrow \mathbf{E}(\mathbb{Z}_{m_1})$ defined by
\[
 (x:y:z) \mapsto  (x \bmod m_2:y\bmod m_2:z\bmod m_2 )
\]
is a group homomorphism.

\section{Predicting the congruential generator on elliptic curve over rings}\label{sec:main}

Suppose we are given an initial segment $x_1,\dots, x_s$ of a sequence $(x_n)$ generated by an elliptic curve generator as non-negative integers. We would like to predict the remainder part of this sequence and specially, to compute the parameters of the generator, namely, the prime $p$, the parameters of the curve $A$, $B$ and the points $G$, $W_0$.

If for two different generators with primes $p$ and $q$, the revealed initial segments coincide, then the same initial segment is generated by an elliptic curve generator over $\mathbb{Z}_{p\cdot q}$.
Clearly, in this case only the generator over the ring $\mathbb{Z}_{p\cdot q}$ is computable (without assuming the easiness of the integer factorization problem) and to recover the private parameters further  revealed elements are needed.

On the other hand, if
the curve $\mathbf{E}$ is defined over $\mathbb{Q}$ with non-negative integers $A,B$ and $G,W_0\in\mathbf{E}(\mathbb{Q})$ are points such that $x(iG\oplus W_0)$ are all integers for $i=1,\dots, s$, then there are infinitely many possible primes $p$ (and generators) exist, namely all large enough primes are suitable.

Thus our aim is to determine the most general elliptic curve generator (possibly over $\mathbb{Q}$ or over a ring $\mathbb{Z}_m$) which generates the same initial segment.

The following theorem shows that if at least seven initial values are revealed, then it can be computed
a curve $\mathbf{E}$ over $\mathbb{Q}$ or over $\mathbb{Z}_m$ with $p\mid m$ and points $G$, $W_0$ such that these revealed values are the initial segment of a sequence generated by an elliptic curve generator with $\mathbf{E}$, $G$ and $W_0$. If more values are revealed, then a better approximation can be given to the generator (i.e. to the prime $p$).

\begin{thm}\label{thm:1}
There is an algorithm such that for given pairwise distinct non-negative integers $x_1,\dots, x_7$ generated by the elliptic curve congruential generator \eqref{eq:def-x} with prime number $p>3$, non-singular elliptic curve $\mathbf{E}=\mathbf{E}_{A,B}$ ($A,B\in\mathbb{F}_p$), and points $W_0,G\in\E(\fp)$ with $W_n\neq \pm G$ for $n=1,\dots, 6$, the algorithm computes a curve $\mathbf{E}_{\widetilde{A},\widetilde{B}}$ over $\mathbb{Q}$ or over $\mathbb{Z}_m$ with $p\mid m$ and a pair $(\widetilde{x_G}, \widetilde{(y_G^2)})$ in polynomial time (polynomial in $\log p$) such that if
\begin{equation}\label{eq:gen}
\widetilde{x}_{n}=2\frac{\widetilde{x}_{n-1}^3+\widetilde{A} \widetilde{x}_{n-1}+\widetilde{B}+\widetilde{\left(y_G^2\right)}}{(\widetilde{x}_{n-1}-\widetilde{x_G})^2}-2(\widetilde{x}_{n-1}+\widetilde{x_G})-\widetilde{x}_{n-2}, \quad n=2,3,\dots,
\end{equation}
then $\widetilde{x}_{n}\equiv x_n \mod p$ whenever $W_n\neq \mathcal{O}$.
\end{thm}

We summarize the steps of computing in Algorithm \ref{alg:1}.

\renewcommand{\algorithmicrequire}{\textbf{Input:}}
\renewcommand{\algorithmicensure}{\textbf{Output:}}
\begin{algorithm}[ht]
\caption{Predicting the EC-LCG}\label{alg:1}

\begin{algorithmic}[1]
\Require non-negative integers $x_1,\dots, x_7$ generated by an elliptic curve congruential generator \eqref{eq:def-x}
\Ensure $\left(\widetilde{x},\widetilde{\left(y^2\right)},\widetilde{A},\widetilde{B},m\right)$ such that $p\mid m$, and if   $\widetilde{x}_n$ ($n=2,3,\dots$) generated by \eqref{eq:gen}, then $\widetilde{x}_n\equiv x_n \mod p$ whenever $W_n\neq \mathcal{O}$.
  \State define the vectors $\c_1,\c_2,\c_3,\c_4$ as the columns of the matrix $C$ in \eqref{eq:C} and $\mathbf{u}$ as \eqref{eq:u}
 \State $m\leftarrow\det(\c_1,\c_2,\c_3,\c_4,\mathbf{u})$.

 \If {$m=0$}
  \State write $\lambda_1\c_1+\lambda_2\c_2+\lambda_3\c_3+\lambda_4\c_4= \mu \mathbf{u}$ with $\lambda_1,\dots, \lambda_4,\mu\in\mathbb{Z}$, $\mu>0$
  \State $m\leftarrow(\lambda_1^2-\lambda_2\mu)/\gcd(\lambda_1^2,\mu)$

  \If {$m=0$}
    \State $ \widetilde{x}\leftarrow\frac{\lambda_1}{\mu}$, $\widetilde{\left(y^2\right)}\leftarrow \left(\frac{\lambda_1}{\mu}\right)^3+ \frac{\lambda_4}{2\mu}+\frac{\lambda_1\lambda_3}{2\mu^2}$,$\widetilde{A}\leftarrow\frac{\lambda_3}{\mu}$, $\widetilde{B}\leftarrow\frac{\lambda_4}{2\mu}-\frac{\lambda_1\lambda_3}{2\mu^2}$
  \EndIf

 \EndIf
\If {$m\neq 0$}
  \While {$\c_1$, $\c_2$, $\c_3$, $\c_4$ are linearly dependent modulo $m$}
  \State write $\lambda_1\c_1+\lambda_2\c_2+\lambda_3\c_3+\lambda_4\c_4\equiv 0 \mod m$ with $(\lambda_1,\lambda_2,\lambda_3,\lambda_4)\neq (0,0,0,0)$,
  \State $m \leftarrow\gcd(m,\lambda_1,\lambda_2,\lambda_4,\lambda_5)$
 \EndWhile
 \State write $\lambda_1\c_1+\lambda_2\c_2+\lambda_3\c_3+\lambda_4\c_4\equiv \mu \mathbf{u}  \mod m$ with  $\mu>0$
  \If {$\gcd(m,\mu)>0$}
  \State $m\leftarrow m/\gcd(m,\mu)$
 \EndIf

    \State write $\lambda_1\c_1+\lambda_2\c_2+\lambda_3\c_3+\lambda_4\c_4\equiv \mathbf{u}  \mod m$
    \If {$\lambda_1^2\not\equiv \lambda_2 \mod m$}
      \State $m\leftarrow\gcd(m,\lambda_1^2- \lambda_2 )$
    \EndIf
    \State $\widetilde{x}\leftarrow\lambda_1 \bmod m$, $\widetilde{\left(y^2\right)}\leftarrow \lambda_1^3+ \frac{\lambda_4+\lambda_1\lambda_3}{2} \bmod m$, $\widetilde{A}\leftarrow\lambda_3\bmod m$,
    \State $\widetilde{B}\leftarrow\frac{\lambda_4-\lambda_1\lambda_3}{2}\bmod m$
 \EndIf
 \State \Return  $\left(\widetilde{x},\widetilde{\left(y^2\right)},\widetilde{A},\widetilde{B},m\right)$
\end{algorithmic}
\end{algorithm}

\begin{remark}
 The assumption that we are given the initial segment of the sequence $(x_n)$ is just a technical simplification. If any segment $x_{k},\dots, x_{k+6}$ of length seven are given, substituting $W_0$ by $W_{k-1}$ we get a generator whose initial values are $x_{k},\dots, x_{k+6}$.
\end{remark}

\begin{remark}
If we increase the number of revealed elements, then  the algorithm can be extended which provides  a better approximation $m$ to $p$ and hence a better approximation to $(x_{n})$, see Section \ref{sec:num}.
In practice, $8$ elements contain enough information to learn the exact value of $p$.
\end{remark}

\begin{proof}[Proof of Theorem \ref{thm:1}]
Let us assume, that the integers $x_1,\dots, x_7$ generated by an elliptic curve congruential generator \eqref{eq:def-x} are given.

By \eqref{eq:def-x}, we can write
\[
 x_{i-1}=x(W_i\oplus (-G)), \quad  x_{i}=x(W_i), \quad x_{i+1}=x(W_i\oplus G), \quad i=2,\dots, 6,
\]
where $-G$ is the (additive) inverse of $G$: $-G=(x_G,-y_G)$.
By the addition law, by the assumption that $W_i\neq \pm G$ ($i=1,\dots, 6$) and by \eqref{eq:weierstrass} we have
\begin{align}\label{eq:x-sum}
 x_{i-1}+x_{i+1}&=\left(\frac{y_i+y}{x_i-x_G}\right)^2-x_i-x_G+\left(\frac{y_i-y}{x_i-x}\right)^2-x_i-x_G \\ \notag
                &=2\frac{y_i^2+y^2}{(x_i-x_G)^2}-2(x_i+x_G)\\
                &=2\frac{x_i^3+Ax_i+B+y_G^2}{(x_i-x_G)^2}-2(x_i+x_G),
                 \quad i=2,\dots, 6 \notag
\end{align}
in $\mathbb{F}_p$. Thus
\begin{multline*}
(x_i-x_G)^2(x_{i-1}+x_{i+1})\equiv 2(x_i^3+Ax_i+B+y_G^2)-2(x_i+x_G)(x_i-x_G)^2 \mod p, \\ i=2,\dots, 6
\end{multline*}
i.e.,
\begin{multline}\label{eq:eqsys}
(2x_i^2+2x_i(x_{i-1}+x_{i+1}))x_G+(2x_i-(x_{i-1}+x_{i+1}))x_G^2+2x_iA+2B +2y_G^2  -2x_G^3\\
\equiv (x_{i-1}+x_{i+1})x_i^2 \mod p, \quad i=2,\dots, 6.
\end{multline}
Put
\begin{equation}\label{eq:C}
C=
\left(\!
\renewcommand{\arraystretch}{1.6}
\begin{array}{c@{\hspace{10pt}}c@{\hspace{10pt}}c@{\hspace{10pt}}c@{\hspace{10pt}}c@{\hspace{10pt}}c}
2x_2^2+2x_2 (x_1+x_3) & 2x_2-(x_1+x_3) & 2x_2 & 2 & 2 & -2\\
2x_3^2+2x_3 (x_2+x_4) & 2x_3-(x_2+x_4) & 2x_3 & 2 & 2 & -2\\
2x_4^2+2x_4 (x_3+x_5) & 2x_4-(x_3+x_5) & 2x_4 & 2 & 2 & -2\\
2x_5^2+2x_5 (x_4+x_6) & 2x_5-(x_4+x_6) & 2x_5 & 2 & 2 & -2\\
2x_6^2+2x_6 (x_5+x_7) & 2x_6-(x_5+x_7) & 2x_6 & 2 & 2 & -2\\
\end{array}\!
\right)\in\mathbb{Q}^{5 \times 6}
\end{equation}
and
\begin{equation}\label{eq:u}
\mathbf{u}= \left(\!
\renewcommand{\arraystretch}{1.6}
\begin{array}{c}
(x_1+x_3)x_2^2\\
(x_2+x_4)x_3^2\\
(x_3+x_5)x_4^2\\
(x_4+x_6)x_5^2\\
(x_5+x_7)x_6^2\\
\end{array}\!
\right)\in\mathbb{Q}^{5}.
\end{equation}

Write $C=(\c_1,\dots, \c_6)$ with $\c_1,\dots, \c_6 \in\mathbb{Z}^5$. Then we have
\begin{lemma}\label{lemma:1}
Having the same assumption as in Theorem \ref{thm:1}, the columns $\c_1$, $\c_2$, $\c_3$, $\c_4$ are linearly independent over $\mathbb{F}_p$.
\end{lemma}

Assuming Lemma \ref{lemma:1}, the matrix $C$ has rank 4 over $\mathbb{F}_p$, and by \eqref{eq:eqsys} the congruence
\begin{equation}\label{eq:C-modp}
 C \cdot \mathbf{e}
 \equiv\mathbf{u}
\mod p
\end{equation}
has the  solution  $\mathbf{e}=(x_G,x_G^2, A,B,y_G^2,x_G^3)^T$.

The algorithm looks for an
integer $m$ such that the congruence \eqref{eq:C-modp} has a solution $\mathbf{e}=(e_1,\dots, e_6)^T$ modulo $m$ with the additional restriction $e_1^2 \equiv e_2 \mod m$. During the algorithm we will always have  $p\mid m$. If $m$ takes a finite value, then we can suppose that $\gcd(m,6)=1$, since $p>3$. Moreover, $a\equiv b \mod m$ for $m=0$, means that $a=b$ as rationals.

Since the congruence \eqref{eq:C-modp} has solutions, $\det(\c_1,\c_2,\c_3,\c_4,\mathbf{u})\equiv 0 \mod p$.
First, assume that $\det(\c_1,\c_2,\c_3,\c_4,\mathbf{u})\neq 0$ and put $m=\det(\c_1,\c_2,\c_3,\c_4,\mathbf{u})$.
If $\c_1$, $\c_2$, $\c_3$, $\c_4$ are linearly dependent modulo $m$, say, $\lambda_1\c_1+\lambda_2\c_2+\lambda_3\c_3+\lambda_4\c_4\equiv 0 \mod m$ with $(\lambda_1,\lambda_2,\lambda_3,\lambda_4)\neq (0,0,0,0)$, then we have $p\mid \gcd(\lambda_1,\allowbreak \lambda_2,\allowbreak \lambda_3,\lambda_4)$. Replacing $m$ to $\gcd(m,\lambda_1,\lambda_2,\lambda_4,\lambda_5)$ we reduced $m$ with maintaining the property $p\mid m$. Iterating this step, we can assume that $\c_1$, $\c_2$, $\c_3$, $\c_4$ are linearly independent modulo $m$.

By the vanishing of the determinant, $\c_1$, $\c_2$, $\c_3$, $\c_4$, $\mathbf{u}$ are linearly dependent modulo $m$: $\lambda_1\c_1+\lambda_2\c_2+\lambda_3\c_3+\lambda_4\c_4\equiv \mu \mathbf{u}  \mod m$. By the independence of $\c_1,\c_2,\c_3,\c_4$, $\mu \not \equiv 0 \mod m$ and if $\mu$ is minimal and positive, the coefficients $\lambda_1$, $\lambda_2$, $\lambda_3$, $\lambda_4$, $\mu$ are unique. If $\gcd(m,\mu)>1$, replacing $m$ to $m/\gcd(m,\mu)$ we can assume $\mu=1$, i.e. $\lambda_1\c_1+\lambda_2\c_2+\lambda_3\c_3+\lambda_4\c_4\equiv \mathbf{u}  \mod m$ with unique $\lambda_1$, $\lambda_2$, $\lambda_3$, $\lambda_4$. Finally, all solutions of the congruence
\begin{equation}\label{eq:eqsys-modm}
  C \cdot \mathbf{e}
 \equiv\mathbf{u}
\mod m
\end{equation}
can be expressed as
\begin{equation}\label{eq:eqsys-gen}
 e_1=\lambda_1, \quad e_2=\lambda_2, \quad e_3=\lambda_3, \quad e_4+e_5-e_6=\lambda_4.
\end{equation}
If $\lambda_1^2\not\equiv \lambda_2 \mod m$, we have to replace $m$ to $\gcd(m,\lambda_1^2- \lambda_2 )$.

Next, consider the case when  $\det(\c_1,\c_2,\c_3,\c_4,\mathbf{u})= 0$ (over $\mathbb{Q}$). Now, the equation
\[
   C \cdot \mathbf{e} = \mathbf{u}
\]
has solutions over $\mathbb{Q}$. By Lemma \ref{lemma:1}, $\c_1,\c_2,\c_3,\c_4$ are linearly independent over $\mathbb{F}_p$ and thus over $\mathbb{Q}$. As before
\begin{equation}\label{eq:solutionQQ}
\lambda_1\c_1+\lambda_2\c_2+\lambda_3\c_3+\lambda_4\c_4= \mu \mathbf{u}, \quad \gcd(\lambda_1,\lambda_2,\lambda_3,\lambda_4,\mu)=1, \quad \mu>0
\end{equation}
with unique integers $\lambda_1,\allowbreak \lambda_2,\lambda_3, \allowbreak \lambda_4,\mu$. Moreover, every  integer solution $(\lambda_1',\allowbreak \lambda_2',\allowbreak \lambda_3',\allowbreak \lambda_4',\mu')$ of \eqref{eq:solutionQQ} has the form $\gamma\cdot (\lambda_1,\lambda_2,\lambda_3,\lambda_4,\mu)$ with some $\gamma\in\mathbb{Z}$.   Since $(x,x^2,A,B+y^2-x^3,1)$ is a solution modulo $p$ we have for $\gcd(p,\mu)=1$, that $(\lambda_1/\mu)^2\equiv \lambda_2/\mu \mod p$. Put $m=(\lambda_1^2-\lambda_2\mu)/\gcd(\lambda_1^2,\mu)$ so $p\mid m$. If $m\neq 0$, write the solutions of \eqref{eq:eqsys-modm} in the same form as \eqref{eq:eqsys-gen}.

In both cases, if $m\neq 0$ write
\[
\widetilde{x_G}=\lambda_1,\quad \widetilde{\left(y_G^2\right)}= \lambda_1^3+ \frac{\lambda_4+\lambda_1\lambda_3}{2},\quad
\widetilde{A}=\lambda_3,\quad \widetilde{B}=\frac{\lambda_4-\lambda_1\lambda_3}{2} \quad
\text{over } \mathbb{Z}_m.
\]
Clearly, if $\widetilde{G}=(\widetilde{x_G},\widetilde{y_G})$ with a $\widetilde{y_G}\in\mathbb{Z}_m[\zeta]/(\zeta^2-\widetilde{\left(y_G^2\right)})$,  $\widetilde{y_G}^2=\widetilde{\left(y_G^2\right)}$, then $\widetilde{G}\in\mathbf{E}_{\widetilde{A},\widetilde{B}}$ over $\mathbb{Z}_m[\zeta]/(\zeta^2-\widetilde{\left(y_G^2\right)})$. Moreover, the vector
\[
\mathbf{e}=(\widetilde{x_G}, \widetilde{x_G}^2,  \widetilde{A},\widetilde{B},\widetilde{\left(y_G^2\right)},\widetilde{x_G}^3 )^T \in\mathbb{Z}_m^6
\]
is a solution of \eqref{eq:eqsys-modm}. If $\widetilde{y}_1\in\mathbb{Z}_m[\xi]/(\xi^2-x_1^3-\widetilde{A}x_1-\widetilde{B})$ is an element such that $\widetilde{W_1}=(x_1,\widetilde{y}_1)$ is on the curve $\mathbf{E}_{\widetilde{A},\widetilde{B}}$, then writing $\widetilde{W_0}=\widetilde{W_1}\oplus(-\widetilde{G})$, the sequence $(\widetilde{x}_n)$  generated by the elliptic curve congruence generator (with $\widetilde{W_0}$, $\widetilde{G}$) satisfies \eqref{eq:gen}
so $\widetilde{x}_1, \widetilde{x}_2,\dots$ are in $\mathbb{Z}_m$ and thus $\widetilde{x}_i=x_i$ ($i=1,2,\dots, 7$) as integers.

Finally, if $m=0$, write
\[
\widetilde{x_G}=\frac{\lambda_1}{\mu},\quad \widetilde{\left(y_G^2\right)}= \left(\frac{\lambda_1}{\mu}\right)^3+ \frac{\lambda_4}{2\mu}+\frac{\lambda_1\lambda_3}{2\mu^2},\quad
\widetilde{A}=\frac{\lambda_3}{\mu},\quad \widetilde{B}=\frac{\lambda_4}{2\mu}-\frac{\lambda_1\lambda_3}{2\mu^2}
\]
over $\mathbb{Q}$.
Then $\widetilde{G}=(\widetilde{x_G},\widetilde{y_G})\in\mathbf{E}_{\widetilde{A},\widetilde{B}}$ over $\mathbb{Q}\left(\sqrt{ \widetilde{\left(y^2\right)} }\right)$ and the integer sequence $x_1\dots, x_7$ is generated by $\widetilde{G}$. The set of possible primes $p$ are those ones which $p\nmid \mu$ and $p>\max\{ x_i: \ i=1,\dots, 7\}$.
\end{proof}

Finally, it is remain to prove Lemma \ref{lemma:1}.

\begin{proof}[Proof of Lemma \ref{lemma:1}]
Clearly, it is enough to show that the vectors
\begin{equation*}
\renewcommand{\arraystretch}{1.6}
\mathbf{v}_1=
\left(\!
\begin{array}{c}
 x_2^2+x_2(x_1+x_3)  \\ x_3^2+x_3(x_2+x_4)\\ x_4^2+x_4(x_3+x_5)\\x_5^2+x_5(x_4+x_6) \\ x_{6}^2+x_{6}(x_{5}+x_{7})
\end{array}\!
\right)\!,
\mathbf{v}_2=
\left(\!
\begin{array}{c}
 x_1+x_3  \\ x_2+x_4\\ x_3+x_5\\x_4+x_6 \\ x_{5}+x_{7}
\end{array}\!
\right)\!,
\mathbf{v}_3=
\left(\!
\begin{array}{c}
 x_2\\ x_3\\ x_4\\x_5 \\x_{6}
\end{array}\!
\right)\!,
\mathbf{v}_4=
\left(\!
\begin{array}{c}
 1  \\ 1\\ 1\\1 \\ 1
\end{array}\!
\right)
\end{equation*}
are linearly independent over $\mathbb{F}_p$.

Suppose to the contrary that there are $\alpha_1,\dots, \alpha_4\in\mathbb{F}_p$, $(\alpha_1,\dots, \alpha_4)\neq (0,\dots, 0)$ such that
\begin{equation*}
 \alpha_1\mathbf{v}_1+\alpha_2\mathbf{v}_2+\alpha_3\mathbf{v}_3+\alpha_4\mathbf{v}_4=0
\end{equation*}
i.e.,
\begin{equation}\label{eq:linComb-i}
\alpha_1( x_i^2+x_i(x_{i-1}+x_{i+1}))+\alpha_2(x_{i-1}+x_{i+1})+\alpha_3 x_i+\alpha_4=0, \quad i=2,\dots, 6.
\end{equation}
Substituting \eqref{eq:x-sum} to \eqref{eq:linComb-i} we have
\begin{multline*}
 \alpha_1\left( x_i^2+2x_i\left( \frac{x_i^3+Ax_i+B+y_G^2}{(x_i-x_G)^2}-(x_i+x_G)  \right)    \right)\\+2\alpha_2\left( \frac{x_i^3+Ax_i+B+y_G^2}{(x_i-x_G)^2}-(x_i+x_G)  \right) +\alpha_3 x_i+\alpha_4=0, \quad i=2,\dots, 6.
\end{multline*}
Clearing the denominator we get
\begin{multline*}
\alpha_1\left( x_i^2(x_i-x_G)^2+2x_i\left(x_i^3+Ax_i+B+y_G^2-(x_i+x_G)(x_i-x_G)^2  \right)    \right)\\
+2\alpha_2\left( x_i^3+Ax_i+B+y_G^2-(x_i+x_G)(x_i-x_G)^2  \right)\\
+\alpha_3 x_i(x_i-x_G)^2+\alpha_4(x_i-x_G)^2=0
\end{multline*}
for $i=2,\dots, 6$, which means that the polynomial $F(X)\in\mathbb{F}_p[X]$
\begin{align*}
 F(X)&=\\
 &\alpha_1\left( X^2(X-x_G)^2+2X\left(X^3+AX+B+y_G^2-(X+x_G)(X-x_G)^2  \right)\right)\\
 +&2\alpha_2\left( X^3+AX+B+y_G^2(X+x_G)(X-x_G)^2  \right) +\alpha_3 X(X-x_G)^2\\
 +& \alpha_4(X-x_G)^2
\end{align*}
has at least five zeros: $x_i$, $i=2,\dots 6$.

Write $F(X)$ into the following form
\begin{multline}\label{eq:F-reduced}
 F(X)=\alpha_1 \left(X^4 + f_4(X)\right)+\alpha_3 \left(X^3 + f_3(X)\right)+(2\alpha_2x+\alpha_4)X^2\\+ (\alpha_2(2x^2+2A)-2\alpha_4x)X+\alpha_2(-2x^3+2B+2y^2)+\alpha_4x^2,
\end{multline}
where $f_3,f_4\in\mathbb{F}_p[X]$ with $\deg f_3<3$  and $\deg f_4<4$.

Since $\deg F\leq 4$, we must have that $F(X)$ is the zero polynomial. In this case we have
\begin{equation}\label{eq:zero-coefficients}
 \alpha_1=0, \quad \alpha_3=0, \quad 2\alpha_2x+\alpha_4=0.
\end{equation}
Then the coefficients of $X$ and 1 in \eqref{eq:F-reduced} are
\[
\alpha_2(3x_G^2+A)=0, \quad \alpha_2(2x_G^3-B-y_G^2)=0.
\]

By \eqref{eq:zero-coefficients} and $(\alpha_1,\dots, \alpha_4)\neq (0,\dots, 0)$ we also have
$\alpha_2,\alpha_4\neq 0$, thus
\begin{equation*}
 3x_G^2+A=0, \quad 2x_G^3-B-y_G^2=0,
\end{equation*}
whence using \eqref{eq:weierstrass} we get
\begin{equation*}
 2y_G^2=x(3x_G^2+A)-(x_G^3+Ax_G+B-y_G^2)-(2x_G^3-B-y_G^2)=0,
\end{equation*}
thus
\begin{equation*}
 y_G=0,\quad 3x_G^2+A=0.
\end{equation*}
Since $G=(x_G,y_G)\in\E(\fp)$, we get that $x_G$ is a multiple root of the right hand side of \eqref{eq:weierstrass}, which contradict that the discriminant of the curve is non-zero.
\end{proof}

\section{Numerical tests}\label{sec:num}
I have implemented the algorithm of Theorem \ref{thm:1} in SAGE. The algorithm have been tested for 1000 random examples of generators with 500-bit primes $p$. For seven revealed sequence elements, the algorithm computed the exact values of the parameters ($p$, $A$, $B$, $G$, $W_0$) in 95,2\% of the cases. In the remainder cases, the algorithm provided a composite integer $m$ and parameters $\widetilde{A}$, $\widetilde{B}$, $\widetilde{G}$, $\widetilde{W_0}$ such that $p\mid m$ and $\widetilde{A}\equiv A$, $\widetilde{B}\equiv B$, $\widetilde{G}\equiv G$, $\widetilde{W_0}\equiv W_0 \mod p$.

If the number of revealed sequence elements increases, then the algorithm can be modified to become more effective. Namely, if there are eight revealed sequence elements, then applying the algorithm for the first and the last seven elements, two approximation $(m_1,A_1,B_1,G_1)$ and $(m_2,A_2,B_2,G_2)$ are provided. Putting $m=\gcd(m_1,m_2, A_1-A_2,B_1-B_2, x(G_1)-x(G_2))$, $m$ is a better approximation of $p$. Modifying the SAGE program in this way, the algorithm was successful in 100\% of the cases.

\section{Final remarks, open questions}
In this paper we showed that the sequence $x_n=x(nG\oplus W_0)$ ($n=1,2,\dots$) is highly predictable if at least seven consecutive elements are revealed. A natural question is whether one can reduce the number of revealed elements.

In the literature, the distribution and the linear complexity profile of the general sequence $f(nG\oplus W_0)$ with $f\in\fp(\E)$ have been also studied (where $\fp(\E)$ is the function field of $\E(\fp)$), see \cite{BD,Clegendre,CGP,CLX,HessShparlinski,HHF,Liu,LWZ,merai-1,merai-2}. The predictability of these sequences could be handle for individual functions $f$, but it is not clear whether there is a universal algorithm for all function $f$ (or at least all function $f$ with small degree).

An other possible question connected to the result is how much information about $x(nG\oplus W_0)$  we really need to recover all the private parameters.
Gutierrez and Ibeas \cite{GutierrezIbeas} showed that when the prime $p$ and $G$ are known, then the sequence $(W_n)$ is predictable even if just an approximation of $W_n$, $W_{n+1}$ are revealed. However, the problem is still open in the case when the prime is unknown.
One can also  consider this problem over arbitrary (not prime) finite fields.
Namely, if the curve is defined over a finite field $\mathbb{F}_{p^s}$ with degree $s>1$, then one can define an integer sequence as $(x_1(nG\oplus W_0),\dots, x_r(nG\oplus W_0))$ $n=1,2,\dots$   where $r\leq s$ and $x_1(P),\dots, x_s(P)$ are the coordinates of $x(P)$ with respect to a fixed basis of $\mathbb{F}_{p^s}$ over $\mathbb{F}_p$. Is it possible to recover the whole sequence from an initial segment, at least when the degree $s$ is fixed?  Clearly, the most interesting case when $p=2$ and the generator builds a binary sequence.

\section*{Acknowledgments}
The author would like to thank Harald Niederreiter for suggesting the problem. The authors would also like to thank the reviewers for their perceptive and useful comments that significantly improved the work.

The author is partially supported by the Austrian Science Fund FWF Project F5511-N26 which is part of the Special Research Program "Quasi-Monte Carlo Methods: Theory and Applications" and by Hungarian National Foundation for Scientific Research,Grant No. K100291.