\section{Introduction}

Semantic segmentation is a vital part of the visual understanding system. It aims to parse images through assigning a class label to each pixel for an image. Currently, semantic segmentation technology has been widely used in various fields such as automatic driving \cite{Cao_2019_CVPR,Yu_2018_ECCV}, and remote sensing image interpreting \cite{marmanis2018classification,lu2019multi-scale}.

\begin{figure*}[t]
\setlength{\abovecaptionskip}{0.cm}
\setlength{\belowcaptionskip}{-0.cm}
\begin{center}
   \includegraphics[width=17cm]{introduction1.pdf}
\end{center}
   \caption{Comparison between self-attention and our proposed hierarchical context modules. The red square represents the current position, and the yellow square represents the remaining positions. The dotted line in (a) indicates weak correlation between pixels. It can be seen that the non-local module models the dense pixel-level correlation between the current position and all other positions indiscriminately. But our method first models the sparse pixel-level context between pixels of the same class. Then, the regional context between different classes is captured through the proposed region context module.}
\label{introduction1}
\end{figure*}

Traditional methods mostly adopt machine learning algorithms to perform image segmentation with various techniques, such as thresholding \cite{yuzhi}, region growing \cite{grow1}, edge detection \cite{edge1, edge2}, clustering \cite{julei1,julei2}, etc. Most successful works are based on hand-crafted features, such as HOG \cite{hog}, SIFT \cite{sift}, etc. However, with the rise of deep learning, traditional methods relying on feature engineering have gradually been replaced by the convolutional neural network (CNN) with adaptive feature learning. Block-based semantic segmentation is an early representative method based on CNN. This method first extracts regular blocks from the image in a sliding window and performs classification using common CNN (such as AlexNet \cite{alexnet}, VGG \cite{VGG}, GoogLeNet \cite{googlenet} and ResNet \cite{ResNet}). The prediction result of the image block is regarded as the class of the center pixel.  For example, Sakrapee et al. \cite{block} exploit CNN for semantic pixel labelling by cropping multi-resolution image blocks. However, the block-based methods have great redundancy due to repeated feature extraction in overlapping regions.

Things changed thoroughly after the emergence of fully convolutional network (FCN). It learns a mapping from pixels to pixels, without extracting image blocks. However, due to the fixed geometric structures, the conventional FCN is inherently limited to local receptive fields. The limitation of insufficient global context information imposes a great adverse effect on its segmentation accuracy. To make up for the above deficiency of FCN, some works obtain global context information from the perspective of multi-scale aggregation. Multiple studies \cite{FCN, UNet, RefineNet, PSPNet} adopt the pooling operation to generate multi-resolution features, which are then up-sampled and aggregated for prediction. In addition, other works \cite{DeepLabV3, DeepLabV3plus} apply dilated convolution with diverse dilated rates to acquire multi-scale contextual information. However, The above method based on multi-scale aggregation adopts a non-adaptive extraction process for all pixels, which cannot meet the requirements of different context dependencies for specific positions.

Recently, some works have focused on using the self-attention mechanism \cite{Attention_is_All_you_Need} to capture global context information for semantic segmentation. OCNet \cite{OCNet} aggregates objects context by computing the correlations of each pixel and all the other pixels. Similarly, DANet\cite{DANet} and Relational Context-aware Network \cite{zhuxiaoxiang} explore dense pixel-level contextual correlations through the self-attention mechanism in both spatial and channel dimensions. However, we found that the correlation between pixels belonging to different classes is usually weak in these methods, which means that these low correlation positions have minimal impact on the feature representation of the current position. Therefore, performing dense pixel-level modeling between these pixels will give rise to enormous redundant computation.

To address the drawback of the self-attention mechanism, we propose a hierarchical context network (HCNet) to model global context information. Specifically, pixel-level correlation is still captured between pixels of the same class with strong correlation, and a unified region-level correlation is modeled for heterogeneous pixels with weak correlations. As illustrated in Figure \ref{overall_structure}, we append two streams named context stream and prior stream at the end of dilated ResNet. The prior stream is designed to provide region partition result to the context stream by the proposed multi-scale guided pre-segmentation. The context stream consists of a pixel context module (PCM) and a region context module (RCM). Concretely, the PCM is first proposed for modeling pixel-level correlation between any two positions within each homogeneous region, as illustrated in Figure \ref{introduction1} (b). Subsequently, instead of performing dense pixel-level modeling between different homogeneous regions in self-attention mechanism, we capture the correlation between the region representations by proposed RCM, as illustrated in Figure \ref{introduction1} (c). The region representation is obtained by the proposed region pooling, and the enhanced region representation is restored to the pixel representation by region unpooling. Finally, we aggregate the output of the above hierarchical context modules to obtain features with global representation. In summary, our main contributions are three-fold:

$\bullet$ In order to improve the heterogeneous pixel redundancy modeling of the self-attention mechanism, We designed a HCNet to efficiently capture global context information for more accurate semantic segmentation.

$\bullet$ A PCM is proposed to learn pixel-level dependencies within each homogeneous region generated by proposed prior pre-segmentation. A RCM is designed to model region-level context between different regions with the help of the proposed region pooling and unpooling. Through aggregating fine-grained pixel context features and coarse-grained region context features, HCNet can harvest multi-granularity representations to more robustly identify multi-scale objects.

$\bullet$ The proposed HCNet achieves leading performance on two authoritative segmentation datasets used for autonomous driving and aerial interpretation, including Cityscapes and ISPRS Vaihingen datasets.

\section{Related work}
\noindent
\textbf{Multi-scale context for segmentation.}
Fully Convolutional Networks (FCNs) \cite{FCN} successfully transform semantic segmentation into a per-pixel labeling task by replacing fully connected layers in DCNN \cite{alexnet, ResNet, VGG, googlenet, Densenet} with convolutional ones. Following that, several FCN-based works have been proposed to capture rich contextual information from the perspective of multi-scale aggregation. RefineNet \cite{RefineNet}, ExFuse \cite{ExFuse} and CCL \cite{CCL} fuse multi-resolution features through an encoder-decoder structure, which achieves the complementation of detail information and semantic information as well as obtains rich multi-scale context. Correspondingly, PSPNet \cite{PSPNet} and Deeplabv3 series \cite{DeepLabV3, DeepLabV3plus} possess abundant contextual information using parallel multi-scale branches with different sizes of pooling kernel or dilated convolution with diverse dilated rates.

\noindent
\textbf{Self-attention for segmentation.}
The Self-attention mechanism, which is first proposed in machine translation \cite{Attention_is_All_you_Need}, has been widely used to re-model feature space according to pixel-level dependencies between each pair of pixels in computer vision. \cite{Non_Local} proposes a self-attention module to model dependencies in the space-time dimension. Due to the outstanding performance on capturing contextual information, the self-attention module has been increasingly applied in various computer vision tasks \cite{Non_Local, DANet, Detection}. OCNet \cite{OCNet} aggregates objects context by computing the similarities of each pixel and all the other pixels, which is essentially equivalent to the self-attention module. Similarly, DANet \cite{DANet} and Relational Context-aware Network \cite{zhuxiaoxiang} explore contextual dependencies through the self-attention module in both spatial and channel dimensions.

Considering that pixel-level similarities between different classes are commonly insignificant, establishing dense pixel-level dependencies leads to massive redundant relationships and high complexity in time and space. Accordingly, our proposed method captures the dependencies between pixels within each homogeneous region and continues to model the correlation between different regions. It can hierarchically capture global context information while effectively reducing the computational complexity.

\noindent
\textbf{Hierarchical structure.}
There are a lot of successful applications of hierarchical structure, such as document classification \cite{Document}, response generation \cite{response} and action recognition \cite{Skeleton}. A hierarchical structure is designed to extract contextual information from both word and sentence levels for document classification and response generation \cite{Document, response}. In addition, for action recognition task, \cite{Skeleton} divides the human skeleton into five parts according to the physical structure of the human body, and then extracts their features respectively, which is fused to produce a final representation of the skeleton sequence at a higher hierarchy.

Our method introduces the hierarchical structure into the semantic segmentation task for the first time, in which we partition the whole feature map into several class-based homogeneous regions and then explore contextual information within and between regions from pixel-level and region-level in a hierarchical structure.

\begin{figure*}[htbp]
\setlength{\abovecaptionskip}{0.cm}
\setlength{\belowcaptionskip}{-0.cm}
\centering
\includegraphics[width=17cm]{overallstructure.pdf}
\caption{The main structure of Hierarchical Context Network, which consists of two streams. The upper context stream is designed to hierarchically capture global context information, and the lower prior stream provides prior information of regions for the context stream.}
\label{overall_structure}
\end{figure*}

\section{Hierarchical Context Network}

\subsection{Overview}
The overall architecture of our proposed network is shown in Figure \ref{overall_structure}, which consists of a context stream and a prior stream. To begin with, an input image is processed by dilated ResNet \cite{dilated_resnet} pre-trained on ImageNet dataset \cite{ImageNet} to produce a feature map $F$ with the spatial size of $H \times W$. Considering the importance of global context, we further introduce two hierarchical context modules on the top of dilated ResNet in the context stream, including PCM and RCM, to hierarchically capture global context information from pixel-level and region-level. Meanwhile, the prior stream is designed to provide region prior information for the context stream, in which we conceived a multi-scale guided pre-segmentation strategy to partition the feature map into several class-based homogeneous regions. Finally, the feature map enhanced by global information is up-sampled to the original resolution, which is then fed into the softmax function to obtain the probability of each pixel belonging to each class. The class with the highest activation probability is considered as the final prediction of the pixel.

\subsection{Multi-scale guided pre-segmentation}

\begin{figure}[t]
\setlength{\abovecaptionskip}{0.cm}
\setlength{\belowcaptionskip}{-0.cm}
\begin{center}
  \includegraphics[width=8cm]{preseg.pdf}
\end{center}
  \caption{Multi-scale guided pre-segmentation module for obtaining homogeneous regions.}
\label{preseg}
\end{figure}

The prior stream aims to partition the input feature map and provide the region partition result for the context stream.
At first, we tried to use superpixel to achieve that. However, most of superpixel segmentation methods perform unsupervised iterative clustering. On the one hand, the iterative process leads to huge computational complexity. On the other hand, these methods cannot guarantee accurate semantic homogeneous regions considering the difference in object appearance and unsupervised optimization process.

Therefore, we propose a multi-scale guided pre-segmentation module, which can flexibly partition features into class-based homogeneous regions according to the supervised guidance of ground truth. As shown in Figure \ref{preseg}, the input feature map $F$ is fed into three parallel dilated convolutions with dilation rates (1, 3, 5). Each convolution output feature map has 64 channels Then, the feature maps of the three branches are aggregated through element-wise addition. Finally, a $1 \times 1$ convolution layers and softmax function is applied to obtain affiliated probability prediction $Q \in \mathbb{R}^{N \times H \times W}$, in which $N$ represents the number of classes. During training, we use prior loss $L_{prior}$ to supervise affiliated probability prediction $Q$.

The proposed multi-scale guided pre-segmentation can generate semantic homogeneous regions by introducing few convolution parameters. In particular, the convolutions with three different dilated rates can integrate multi-scale features to enhance the sensitivity of proposed pre-segmentation module to multi-scale objects. Moreover, the auxiliary supervision $L_{prior}$ can directly transfer the gradient to the shallower layer while accelerating the network training process.

\subsection{Context stream}

\subsubsection{Pixel Context module}

\begin{figure}[t]
\setlength{\abovecaptionskip}{0.cm}
\setlength{\belowcaptionskip}{-0.cm}
\begin{center}
   \includegraphics[width=8cm]{intra.pdf}
\end{center}
   \caption{llustration of the proposed Pixel Context module. The red arrow represents the prior information guidance.}
\label{intra}
\end{figure}

Establishing pixel-level dependencies can capture rich contextual information to enhance the representation of features. Different from \cite{DANet} \cite{Non_Local} modeling dense pixel-level dependencies on the entire feature map, we introduce a relatively sparse PCM to establish pixel-level dependencies within each homogeneous region.

As illustrated in Figure \ref{intra}, taking the affiliated probability prediction $Q \in \mathbb{R}^{N \times H \times W}$ obtained by multi-scale guided pre-segmentation, we first perform $arg max$ to get explicit region boundary prediction $T \in \mathbb{R}^{H \times W}$. And then divide the input feature map $X$ into $N$ several homogeneous regions $\{B_i|1,2,...,N\}$. For each region $B_i \in \mathbb{R}^{C \times K_{i}}$, in which $K_i$ represents the number of pixels of the $i^{th}$
region, we capture pixel-level dependencies using Position-independent Attention Module (PiAM, detailed description below) to generate $B_{i}^{'} \in \mathbb{R}^{C \times K_{i}}$,  and $\sum_{i=1}^{N} K_{i}=H \times W$. Finally, we aggregate features $\{B_{i}^{'}|1,2,...,N\}$ to reconstruct a new feature map $X^{'} \in \mathbb{R}^{C \times H \times W}$ according to the explicit region boundary prediction $T$.

Quantitatively, given the number of regions $N$, the time complexity and spatial complexity of PCM are both $O(\sum_{i=1}^{N} K_{i}^{2}C^2)$. In particular, the complexity takes the minimum value $O(\frac{H^{2}W^{2}C^{2}}{N})$ with $K_{i}=\frac{HW}{N}$ for each block. With the help of explicit region boundary prediction $T$, our PCM models pixel-level context within the class-based homogeneous region. On the one hand, it can aggregate the features of strongly associated positions in a more sparse way to enhance the pixel representation. On the other hand, ignoring the feature of weakly associated positions at the pixel-level can improve the redundancy of self-attention while not affecting the model performance.

\subsubsection{Region Context Module}

\begin{figure}[t]
\setlength{\abovecaptionskip}{0.cm}
\setlength{\belowcaptionskip}{-0.cm}
\begin{center}
   \includegraphics[width=8cm]{inter.pdf}
\end{center}
   \caption{llustration of the proposed Region Context Module.}
\label{inter}
\end{figure}

The PCM only obtains the context between pixels within each homogeneous region. This section further proposes an RCM to capture the context between different regions. By combining the proposed hierarchical PCM and RCM, the global context can be completely constructed while avoiding redundant self-attention methods. As shown in Figure \ref{inter}, our RCM mainly includes region pooling, region-level attention and region unpooling.

The purpose of region pooling is to achieve the scale conversion from fine-grained pixel representation to representative region representation.
Considering that the explicit region boundary prediction $T$ cannot provide sufficient information about the relationship between pixels and each superpixels, we adopt the affiliated probability prediction $Q$ as the mapping index from pixel to region. To be more specific, given the input feature map $X^{'} \in \mathbb{R}^{C \times H \times W}$ and affiliated probability prediction $Q \in \mathbb{R}^{N \times H \times W}$, we reshape them to $\mathbb{R}^{C \times HW}$ and $\mathbb{R}^{N \times HW}$, respectively. Then the
region representation $R \in \mathbb{R}^{C \times N}$ can be calculated as follows:
\begin{equation}
R_{i,j}= \frac{X^{'}_{i,:} Q_{j,:}^{T}}{ \sum_{k=0}^{HW}Q_{j,k}},
\end{equation}
\noindent where $R_{i,j}$ represents the feature of $i^{th}$ channel in the $j^{th}$ region. After getting the region representation, we apply the proposed PiAM to adaptively capture the region correlations and enhance region representations. Specifically, we feed the region representations $R$ into PiAM to obtain new region features $R^{'} \in \mathbb{R}^{C \times N}$. Finally, the region unpooling directly performs matrix multiplication between region features $R^{'}$ and affiliated probability prediction $Q \in \mathbb{R}^{N \times H \times W}$ to recover pixel representation $X^{''}\in \mathbb{R}^{C \times H \times W}$ as follows:
\begin{equation}
X_{i,j,k}^{''}=R_{i,:}^{'} Q_{:,j,k},
\end{equation}
\noindent where $X_{ijk}^{''}$ represents the feature of $i^{th}$ channel in row $j$ and column $k$ of feature map $X^{''}$.

Previous pooling operations usually aggregate features of regular regions indiscriminately, while our region pooling selectively aggregates pixel features according to the pixel-region affiliated probability. Therefore, it can effectively deal with irregular region pooling, and the resulting coarse-grained region features are more representative. Following that, region-level attention can adaptively capture region-level dependencies and update region  representations accordingly. The most important part is that the region unpooling can effectively recover the pixel representation with more accurate details through pixel-region affiliated probability.

\subsubsection{Position-independent attention module}

\begin{figure}[t]
\setlength{\abovecaptionskip}{0.cm}
\setlength{\belowcaptionskip}{-0.cm}
\begin{center}
   \includegraphics[width=8cm]{piam.pdf}
\end{center}
   \caption{llustration of the proposed Position-independent attention module. $\oplus$ represents the element-wise sum operation. And $\otimes$ means matrix multiplication operation.}
\label{piam}
\end{figure}

The self-attention mechanism has a talent for capturing the internal correlation of features, which is then used to update the original features. However, the previous self-attention module is usually designed for regular feature maps, and the proposed PCM and RCM put forward requirements for modeling the correlation between irregular feature set. Here we proposes a PiAM for feature correlation modeling and enhancement of irregular feature set.

As the structure shown in Figure \ref{piam}, given an input feature set $B \in \mathbb{R}^{C \times K}$, where $C$ and $K$ represent the number of feature channels and set length respectively, we first apply two different convolution layers to generate two feature maps $O$ and $P$, in which $\{O,P\} \in \mathbb{R}^{\frac{C}{4} \times K}$. Different from the squared difference of Euclidean distance, we calculate the correlation coefficient between any elements of the feature set through matrix multiplication operation:
\begin{equation}
A_{i,j} = O_{:, i} \times P_{:, j},
\end{equation}
\noindent where $A \in \mathbb{R}^{K \times K}$ represents the correlation coefficient matrix, and $A_{i, j}$ represents the correlation coefficient between elements $i$ and $j$. Subsequently, we perform normalizing rank aggregation (NRA) on matrix $A$:
\begin{equation}
A_{i,j} = \frac{A_{i,j}}{\sum{A_{i,:}}},
\end{equation}

\noindent It ensures that the enhanced feature statistics will not change significantly, because the sum of the weight of different feature elements is 1. The feature set will then be updated through matrix multiplication between $B$ and $A$ to enhance feature representation. At last, the updated  feature set is multiplied by a scale parameter $\alpha$ and performs a element-wise sum operation with the original feature set $B$ to acquire the final output $B' \in \mathbb{R}^{C \times K}$, where $\alpha$ is initialized to 0 and can be used to confirm the stability at the beginning of training.

In general, our proposed PiAM can perform feature similarity measurement and enhancement on irregular sets. Compared with the previous attention module for regular feature maps, it can effectively deal with irregular feature sets. Moreover, by ignoring the two square terms in the Euclidean distance expansion term and directly calculating the correlation coefficient matrix, the calculation cost can be significantly reduced.

\subsection{Loss function}
Considering the large variation in the number of pixels in each class in the training set, we adopt weighted cross-entropy loss function $L_{context}$ to train proposed model:

\begin{equation}
L_{context}(y_i, p_i)=\sum_{i=1}^{N}{-w_iy_ilog(p_i)},
\end{equation}

\noindent where $y_i$ represent the ground truth of current pixel and $p_i$ is the prediction result by softmax. $w_i$ represents the weight of the $i_{th}$ class, which is calculated through the median frequency balance \cite{median}:

\begin{equation}
w_i = f_{median}/f_i,
\end{equation}

\begin{equation}
f_i = n_i / \sum_{i=1}^{N}{n_i},
\end{equation}

\noindent where $f_{median}$ is the median of all these frequencies.
In addition, we introduce an auxiliary loss $L_{prior}$ to the output of multi-scale guided pre-segmentation module in Figure \ref{overall_structure} based on above weighted cross-entropy loss. And the total loss is denoted as:
\begin{equation}
L = L_{context} + \lambda L_{prior}
\end{equation}
\noindent in which $\lambda$ is a hyperparameter used to control the weight between $L_{context}$ and $L_{prior}$.

\section{Experiments}
To evaluate our proposed method, we conduct extensive experiments on Cityscapes dataset \cite{cityscapes} and ISPRS Vaihingen dataset \cite{vaihingen}, which differ greatly in the spatial distribution and the scale of objects.
The spatial distribution is related to the imaging perspective, the former is the front view and the latter is the top view. The scale variation depends on the distance between the object and the camera and the size of the object itself, where the distance between the objects and the camera in aerial image is approximately equal. For instance, the size of the car varies greatly in autopilot images, but roughly the same in aerial images, as shown in Figure \ref{contrast}. And the most important is that they are the authoritative benchmarks in the field of autonomous driving and remote sensing for semantic segmentation, respectively.

\begin{figure}[t]
\setlength{\abovecaptionskip}{0.cm}
\setlength{\belowcaptionskip}{-0.cm}
\begin{center}
   \includegraphics[width=8cm]{contrast.pdf}
\end{center}
   \caption{Illustration of the difference between autopilot images and aerial images. Due to the different imaging views (front-view and top-view), there are great differences in the distribution and scale of objects.}
\label{contrast}
\end{figure}

\subsection{Experimental Setup}
\subsubsection{Datasets}
\noindent
\textbf{Cityscapes.}
This is a large-scale dataset used for semantic urban scene understanding, which contains 5,000 images with fine annotations and 20,000 images with coarse annotations. This dataset is collected from 50 different cities and includes a total of 30 classes, 19 of which are used for actual training and validation. It is noted that in our experiments, we only use 5,000 images with fine annotations as our dataset, which is divided into 2,975, 500 and 1,525 images for training, validation and online testing.

\noindent
\textbf{ISPRS Vaihingen.}
This dataset consists of 33 airborne tiles of Vaihingen, whose size is about 2500 $\times$ 2000 pixels. Each of them contains a high-resolution TOP (True Ortho Photo) tile and corresponding DSM (Digital Surface Model) and nDSM (normalized Digital Surface Model) data with a GSD (Ground Sampling Distance) of 9 $cm$. The TOP file contains three bands corresponding to IR (near-infrared), R (red) and G (green) bands respectively. Among these images, 16 tiles are used for training, in which all pixels are classified as impervious surface, building, low vegetation, tree, car and background. The remaining 17 tiles are withheld for testing. Note that we only adopt IRRG images without DSM and nDSM data during the process of training.

\subsubsection{Implementation Details}
We implement our method in Pytorch. Following \cite{PSPNet, DeepLabV3}, we initialize the learning rate to 0.01 and adopt the poly learning rate policy whose learning rate is updated by $\left ( 1-\frac{iter}{total\_iter} \right )^{0.9}$ after each iteration. For the optimizer, we use stochastic gradient descent (SGD) with weight decay 0.0005 and momentum 0.9. To ensure the stability of parameters in normalization layers, Synchronized BN \cite{syncbn} is adopted to collect the statistics of batch normalization on the whole mini-batch. Specifically, all experiments are trained for 200 epochs with batch size 8 on 4 Tesla V100 GPUs with 16G memory per GPU. To avoid overfitting, we employ common data augmentation strategies, including random scaling in the range of [0.5, 2], random horizontal flipping, and random cropping. In particular, rotating at 90$^\circ$ interval is employed for ISPRS Vaihingen dataset to simulate the changes in flight direction. We set the crop size to 512 $\times$ 1024 for Cityscapes dataset and 768 $\times$ 768 for ISPRS Vaihingen dataset. As for the loss function, the weight $\lambda$ of the prior loss for multi-scale guided pre-segmentation is set to 0.8.

\subsubsection{Standard Pixel-wise Evaluation Metrics}
\noindent
\textbf{Cityscapes.}
To assess performance, Cityscapes benchmark relies on the standard Jaccard Index, commonly known as intersection-over-union (IoU):
\begin{equation}
IoU = TP ⁄ (TP+FP+FN),
\end{equation}
\noindent where TP, FP, and FN are the numbers of true positive, false positive, and false negative pixels, respectively, determined over the whole test set.

~\\
\noindent
\textbf{ISPRS Vaihingen.}
Following the evaluation metrics of ISPRS 2D Semantic Labeling Challenge \cite{vaihingen}, per-class $F_{1}$ score and overall accuracy (OA) are adopted to evaluate the performances of our proposed model. F1 score and OA are defined as follows:
\begin{equation}
\footnotesize
F_{1}=\frac{2TP}{2TP+FP+FN},~~~~
OA=\frac{TP+TN}{TP+FP+TN+FN},
\end{equation}

\subsubsection{Scale-sensitive IoU (S-IoU)}

\begin{figure}[t]
\setlength{\abovecaptionskip}{0.cm}
\setlength{\belowcaptionskip}{-0.cm}
\begin{center}
   \includegraphics[width=8cm]{siouyl.pdf}
\end{center}
   \caption{Illustration of S-IoU used for object evaluation at different scales.}
\label{siouyl}
\end{figure}

However, the above standard evaluation metrics regard pixel as the basic evaluation elements, and then count the pixel-wise performance of each class separately. This method can only reflect the performance of the model in each class, and cannot evaluate the performance of the model against multi-scale objects. To evaluate the performance of models at various scale objects, previous works \cite{gao2019focusnet:, gscnn} usually qualitatively divide objects of the same class into specific scales, for example, it was intuitively believed that buses belong to objects with large size while bicycles belong to small objects. Since the size of the objects varies greatly with the distance to camera in natural images, these methods cannot quantitatively reflect the performance of the model for multi-scale objects.

Thus, we propose an evaluation metric called scale-sensitive IoU (S-IoU) for quantitative evaluation. It regards each instance object as the basic evaluation unit with the help of the ground truth of the instance segmentation task.

Specifically, for each instance object, we first calculate its scale (area) and then match its mask in the prediction result. Subsequently, the intersection-over-union between the predicted mask and the label mask is calculated:
\begin{equation}
S-IoU=\frac{area(label) \cap area(prediction)}{area(label) \cup area(prediction)},
\end{equation}

\subsection{ABLATION study}

\begin{figure*}[htbp]
\centering
\includegraphics[width=17cm]{context_visual.pdf}
\caption{Visualization results of context modules on Cityscapes val set. For each example, we show the input image, ground truth, pre-segmentation result, final segmentation result and intra-class similarity map correspondence to the markers in the input image.}
\label{context_visual}
\end{figure*}

\begin{figure}
\centering
\includegraphics[width=8cm]{inter_class.pdf}
\caption{Visualization results of class-based region dependencies.}
\label{inter_class}
\end{figure}

In our network, PCM and RCM are used to hierarchically capture global context information from pixel-level and region-level. In this section, we show the results comparing to our baseline dilated ResNet-101 in terms of standard mIoU.

\begin{table}
\linespread{0.8}
\begin{center}
\setlength{\tabcolsep}{5pt}
\begin{tabular}{lcccc}
\toprule[0.5pt]
Method  & PCM & RCM& mIoU (\%) \\

\hline \hline

Baseline  &         &      &   75.37 \\
HCNet  &        \checkmark &         &      78.59 \\
HCNet  &         &       \checkmark &      77.26 \\
HCNet  &        \checkmark &        \checkmark &      79.86 \\
\bottomrule[0.5pt]
\end{tabular}
\caption{Ablation study for our proposed hierarchical context modules on Cityscapes val set. For fair comparison, both baseline and HCNet use ResNet-101 as the backbone.}
\label{Ablation}
\end{center}
\end{table}

In Table \ref{Ablation}, we evaluate the effectiveness of our proposed context modules on the Cityscapes validation set. It can be seen that employing the PCM individually yields a result of 78.59\% in mIoU, which outperforms baseline by 3.22\%. Meanwhile, HCNet with RCM individually realizes 77.26\% in mIoU, which results in a 1.89\% improvement compared to baseline.
It can also be seen that the RCM has achieved fewer improvements than the PCM. This is because only capturing region-level context information makes the model lack pixel-level detailed information.
After integrating pixel-level and region-level context, the performance of our hierarchical context network improves to 79.86\%  as expected.

To better understand the context modules visually, we provide visualization results for our hierarchical context modules. As shown in Figure \ref{context_visual}, we also visualize the pre-segmentation and final segmentation results in columns 3 and 4, with the yellow dashed boxes marking the improved areas. For PCM, we select one marker on each image and visualize the similarities with all the other positions belonging to the same homogeneous region as a similarity map of size $H \times W$ in column 5. Specially, the similarities between positions belonging to different regions are filled with zero according to the pre-segmentation results. For example, the first marker in the image (1) is pre-segmented as a car and the corresponding highlighted areas on the similarity map also belong to the car. Similarly, the second marker is pre-segmented into person, and its corresponding similarity map is only responded in areas of the corresponding region. For RCM, we can obtain the similarities of a certain class-based region with all class-based regions with shape of $1 \times N$. Corresponding to the two markers in Figure \ref{context_visual}, we visualize the correlations of car and person as a histogram in Figure \ref{inter_class}. It can be seen that these are specific correlations between different class-based regions, which can further enhance the discriminability of features.

We further exhibit the comparison of the increased FLOPs and GPU memory of HCNet and self-attention compared to baseline (dilated ResNet-101). As shown in Table \ref{flops}, HCNet achieves a 1.03\% improvement in mIoU compared to self-attention, but significantly reduces FLOPs and GPU memory by about 40\% and 85\%, demonstrating that our approach can capture global context information in a more sparse and effective way.

\begin{table}
\linespread{1.25}
\begin{center}
\setlength{\tabcolsep}{5pt}
\begin{tabular}{lrrc}
\toprule[0.5pt]
Method &  GFLOPs ($\blacktriangle$) & Memory(M $\blacktriangle$) & mIoU \% \\
\hline\hline
Baseline & 0.00  & 0.00   & 75.37\\
self-attention & 121.62   &  471.04 & 78.83\\
HCNet &71.63 &  70.08 & 79.86 \\
\bottomrule[0.5pt]
\end{tabular}
\caption{Comparison of self-attention and HCNet. FLOPs and Memory usage are estimated for an input of $1 \times 3 \times 512\times 1024$.}
\label{flops}
\end{center}
\end{table}

\subsection{Scale-sensitive Evaluation}

\begin{table*}[t]
\linespread{0.8}
\begin{center}
\setlength{\tabcolsep}{5pt}
\begin{tabular}{lcccc}
\toprule[0.5pt]
Method   & \begin{tabular}[c]{@{}c@{}}mS-IoU\\ area $\in$ {[}0, 2500) \end{tabular} & \begin{tabular}[c]{@{}c@{}}mS-IoU\\ area $\in$ {[}2500, 62500) \end{tabular} & \begin{tabular}[c]{@{}c@{}}mS-IoU\\ area $\in$ {[}62500, $\infty$)\end{tabular} & mS-IoU \\
\hline\hline
Baseline & 49.32  & 74.51   & 87.43 & 62.41\\
PSPNet & 51.01 & 75.26 & 89.98 & 65.07 \\
DANet & 53.47&	75.29&	89.06&	66.43\\
HCNet &54.17&	75.80&	90.82&	67.78 \\
\bottomrule[0.5pt]
\end{tabular}
\caption{Comparison in terms of S-IoU vs state-of-the-art models on the Cityscapes val set. It can be seen that HCNet achieves the best interval and overall performance compared to the other three methods.}
\label{siou}
\end{center}
\end{table*}

In Table \ref{siou}, we compare the performance of our HCNet with Baseline, DANet and PSPNet using our proposed S-IoU. Specifically, we first calculate the scale (measured by area) and S-IoU of all instance objects. Whereafter, all the objects are divided into three scale intervals by area (number of pixels) as shown in Table \ref{siou}. For the sake of comparison, we calculate the interval evaluation metric mS-IoU by averaging the S-IoU of all objects in each interval. Similarly, the overall evaluation metric can also be obtained by averaging the S-IoU of all objects. It can be seen that HCNet achieves the best interval and overall performance compared to the other three methods. This confirms that our proposed hierarchical context modules are extremely effective  for identifying multi-scale objects.

\begin{table*}[t]
\setlength{\abovecaptionskip}{3pt}
\setlength{\belowcaptionskip}{-3pt}
\linespread{0.8}
\begin{center}
\setlength{\tabcolsep}{3pt}

\small
\begin{tabular}{l|c|ccccccccccccccccccc}
Methods	& \rotatebox{90}{mIoU}	& \rotatebox{90}{road}	& \rotatebox{90}{sidewalk}	& \rotatebox{90}{building}	&\rotatebox{90}{wall}	&\rotatebox{90}{fence}	&\rotatebox{90}{pole}	&\rotatebox{90}{traffic light}	&\rotatebox{90}{traffic sign}	&\rotatebox{90}{vegetation}	&\rotatebox{90}{terrain}&	\rotatebox{90}{sky}	&\rotatebox{90}{person}	&\rotatebox{90}{rider}	&\rotatebox{90}{car}	& \rotatebox{90}{truck}&	\rotatebox{90}{bus}	&\rotatebox{90}{train}&	\rotatebox{90}{motorcycle}	&\rotatebox{90}{bicycle} \\

\hline\hline

DeepLab-v2 \cite{deeplabv2:}	& 70.4 &	97.9 &	81.3 &	90.3 &	48.8 &	47.4 &	49.6 &	57.9 	&67.3 &	91.9 &	69.4& 	94.2 &	79.8 &	59.8 	&93.7 &	56.5 	&67.5 &	57.5 &	57.7 &	68.8\\

RefineNet \cite{RefineNet}	&73.6 	&98.2 	&83.3 	&91.3 	&47.8 	&50.4 	&56.1 	&66.9 	&71.3 	&92.3 	&70.3& 	94.8 	&80.9 	&63.3 	&94.5 	&64.6 	&76.1 	&64.3 	&62.2 	&70.0 \\

DenseASPP \cite{DenseASPP}	&80.6 	&98.7 	&87.1 	&93.4 	&60.7 	&62.7 	&65.6 	&74.6 	&78.5 	&93.6 	&72.5& 	95.4 	&86.2 	&71.9 	&96.0 	&78.0 	&90.3 	&80.7 	&69.7 	&76.8 \\

BFP \cite{bfp} &81.4 &98.7 &87.0 &93.5 &59.8 &63.4 &68.9 &76.8 &80.9 &93.7 &72.8 &95.5 &87.0 &72.1 &96.0 &77.6 &89.0 &86.9 &69.2 &77.6 \\

DANet \cite{DANet}	&81.5 	&98.6 	&86.1 	&93.5 	&56.1 	&63.3 	&69.7 	&77.3 	&81.3 	&93.9 	&72.9 	&95.7 	&87.3 	&72.9 	&96.2 	&76.8 	&89.4 	&86.5 	&72.2 	&78.2 \\

ACFNet \cite{acfnet}  &81.8 &98.7 &87.1 &93.9 &60.2 &63.9 &71.1 &78.6 &81.5 &94.0 &72.9 &95.9 &88.1 &74.1 &96.5 &76.6 &89.3 &81.5 &72.1 &79.2 \\
CCNet \cite{PSANet}	&81.9 	&-	&-	&-	&-	&-	&-	&-	&-	&-	&-	&-	&-	&-	&-	&-	&-	&-	&-	&- \\

HCNet	&82.8 	&98.6 	&87.5 	&94.1 	&61.4 	&64.9 	&73.0 	&78.8 	&82.1 	&93.9 	&74.0 	&95.7 	&88.3 	&73.6 	&96.2 	&80.0 	&91.8 	&87.6 	&71.6 	&79.7 \\

\bottomrule[0.5pt]
\end{tabular}
\end{center}
\caption{Comparison with state-of-the-art on Cityscapes test set. For fair comparison, all methods only adopt data with fine annotations for training. HCNet outperforms existing approaches and achieves 82.8\% in Mean IoU.}
\label{city_table}
\end{table*}

\subsection{Benchmark Evaluation}
To get optimal performance on the benchmark test set, we use our best model (i.e., HCNet with multi-scale guided pre-segmentation and hierarchical context modules). Additionally, we use multi-scale inference schemes with scales 0.5, 1.0, and 2.0.

\subsubsection{Cityscapes Benchmark}
In Table \ref{city_table}, we compare against published methods on the Cityscapes test set without using the coarse data. Among these methods, DeepLab-v2, RefineNet and DenseASPP, which utilize multi-scale context aggregation, have not achieved leading performance due to the lack of adaptive feature aggregation for each position. However, DANet, ACFNet, and CCNet using self-attention mechanisms generally outperform multi-scale methods, which reflects the necessity of adaptive modeling global context information. By combining the PCM and the RCM, HCNet can model global context information in a sparse manner and effectively capture multi-granularity features. Experimental results show that Our model achieves the best performance 82.8\%, which is extremely competitive with recent state-of-the-art models. It is also important to stress that our model outperforms very strong to other state-of-the-art models in classes with large object size (such as bus, train, and wall) and small object size (such as pole and traffic light).

\subsubsection{ISPRS Vaihingen Benchmark}
We carry out experiments on ISPRS Vaihingen Benchmark to further evaluate the effectiveness of our method. Table \ref{vaihingen_table} shows the quantitative results of our HCNet on ISPRS Vaihingen test set. ADL\_3 \cite{ADL_3}, which adopts CNN and uses some post-processing schemes such as conditional random field for classification, obtains an overall accuracy of 88.0\%. DST\_2  \cite{DST_2} using the FCN variant and ONE\_7 \cite{ONE_7} using the SegNet variant are limited by the local receptive field, achieving overall accuracy of 89.1\% and 89.8\% respectively. The most powerful competitors are GSN \cite{GSN} using gated convolution and DLR\_10 \cite{DLR_10} combined with boundary detection. Although we did not use auxiliary techniques such as conditional random fields and boundary detection, HCNet ranks $1^{st}$ both in per-class $F_{1}$ score and overall accuracy, compared with other state-of-the-art methods. In particular, our approach achieves 88.6\% in car class (average size $38\times 18$ pixels), outperforming the previous best model by a large margin. This is because our proposed HCNet can simultaneously model global context and capture multi-granularity features.

\begin{table}
\setlength{\abovecaptionskip}{3pt}
\setlength{\belowcaptionskip}{15pt}
\linespread{0.8}
\begin{center}
\setlength{\tabcolsep}{3pt}
\small
\begin{tabular}{lcccccc}
\toprule[0.5pt]

Method & Imp surf  & Building & Low veg	&Tree& 	Car & OA  \\
\hline\hline
 ADL\_3   &  89.5 &93.2 &82.3   &88.2       &63.3 & 88.0\\
DST\_2  & 90.5    &93.7    &83.4   & 89.2 &  72.6    & 89.1\\
ONE\_7  &91.0    & 94.5  & 84.4  & 89.9 & 77.8     & 89.8\\
DLR\_10  & 92.3     & 95.2  & 84.1 & 90.0 & 79.3     & 90.3\\
GSN  & 92.2     & 95.1  & 83.7 & 89.9 &82.4      & 90.3\\
HCNet  & 93.1     & 96.1  & 84.9 & 90.1 & 88.6     & 91.4\\
\bottomrule[0.5pt]
\end{tabular}
\end{center}
\caption{Quantitative comparisons between our method and other related methods on ISPRS Vaihingen test set.}
\label{vaihingen_table}
\end{table}

\begin{figure*}[htbp]
\centering
\includegraphics[width=17cm]{city_visual.pdf}
\caption{Visualization results of HCNet on Cityscapes val set. The yellow dashed box marks the challenging areas. In addition, the black area in the ground truth is the unlabeled class, which is not included in the final evaluation.}
\label{city_visual}
\end{figure*}

\begin{figure}
\setlength{\abovecaptionskip}{3pt}
\setlength{\belowcaptionskip}{-3pt}
\centering
\includegraphics[width=8cm]{vaihingen.pdf}
 \caption{Visual results on ISPRS Vaihingen val set, where the impervious surface, building, tree, low vegetation and car are marked as white, blue, green, cyan, and yellow respectively.}
\label{vaihingen}
\end{figure}

\subsection{qualitative results}
In Figure \ref{city_visual}, we provide the qualitative comparison between the results of baseline and our HCNet on Cityscapes validation set. As shown in the yellow dashed boxes of the first example, baseline cannot identify truck considering of local receptive field, while our HCNet makes relatively accurate predictions due to the rich contextual information capturing through the hierarchical context modules. Furthermore, the scene in the second line has a highly heterogeneous appearance, which leads to a misjudgment of baseline. Thanks to the global context information provided by PCM and RCM, HCNet is able to effectively enhance the discriminability of features for accurate reasoning and prediction. In particular, in the third line, we can see that our model also has good performance for small objects.

Figure \ref{vaihingen} illustrates a few examples of segmentation results on ISPRS Vaihingen test set. It can be seen that HCNet can accurately segment objects with great differences in appearance (such as buildings) and small-scale objects (such as cars), which is benefited from our proposed context modules that can effectively capture contextual dependencies and multi-scale features.

\section{Conclusion}
In this paper, we propose a hierarchical context network (HCNet) for semantic segmentation, which can capture global context information more sparsely than the self-attention mechanism. Specifically, based on the region partition prior generated by proposed multi-scale guided pre-segmentation, PCM and RCM are designed to hierarchically model global context dependencies from pixels-level and region-level. Meanwhile, through aggregating fine-grained pixel context features and coarse-grained region context features, HCNet can harvest multi-granularity representations to more robustly identify multi-scale objects. Extensive experiments conducted on Cityscapes and ISPRS Vaihingen dataset prove the effectiveness of our proposed modules. Furthermore, our method also achieves the optimal multi-scale performance, which is evaluated through our proposed scale-sensitive IoU.

{\small
\bibliographystyle{ieee}


\section{Introduction}

Let $U=\{0\leq u_1<u_2<...<u_k<n \}$ be a subset of size $k\geq2$ of a circle ${\mathbb T}_n={\mathbb R}/n\mathbb{Z}$ or of a cyclic group $\mathbb{Z}_n=\mathbb{Z}/n\mathbb{Z}$. Denote by $s_i=u_{i+1}-u_i\in{\mathbb R}^+$, $i=1,...,k-1$, the consecutive distances between elements of $U$, and write $s_k=u_1+n-u_k\in{\mathbb R}^+$. Following Arnold \cite[\S9]{Arn}, we define the stochasticity parameter $S(U)$ of the set $U$ by
$$S(U)=\sum_{i=1}^{k}s_i^2.
$$
While we study this parameter only in the discrete case of $\mathbb{Z}_n$, we would like to  focus for a moment on the continuous case of ${\mathbb T}_n$, which is useful for some heuritics in what follows. Since $\frac{n}{k} =\frac{1}{k}\sum_{i=1}^ks_i\leq \left(\frac{1}{k}\sum_{i=1}^ks_i^2\right)^{1/2}$ and $\sum_{i=1}^ks_i^2 < \left(\sum_{i=1}^ks_i\right)^2=n^2$, we see that
$$ \min_{|U|=k} S(U)= \frac{n^2}{k}
$$
(the minimum is attained when all $s_i$ are equal to $n/k$) and
$$\sup_{|U|=k} S(U)=n^2
$$
($S(U)$ is close to this value when $U$ is contained in a small interval). Now let $U$ be a random subset of ${\mathbb T}_n$ of size $k$ (we identify it with a random point of $\{(s_1,...,s_k)\in[0,1]^k: \sum_{i=1}^ks_i=n\}$ with the normalized Lebesgue measure).
It is not hard to show that $\P\left(s_i>t\right)=\P\left(s_1>t\right) = \left(\frac{n-t}{n}\right)^{k-1}$ for each $i$ and $t\in(0,n)$, and we thus have
${\mathbb E} s_i=n/k$ and ${\mathbb E} s_i^2=\frac{2n^2}{k(k+1)}$. It follows that
\begin{equation}\label{ES(U)}
{\mathbb E} S(U)=\frac{2n^2}{k+1}
\end{equation}
and
\begin{equation}\label{rand}
\P(s_i>t{\mathbb E} s_i) = (1-t/k)^{k-1} =  e^{-t} (1+o(1)), \quad k\to\infty,
\end{equation}
uniformly for $0\leq t \leq t_0$ for any fixed $t_0>0$, and so for large $k$ the normalized gaps $s_i/{\mathbb E} s_i$ have the exponential distribution with parameter $1$.

Now let $U$ be a subset of a cyclic group $\mathbb{Z}_n$ of a fixed size $k$. Again, it is easy to show that $S(U)$ is minimal when the points of $U$ are equidistributed, and is maximal when $U$ is an interval; thus, too small or too large value of $S(U)$ points out nonrandom behaviour of $U$. Having noted these extreme cases,
we will be interested in the situations where  $S(U)$ is less or greater than the average $s(k,\mathbb{Z}_n)$ of the stochasticity parameter taken over all subsets of $\mathbb{Z}_n$ of size $k$, which can be interpreted as the stochasticity parameter of a random set of size $k$. In what follows we write $s(k)$ instead of $s(k,\mathbb{Z}_n)$ for the brevity, since the ambient group will always be clear from context; we also mention here the discrete analog of (\ref{ES(U)}):
\begin{equation}\label{1.3}
s(k)=\frac{n(2n-k+1)}{k+1};
\end{equation}
note that this agrees with (\ref{ES(U)}) in the sense that  $s(k)\sim \frac{2n^2}{k+1}$ whenever $k=o(n)$.

In this work we study the set $R_M=\{a^2: a\in\mathbb{Z}_M \}$ of quadratic residues modulo $M$. It has multiplicative structure, and, as one may expect, there are many evidences for its random behaviour as an additive set: say, in the case of a prime modulo $p$, the probability that a random element of $\mathbb{Z}_p$ is a quadratic residue is close to $1/2$, and it can be shown (see Lemma \ref{lem2.2}) that for any arbitrary disjoint subsets $C_1$ and $C_2$ of $\mathbb{Z}_p$ the events $x+c_1\in R_p$ and $x+c_2\notin R_p$, $c_1\in C_1$, $c_2\in C_2$, are mutually independent, provided that $|C_1|+|C_2|\ll \log p$ (in particular, $R_p$ contains any additive configurations of length $\ll\log p$). The last phenomenon also comes up in the case of square-free modulo $m$; see Lemma \ref{lem2.3}. In the general case, the main evidence for random behaviour of $R_M$ is the following. Write $R_M=\{0=r_1<r_2<...<r_{|R_M|}\}$, set $r_{|R_M|+1}=M$, and take an index $i$ randomly and uniformly in $\{1,...,|R_M|\}$; then on average we have ${\mathbb E}(r_{i+1}-r_i)=\frac{M}{|R_M|}$. The limit distribution of the distances $r_{i+1}-r_i$ was found by Kurlberg and Rudnick  (see \cite{KR} and \cite{Kur}); they proved that
\begin{equation}\label{KR}
\P\left(r_{i+1}-r_i>u\frac{M}{|R_M|}\right)=e^{-u}(1+o(1)), \qquad  \o(M)\to\infty,
\end{equation}
uniformly in the range $0\leq u\leq u_0$ for any fixed $u_0$, which is exactly the distribution (\ref{rand}) of gaps for a random set. This result supports the conjecture that
$$
S(R_M)=\frac{M^2}{|R_M|}{\mathbb E}\left(\frac{r_{i+1}-r_i}{M/|R_M|}\right)^2 \sim \frac{M^2}{|R_M|}\int_0^{\infty}x^2e^{-x}dx=\frac{2M^2}{|R_M|}, \quad \o(M)\to\infty;
$$
however, to make this rigorous we need good upper bounds for the contribution of large gaps between residues. On the other hand, if $\o(M)\to\infty$, then $|R_M|\to\infty$ and $M/|R_M|\to\infty$, and we have by (\ref{1.3})
$$s(|R_M|)=M\frac{2M-|R_M|+1}{|R_M|+1}\sim\frac{2M^2}{|R_M|}, \quad \o(M)\to\infty,
$$
as we should expect. The quantity $s(|R_M|)$ turns out to be a better approximation for $S(R_M)$ than $\frac{2M^2}{|R_M|}$, since it also covers the case of a prime modulo $M=p$: Garaev, Konyagin, and Malykhin \cite{GKM} proved that
\begin{equation}\label{GKM}
S(R_p)= s(|R_p|)(1+o(1)), \quad p\to\infty,
\end{equation}
(and in fact they showed the same for all large multiplicative subgroups of $\mathbb{Z}_p^*$), but it is easy to compute that $s(|R_p|)=p\frac{2p-(p+1)/2+1}{(p+1)/2}\sim 3p$, whereas $\frac{2p^2}{|R_p|}=\frac{2p^2}{(p+1)/2}\sim 4p$.

One may wonder whether (\ref{GKM}) is still true if we consider all positive integers instead of primes. We show that this is false.

\begin{theorem}\label{th1.1}
	We have
	$$\varliminf_{M\to\infty}\frac{S(R_M)}{s(|R_M|)}<1<\varlimsup_{M\to\infty}\frac{S(R_M)}{s(|R_M|)}\,.
	$$
\end{theorem}

Note that the trivial lower bound $S(R_M)\geq\frac{M^2}{|R_M|}$ (obtained from Cauchy-Schwarz), together with (\ref{1.3}), implies that $\varliminf_{M\to\infty}\frac{S(R_M)}{s(|R_M|)}\geq0.5$; on the other hand, we do not know if $\varlimsup_{M\to\infty}\frac{S(R_M)}{s(|R_M|)}$ is finite and just mention here a recent work of Aryan \cite{Ary}, where it was shown that
$$S(R_M)\ll M2^{\o(M)}\log M \prod_{p|M}\left(1+\frac{1}{\sqrt{p}}\right)\left(1-\frac1p\right)
$$
for square-free $M$.

Our second result is concerned with a more delicate comparison $S(R_M)$ with $s(|R_M|)$.

\begin{theorem}\label{th1.2}
	The set $\{ M\in \mathbb{N}: S(R_M)<s(|R_M|) \}$
	has positive density.
\end{theorem}
(Here and in what follows we consider the lower asymptotic density of a set $A\subseteq\mathbb{N}$, which is defined as $\varliminf_{N\to\infty}\frac{|A\cap\{1,\ldots,N\}|}{N}$.)

\bigskip

Both Theorems \ref{th1.1} and \ref{th1.2} are obtained by the same method. To prove Theorem \ref{th1.1} (which is technically easier), we consider moduli of the form $M=Ap$, where $A$ is fixed and $p$ is a large prime; it turns out that the situation highly depends on this additional factor.

\begin{theorem}\label{th1.3}
	Let $M=Ap$, where $A$ is square-free, $A\leq 0.1p^{0.1}$, and $(A,p)=1$. Then
	$$S(R_M)= 2f_A(0.5)p + O(A^3p^{0.95}),
	$$
	where $f_A$ is a rational function determined by $A$.
\end{theorem}

The coefficients of the function $f_A$ depend on the distances between quadratic residues modulo $A$ and can be written down explicitly; see Lemma \ref{lem5.3}.

Using (\ref{1.3}) for the same $M$, we find
\begin{equation*}\label{rand2}
s(|R_M|)=\left(\frac{4A^2}{|R_A|}-A\right)p+O(A^2|R_A|^{-1});
\end{equation*}
thus, to deduce Theorem \ref{th1.1} from Theorem \ref{th1.3}, it is enough to compare $2f_A(0.5)$ and $\frac{4A^2}{|R_A|}-A$ and find two appropriate values of $A$ which give us both possible inequalities. Numerical computations\footnote{See https://github.com/fractalon/arXiv2010.04982\_appendix. The author is grateful to Konstantin Olmezov for this code.} show that $2f_A(0.5)<\frac{4A^2}{|R_A|}-A$
for all $3\leq A\leq 100$, $A\neq89$, and
$2f_A(0.5)>\frac{4A^2}{|R_A|}-A$ for $A=89$. It implies Theorem \ref{th1.1}.

Now we turn to Theorem \ref{th1.2}. To prove it, we develop the previous idea and allow $M$ to have not one but many large prime divisors, and also vary the parameter $A$ (keeping it relatively small). Let us be more precise. Let $C_0>2$ be a large absolute constant; we denote by $\o(M)$ the number of prime divisors of $M$ and by $P^-(n)$ the least prime divisor of $n$. We consider the set $\O$ of sufficiently large positive integers $M$ such that

\begin{itemize}

	\item[(i)] $M$ is square-free;

	\item[(ii)] $0.9\log\log M < \o(M) < 1.1\log\log M$;

	\item[(iii)] $M$ has the representation of the form $M=Am$, where $3\leq A \leq (\log M)^{0.1}$ and $P^-(m)>(\log M)^{C_0}$ (note that this representation is unique).

\end{itemize}

We will show that $S(R_M)<s(|R_M|)$ for $M\in\O$. In fact, one can avoid the condition (i) here, but it simplifies our argument and added for convenience. Further, due to Hardy and Ramanujan \cite{HR}, (ii) holds for almost all integers (that is, for a set of density $1$), and thus we can assume it without affecting the density of $\O$. Finally, the third condition (iii) is crucial for us: we will need it for making some error terms small, and anyway it occurs in different steps of our proof.

We deduce Theorem \ref{th1.2} from the following two results.

\begin{theorem}\label{th1.4}
The density of the set $\O$ is equal to $e^{-\gamma}\frac{3}{5\pi^2C_0}$, where $\gamma$ is the Euler constant.
\end{theorem}

Denote by $\delta$ the density of the set $\{M\in\mathbb{N}: S(R_M)<s(|R_M|)\}$. One could extract the exact value of $C_0$ from our argument and thus obtain a lower bound for $\delta$, but this bound would be very poor. It is possible (using some additional ideas we do not present in this work) to get some better bounds, but it seems to be very hard to show that $\delta\geq 0.1$ using our method. Thus we did not try to optimize our lower bounds for $\d$, and keep the paper relatively short.

\begin{theorem}\label{th1.5}
For $M\in \O$ we have
	\begin{equation*}
	S(R_M)= \frac{2M^2}{|R_M|} - MA|R_A|^{-1}+O(M(\log M)^{-0.4}).
	\end{equation*}
\end{theorem}

On the other hand, (\ref{1.3}) implies $M^{-1}s(k)=\frac{2M-k+1}{k+1}=\frac{2(M+1)}{k+1}-1=
\frac{2M}{k}-1+O\left(\frac{M}{k^2}\right).
$
Applying this to $M\in\O$, we obtain
\begin{equation}\label{1.6}
s(|R_M|)=\frac{2M^2}{|R_M|}-M+O((\log M)^2).
\end{equation}
Since $|R_A|<A$ for any $A\geq3$, Theorem \ref{th1.2} follows from Theorems \ref{th1.4} and \ref{th1.5} and the identity (\ref{1.6}). Besides, we immediately obtain an analog of (\ref{GKM}).

\begin{cor} We have
	$$S(R_M)=s(|R_M|)(1+o(1)), \quad M\to\infty, \, M\in \O.
	$$
\end{cor}

Let $K_l$ be the number of gaps of length $l$ in $R_M$, that is, the number of elements $x\in R_M$ such that $x+l\in R_M$ but $x+j\notin R_M$ for $1\leq j\leq l-1$. To deal with $S(R_M)$, we write it as $\sum_{l\geq1}K_ll^2$ and, choosing a parameter $D$, divide the range of summation into small ($l\leq D$) and large ($l>D$) gaps. Note that for computing the contribution of small gaps in the setting of Theorem \ref{th1.5} one may use the limit distribution (\ref{KR}) of the distances between quadratic residues. Nevertheless, in order to be able to compute the second term in the asymptotics for $S(R_M)$, one should avoid $o(1)$-terms which occur in (\ref{KR}); so we use another technique (and consider special moduli). Finally, we note that our method seems to be insufficient to prove that the inequality $S(R_M)>s(|R_M|)$ holds for a large set of $M$; nevertheless, we believe that the set of such $M$ also has positive density.

In Sections \ref{Sec2}-\ref{Sec5} we will be simultaneously proving Theorems \ref{th1.3} and \ref{th1.5}. Section \ref{Sec2} is devoted to general results concerning additive configurations related to quadratic residues. In Section \ref{Sec3} we find the asymptotics for the contribution of small gaps into $S(R_M)$, and in Section \ref{Sec4} we show that the contribution of large gaps is asymptotically small. In Section \ref{Sec5} we finish the proofs of both Theorems \ref{th1.3} and \ref{th1.5}. Further, in Section \ref{Sec6} we prove Theorem \ref{th1.4}. Finally, in Section \ref{Sec7} we prove for completeness the identity (\ref{1.3}).

\section{Additive configurations related to quadratic residues}\label{Sec2}

We begin with a classical estimate for complete character sums.

\begin{lem}\label{lem2.1}
	Let $p$ be a prime, $\chi_1,\ldots,\chi_r$ be characters modulo $p$ with some $\chi_i\neq\chi_0$, where $\chi_0$ is the principal character, and $a_1,\ldots,a_r$ be distinct elements of $\mathbb{Z}_p$. Then
	$$\left|\sum_{x\in\mathbb{Z}_p}\chi_1(x+a_1)\chi_2(x+a_2)\ldots\chi_r(x+a_r)\right|\leq r\sqrt p.
	$$
\end{lem}

\begin{proof}
See \cite{Joh}, Lemma 1.
\end{proof}

To prove our Theorems \ref{th1.3} and \ref{th1.5}, we will need estimates for the number of some additive configurations related to quadratic residues. Let $M$ be a positive integer; for disjoint subsets $C_1$ and $C_2$ of $\mathbb{Z}_M$ we define the sets
\[
R_M(C_1)=\{x\in\mathbb{Z}_M: x+c_1\in R_M \mbox{ for all $c_1\in C_1$} \},
\]
\[
N_M(C_2)=\{x\in\mathbb{Z}_M: x+c_2\notin R_M   \mbox{ for all $c_2\in C_2$} \},
\]
and
\begin{equation*
RN_M(C_1,C_2)=R_M(C_1)\cap N_M(C_2)
\end{equation*}
(we do not exclude the cases where $C_1$ or $C_2$ is empty). The following two lemmas provide bounds for these quantities in the cases of prime and square-free moduli. For $x\in\mathbb{Z}_M$, we denote by $R_M(x)$ the indicator $\mathbb{I}(x\in R_M)$; $\theta$ is used for a number bounded in magnitude by $1$, which can vary from line to line.

\begin{lem}\label{lem2.2}
	Let $p$ be a prime. Then
	\begin{equation*}
	|RN_p(C_1,C_2)|= p2^{-(|C_1|+|C_2|)}+0.5\theta(|C_1|+|C_2)|(p^{1/2}+1).
	\end{equation*}
\end{lem}

\begin{proof} We use the method from \cite{Dav}. Let $\chi$ be the Legendre symbol modulo $p$ and $1_p(x)$ be the indicator $\mathbb{I}(x=0)$. Then $R_p(x)=\frac12(1+\chi(x)+1_p(x))$
and $1-R_p(x)=\frac12(1-\chi(x)-1_p(x))$. Since the functions $R_p(x)$, $1-R_p(x)$, and $\frac{1\pm\chi(x+c_i)}{2}$ are bounded in magnitude by $1$, we have
\begin{multline*}
|RN_p(C_1,C_2)|=\sum_{x\in\mathbb{Z}_p}\prod_{c_1\in C_1}R_p(x+c_1)\prod_{c_2\in C_2}(1-R_p(x+c_2))=\\
\sum_{x\in\mathbb{Z}_p}\prod_{c_1\in C_1}\frac{1+\chi(x+c_1)+1_p(x+c_1)}{2}\prod_{c_2\in C_2}\frac{1-\chi(x+c_2)-1_p(x+c_2)}{2}=\\
2^{-(|C_1|+|C_2|)}\sum_{x\in\mathbb{Z}_p}\prod_{c_1\in C_1}(1+\chi(x+c_1))\prod_{c_2\in C_2}(1-\chi(x+c_2))+0.5\theta(|C_1|+|C_2|).
\end{multline*}
Denote $k=|C_1|+|C_2|$. Expanding the product in the last expression, we obtain the term $p$ and $2^k-1$ character sums; there are exactly ${k\choose r}$ of them of the type $\sum_{x\in\mathbb{Z}_p}\prod_{j=1}^r\chi(x+a_j)$, each of which is bounded by $rp^{1/2}$ by Lemma \ref{lem2.1}. Since $\sum_{r=1}^k{k\choose r}r=k2^{k-1}$, we get
\begin{equation*}
2^k|RN_p(C_1,C_2)|=p+0.5\theta k2^k(p^{1/2}+1),
\end{equation*}
and the claim follows.
\end{proof}

\begin{lem}\label{lem2.3}
Let $m$ be a square-free positive integer, $p_1=P^-(m)$ and $t=\o(m)$, and let $C_1$ and $C_2$ be disjoint subsets of non-negative integers less than $p_1$, $|C_1|=s_1$, $|C_2|=s_2$. Then for all $r\in\mathbb{N}$ we have
	\begin{equation*}
	|RN_m(C_1,C_2)|=m2^{-s_1t}(1-2^{-t})^{s_2}+O(m2^{-s_1t}\left(E_1+E_2\right)),
	\end{equation*}
	where
$$E_1= \sum_{k=0}^{\min\{r,s_2\}}{s_2\choose k}2^{-kt}E(s_1+k),
$$
$$E(\nu)=\sum_{d|m,d>1}\frac{\left(\nu2^{\nu}\right)^{\o(d)}}{d^{1/2}}, \quad 	E_2=\sum_{k=r}^{s_2}{s_2\choose k}2^{-kt}
$$
(we adopt the convention that $E_2=0$ for $r>s_2$).
\end{lem}

Note that the probability that a random element $x\in \mathbb{Z}_m$ is a quadratic residue is about $2^{-t}$. The main term $m2^{-s_1t}(1-2^{-t})^{s_2}$ in this lemma agrees with this, also telling us that the events $x+c_1\in R_m$, $c_1\in C_1$, and $x+c_2\notin R_m$, $c_2\in C_2$, are nearly independent, provided that the error terms are small.

\begin{proof} We may assume $r\leq s_2+1$. For a subset $C\subseteq C_2$, we write for the brevity
$$R_C:=R_m(C_1\cup C)=\{x\in\mathbb{Z}_m: x+c\in R_m \mbox{ for all $c\in C_1\cup C$}\}.
$$
Note that $|RN_m(C_1,C_2)|=\left|R_{\varnothing}\setminus\bigcup_{c\in C_2}R_{\{c\}}\right|=|K_{\varnothing}|-|\bigcup_{c\in C_2}R_{\{c \}}|$.
Then by the inclusion-exclusion principle
\begin{equation}\label{2.2}
|RN_m(C_1,C_2)|=|R_{\varnothing}|-\sum_{|C|=1}|R_{C}|+\sum_{|C|=2}|R_C|-\ldots+(-1)^{s_2}|R_{C_2}|,
\end{equation}
and by the Bonferroni inequalities (see, for example, \cite{TV}, Ch.1, Exercise 1.1.3) we have
\begin{equation}\label{2.3}
|RN_m(C_1,C_2)|=|R_{\varnothing}|-\sum_{|C|=1}|R_C|+\sum_{|C|=2}|R_C|-\ldots+(-1)^{r-1}\sum_{|C|=r-1}|R_C|+O\left(\sum_{|C|=r}|R_C|\right)
\end{equation}
for all $r\leq s_2$.

Now we fix $C\subseteq C_2$ and work with $|R_C|$. By multiplicity and the assumption that elements of $C_1\bigsqcup C_2$ are less than $p_1$ (and, hence, are distinct modulo $p$ for any $p|m$), we can rewrite it as
\begin{equation*}
|R_C|=\prod_{p|m}|R_p(C_1\cup C)|.
\end{equation*}
Let $|C|=k$. By Lemma \ref{lem2.2}, $|R_p(C_1\cup C)|=p2^{-(s_1+k)}\left(1+\theta(s_1+k)2^{s_1+k}p^{-1/2}\right)$ for some $|\theta|\leq1$. Thus, we see that
\begin{equation}\label{2.4}
|R_C|=m2^{-(s_1+k)t}+O\left(m2^{-(s_1+k)t}E(s_1+k)\right),
\end{equation}
with $E(s_1+k)$ defined in the statement of the lemma.

We first consider the case $r=s_2+1$. By (\ref{2.2}) and (\ref{2.4}),
\begin{multline*}
2^{s_1t}m^{-1}|RN_m(C_1,C_2)|=1-{s_2\choose 1}2^{-t}+{s_2\choose 2}2^{-2t}+\ldots  +(-1)^{s_2}2^{-s_2t}+\\+O\left(\sum_{k=0}^{s_2}{s_2\choose k}2^{-kt}E(s_1+k)\right)=
(1-2^{-t})^{s_2}+O(E_1),
\end{multline*}
and the claim follows.

Now suppose that $r\leq s_2$. From (\ref{2.3}) and (\ref{2.4}) we obtain
\begin{multline*}
2^{s_1t}m^{-1}RN_m(C_1,C_2)=1-{s_2\choose 1}2^{-t}+{s_2\choose 2}2^{-2t}+\ldots +(-1)^{r-1}{s_2\choose {r-1}}2^{-(r-1)t}+\\
+O\left(\sum_{k=0}^{r}{s_2\choose k}2^{-kt}E(s_1+k)+{s_2\choose {r}}2^{-rt}\right)=
(1-2^{-t})^{s_2}+E,
\end{multline*}
where
\[
E\ll E_1+{s_2\choose {r}}2^{-rt}+\left|\sum_{k=r}^{s_2}(-1)^k{s_2 \choose k}2^{-kt}\right| \ll E_1+\sum_{k=r}^{s_2}{s_2 \choose k}2^{-kt}=E_1+E_2,
\]
as desired. This concludes the proof.
\end{proof}

\begin{cor}\label{cor2.4}
Let $m$ be a square-free positive integer with $p_1=P^-(m)$ large enough, $t=\o(m)$, and let $C_1=\{0,l\}$ and $C_2$ be disjoint subsets of non-negative integers less than $p_1$. Suppose that $s=|C_2|\leq 0.13\cdot2^t\log p_1$ and $t\leq p_1^{0.01}$. Then
\[ |RN_m(C_1,C_2)|=m2^{-2t}(1-2^{-t})^s+O(m2^{-2t}(1+2^{1-t})^sp_1^{-0.48}),
\]
\end{cor}

\begin{proof} We set $r=[\frac{0.49\log p_1}{\log2}]$ in the previous lemma and estimate the quantities $E(k+2)$ for $0\leq k\leq r$.
Since then $t(k+2)2^{k+2}<p_1^{1/2}$, we have
\begin{multline*}
E(k+2)=\sum_{d|m,d>1}\frac{((k+2)2^{k+2})^{\o(d)}}{d^{1/2}} \leq \sum_{l=1}^t\left(\frac{(k+2)2^{k+2}}{p_1}\right)^l{t\choose l} = \\ \left(1+(k+2)2^{k+2}p_1^{-1/2}\right)^t-1 \leq
\exp(t(k+2)2^{k+2}p_1^{-1/2})-1 \ll tr2^kp_1^{-1/2}.
\end{multline*}
We thus see that
$$E_1=\sum_{k=0}^{\min\{r,s\}}{s \choose k}2^{-kt}E(k+2) \ll trp_1^{-1/2}\sum_{k=0}^s{s\choose k}2^{(1-t)k} = (1+2^{1-t})^strp_1^{-1/2} \ll (1+2^{1-t})^sp_1^{-0.48}.
$$
Further,
$${s \choose k}2^{-kt}\leq \frac{s^k}{k!2^{kt}} \ll   \left(\frac{se}{k2^t}\right)^k.
$$
Besides, we have $se\leq 0.13e2^t\log p_1 < 0.5r2^t$, and therefore
$$E_2=\sum_{k=r}^s{s \choose k}2^{-kt}  \ll \sum_{k=r}^{\infty}\left(\frac{se}{r2^t}\right)^k \ll 2^{-r} \ll  p_1^{-0.49}.
$$
The claim now follows from Lemma \ref{lem2.3}.
\end{proof}

\begin{cor}\label{cor2.5}
	Let $m$ be a square-free positive integer with $p_1=P^-(m)$ large enough, $t=\o(m)$, and let $C$ be a subset of non-negative integers less than $p_1$. Suppose that $|C|\geq \a2^t\log p_1$ for some fixed $\a\in(0,0.13)$ and $t\leq p_1^{0.01}$. Then
	\[ |N_m(C)| \ll mp_1^{-\a}.
	\]
\end{cor}

\begin{proof} Let $s=|C|$; since $N_m(C')\leq N_m(C'')$ whenever $C'\supseteq C''$, we may assume that $s=[\a2^t\log p_1]+1$. Arguing as in the proof of the previous lemma and using the inequality $(1+x)\leq e^x$, $x\in{\mathbb R}$, we get from Lemma \ref{lem2.3}
$$m^{-1}|N_m(C)| \ll (1-2^{-t})^s+(1+2^{1-t})^sp_1^{-0.48}\ll p_1^{-\a},
$$
as desired.
\end{proof}

\section{Contribution of small gaps}\label{Sec3}

For $2\leq j\leq A$ and $i\in \mathbb{Z}_A$, we set $I(i,j)=\{i+1,...,i+j-1\}\pmod{A}$ and
\begin{equation}\label{3.1}
r(i,j)=\#\big(I(i,j)\cap R_A\big);
\end{equation}
also we adopt the convention that $r(i,1)=0$. Obviously, $0\leq r(i,j)\leq |R_A|$.

\begin{lem}\label{lem3.1}
Let $D\geq1$, $(A,m)=1$, $M=Am$, and $D_0=D|R_A|A^{-1}+|R_A|$. Suppose that for any $s\leq D_0$ we have an estimate of the type
$$\sum_{u\in\mathbb{Z}_{m}} R_m(u)R_m(u+l)\prod_{j=1}^{s}(1-R_m(u+\nu_j))=m2^{-2t}(1-2^{-t})^s+O(L(m,s))
$$
uniformly on distinct positive integers $\nu_1,\ldots,\nu_s,l$ less than $D$.
Then
\begin{equation}\label{3.2}
\sum_{l\leq D} K_ll^2 =
m2^{-2t}\sum_{j=1}^{A}\sum_{i\in R_A\cap(R_A-j)}\sum_{Ak+j\leq D}y_t^{|R_A|k+r(i,j)}(Ak+j)^2+O(E(M,D_0)),
\end{equation}
where
$$E(M,D_0)=A^3|R_A|^{-1}\sum_{s\leq D_0}L(m,s)s^2 + A^3|R_A|\sum_{s=0}^{|R_A|}L(m,s).
$$
\end{lem}

\begin{proof} Since $(A,m)=1$, we have $R_M(x)=R_A(x)R_m(x)$, and hence
\begin{equation}\label{3.3}
R_M(x)=\begin{cases}
R_m(x), \, x\in R_A,\\
0, \, \mbox{otherwise.}
\end{cases}
\end{equation}
Fix $l\leq D$. We can write
\begin{equation*}
K_l=\sum_{x\in\mathbb{Z}_M}R_M(x)R_M(x+l)\prod_{\nu=1}^{l-1}(1-R_M(x+\nu))=\sum_{i=0}^{A-1}K_l(i),
\end{equation*}
where
\begin{equation*}
K_l(i)=\sum_{x\in\mathbb{Z}_M: x\equiv i\pmod{A}}R_M(x)R_M(x+l)\prod_{\nu=1}^{l-1}(1-R_M(x+\nu)).
\end{equation*}
Further, using (\ref{3.3}) we get
\begin{equation}\label{3.4}
K_l=\sum_{i\in R_A\cap(R_A-l)}K_l(i).
\end{equation}
Now fix  $i\in R_A\cap(R_A-l)$. Let $\nu_1,\ldots,\nu_s$ be the numbers $\nu\in\{1,\ldots,l-1\}$ such that $(\nu+i)\pmod{A}\in R_A$. We write $l=Ak+j$ for some $k\geq0$ and $1\leq j\leq A$; then $k=\lfloor\frac{l}{A}\rfloor$ for $1\leq j\leq A-1$ and $k=\frac{l}{A}-1$ otherwise. Therefore, $\frac{l}{A}-1\leq k\leq \frac{l}{A}$ and
\begin{equation}\label{3.5}
s=s(l,i,A)=|R_A|k+r(i,j) = l|R_A|A^{-1}+\theta'|R_A|
\end{equation}
for some $\theta'\in[-1,1]$. In particular, $s(l,i,A)\leq D|R_A|A^{-1}+|R_A|=D_0$ for any $l\leq D$. Taking into account (\ref{3.3}) and recalling that $(A,m)=1$, we obtain from the assumption of the lemma
\begin{multline*}
K_l(i)=\sum_{x\in\mathbb{Z}_M: x\equiv i\pmod{A}}R_m(x)R_m(x+l)\prod_{j=1}^{s}(1-R_m(x+\nu_j))=\\
\sum_{u\in\mathbb{Z}_{m}} R_m(Au+i)R_m(Au+i+l)\prod_{j=1}^{s}(1-R_m(Au+i+\nu_j))=\\
\sum_{u\in\mathbb{Z}_{m}} R_m(u)R_m(u+l)\prod_{j=1}^{s}(1-R_m(u+\nu_j))=m2^{-2t}(1-2^{-t})^{s(l,i,A)}+O(L(m,s(l,i,A))).
\end{multline*}
Since $l\equiv j\pmod{A}$, from here, (\ref{3.4}), and (\ref{3.5}) we have
\begin{equation*}\label{N_l}
K_l=m2^{-2t}\sum_{i\in R_A\cap(R_A-j)}y_t^{|R_A|k+r(i,j)}+O\left(\sum_{i\in R_A\cap (R_A-l)}L(m,s(l,i,A))\right),
\end{equation*}
where we write $y_t=1-2^{-t}$ for the brevity. Summing this estimate over $l\leq D$, we get
\begin{equation*}
\sum_{l\leq D} K_ll^2 =
m2^{-2t}\sum_{j=1}^{A}\sum_{i\in R_A\cap(R_A-j)}\sum_{Ak+j\leq D}y_t^{|R_A|k+r(i,j)}(Ak+j)^2+O(E(M,D_0)),
\end{equation*}
where
$$
E(M,D_0)=\sum_{l\leq D}\sum_{i\in R_A\cap (R_A-l)}L(m,s(l,i,A))l^2.
$$
Note that (\ref{3.5}) implies $l\leq sA|R_A|^{-1}+A$. Further, if $l$ runs from $Ak+1$ to $Ak+A$ and $i$ runs over $R_A\cap (R_A-l)$, then $s(l,i,A)$ takes values from  $|R_A|k$ to $|R_A|k+|R_A|$. Hence, for any $s\geq0$ there are at most $O(A|R_A|)$ pairs $(l,i)$ with $s(A,l,i)=s$. Thus,
\begin{multline*}
E(M,D_0)\ll \sum_{s=0}^{D_0}L(m,s)A|R_A|\left(sA|R_A|^{-1}+A\right)^2 \ll \\ A^3|R_A|^{-1}\sum_{s\leq D_0}L(m,s)s^2 + A^3|R_A|\sum_{s=0}^{|R_A|}L(m,s).
\end{multline*}
This concludes the proof.
\end{proof}

Let $\varepsilon_0\in(0,1)$ be a parameter to be chosen later.

\begin{lem}\label{lem3.2}
	Let $M=Ap$ with prime $p$ and $(A,p)=1$. Suppose that $A\leq 0.1p^{\varepsilon_0}$ and let $D=A|R_A|^{-1}p^{\varepsilon_0}$. Then
	$$E(M,D_0)\ll A^3|R_A|^{-1}p^{1/2+4\varepsilon_0}.
	$$
\end{lem}

\begin{proof} Lemma \ref{lem2.2} gives us the assumption of Lemma \ref{lem3.1} with $L(p,s)=(s+2)p^{1/2}$; also we have $D_0=p^{\varepsilon_0}+|R_A|$. Thus,
\begin{multline*}
E(M,D_0)\ll A^3|R_A|^{-1}p^{1/2}\sum_{s\leq D_0}s^3+A^3|R_A|p^{1/2}\sum_{s\leq |R_A|}s\ll A^3|R_A|^{-1}p^{1/2}\left(D_0^4+|R_A|^4\right) \ll \\
A^3|R_A|^{-1}p^{1/2}\left(p^{4\varepsilon_0}+|R_A|^4\right) \ll A^3|R_A|^{-1}p^{1/2+4\varepsilon_0},
\end{multline*}
as desired.
\end{proof}

Let $M\in\O$ and $M=Am$. In what follows we will often use the fact that the condition (iii) in the definition of $\O$ implies $A\leq 2^t$, where $t=\o(m)$.

\begin{lem}\label{lem3.3}
Let $M\in \O$ and $D=0.13 A|R_A|^{-1}2^t\log p_1$. Then
$$E(M,D_0)\ll MA^22^tp_1^{-0.21}.
$$
\end{lem}

\begin{proof} Note that $M\in\O$ implies $t\leq p_1^{0.01}$. Then Corollary \ref{cor2.4} gives us the assumption of Lemma \ref{lem3.1} with
$$L(m,s) = m2^{-2t}(1+2^{1-t})^sp_1^{-0.48};
$$
set $D_0=0.13\cdot2^t\log p_1+|R_A|$. Also $M\in\O$ gives us $|R_A|\leq A\leq 2^t$. Then
\begin{equation*}
m^{-1}\sum_{s=0}^{|R_A|}L(m,s) \ll p_1^{-0.48}2^{-t}\exp(O(|R_A|2^{-t})) \ll p_1^{-0.48}2^{-t}.
\end{equation*}
Further,
\begin{multline*}
m^{-1}\sum_{s\leq D_0}L(m,s)s^2\ll 2^{-2t}p_1^{-0.48}\sum_{s\leq D_0}s^2(1+2^{1-t})^{s}\ll 2^{-2t}p_1^{-0.48}D_0^2\sum_{s\leq D_0}(1+2^{1-t})^s\ll\\
2^{-t}p_1^{-0.48}D_0^2(1+2^{1-t})^{D_0} \ll 2^{t}p_1^{-0.21}.
\end{multline*}
Then
\[
E(M,D_0)\ll mA^3|R_A|^{-1}\left(2^tp_1^{-0.21}+|R_A|^22^{-t}p_1^{-0.48} \right) \ll MA^22^tp_1^{-0.21}.
\]
The claim follows.
\end{proof}

To work with the main term in (\ref{3.2}), we introduce the function
\begin{equation}\label{3.6}
g_A(y)=\sum_{j=1}^{A}\sum_{i\in R_A\cap(R_A-j)}y^{r(i,j)}\sum_{k\geq 0}y^{|R_A|k}(Ak+j)^2, \quad y\in(0,1);
\end{equation}
in the next section we will simplify this expression. Note that $g_A$ does not depend on $m$.

Now we need an estimate for the ``tail'' of $g_A$. We set
\begin{equation*}
G_A(y,D)=\sum_{j=1}^{A}\sum_{i\in R_A\cap(R_A-j)}y^{r(i,j)}\sum_{Ak+j>D}y^{|R_A|k}(Ak+j)^2.
\end{equation*}

\begin{lem}\label{lem3.4}
	For any $N\geq2$ and $x\in (0,1)$, we have
	\begin{equation*}
	\sum_{n\geq N}n^2x^{n}\ll  x^N(1-x)^{-1}(N^2+(1-x)^{-2}).
	\end{equation*}
\end{lem}

\begin{proof} We may assume that $N$ is a positive integer. For $x\in(0,1)$, we have $\sum_{n\geq N}x^n=x^N(1-x)^{-1}$;
by differentiating it twice, we obtain
\begin{multline*}
\sum_{n\geq N}n(n-1)x^{n-2}=N(N-1)x^{N-2}(1-x)^{-1}+2Nx^{N-1}(1-x)^{-2}+2x^N(1-x)^{-3} \\
\ll x^{N-2}(1-x)^{-1}\left(N^2+(1-x)^{-2}\right).
\end{multline*}
Since $n^2\ll n(n-1)$ for $n\geq2$, we have
$$	\sum_{n\geq N}n^2x^{n} \ll x^2\sum_{n\geq N}n(n-1)x^{n-2} \ll x^N(1-x)^{-1}\left(N^2+(1-x)^{-2}\right),
$$
and the claim follows.
\end{proof}

\begin{lem}\label{lem3.5}
	Let $M=Ap$ with prime $p$ and $(A,p)=1$. Suppose that $A\leq 0.1p^{\varepsilon_0}$ and let $D=A|R_A|^{-1}p^{\varepsilon_0}$. Then
	$$G_A(0.5,D)\ll 2^{-0.8p^{\varepsilon_0}}.
	$$
\end{lem}

\begin{proof} First, we note that
\[
\sum_{j=1}^A\sum_{i\in R_A\cap(R_A-j)}1\leq \sum_{i=1}^AR_A(i)\sum_{j=1}^AR_A(i+j)=|R_A|^2.
\]
Hence,
\begin{multline*}
G_A(0.5,D)\leq \sum_{j=1}^A\sum_{i\in R_A\cap(R_A-j)}\max_{1\leq j\leq A}\sum_{Ak+j>D}2^{-|R_A|k}(Ak+j)^2 \leq \\
|R_A|^2\sum_{Ak+A>D}2^{-|R_A|k}(Ak+A)^2 \ll A^2|R_A|^2\sum_{k>|R_A|^{-1}p^{\varepsilon_0}-1}k^22^{-|R_A|k}.
\end{multline*}
Using the previous lemma, we get from here
\begin{multline*}
G_A(0.5,D)\ll A^2|R_A|^22^{-|R_A|(p^{\varepsilon_0}|R_A|^{-1}-1)}\left(|R_A|^{-2}p^{2\varepsilon_0}+1\right) \ll A^2p^{2\varepsilon_0}2^{-0.9p^{\varepsilon_0}} \ll 2^{-0.8p^{\varepsilon_0}}.
\end{multline*}
This concludes the proof.
\end{proof}

\begin{lem}\label{lem3.6}
	Let $M\in \O$ and $D=0.13A|R_A|^{-1}2^t\log p_1$. Then
	\begin{equation*}
	G_A(y_t,D) \ll A^42^{3t}p_1^{-0.12}.
	\end{equation*}
\end{lem}

\begin{proof} As before,
\begin{multline*}
G_A(y_t,D)\leq \sum_{j=1}^A\sum_{i\in R_A\cap(R_A-j)}\max_{1\leq j\leq A}\sum_{Ak+j>D}y_t^{|R_A|k}(Ak+j)^2 \leq \\
|R_A|^2\sum_{Ak+A>D}y_t^{|R_A|k}(Ak+A)^2 \ll A^2|R_A|^2\sum_{k>D/A-1}k^2y_t^{|R_A|k}.
\end{multline*}
Further, $y_t^{|R_A|}\leq y_t = 1-2^{-t}$, and hence $(1-y_t^{|R_A|})^{-1}\leq 2^t$; also $1/2\leq y_t<1$. Using Lemma \ref{lem3.4} and the inequality $1-x\leq e^{-x}$, we see that
\begin{multline*}
G_A(y_t,D)\ll A^2|R_A|^2 y_t^{|R_A|D/A}2^t\left( D^2A^{-2}+2^{2t}\right) \ll \\ A^2|R_A|^22^{3t}(|R_A|^{-2}(\log p_1)^2+1)p_1^{-0.13} \ll A^42^{3t}p_1^{-0.12}.
\end{multline*}
This concludes the proof.
\end{proof}

Combining the results of this section, we obtain the following two corollaries.

\begin{cor}\label{cor3.7}
	Let $M=Ap$ with prime $p$ and $(A,p)=1$. Suppose that $A\leq 0.1p^{\varepsilon_0}$ and let $D=A|R_A|^{-1}p^{\varepsilon_0}$. Then
$$\sum_{l\leq D}K_ll^2 = 0.25g_A(0.5)p+O(A^3p^{1/2+4\varepsilon_0}).
$$
\end{cor}

\begin{proof} Set $D_0=p^{\varepsilon_0}+|R_A|$. By Lemmas \ref{lem3.1}, \ref{lem3.2} and \ref{lem3.5}, we have
\begin{multline*}
\sum_{l\leq D}K_ll^2 = 0.25p\left(g_A(0.5)+G_A(0.5,D)\right)+O(E(M,D_0)) = \\ 0.25g_A(0.5)p+O(A^3|R_A|^{-1}p^{1/2+4\varepsilon_0}+p2^{-0.8p^{\varepsilon_0}}).
\end{multline*}
The claim follows.
\end{proof}

\begin{cor}\label{cor3.8}
Let $M\in \O$ and $D=0.13 A|R_A|^{-1}2^t\log p_1$. Then
$$\sum_{l\leq D}K_ll^2 = m2^{-2t}g_A(y_t) + O(MA^32^tp_1^{-0.12}).
$$
\end{cor}

\begin{proof}
Set $D_0=0.13\cdot2^t\log p_1+|R_A|$. By Lemmas \ref{lem3.1}, \ref{lem3.3} and \ref{lem3.6}, we have
\begin{multline*}
\sum_{l\leq D}K_ll^2 = m2^{-2t}\left(g_A(y_t)+G_A(y_t,D)\right)+O(E(M,D_0)) = \\ m2^{-2t}g_A(y_t)+O(MA^32^tp_1^{-0.12}+MA^22^tp_1^{-0.21}).
\end{multline*}
The claim follows.
\end{proof}

\section{Contribution of large gaps}\label{Sec4}

To bound the contribution of large gaps, we will need Burgess' estimates for character sums.

\begin{lem}\label{lem4.1}
	Let $k$ be square-free, $\chi$ be a nonprincipal character modulo $k$, and $\varepsilon>0$, $r\in\mathbb{N}$. Then for any $x\in\mathbb{Z}_k$ and $h\in \mathbb{N}$ we have
	$$\left|\sum_{i=0}^{h-1}\chi(x+i)\right|\ll_{r,\varepsilon}h^{1-1/r}k^{(r+1)/4r^2+\varepsilon}.
	$$
\end{lem}

\begin{proof} See \cite{burg1},  Theorem 2.
\end{proof}

\begin{lem}\label{lem4.2}
	Let $k$ be square-free, $\chi$ be a primitive character modulo $k$, and $\varepsilon>0$, $r\in\mathbb{N}$. Then
	$$\sum_{x\in\mathbb{Z}_k}\left|\sum_{i=0}^{h-1}\chi(x+i)\right|^{2r} \ll_{r,\varepsilon} kh^r+k^{1/2+\varepsilon}h^{2r}.
	$$
\end{lem}

\begin{proof} See \cite{burg2}, Lemma 8.
\end{proof}

Let $H\in\mathbb{N}$. For the brevity, we denote by $N_M(H)=N_M(\{0,1,...,H-1\})$ the set of $x\in\mathbb{Z}_M$ such that all elements $x,x+1,\ldots,x+H-1$ are quadratic non-residues modulo $M$.

\begin{lem} \label{lem4.3}
	Let $M$ be a square-free positive integer. Suppose that $h\geq (2C\tau(M))^{50}$,
where $C$ is an absolute constant. Then
	$$|N_M(h)| \ll M\tau(M)^4h^{-1.47}.
	$$
\end{lem}

\begin{proof} Let $\L$ be the set of nontrivial real characters modulo $M$, and $\chi_p$ be the Legendre symbol modulo $p$. For any quadratic nonresidue $x\in\mathbb{Z}_M$ we have
	$$
	\sum_{\chi\in\L}\chi(x)=\prod_{p|M}(1+\chi_p(x))-1=-1;
	$$
	therefore, for any $x\in N_M(h)$
	\begin{equation}\label{4.1}
	\left|\sum_{\chi\in\L}\sum_{i=0}^{h-1}\chi(x+i)\right|=h.
	\end{equation}
	Denote by $k(\chi)$ the conductor of $\chi\in\L$ and set
	$$\L_1=\{\chi\in\L: k(\chi)\leq h^3\}, \quad \L_2=\{\chi\in L: k(\chi)>h^3\}.
	$$
	We want to show that
	\begin{equation}\label{4.2}
	\sum_{\chi\in\L_1}\left|\sum_{i=0}^{h-1}\chi(x+i)\right|\leq h/2.
	\end{equation}
	For any $\chi\in\L_1$ Lemma \ref{lem4.1} implies
	$$\left|\sum_{i=0}^{h-1}\chi(x+i)\right|\ll_{r,\varepsilon} h^{1-1/r}h^{3(1/4r+1/4r^2+\varepsilon)}= h^{1-(1-3/r)/(4r)+3\varepsilon}.
	$$
	Taking $r=6$ and $\varepsilon=0.0001$, we obtain
	$$\sum_{\chi\in\L_1}\left|\sum_{i=0}^{h-1}\chi(x+i)\right|\leq Ch^{1-1/48+0.0003}|\L_1|\leq Ch^{0.98}|\L_1|
	$$
	for some absolute constant $C>0$, and (\ref{4.2}) follows from $2C|\L_1|\leq 2C|\L|=2C\tau(M)\leq h^{0.02}$.

From (\ref{4.1}) and (\ref{4.2}) we see that
	$$\sum_{\chi\in\L_2}\left|\sum_{i=0}^{h-1}\chi(x+i)\right|\geq h/2,
	$$
	for all $x\in N_M(h)$, and thus
	\begin{equation}\label{4.3}
	\sum_{x\in Y(h)}\sum_{\chi\in\L_2}\left|\sum_{i=0}^{h-1}\chi(x+i)\right|\geq |Y(h)|h/2.
	\end{equation}
	On the other hand, for any $\chi\in\L_2$ we have from Lemma \ref{lem4.2} (with $\varepsilon=0.01$)
	\begin{equation*}
	\sum_{x\in Y(h)}\left|\sum_{i=0}^{h-1}\chi(x+i)\right|^{2r}\leq \sum_{x\in\mathbb{Z}_M}\left|\sum_{i=0}^{h-1}\chi(x+i)\right|^{2r}\ll \frac{M}{k_{\chi}} \sum_{x\in \mathbb{Z}_{k_{\chi}}}\left|\sum_{i=0}^{h-1}\chi(x+i)\right|^{2r}\ll_r Mh^r+Mh^{2r}k_{\chi}^{-0.49}.
	\end{equation*}
	Setting $r=2$, we have for all $\chi\in\L_2$
	$$\sum_{x\in Y(h)}\left|\sum_{i=0}^{h-1}\chi(x+i)\right|^{4} \ll Mh^{2.53}.
	$$
	Using H\"older's inequality, we find from (\ref{4.3})
	$$
	|Y(h)|h\ll \sum_{\chi\in\L_2}\sum_{x\in Y(h)}\left|\sum_{i=0}^{h-1}\chi(x+i)\right| \ll \tau(M)|Y(h)|^{3/4}(Mh^{2.53})^{1/4},
	$$
	and hence,
	$$
	|Y(h)|\ll M\tau(M)^4h^{-1.47}.
	$$
	The claim follows.
\end{proof}

\begin{rem}\label{rem4.4}
	One can deduce from this argument that the size of maximal gap in the set $R_M$ is at most $O(M^{1/3})$, but we do not need it here.
\end{rem}

Let $R_M=\{0=r_1<r_2<\ldots<r_{|R_M|}\}$ be the set of quadratic residues modulo $M$.

\begin{lem}\label{lem4.5}
	We have
	$$
	\sum_{2^{\nu+1}\leq r_i-r_{i-1} < 2^{\nu+2} }(r_i-r_{i-1})^2 \leq |N_M(2^\nu)|2^{\nu+4}.
	$$
\end{lem}

\begin{proof} If $r_i-r_{i-1}\geq 2^{\nu+1}$, then for all $1\leq k\leq 2^{\nu}$ the elements $r_{i-1}+k$ belong to $N_M(2^{\nu})$. Therefore the number of such $i$ is at most $2^{-\nu}|N_M(2^{\nu})|$. The claim follows.
\end{proof}

\begin{lem}\label{lem4.6}
Let $M=Ap$, where $(A,p)=1$, $p$ is large enough, and $A\leq p^{\varepsilon_0}$. Then for $D=A|R_A|^{-1}p^{\varepsilon_0}$ we have
$$\sum_{l>D}K_ll^2 \ll A^2p^{1-0.47\varepsilon_0}.
$$
\end{lem}

\begin{proof} Using the previous lemma, we obtain
\[
\sum_{l>D}K_ll^2 \leq \sum_{\nu\,:\, 2^{\nu}>D}\sum_{2^{\nu-1}\leq u_j-u_{j-1}<2^{\nu}}(u_j-u_{j-1})^2 \ll \sum_{\nu\,:\, 2^{\nu}>D}|N_M(2^{\nu-2})|2^{\nu}.
\]
Note that since $p$ is large enough, then for some absolute constant $C_1$ we have $(4C\tau(A))^{50}\leq C_1A^{1/2} \leq 0.25p^{\varepsilon_0} \leq 0.25D$; thus $2^{\nu}>D$ implies $2^{\nu-2}>(2C\tau(M))^{50}$. Then Lemma \ref{lem4.3} gives us
\[
\sum_{l>D}K_ll^2\ll 2^{4\o(M)}M\sum_{\nu\,:\, 2^{\nu+2}>D}2^{-0.47\nu} \ll \tau(A)^4ApD^{-0.47} \ll A^2p^{1-0.47\varepsilon_0},
\]
as desired.
\end{proof}

\begin{cor}\label{cor4.7}
	Let $M=Ap$ with large enough prime $p$ and $(A,p)=1$. Suppose that $A\leq 0.1p^{0.1}$. Then
	$$S(R_M) = 0.25g_A(0.5)p+O(A^3p^{0.95}).
	$$
\end{cor}

\begin{proof} Set $D=A|R_A|^{-1}p^{\varepsilon_0}$ with some $\varepsilon_0\geq0.1$. By Corollary \ref{cor3.7} and Lemma \ref{lem4.6}, we have
	$$S(R_M)=\sum_{l\leq D}K_ll^2+\sum_{l>D}K_ll^2 = 0.25g_A(0.5)p+O(A^3p^{1/2+4\varepsilon_0}+A^2p^{1-0.47\varepsilon_0})	.
	$$
	Now we set $\varepsilon_0=\frac{1}{2\cdot4.47}=0.11...$. Then $1-0.47\varepsilon_0=1/2+4\varepsilon_0<0.95$, and the claim follows.
\end{proof}

Now we work with $M\in\O$.

\begin{lem}\label{lem4.8}
Let $M\in\O$, $D=0.13\frac{A}{|R_A|}2^t\log p_1$, and $D'<p_1$. Then
$$\sum_{D<l\leq D'}K_ll^2 \ll mD'p_1^{-0.03}.
$$
\end{lem}

\begin{proof}
Lemma \ref{lem4.5} gives us
\begin{equation}\label{4.4}
\sum_{D<l\leq D'} K_ll^2 \leq \sum_{\nu\,:\, D<2^{\nu}\leq 2D'}\sum_{2^{\nu-1}\leq u_j-u_{j-1}<2^{\nu}}(u_j-u_{j-1})^2 \ll \sum_{\nu\,:\, D<2^{\nu}\leq 2D'}|N_M(2^{\nu-2})|2^{\nu}.
\end{equation}
Recall that $M\in\O$ implies $A\leq 2^t<D$. Suppose that $2^{\nu}>D$ and denote $h=2^{\nu-2}$. Then in the interval $\{0,...,2^{\nu-2}-1\}$ there are at least
$$|R_A|[hA^{-1}]\geq |R_A|A^{-1}h-|R_A|\geq [0.03\cdot2^t\log p_1]=:s
$$
quadratic residues modulo $A$; denote the set of them by $C$. By Corollary \ref{cor2.5}
$$N_M(h)\leq N_m(C) \ll mp_1^{-0.03},
$$
and from (\ref{4.4}) we have
\[
\sum_{D<l\leq D'} K_ll^2 \ll D'mp_1^{-0.03}.
\]
This completes the proof.
\end{proof}

\begin{lem}\label{lem4.9}
Let $M\in\O$. Then for $D'=2^{51\o(M)}$ we have
\[
\sum_{l>D'}K_ll^2\ll M2^{-19\o(M)}.
\]
\end{lem}

\begin{proof}
Lemma \ref{lem4.5} gives us
\[
\sum_{l>D'}K_ll^2 \leq \sum_{\nu\,:\, 2^{\nu}>D'}\sum_{2^{\nu-1}\leq u_j-u_{j-1}<2^{\nu}}(u_j-u_{j-1})^2 \ll \sum_{\nu\,:\, 2^{\nu}>D'}|N_M(2^{\nu-2})|2^{\nu}.
\]
Note $2^{\nu}>D=\tau(M)^{51}$, together with $M\in\O$, implies $2^{\nu-2}>(2C\tau(M))^{50}$. Then we apply Lemma \ref{lem4.3} and obtain
\[
\sum_{l>D'}K_ll^2 \ll M2^{4\o(M)}\sum_{\nu: 2^{\nu}>D'}2^{-0.47} \ll M2^{-19\o(M)}.
\]
The claim follows.
\end{proof}

\begin{cor}\label{cor4.10}
	Let $M\in \O$. Then
	$$S(R_M) = m2^{-2t}g_A(y_t) + O(M(\log M)^{-1}).
	$$
\end{cor}

\begin{proof}
Set $D=0.13 A|R_A|^{-1}2^t\log p_1$ and $D'=2^{51\o(M)}$. Combining Corollary \ref{cor3.8}, Lemma \ref{lem4.8} and Lemma \ref{lem4.9}, we have
\begin{multline*}
S(R_M)=\sum_{l\leq D}K_ll^2+\sum_{l<l\leq D'}K_ll^2 + \sum_{l>D'}K_ll^2=\\ m2^{-2t}g_A(y_t) + O(MA^32^tp_1^{-0.13}+m2^{51\o(M)}p_1^{-0.03}+M2^{-19\o(M)}).
\end{multline*}
Recall that $M\in\O$ implies $A\leq 2^t$, $p_1>(\log M)^{C_0}$, and $\o(M)/\log\log M\in[0.9,1.1]$. Then taking $C_0$ large enough concludes the proof.
\end{proof}

\section{The end of the proofs of Theorems \ref{th1.3} and \ref{th1.5}}\label{Sec5}

Having proved Corollaries \ref{cor4.7} and \ref{cor4.10}, we found the assumed asymptotics for $S(R_M)$ in terms of the function $g_A$. Now we are aiming to extract further information about it. Recall that for any $y\in(0,1)$
\begin{equation}\label{5.1}
g_A(y)=\sum_{j=1}^{A}\sum_{i\in R_A\cap(R_A-j)}y^{r(i,j)}\sum_{k\geq 0}y^{|R_A|k}(Ak+j)^2;
\end{equation}
We begin with the following simple lemma, which allows us to simplify this expression.

\begin{lem}\label{lem5.1}
For all $x\in(0,1)$ we have
\begin{equation*}
\sum_{n=0}^{\infty}x^n=\frac{1}{1-x}, \qquad \sum_{n=1}^{\infty}nx^{n}=\frac{x}{(1-x)^2},
\qquad \sum_{n=1}^{\infty}n^2x^{n}=\frac{x(1+x)}{(1-x)^3}.
\end{equation*}
\end{lem}

\begin{proof} This can be obtained by standard direct calculations.
\end{proof}

We apply this lemma to the inner sum in the definition (\ref{5.1}) of the function $g_A$. We have
\begin{equation*}
\sum_{k=0}^{\infty}y^{|R_A|k}=\frac{1}{1-y^{|R_A|}}, \qquad  \sum_{k=0}^{\infty}ky^{|R_A|k}=\frac{y^{|R_A|}}{(1-y^{|R_A|})^2}, \qquad
\sum_{k=0}^{\infty}k^2y^{|R_A|k}=\frac{y^{|R_A|}(1+y^{|R_A|})}{(1-y^{|R_A|})^3},
\end{equation*}
and hence
\begin{equation}\label{5.2}
g_A(y)=\frac{A^2y^{|R_A|}(1+y^{|R_A|})P_0(y)}{(1-y^{|R_A|})^3}+\frac{2Ay^{|R_A|}P_1(y)}{(1-y^{|R_A|})^2} +\frac{P_2(y)}{1-y^{|R_A|}},
\end{equation}
where
$$P_0(y)=\sum_{i\in R_A}\sum_{j\in R_A-i}y^{r(i,j)},
$$
$$P_1(y)=\sum_{i\in R_A}\sum_{j\in R_A-i}jy^{r(i,j)},
$$
$$P_2(y)=\sum_{i\in R_A}\sum_{j\in R_A-i}j^2y^{r(i,j)}.
$$

Let $s_1,\ldots,s_{|R_A|}$ be the distances between consecutive quadratic residues modulo $A$; here and further we consider the indeces of $s_i$ modulo $|R_A|$. For $k=0,...,|R_A|-1$ we set $\a_k=\sum_{i=1}^{|R_A|}\left(\sum_{l=i}^{i+k}s_l\right)^2$. In particular, $\a_0=\sum_{i=1}^{|R_A|}s_i^2=S(R_A)$ and $\a_{|R_A|-1}=A^2|R_A|$.

Now we compute the polynomials $P_0$, $P_1$, and $P_2$.

\begin{lem}\label{lem5.2}
We have $P_0(y)=|R_A|\sum_{k=0}^{|R_A|-1}y^k$, $P_1(y)=A\sum_{k=0}^{|R_A|-1}(k+1)y^k$, and $P_2(y)=\sum_{k=0}^{|R_A|-1}\a_ky^k$.
\end{lem}

\begin{proof} Recall that $r(i,j)=\#((i,i+j)\cap R_A)$  (see (\ref{3.1})). Clearly, for any fixed $i\in R_A$, the function $r(i,j)$ takes values $0,1,\ldots,|R_A|-1$ exactly once. Let $s\in\{0,1,2\}$, and denote by $\{i_{\nu} \}_{\nu\in\mathbb{Z}_{|R_A|}}$ the set $R_A$. We have
\begin{multline*}
P_s(y)=\sum_{i\in R_A}\sum_{j\in R_A-i}j^sy^{r(i,j)}=
\sum_{k=0}^{|R_A|-1}y^k\sum_{i,i+j\in R_A}j^s\mathbb{I}(\#((i,i+j)\cap R_A)=k))=\\
\sum_{k=0}^{|R_A|-1}y^k\sum_{\nu=1}^{|R_A|}(i_{\nu+k+1}-i_{\nu})^s=
\sum_{k=0}^{|R_A|-1}y^k\sum_{i=1}^{|R_A|}\left(\sum_{l=i}^{i+k}s_l\right)^s.
\end{multline*}
Now if $s=0$, then we obtain $P_0(y)=|R_A|\sum_{k=0}^{|R_A|-1}y^k$. If $s=1$, then $\sum_{i=1}^{|R_A|}\sum_{l=i}^{i+k}s_l=\sum_{l=1}^{|R_A|}s_l(k+1)=(k+1)A$, and thus $P_1(y)=A\sum_{k=0}^{|R_A|-1}(k+1)y^k$. Finally, if $s=2$, then $P_2(y)=\sum_{k=0}^{|R_A|-1}\a_ky^k$ by the definition of the numbers $\a_k$. This completes the proof.
\end{proof}

Lemma \ref{lem5.2} gives us
$$(1-y)P_0(y)=|R_A|(1-y^{|R_A|})
$$
and
$$
(1-y)^2P_1(y)=A(1-y)\left(\sum_{k=0}^{|R_A|-1}y^k-|R_A|y^{|R_A|}\right)=
A(1-y^{|R_A|})-A|R_A|(1-y)y^{|R_A|}.
$$
Hence,
\begin{equation}\label{5.3}
\frac{A^2y^{|R_A|}(1+y^{|R_A|})(1-y)^3P_0(y)}{(1-y^{|R_A|})^3}=
\frac{A^2|R_A|y^{|R_A|}(1+y^{|R_A|})(1-y)^2}{(1-y^{|R_A|})^2}.
\end{equation}
and
\begin{equation}\label{5.4}
\frac{2Ay^{|R_A|}(1-y)^3P_1(y)}{(1-y^{|R_A|})^2}=
\frac{2A^2(1-y)y^{|R_A|}}{1-y^{|R_A|}}-
\frac{2A^2|R_A|y^{2|R_A|}(1-y)^2}{(1-y^{|R_A|})^2}.
\end{equation}
Now we compute the function $f_A(y)=(1-y)^3g_A(y)$. From (\ref{5.2}), (\ref{5.3}), and  (\ref{5.4}) we have
\begin{multline}\label{5.5}
f_A(y)=
\frac{A^2|R_A|(1-y)^2y^{|R_A|}}{1-y^{|R_A|}}+
\frac{2A^2(1-y)y^{|R_A|}}{1-y^{|R_A|}}+
\frac{(1-y)^3P_2(y)}{1-y^{|R_A|}}=\\
\frac{A^2y^{|R_A|}(2+|R_A|(1-y))+(1-y)^2P_2(y)}{1+y+\ldots+y^{|R_A|-1}}.
\end{multline}

We continue working with the function $f_A$.

\begin{lem}\label{lem5.3}
	We have
	$$f_A(y)=\frac{F(y)}{Q(y)},
	$$
where $Q(y)=Q_A(y)=1+y+\ldots+y^{|R_A|-1}$ and $F(y)=F_A(y)=\sum_{k=0}^{|R_A|}\b_ky^k$ be the reciprocal polynomial with the coefficients
$\b_0=\b_{|R_A|}=\sum_is_i^2$
and
$\b_k=2\sum_{i=1}^{|R_A|}s_is_{i+k}$
for $0<k<|R_A|$.

\end{lem}

\begin{proof} From (\ref{5.5}) and Lemma \ref{lem5.2} we see that
\begin{equation*}
F(y)=(1-2y+y^2)\sum_{k=0}^{|R_A|-1}\a_ky^k+A^2(|R_A|+2)y^{|R_A|}-A^2|R_A|y^{|R_A|+1}.
\end{equation*}
Now we write $F(y)=\sum_{k=0}^{|R_A|+1}\gamma_ky^k$ and compute $\gamma_k$. Firstly,
$$\gamma_0=\a_0=\sum_{i=1}^{|R_A|}s_i^2,
$$
and
$$\gamma_1=\a_1-2\a_0=\sum_{i=1}^{|R_A|}(s_i+s_{i+1})^2-2\sum_{i=1}^{|R_A|}s_i^2=2\sum_{i}s_is_{i+1}.
$$
Further, observe that
\begin{equation*} \a_k=\sum_i(s_i+\ldots+s_{i+k})^2=
(k+1)\sum_is_i^2+2\sum_{l=1}^k(k+1-l)\sum_{i}s_is_{i+l};
\end{equation*}
then for $2\leq k\leq |R_A|-1$ we have
$$\gamma_k=\a_{k-2}-2\a_{k-1}+\a_k=2\sum_is_is_{i+k}.
$$
It remains to calculate $\gamma_{|R_A|}$ and $\gamma_{|R_A|+1}$. Recall that $\a_{|R_A|-1}=A^2|R_A|$. Then we have
$$
\gamma_{|R_A|+1}=\a_{|R_A|-1}-A^2|R_A|=0
$$
and
$$\gamma_{|R_A|}=-2\a_{|R_A|-1}+\a_{|R_A|-2}+(|R_A|+2)A^2=\a_{|R_A|-2}+2A^2-A^2|R_A|.
$$
Finally,
\begin{multline*}
\a_{|R_A|-2}=\sum_{i}\left(\sum_{l=i}^{i+|R_A|-2}s_{i+l}\right)^2=\sum_{i=1}^{|R_A|}(A-s_i)^2=\\
A^2|R_A|-2A\sum_is_i+\sum_is_i^2=A^2|R_A|-2A^2+\a_0,
\end{multline*}
and therefore $\gamma_{|R_A|}=\gamma_0=\sum_is_i^2$. Thus $F(y)=\sum_{k=0}^{|R_A|}\b_ky^k$, as desired.
\end{proof}

Since $f_A(y)=(1-y)^3g_A(y)$, we have $g_A(0.5)=8f_A(0.5)$ and $g_A(y_t)=2^{3t}f_A(y_t)$ for $y_t=1-2^{-t}$. Now Theorem \ref{th1.3} follows immediately from Corollary \ref{cor4.7}; also we can rewrite Corollary \ref{cor4.10} as
\begin{equation}\label{5.6}
S(R_M)=m2^tf_A(y_t)+O(M(\log M)^{-1})
\end{equation}
for $M\in\O$. To finish the proof of Theorem \ref{th1.5}, we expand the function $f_A$ to the Taylor series at the point $y=1$.

We see from Lemma \ref{lem5.3} that $f_A$ obeys the functional equation
$$f_A(y)=yf_A\left(\frac1y\right).
$$
We then obtain
\begin{equation}\label{5.7}
f'_A\left(\frac1y\right)=f_A(y)-yf_A'(y)
\end{equation}
and
\begin{equation}\label{5.8}
f_A''\left(\frac1y\right)=y^3f''_A(y).
\end{equation}
Also from (\ref{5.1}) we have
\begin{equation}\label{5.9}
f_A(1)=\frac{2A^2}{|R_A|},
\end{equation}
and (\ref{5.7}) immediately implies $f'_A(1)=0.5f_A(1)=A^2|R_A|^{-1}$. Then by the Taylor expansion
\begin{equation}\label{5.10}
m2^tf_A(y_t)=m2^{t+1}A^2|R_A|^{-1}-mA^2|R_A|^{-1}+\frac12f''_A(\theta_{A,t})m2^{-t}
\end{equation}
for some $\theta_{A,t}\in(1-2^{-t},1)$, and it remains to estimate the second derivative of $f_A$.

\begin{lem}\label{lem5.4}
	Suppose that $|R_A|\leq 2^{t}$. Then
	$$
	f''_A(y)\ll A^2|R_A|$$
	for $y\in(y_t,1)$.
\end{lem}

\begin{proof} Due to (\ref{5.8}), it suffices to prove the lemma for $1\leq y\leq 1+O(2^{-t})$ instead of $y\in(y_t,1)$. Firstly,
\begin{equation}\label{5.11}
\left(\frac{F(y)}{Q(y)}\right)''=\frac{F''(y)Q^2(y)-2F'(y)Q'(y)Q(y)-F(y)Q''(y)Q(y)+2F(y)Q'^2(y)}{Q^3(y)}.
\end{equation}
Further, for $1\leq y\leq 1+O(2^{-t})$ by the assumption on $|R_A|$ we have $y^k \asymp 1$ uniformly for $0\leq k\leq |R_A|$ and thus
$$Q(y)\asymp |R_A|, \qquad Q'(y)\asymp |R_A|^2, \qquad Q''(y) \asymp |R_A|^3.
$$
Taking into account (\ref{5.9}), we find
$$F(y)\ll \sum_{k=0}^{|R_A|}\b_k=F(1)=f_A(1)Q(1)=2A^2,
$$
$$F'(y)=\sum_{k=1}^{|R_A|}k\b_ky^{k-1}\ll A^2|R_A|, \qquad F''(y)\ll A^2|R_A|^2.
$$
Thus for $1\leq y\leq 1+O(2^{-t})$ we have from (\ref{5.11})
$$f_A''(y) \ll \frac{A^2|R_A|^4}{|R_A|^3}=A^2|R_A|,
$$
and the claim follows.
\end{proof}

Now combining (\ref{5.6}), (\ref{5.10}) and the bound from Lemma \ref{lem5.4}, we get
\begin{equation}\label{5.12}
S(R_M) = M2^{t+1}A|R_A|^{-1}-MA|R_A|^{-1} + O(M(\log M)^{-1}+MA|R_A|2^{-t}).
\end{equation}
Recalling that $M\in\O$ implies $A\leq (\log M)^{0.1}$ and $2^{-t}\leq2^{-(0.9+o(1))\log\log M}$, we get
\begin{equation}\label{5.13}
A|R_A|2^{-t}\leq (\log M)^{-0.4}.
\end{equation}
Further, since $\prod_{i=1}^t(1+p_i^{-1})\leq\exp(\sum_{i=1}^tp_i^{-1})=1+O(tp_1^{-1}),
$
we have
\begin{equation*}
|R_M|=|R_A|\prod_{i=1}^t\frac{p_i+1}{2}=|R_A|m2^{-t}(1+O(tp_1^{-1})),
\end{equation*}
and hence
\begin{equation*}
\frac{M}{|R_M|}=A|R_A|^{-1}2^t(1+O(tp_1^{-1}))=A|R_A|^{-1}2^t+O((\log M)^{1-C_0})
\end{equation*}
and
\begin{equation}\label{5.14}
\frac{2M}{|R_M|^2}=M2^{t+1}A|R_A|^{-1}+O(M(\log M)^{1-C_0}).
\end{equation}
Putting (\ref{5.12}), (\ref{5.13}), and (\ref{5.14}) together and recalling that $C_0>2$, we conclude the proof of Theorem \ref{th1.5}.

\section{Proof of Theorem \ref{th1.4}} \label{Sec6}

Let $X>0$ be a large number and $\O(X)=\{M\leq X: M\in\O \}$. We first prove the lower bound for $|\O(X)|$. To do this, we consider the set $\O_1(X)$ of sufficiently large positive integers $\sqrt{X}<M\leq X$ such that

\begin{itemize}

	\item[(i)] $M$ is square-free;

	\item[(ii)] $0.9\log\log X < \o(M) < 1.05\log\log X$;

	\item[(iii)] $M$ has the representation $M=Am$, where $P^-(m)>(\log X)^{C_0}$ and \\ $3\leq A\leq0.5(\log X)^{0.1}$.

\end{itemize}
Clearly, $\O_1(X)\subseteq \O(X)$. We are going to show that $|\O_1(X)|= (\eta+o(1))X$, where $\eta=e^{-\gamma}\frac{3}{5\pi^2C_0}$ and $\gamma$ is the Euler constant. It is well-known that the number of $M\leq X$ violating (ii) is $O(X(\log X)^{-c})$ for some absolute constant $c>0$ (see, for example, \cite{Ten}, Exercise 04). Further, let $S(X)$ be the set of square-free $m\leq X$ such that $P^-(m)>(\log X)^{C_0}$. We are aiming to prove that
\begin{equation}\label{6.1}
|S(X)| = \frac{e^{-\gamma}X}{C_0\log\log X}+O\left(\frac{X}{(\log\log X)^2}\right).
\end{equation}
Set $z=(\log X)^{C_0}$ and $P(z)=\prod_{p\leq z}p$. Obviuosly $S(X) = S_1\setminus S_2$, where
$$S_1=\{m\leq X: (m,P(z))=1\}
$$
and
$$S_2=\{m\leq X: \exists p>z \mbox{ with } p^2|m \}.
$$
Standard sieve methods (see, for example, \cite{Mur}, Exercises 9.1.1 and 9.1.9) yield  that
\begin{equation*}
|S_1| = \frac{e^{-\gamma}X}{\log z}+O\left(\frac{X}{\log^2 z}\right) = \frac{e^{-\gamma}X}{C_0\log\log X}+O\left(\frac{X}{(\log\log X)^2}\right);
\end{equation*}
also we have
\begin{equation*}
|S_2|\leq \sum_{p>z}\left[\frac{X}{p^2}\right] \leq Xz^{-1} = X(\log X)^{-C_0}.
\end{equation*}
These two estimates imply (\ref{6.1}).

Now for any $A\leq 0.5(\log X)^{0.1}$ we have $\log\log (X/A)=\log\log X+O\left(\frac{\log\log X}{\log X}\right)$ and thus
$$|S(X/A)|=\frac{e^{-\gamma}X}{C_0A\log\log X}+O\left(\frac{X}{A(\log\log X)^2}\right).
$$
Also it is well-known that $\sum_{A\leq y}\mu^2(A)=\frac{6}{\pi^2}y+O(\sqrt{y})$; from here we have $\sum_{3\leq A\leq y}\frac{\mu^2(A)}{A}=\frac{6}{\pi^2}\log y+O(1)$. Thus, using (\ref{6.1}), we obtain
\begin{multline*}
|\O_1(X)| = \sum_{3\leq A\leq 0.5(\log X)^{0.1}} \mu^2(A)|S(X/A)|+O\left(\frac{X}{(\log X)^c}\right) = \\
\frac{e^{-\gamma}X}{C_0\log\log X}\sum_{3\leq A\leq0.5(\log X)^{0.1}}\frac{\mu^2(A)}{A}+o(X)=(\eta+o(1))X.
\end{multline*}

To prove the upper bound $|\O(X)|$, we consider the set $\O_2(X)$ of positive integers $\sqrt{X}<M\leq X$ such that

\begin{itemize}

	\item[(i')] $M$ is square-free;

	\item[(ii')] $0.85\log\log X < \o(M) < 1.1\log\log X$;

	\item[(iii')] $M$ has the representation $M=Am$, where $P^-(m)>(0.5\log X)^{C_0}$ and \\ $3\leq A\leq(\log X)^{0.1}$.

\end{itemize}

Arguing as before, we see that $|\O_2(X)|=(\eta+o(1))X$. Finally, any $M\in\O(X)$ with $M>\sqrt{X}$ belongs to $\O_2(X)$, and thus $|\O(X)|\leq(\eta+o(1))X$.

This completes the proof of Theorem \ref{th1.4}.

\section{Proof of the identity (\ref{1.3})} \label{Sec7}

We have
\begin{equation}\label{7.1}
{M \choose k}s(k)=\sum_{|U|=k}S(U)=\frac{1}{k}\sum_{x\in\mathbb{Z}_M}\sum_{\substack{|U|=k\\mathbf U\ni x}}S(U)=\frac{M}{k}\sum_{\substack{|U|=k\\mathbf U\ni 0}}S(U).
\end{equation}
Set $\ensuremath{\mathcal{A}}_{M,k}=\left\{\overline{x}=(x_1,\ldots,x_k)\in \mathbb{N}^k : \sum_{i=1}^kx_i=M\right\}$; then $|\ensuremath{\mathcal{A}}_{M,k}|={M-1 \choose k-1}$. Note that a set $U=\{0=u_1<u_2<\ldots<u_k<M\}\subseteq \mathbb{Z}_M$ containing $0$ is uniquely determined by a $k$-tuple $\overline{x}\in \ensuremath{\mathcal{A}}_{M,k}$ of the distances between its elements. Then by (\ref{7.1}) and the symmetry
$$
{M-1 \choose k-1}s(k)=\sum_{\substack{|U|=k\\mathbf U\ni 0}}S(U) = \sum_{\overline{x}\in\ensuremath{\mathcal{A}}_{M,k}}\sum_{i=1}^k{x_i^2}=k\sum_{\overline{x}\in\ensuremath{\mathcal{A}}_{M,k}}{x_1^2}.
$$
Further, the number of $\overline{x}\in\ensuremath{\mathcal{A}}_{M,k}$ with $x_1=a$ (here $a=1,...,M-k+1$) is equal to $|\ensuremath{\mathcal{A}}_{M-a,k-1}|={M-a-1 \choose k-2}$. Thus
\begin{equation}\label{s(k)2}
{M-1 \choose k-1}s(k)=k\sum_{\nu=1}^{M-k+1}a^2{M-1-a \choose k-2}.
\end{equation}
It remains to compute the latter sum. By induction on $v$ it is easy to prove the identities (first (\ref{7.3}), then using it, (\ref{7.4}), and then (\ref{7.5}))
\begin{equation}\label{7.3}
\sum_{a=0}^v{u+v-a \choose u} = {u+v+1 \choose u+1},
\end{equation}
\begin{equation}\label{7.4}
\sum_{a=1}^va{u+v-a \choose u} = {u+v+1 \choose u+2},
\end{equation}
\begin{equation}\label{7.5}
\sum_{a=1}^va^2{u+v-a \choose u} = 2{u+v+1 \choose u+3}+{u+v+1 \choose u+2}.
\end{equation}
Using (\ref{7.5}) for $u=k-2$ and $v=M-k+1$, we get from (\ref{s(k)2})
$$s(k)=k\frac{2{M \choose k+1}+{M \choose k}}{{M-1 \choose k-1}}=M\frac{2M-k+1}{k+1}\,,
$$
as desired.

\bigskip

\textbf{Acknowledgements.} This work is supported by the Russian Science Foundation under grant 19--11--00001. The author is a Young Russian Mathematics award winner and would like to thank its sponsors and jury. The author also would like to thank Sergei Konyagin and Yuri Malykhin for helpful discussions, and Pavel Sobyanin and Konstantin Olmezov for computational work concerned with the values $f_A(0.5)$ for small $A$.