# Reconstructing Speech from Real-Time Articulatory MRI Using Neural Vocoders

Yide Yu  Institute of Informatics
University of Szeged
Szeged, Hungary
<EMAIL_ADDRESS>Amin Honarmandi Shandiz  Institute of Informatics
University of Szeged
Szeged, Hungary
<EMAIL_ADDRESS>László Tóth  Institute of Informatics
University of Szeged
Szeged, Hungary
<EMAIL_ADDRESS>
###### Abstract

Several approaches exist for the recording of articulatory movements, such as
eletromagnetic and permanent magnetic articulagraphy, ultrasound tongue
imaging and surface electromyography. Although magnetic resonance imaging
(MRI) is more costly than the above approaches, the recent developments in
this area now allow the recording of real-time MRI videos of the articulators
with an acceptable resolution. Here, we experiment with the reconstruction of
the speech signal from a real-time MRI recording using deep neural networks.
Instead of estimating speech directly, our networks are trained to output a
spectral vector, from which we reconstruct the speech signal using the
WaveGlow neural vocoder. We compare the performance of three deep neural
architectures for the estimation task, combining convolutional (CNN) and
recurrence-based (LSTM) neural layers. Besides the mean absolute error (MAE)
of our networks, we also evaluate our models by comparing the speech signals
obtained using several objective speech quality metrics like the mean cepstral
distortion (MCD), Short-Time Objective Intelligibility (STOI), Perceptual
Evaluation of Speech Quality (PESQ) and Signal-to-Distortion Ratio (SDR). The
results indicate that our approach can successfully reconstruct the gross
spectral shape, but more improvements are needed to reproduce the fine
spectral details.

###### Index Terms:

Real-Time MRI, articulatory-to-acoustic mapping, deep learning

##  I Introduction

Human speech production is a complex process that requires precisely
coordinated movements from the respiratory organs, the larynx and the
articulators. Theoretically, the configuration of the articulatory organs
determines the speech signal that is being produced. Articulatory-to-acoustic
mapping seeks to determine whether it is possible to estimate the speech
signal if we know the physical positions of these organs. During the last
decade, several types of devices have been proposed to record the movement of
the articulators. The most important of these are ultrasound tongue imaging
(UTI) [1, 2, 3, 4, 5, 6], electromagnetic articulography (EMA) [7, 8, 9, 10],
permanent magnetic articulography (PMA) [11, 12], and surface electromyography
(sEMG) [13, 14, 15, 16, 17]. While articulatory-to-acoustic mapping can also
be considered as a theoretical problem, it has an important application, which
is the creation of silent speech interfaces [18]. Silent speech interfaces
(SSI) seek to generate speech from the articulatory movement without actually
having access to the speech signal, with the aim of aiding the communication
of speaking impaired people who can move their articulators, but cannot
actually produce speech. SSIs may also be applied in situations where the
speech signal cannot be heard, as in extremely noisy industrial environments
and certain military applications.

Currently, the devices listed above typically record noisy and low resolution
signals. Magnetic resonance imaging (MRI) may offer an alternative solution,
as it can produce high resolution images. Moreover, while the other methods
capture only the lingual, labial and jaw motions, MRI also covers the
pharyngeal and nasal regions, which are not reachable by the other methods
[19]. Of course, MRI is costly and it requires a huge machine, so currently
this technology cannot provide a basis for wearable silent speech interfaces,
but treating the articulatory-to-acoustic mapping problem with MRI as input is
a very interesting theoretical research topic in itself. The current
developments of MRI technology allow us to create real-time recordings of the
speech production process at rates of about 20-50 frames per second. While
this is lower than the sampling rate of other tools such as ultrasound tongue
imaging, the spatial resolution of MRI is typically much better [20].
Therefore, compared to the other tools, it provides complementary information,
and thus examining the applicability of MRI recordings for articulatory-to-
acoustic conversion may lead to some intriguing conclusions.

Several previous studies have applied articulatory MRI for speech-related
tasks. The majority of these tried to perform speech recognition using the MRI
as input [21, 22, 23]. Interestingly, some authors tried to estimate the MRI
from the speech signal, which is the inverse of the problem we study here
[24]. However, relatively few attempts have been made to synthesize speech
signals based on the MRI [25, 26]. Here, we extend the recent study by Tamás
G. Csapó [26]. Similar to his approach, we use the MRI recordings of the USC-
TIMIT data set as input [27]. However, while he used a simple conventional
vocoder for the speech synthesis process, here we apply more recent neural
vocoders for this task. As the name suggests, these are based on deep learning
models, and they have been reported to produce higher quality speech [28].
Hence, the we approach we present here is purely neural. In the next section
we shall introduce the concept of our articulatory-to-acoustic mapping
framework, with a special focus on the application of neural vocoders.

##  II MRI-based Articulatory-to-Acoustic Conversion

###  II-A Real-Time Articulatory MRI

As the input MRI data we used the freely available USC-TIMIT data set [27].
This data set contains synchronized speech and real-time MRI recordings from
American English speakers. The subjects lying in an MRI device were asked to
read 460 sentences from the MOCHA-TIMIT database. The MRI data was captured
using an 1.5 Tesla Sigma Excite HD MRI scanner. The resolution of the recorded
images is 64x64 pixels, and their orientation is adjusted to the mid-saggital
plane. The time resolution of the recordings was 23 frames per second. The
speech signal was recorded simultaneously with the MRI video inside the MRI
scanner at a sampling rate of 20 kHz. Due to the noisy operation of the
scanner, the recorded speech signals were very noisy, so they had to be post-
processed using noise cancellation algorithms.

Fig 1 shows two example images from a male and a female speaker from the
database. As can be seen, the images have significant inter-speaker
differences, so in the experiments we created speaker-dependent models, that
is, separate models for each speaker. The database contains samples from five
male and five female subjects, and we worked with the data of two males and
two females. Since Csapó reported earlier that the data of speaker ’M1’ has
significant misalignment problems [26], we worked with speakers ’M2’, ’M3’,
’F2’ and ’F3’.

Figure 1: Sample MRI images from two speakers.

###  II-B Speech Generation using Neural Vocoders

As shown in Fig 2, our goal is to convert a sequence of MRI images into a
speech signal. This task can be formalized as a sequence-to-sequence mapping
problem. For such tasks, various sophisticated deep neural network (DNN)
structures have been proposed recently that do not even require aligned data.
However, our MRI and speech sequences were synchronized, allowing us to use
much simpler architectures that estimate the mapping in a pairwise manner,
that is, give one output vector for each element of the input sequence.

Figure 2: Schematic diagram of the MRI-to-speech conversion process applied here. Out goal is to find the optimal structure and parameters for the first network of the processing chain. TABLE I: The layers of the 3 network architectures, along with their main parameters. To save space, some details such as dropout layers and regularization parameters are not shown. 2D-CNN+BiLSTM | 3D-CNN | 3D-CNN+BiLSTM
---|---|---
TimeDistributed(Conv2D(30,(13,13),strides=(2,2)) | Conv3D(30, (5, 13, 13), strides=(sts, 2, 2) | Conv3D(30, (5, 13, 13), strides=(sts, 2, 2)
TimeDistributed(Conv2D(60,(13,13),strides=(2,2)) | Conv3D(60, (1, 13, 13), strides=(1, 2, 2) | Conv3D(60, (1, 13, 13), strides=(1, 2, 2)
TimeDistributed(MaxPooling2D((2,2))) | MaxPooling3D((1, 2, 2)) | MaxPooling3D((1, 2, 2))
TimeDistributed(Conv2D(90,(13,13),strides=(1, 1)) | Conv3D(90, (1, 13, 13), strides=(1, 1, 1) | Conv3D(90, (1, 13, 13), strides=(1, 1, 1)
TimeDistributed(Conv2D(85,(13,13),strides=(2,2)) | Conv3D(120, (5, 3, 3), strides=(1, 2, 2) | Conv3D(120, (5, 3, 3), strides=(1, 2, 2)
TimeDistributed(MaxPooling2D((2,2)) | MaxPooling3D((1, 2, 2)) | MaxPooling3D((1, 2, 2))
TimeDistributed(Flatten()) | Flatten() | Reshape((5, 480))
Bidirectional(LSTM(320) | Dense(1000) | Bidirectional(LSTM(370)
Dense(80, activation=’linear’) | Dense(80, activation=’linear’) | Dense(80, activation=’linear’)

Thus, the input of our network is an MRI image, or a short sequence of
consecutive images. However, the optimal choice for the output or target
vector is not trivial. Creating DNNs that output speech signals directly is
feasible [29, 30]. However, such networks generate tens of thousands of speech
samples per second, and, consequently, these models are enormous and their
training requires a huge amount of speech data. This is why we applied an
alternative, indirect approach here, which is depicted in Fig 2. This approach
was motivated by neural speech synthesis, where DNNs are applied in two steps
[30]. First, there is a network that estimates a spectral representation from
the text to be synthesized. Then, there is a second network that generates the
speech signal from the spectral representation. Adapting this approach to our
task requires modifications only in the first network, as our input is an MRI
video and not a text. However, the second step is the same, so we can borrow
large, pre-trained networks for the second task from text-to-speech synthesis.
Doing so, we had to create and train only the first network, which had to
estimate a dense spectral representation instead of the speech waveform
itself. As for the second, the speech generation task, several neural vocoders
are available [28], and we chose to use the WaveGlow model [30], as it worked
well in our earlier study [31]. WaveGlow requires a sequence of 80-dimensional
mel-scaled spectral vectors as the input, so our task was to create a network
that can estimate such a spectral vector for each frame of the MRI video. As
the default frame rate of WaveGlow (86 fps) is almost 4 times higher than the
frame rate of our MRI video (23 fps), we estimated the missing spectral
vectors by applying interpolation before the actual speech synthesis step.

##  III DNN Architectures for MRI-to-Spectrogram Conversion

To prepare the MRI images for DNN training, the pixel intensities of each
image were min-max normalized to the $[-1,1]$ range, and the speech recordings
were resampled at 22050 Hz, as this is the sampling rate required by WaveGlow.
The speech signals were then converted to a mel-spectral representation, and
the resulting 80-dimensional mel-spectral vectors were standardized before
using them as the DNN training targets. From the 92 sentences of each subject,
4 were used for validation and 2 for testing.

Formally, our networks has to map each MRI image to a spectral vector.
However, using several consecutive input frames instead of a single frame can
significantly improve the results [32, 6]. Hence, the input for all our
network configurations was a 3D array, treating time as the the third axis
besides the two spacial axes of the images. Table I summarizes the main
parameters of our experimental network configurations, and we also give a
brief description of each below. Our networks were implemented in Keras. We
applied the Swish activation function [33] in the hidden layers. The input
window of the network contained 13 consecutive video frames. The number of
trainable parameters was approximately the same for all 3 configurations. For
the loss function we applied the mean absolute error (MAE), as it was reported
to give slightly better results than the more conventional mean-squared error
for speech-related tasks [34].

2D-CNN+BiLSTM: When working with a sequence of images, a popular technique is
to process each image using a convolutional neural network (CNN), and then
combine the results along the time axis using recurrent neural structures such
as the long short-term memory (LSTM) layer [35]. Thus, our first network
combined 2D-CNN layers that process each image with an LSTM layer on top to
fuse the information along the time axis. As shown in Table I, the 2D
convolutional and max-pooling layers are applied to each image of the input
sequence using the TimeDistributed() function of Keras. Then their outputs are
combined along the time axis using an LSTM layer, and the actual regression
task with respect to the spectral target vectors is performed by the topmost
linear layer.

3D-CNN: Several authors argued recently that good video classification is also
achievable using purely convolutional structures by extending the convolution
to the temporal axis [36, 6, 32]. This is why a 3D convolutional model served
as our second model (see Table I). Its lowest layer processes the input in
5-frame blocks with a hop size controlled by the $sts$ parameter. This
parameter allows us to analyze input blocks that are placed at bigger time
intervals. In an earlier study we processed ultrasound videos which had a much
larger frame rate, and the optimal value for $sts$ was found to be 5 [32].
Here, we got the best performance with $sts=3$, which is reasonable as the
frame rate was much lower. The subsequent Conv3D layers essentially process
each input block separately (their filter size along the time axis being set
to 1), and the uppermost Conv3D layer performs the fusion along time. This
structure was motivated by the findings of Tran et al. [36]. Finally, the
regression is performed by a Dense layer.

3D-CNN+BiLSTM: A drawback of the 2D-CNN+BiLSTM model is that it cannot skip
input frames, while the drawback with the 3D-CNN model is that recurrent
layers may be more effective in fusing information obtained at several points
along the time axis. Hence, we also experimented with a third model that
combines the advantages of the two previous architectures. As Table I shows,
this model retained the basic architecture of the 3D-CNN network, but we
replaced the uppermost Dense layer with an LSTM layer.

##  IV Results and Discussion

TABLE II: Mean absolute error values for the validation and test sets.

| Mean absolute error (validation / test)
---|---
speaker | 2D-CNN+BiLSTM |  3D-CNN | 3D-CNN+BiLSTM
‘F2’ | 0.26 / 0.26 | 0.28 / 0.28 | 0.26 / 0.26
‘F3’ | 0.48 / 0.40 | 0.48 / 0.40 | 0.45 / 0.40
‘M2’ | 0.37 / 0.33 | 0.36 / 0.33 | 0.36 / 0.32
‘M3’ | 0.30 / 0.32 | 0.31 / 0.32 | 0.29 / 0.32
avg. | 0.35 / 0.33 | 0.36 / 0.33 | 0.34 / 0.33

TABLE III: MCD scores for the test set.

| Mel-Cepstral Distortion (dB)
---|---
speaker | 2D-CNN+BiLSTM |  3D-CNN | 3D-CNN+BiLSTM
‘F2’ | 5.87 | 6.02 | 5.84
‘F3’ | 5.75 | 5.74 | 5.84
‘M2’ | 5.54 | 5.47 | 5.40
‘M3’ | 5.06 | 5.09 | 4.95
average | 5.56 | 5.58 | 5.51

TABLE IV: Objective speech quality scores for the test set. speaker | STOI | PESQ | SDR
---|---|---|---
’F2’ | 0.51 | 1.67 | -22.7
’F3’ | 0.27 | 1.86 | -25.7
’M2’ | 0.52 | 1.65 | -18.9
’M3’ | 0.41 | 1.81 | -22.4
average | 0.43 | 1.75 | -22.4
Figure 3: Spectrograms of a sample sentence from speaker M3. Left: original,
right: reconstructed from MRI.

Table III shows the MAE values obtained by the three network configurations on
the validation and test sets of each speaker. The results indicate that there
are large interpersonal differences compared to the average shown in the
bottom row. However, the three networks architectures produced very similar
results. Based on the average score on the validation set, the 3D-CNN+BiLSTM
model seems to be slightly better than the 2D-CNN+BiLSTM network, and the
latter is slightly better than the 3D-CNN, but the differences are negligible,
and the average loss values are equivalent on the test set.

To compare not only the DNN loss values but also the synthesized speech
signals, we calculated the mel-cepstral distortion (MCD) of the test files.
MCD is frequently used as an objective measure of speech quality in speech
synthesis [37]. The MCD values reported in Table III are not in complete
accord with the MAE loss values obtained for the various speakers. This
reflects the known fact that the simple DNN loss functions such as the MAE or
MSE do not necessarily coincide with human perception [38]. However, the
average MCD values have the same tendency as the MAE values in the sense that
the 3D-CNN+BiLSTM model seems to be slightly better than the other two
networks, but the difference is minimal.

To compare our results with those of the literature, Csapó performed a similar
experiment using a conventional MGLSA vocoder, and he reported MCD scores
around 4.5 for speakers ’F2’ and ’M2’ [26]. However, he attempted to estimate
only the spectral envelope, and he used the residual component of the original
speech signal during the synthesis. In contrast, our approach attempts to
reconstruct the full spectrum from the MRI only, including the fine spectral
details such as the pitch information. To illustrate the performance of our
solution in this respect, we give an example of the narrow-band spectrogram
(created from the output of the 3D-CNN+BiLSTM net) for speaker ’M3’ in Fig. 3.
Comparing the estimated signal with the original, we see that our network is
quite successful in reconstructing the rough spectral shape, but it fails with
the fine spectral details. The horizontal stripes that reflect the fundamental
frequency and its harmonics are preserved only in the lower spectral region.
This is probably due the the mel-scale used by WaveGlow, as it represents the
higher frequencies with a lower resolution. There is considerable smearing
present along the time axis as well. For example. the plosive burst at 2
seconds in the original spectrogram is missing in the reconstructed signal. A
further analysis is required to see whether the relatively low frame rate of
the MRI is responsible for this.

We evaluated additional objective metrics to assess the quality and
intelligibility of the synthesized signals (see Table IV). The Short-Time
Objective Intelligibility (STOI) metric [39] returns values between -1 and 1,
so the average score of 0.43 we attained is fair, but not that good. The
Perceptual Evaluation of Speech Quality (PESQ) metric [40] returns values on
the mean opinion scale between 1 and 5. On this scale, our average score of
1.75 lies between ”bad” and ”poor”. Lastly, the Signal-to-Distortion Ratio
(SDR) metric [41] behaves much like the signal-to-noise ratio, so the negative
values we got indicate that there is a large rate of distortion in our
signals. However, we emphasize that the original speech recordings were also
of low quality, as they were recorded in an MRI device, and they were post-
processed by noise cancellation methods.

##  V Conclusions

Here, we investigated the feasibility of performing articulatory-to-acoustic
mapping using neural vocoders such as WaveGlow. Compared to earlier attempts,
the approach we proposed estimates the whole spectral content of the signal
and not just the spectral envelope. Our experiments showed that our method is
able to reconstruct the gross spectral shape and also some fine details such
as the lowest pitch harmonics, but the spectrograms we got were quite blurred,
and hence the resulting speech signals are at the lower end of the scale,
according to several objective speech quality metrics. In the future we plan
to incorporate perceptually motivated loss functions in the training process,
as this would allow the direct optimization of the output speech quality [38].

## Acknowledgments

This study was supported by grant NKFIH-1279-2/2020 of the Ministry for
Innovation and Technology, Hungary, and by the Ministry of Innovation and the
National Research, Development and Innovation Office by grant FK 124584 and
the framework of the Artificial Intelligence National Laboratory Programme.
The GPU card used for the computations was donated by the NVIDIA Corporation.
The authors are grateful for the comments and source code of Tamás Gábor
Csapó.

## References

  * [1] T. Hueber, E. Benaroya, B. Denby, and G. Chollet,  “Statistical mapping between articulatory and acoustic data for an ultrasound-based silent speech interface,”  in Proc. Interspeech, 2011, pp. 593–596.
  * [2] A. Jaumard-Hakoun, K. Xu, C. Leboullenger, P. Roussel-Ragot, and B. Denby,  “An articulatory-based singing voice synthesis using tongue and lips imaging,”  in Proc. Interspeech, 2016, pp. 1467–1471.
  * [3] T. G. Csapó, T. Grósz, G. Gosztolya, L. Tóth, and A. Markó,  “DNN-based ultrasound-to-speech conversion for a silent speech interface,”  in Proc. Interspeech, 2017, pp. 3672–3676.
  * [4] K. Xu, P. Roussel, T. G. Csapó, and B. Denby,  “Convolutional neural network-based automatic classification of midsagittal tongue gestural targets using b-mode ultrasound images,”  The Journal of the Acoustical Society of America, vol. 141, no. 6, pp. EL531–EL537, 2017.
  * [5] E. Tatulli and T. Hueber,  “Feature extraction using multimodal convolutional neural networks for visual speech recognition,”  in Proceedings of ICASSP, 2017, pp. 2971–2975.
  * [6] P. Saha, Y. Liu, B. Gick, and S. Fels,  “Ultra2speech – a deep learning framework for formant frequency estimation and tracking from ultrasound tongue images,”  in Proc. MICCAI, 2020, pp. 473–482.
  * [7] J. Wang, A. Samal, J.R. Green, and F. Rudzicz,  “Sentence recognition from articulatory movements for silent speech interfaces,”  in Proc. ICASSP, 2012, pp. 4985–4988.
  * [8] W. Jun, S. Ashok, and G. Jordan,  “Preliminary test of a real-time, interactive silent speech interface based on electromagnetic articulograph,”  in Proc. SLPAT, 2014, pp. 38–45.
  * [9] M. Kim, B.g Cao, T. Mau, and J. Wang,  “Speaker-Independent Silent Speech Recognition From Flesh-Point Articulatory Movements Using an LSTM Neural Network,”  IEEE/ACM Trans. ASLP, vol. 25, no. 12, pp. 2323–2336, 2017.
  * [10] F. Taguchi and T. Kaburagi,  “Articulatory-to-speech conversion using bi-directional long short-term memory,”  in Interspeech, 2018, pp. 2499–2503.
  * [11] M. J. Fagan, S. R. Ell, J. M. Gilbert, E. Sarrazin, and P. M. Chapman,  “Development of a (silent) speech recognition system for patients following laryngectomy,”  Medical Engineering and Physics, vol. 30, no. 4, pp. 419–425, 2008\.
  * [12] J. A. Gonzalez, R. K. Moore, J. M. Gilbert, L.A. Cheah, S. Ell, and J. Bai,  “A silent speech system based on permanent magnet articulography and direct synthesis,”  Computer, Speech and Language, vol. 39, pp. 67–87, 2016.
  * [13] K. Nakamura, M. Janke, M. Wand, and T. Schultz,  “Estimation of fundamental frequency from surface electromyographic data: EMG-to-F0,”  in Proc. ICASSP, 2011, pp. 573–576.
  * [14] Y. Deng, J. T. Heaton, and G. S. Meltzner,  “Towards a practical silent speech recognition system,”  in Proc. Interspeech, 2014, pp. 1164–1168.
  * [15] L. Diener, M. Janke, and T. Schultz,  “Direct conversion from facial myoelectric signals to speech using Deep Neural Networks,”  in Proc. IJCNN, 2015, pp. 1–7.
  * [16] G. S. Meltzner, J. T. Heaton, Y. Deng, G. De Luca, S. H. Roy, and J. C. Kline,  “Silent speech recognition as an alternative communication device for persons with laryngectomy,”  IEEE/ACM Trans. ASLP, vol. 25, no. 12, pp. 2386–2398, 2017.
  * [17] M. Wand, T. Schultz, and J. Schmidhuber,  “Domain-adversarial training for session independent EMG-based speech recognition,”  in Proc. Interspeech, 2018, pp. 3167–3171.
  * [18] B. Denby, T. Schultz, K. Honda, T. Hueber, J. M. Gilbert, and J. S. Brumberg,  “Silent speech interfaces,”  Speech Communication, vol. 52, no. 4, pp. 270–287, 2010.
  * [19] V. Ramanarayanan, S. Tilsen, M. Proctor, J. Töger, L. Godstein, K.S. nazak, and S. Narayanan,  “Real-time control of an articulatory-based speech synthesizer for brain computer interfaces,”  Computer Speech and Language, vol. 52, pp. 1–22, 2018.
  * [20] A Toutis, D . Byrd, L. Goldstein, and S. Narayanan,  “Advances in vocal tract imaging and analysis,”  in The Routledge Handbook of Phonetics, 2019, pp. 34–50.
  * [21] K.I. Douros, A. Katsamanis, and P. Maragos,  “Multi-view audio-articulatory features for phonetic recognition on RTMRI-TIMIT database,”  in Proc. ICASSP, 2018, pp. 5514–5518.
  * [22] A. Katsamanis, E. Bresch, V. Ramanarayanan, and S. Narayanan,  “Validating rt-MRI based articulatory representations via articulatory recognition,”  in Proc. Interspeech, 2011, pp. 2841–2844.
  * [23] P. Saha, P. Srungarapu, and S. Fels,  “Towards automatic speech identification from vocal tract shape dynamics in real-time MRI,”  in Proc. Interspeech, 2018, pp. 1249–1253.
  * [24] H. Li, J. Tao, M. Yang, and B. Liu,  “Estimate articulatory mri series from acoustic signal using deep architecture,”  in Proc. ICASSP, 2015, pp. 4854–4858.
  * [25] A. Toutios, T. Sorensen, K. Somandepalli, R. Alexander, and S. S. Narayanan,  “Articulatory synthesis based on real-time magnetic resonance imaging data,”  in Proc. Interspeech, 2016, pp. 1492–1496.
  * [26] T. G. Csapó,  “Speaker-dependent articulatory-to-acoustic mapping using real-time MRI of the vocal tract,”  in Proc. Interspeech, 2020, pp. 2722–2726.
  * [27] S. Narayanan, A. Toutios, V. Ramanarayanan, A. Lammert, J. Kim, S. Lee, K. Nayak, Y-C. Kim, Y. Zhu, L. Goldstein, D. Byrd, E. Bresch, P. Ghosh, Katsamanis A., and M. Proctor,  “Real-time magnetic resonance imaging and electromagnetic articulography database for speech production research (TC),”  Journal of the Acoustical Society of America, vol. 136, pp. 1307–1311, 2014.
  * [28] P. Govalkar, J. Fisher, F. Zalkov, and C. Dittmar,  “A comparison of recent neural vocoders for speech signal reconstruction,”  in Proc. ISCA Speech Synthesis Workshop, 2019.
  * [29] K. Kumar, R. K. Kumar, T. de Boissiere, L. Gestin, W. Z. Teoh, J. Sotelo, A. de Brébisson, Y. Bengio, and A. C. Courville,  “MelGAN: Generative adversarial networks for conditional waveform synthesis,”  in Advances in Neural Information Processing Systems, 2019, vol. 32, pp. 14910–14921.
  * [30] R. Prenger, R. Valle, and B. Catanzaro,  “Waveglow: A flow-based generative network for speech synthesis,”  in Proc. ICASSP, 2019, pp. 3617–3621.
  * [31] T. G. Csapó, C. Zainkó, L. Tóth, G. Gosztolya, and A. Markó,  “Ultrasound-Based Articulatory-to-Acoustic Mapping with WaveGlow Speech Synthesis,”  in Proc. Interspeech 2020, 2020, pp. 2727–2731.
  * [32] L. Tóth and A. H. Shandiz,  “3D convolutional neural networks for ultrasound-based silent speech interfaces,”  in Proc. ICAISC. Springer, 2020, pp. 159–169.
  * [33] P. Ramachandran, B. Zoph, and Q. V. Le,  “Swish: a Self-Gated Activation Function,”  ArXiv e-prints 1710.05941, 2017.
  * [34] A. Pandey and D. Wang,  “On adversarial training and loss functions for speech enhancement,”  in Proc. ICASSP, 2018, pp. 5414–5418.
  * [35] S. Hochreiter and J. Schmidhuber,  “Long short-term memory,”  Neural computation, vol. 9, no. 8, pp. 1735–1780, 1997.
  * [36] D. Tran, H. Wang, L. Torresani, J. Ray, Y. LeCun, and M. Paluri,  “A closer look at spatiotemporal convolutions for action recognition,”  in Proc. CVPR, 2018.
  * [37] R.F. Kubichek,  “Mel-cepstral distance measure for objective speech quality assessment,”  in Proc. ICASSP, 1993, pp. 125–128.
  * [38] M. Kolbæk, Z.-H. Tan, S.H. Jensen, and J. Jensen,  “On loss functions for supervised monaural time-domain speech enhancement,”  IEEE/ACM Trans. ASLP, vol. 28, pp. 825–838, 2020.
  * [39] C.H. Taal, R.C. Hendriks, R. Heusdens, and J. Jensen,  “An algorithm for intelligibility prediction of time–frequency weighted noisy speech,”  IEEE Trans. ASLP, vol. 19, no. 7, pp. 2125–2136, 2011.
  * [40] J. Martín-Doñas, A. Gomez, Gonzalez L.J., and A. Peinado,  “A deep learning loss function based on the perceptual evaluation of the speech quality,”  IEEE Signal Processing Letters, vol. 25, no. 11, pp. 1680 – 1684, 2018.
  * [41] J. Le Roux, S. Wisdom, H. Erdogan, and J.R. Hershey,  “SDR - half-baked or well done?,”  in Proc. ICASSP, 2019.


11institutetext: Computer Science and Engineering, Michigan State University
11email<EMAIL_ADDRESS>
# Multiplication and Modulo are Lattice Linear

Arya Tanmay Gupta  Sandeep S Kulkarni

###### Abstract

In this paper, we analyze lattice linearity of multiplication and modulo
operations. We demonstrate that these operations are lattice linear and the
parallel processing algorithms that we study for both these operations are
able to exploit the lattice linearity of their respective problems. This
implies that these algorithms can be implemented in asynchronous environments,
where the nodes are allowed to read old information from each other and are
still guaranteed to converge within the same time complexity. These algorithms
also exhibit properties similar to snap-stabilization, i.e., starting from an
arbitrary state, the system follows the trace strictly according to its
specification.

###### Keywords:

lattice linear modulo multiplication self-stabilization asynchrony

##  1 Introduction

Development of parallel processing algorithms to solve problems is
increasingly gaining interest. This is because computing machines are
manufactured with multiprocessor chips as we face a bound on the rate of
architectural development of individual microprocessors. Such algorithms, in
general, require synchronization among their processing nodes, otherwise the
nodes perform executions based on inconsistent values possibly resulting in
substantial delay or halt in progress.

In this context, lattice theory has provided with very useful concepts. In
lattice linear systems, a partial order is induced in the state space, which
allows the nodes to perform executions asynchronously. Lattice linearity was
utilized in modelling the problems (which naturally exhibit lattice linearity,
called lattice linear problems) [5] and in developing algorithms (called
lattice linear algorithms) which impose a lattice structure in problems (which
are not lattice linear) [7]. This allows the algorithms that traverse such a
state space can allow the nodes to perform executions asynchronously while
preserving correctness.

Lattice linearity allows inducing single or multiple lattices among the global
states. If self-stabilization is required, then for every induced lattice, its
supremum must be an optimal state. This way, it can be ensured that the system
can traverse to an optimal state from an arbitrary state. Eventually lattice
linear algorithms [7] were developed for problems that require self-
stabilization, such as minimal dominating set, minimal vertex cover, maximal
independent set and graph colouring. These algorithms induce multiple disjoint
lattices in a subset of the state space. Fully lattice linear algorithms were
developed for the minimal dominating set [8]; these algorithms induce multiple
lattices in the entire state space. On the other hand, lattice linear problems
form a single lattice among the states. Many such problems do not allow self-
stabilization [5, 6, 4].

While developing new algorithms for lattice linear and non-lattice linear
problems is very interesting, it is also worthwhile to study if algorithms
already present in the literature exploit lattice linearity of problems or if
lattice linearity is present in existing algorithms. For example, lattice
linearity was found to be exploited by Johnson’s algorithm for shortest paths
in graphs [5] and by Gale’s top trading cycle algorithm for housing market
[6].

In this paper, we focus on parallel implementations of two fundamental
operations in mathematics: multiplication and modulo. Specifically, we study
algorithms for $n\times m$ and $n\mod m$ where $n$ and $m$ are large integers
represented as binary strings. A crucial observation is that these operations
are lattice linear. Also, the algorithms that we study in this paper are
capable of computing the correct answer even if the nodes are initialized in
an arbitrary state. These properties are present in these operations as
opposed to many other lattice linear problems where the lattice structure does
not allow self-stabilization.

###  1.1 Contributions of the paper

  * •

We study the lattice linearity of modulo and multiplication operations. We
analyze some algorithms which exploit lattice linearity of these operations.

  * •

The lattice linear versions of the algorithms exhibit properties similar to
snap-stabilizing algorithms present in [1], i.e., starting from an arbitrary
state, the system follows the trace strictly according to its specification.

  * •

The algorithms present in this paper induce a tree among the nodes. They can
be transformed into algorithms that require the nodes with number equal to the
nodes present in the lowest level. These transformed algorithms are also able
to exploit lattice linearity of these operations and converge within the same
time complexity bounds.

###  1.2 Organization of the paper

In Section 2, we describe the definitions and notations that we use in this
paper. In Section 3, we discuss the related work. We study lattice linearity
of multiplication operation in Section 4. We discuss some classic sequential
algorithms for modulo in Section 5. Then, in Section 6, we study the lattice
linearity modulo operation. Finally, we conclude in Section 7.

##  2 Preliminaries

This paper focuses on multiplication and modulo operations, where the operands
are $n$ and $m$. In the computation of $n\times m$ or $n\mod m$, $n$ and $m$
are the values of these numbers respectively, and $|n|$ and $|m|$ are the
length of the bitstrings required to represent $n$ and $m$ respectively. If
$n$ is a bitstring, then $n[i]$ is the $i^{th}$ bit of $n$ (indices starting
from 1), and $n[i:j]$ is equal to the bitstring from $i^{th}$ to $j^{th}$ bit
of $n$. For simplicity, we stipulate that $n$ and $m$ are of lengths in some
power of $2$. Since size of $n$ and $m$ may be substantially different, we
provide complexity results that are of the form $O(f(n,m))$ in all cases.

One of the ways to compute modulo, that we discuss in Section 5, is by using a
precomputed finite state automaton. A finite state automaton [11, p. 1-100]
$M$ is of the form $M=(Q$, $\Sigma$, $\delta$, $q_{0})$, where $Q$ is the set
of states of $M$, $\Sigma$ is the alphabet set of $M$, $\delta$ is the
transition function, and $q_{0}$ is the initial state of $M$.

We assume that we have access to the following string operations: (1)
$append(a$, $b)$, appends $b$ to the end of $a$ in $O(1)$ time, (2) $a$ >>
$k$, deletes rightmost $k$ bits of $a$ in $O(k)$ time, and (3) $a$ << $k$,
appends $k$ zeros to the right of $a$ in $O(k)$ time.

###  2.1 Modeling distributed programs

$n\times m$ or $n\mod m$ are typically thought of as arithmetic operations.
However, when $n$ and $m$ are large, we view them as algorithms. In this
paper, we view them as parallel/distributed algorithms where the nodes
collectively perform computations to converge to the final output. Next, we
provide the relevant definitions of a parallel/distributed program that we
utilize in this paper.

The parallel/distributed program consists of nodes where each node is
associated with a set of variables. A global state, say $s$, is obtained by
assigning each variable of each node a value from its respective domain. $s$
is represented as a vector such that $s[i]$ denotes the variables of node $i$,
where $s[i]$ itself is a vector of the variables of node $i$. A state space
$S$ is the set of all possible global states that a given system can obtain
under a given algorithm.

Each node is associated with actions. Each action at node $i$ checks the
values of other nodes and updates its own variables. An action is of the form
$g\longrightarrow a_{c}$, where $g$ is the guard (a Boolean predicate
involving variables of $i$ and other nodes) and $a_{c}$ is an instruction that
updates the variables of $i$. A trace is a sequence of global states $\langle
s_{0},s_{1},\cdots\rangle$ where $\forall j\geq 1,s_{j}$ is obtained by
executing some action in $s_{j-1}$.

###  2.2 Execution without synchronization.

Typically, we view the computation of an algorithm as a sequence of global
states $\langle s_{0},s_{1},\cdots\rangle$ where $s_{\ell+1}(\forall l\geq 0)$
is obtained by executing some action by one or more nodes in $s_{\ell}$. Under
proper synchronization, node $i$ evaluates its guards on the current local
states of its neighbours.

For the sake of discussion, assume that only node $i$ is executing in state
$s_{\ell}$. To understand how execution without synchronization works, let
$x(s)$ denote the value of some variable $x$ in state $s$, and let
$\forall\ell\geq 0,s_{\ell+1}$ be the state obtained by executing
$g\longrightarrow a_{c}$ (of process $i$) from $s_{\ell}$. If $g$ is evaluated
without synchronization, then node $i$ may read old values of some variables.
Let the computation prefix until $s_{\ell}$ be $\langle
s_{0},s_{1},\cdots,s_{\ell}\rangle$. Then, the state may appear to $i$ as
$s^{\prime}\neq s_{\ell}$, where some values are old. Specifically, the state
$s_{\ell}$ will appear as $s^{\prime}$ where
$x(s^{\prime})\in\\{x(s_{0}),x(s_{1}),\cdots,x(s_{\ell})\\}.$ If $g$ evaluates
to false, then the next state will stay to be $s_{\ell}$, which may be
incorrect. If, otherwise, $g$ evaluates to true then node $i$ will execute
$a_{c}$, and the resulting state $s^{\prime\prime}$ may not be the same as
$s_{\ell+1}$, where (1) the state of $i$ will be $s^{\prime\prime}[i]$, i.e.,
the value of node $i$ in state $s^{\prime\prime}$, and (2) $\forall j\neq i$,
the state of $j$ will be unchanged, i.e., $s_{\ell+1}[j]=s_{\ell}[j]$.

###  2.3 Embedding a $<$-lattice in global states.

Next, we discuss the structure of a partial order in the state space which,
under proper constraints, allows an algorithm to converge to an optimal state.
To describe the embedding, let that $s$ denotes a global state, and let that
$s[i]$ denotes the state of node $i$ in $s$. First, we define a total order
$<_{l}$; all local states of a node $i$ are totally ordered under $<_{l}$.
Using $<_{l}$, we define a partial order $<_{g}$ among global states as
follows:

We say that $s<_{g}s^{\prime}$ iff $(\forall i:s[i]=s^{\prime}[i]\lor
s[i]<_{l}s^{\prime}[i])\land(\exists i:s[i]<_{l}s^{\prime}[i])$. Also,
$s=s^{\prime}$ iff $\forall i:s[i]=s^{\prime}[i]$. For brevity, we use $<$ to
denote $<_{l}$ and $<_{g}$: $<$ corresponds to $<_{l}$ while comparing local
states, and $<$ corresponds to $<_{g}$ while comparing global states. We use
the symbol ‘$>$’ which is the opposite of ‘$<$’, i.e., $s>s^{\prime}$ iff
$s^{\prime}<s$. Similarly, we use symbols ‘$\leq$’ and ‘$\geq$’; e.g., $s\leq
s^{\prime}$ iff $s=s^{\prime}\vee s<s^{\prime}$. We call the lattice, formed
from such partial order, a $<$-lattice.

###### Definition 1

$<$-lattice. Given a total relation $<_{l}$ that orders the values of $s[i]$
(the local state of node $i$ in state $s$), the $<$-lattice corresponding to
$<_{l}$ is defined by the following partial order: $s<s^{\prime}$ iff
$(\forall i\ \ s[i]\leq_{l}s^{\prime}[i])\wedge(\exists i\ \
s[i]<_{l}s^{\prime}[i])$.

In the $<$-lattice discussed above, we can define the meet and join of two
states in the standard way: the meet (respectively, join), of two states
$s_{1}$ and $s_{2}$ is a state $s_{3}$ where $\forall i,s_{3}[i]$ is equal to
$min(s_{1}[i],s_{2}[i])$ (respectively, $max(s_{1}[i],s_{2}[i])$).

Remark: The $<$-lattice alone is insufficient for permitting asynchronous
execution. Additional constraints are identified in the following part of this
section.

By varying $<_{l}$ that identifies a total order among the states of a node,
one can obtain different lattices. A $<$-lattice, embedded in the state space,
is useful for permitting the algorithm to execute asynchronously. Under proper
constraints on how the lattice is formed, convergence is ensured.

###  2.4 Lattice Linear problems

Next, we discuss lattice linear problems, i.e., the problems where the problem
statement creates the lattice structure automatically. Such problems can be
represented by a predicate that induces a lattice among the states in $S$.

In lattice linear problems, a problem $P$ can be represented by a predicate
$\mathcal{P}$ such that for any node $i$, if it is violating $\mathcal{P}$ in
some state $s$, then it must change its state, otherwise the system will not
satisfy $\mathcal{P}$. Let $\mathcal{P}(s)$ be true iff state $s$ satisfies
$\mathcal{P}$. A node violating $\mathcal{P}$ in $s$ is called a forbidden
node. Formally,

###### Definition 2

[5] Forbidden (local state)
$\textsc{Forbidden}(i,s,\mathcal{P})\equiv\lnot\mathcal{P}(s)\land(\forall
s^{\prime}>s:s^{\prime}[i]=s[i]\implies\lnot\mathcal{P}(s^{\prime}))$.

Definition 2 implies that if a node $i$ is forbidden in some state $s$, then
in any state $s^{\prime}$ such that $s^{\prime}>s$, if the state of $i$
remains the same, then the algorithm will not converge. We define lattice
linear problems as follows.

###### Definition 3

Lattice linear problems. A problem $P$ is lattice linear iff there exists a
predicate $\mathcal{P}$ and an $<$-lattice such that

  * •

$P$ requires that we reach a state where $\mathcal{P}$ is true, and

  * •

$\mathcal{P}$ is lattice linear with respect to the $<$-lattice induced in $S$
[5], i.e., $\forall s:\neg\mathcal{P}(s)\Rightarrow\exists
i:\textsc{Forbidden}(i,s,\mathcal{P})$.

Remark: We study lattice structures that can be induced in multiplication and
modulo: we show that multiplication and modulo are lattice linear problems.
All the lattice structures that we study in this paper allow self-
stabilization: the supremum of the lattice induced in the state space is the
optimal state.

###  2.5 Role of Lattice Linearity in Concurrent Execution

Suppose $\textsc{Forbidden}(i,s,\mathcal{P})$ is false where $s^{\prime}$ is
the state of the system observed by node $i$. However, the actual state of the
system is $s$, as other nodes have changed their states and node $i$ is not
aware of it. In this case also, by definition of
$\textsc{Forbidden}(i,s^{\prime},\mathcal{P})$, node $i$ is still permitted to
execute this action.

In a typical program, checking that a guard was true in $s^{\prime}$ does not
provide any guarantee in determining it is still true in $s$. However, if a
$<$-lattice is induced, then node $i$ needs to only ensure that it does not
execute unless it reads the information from all concerned nodes since its
last execution. This means that all the nodes can execute in asynchrony.

##  3 Related Work

Lattice linear problems: The notion of representing problems through a
predicate under which the states form a lattice was introduced in [5]. We call
the problems for which such a representation is possible lattice linear
problems. Lattice linear problems are studied in [5, 6, 4], where the authors
have studied lattice linearity in, respectively, housing market problem and
several dynamic programming problems. Many of these problems are not self-
stabilizing.

Snap-stabilization: The notion of snap-stabilization was introduced in [1].

###### Definition 4

Snap-stabilizing algorithms. An algorithm is snap-stabilizing iff starting
from any arbitrary state in the state space, it makes the system follow a
trace as per the specification of that system.

For example, the fully lattice linear algorithms [8] follow the trace strictly
according to the required specification for any system that initializes in an
arbitrary feasible state, even if the nodes run asynchronously.

The algorithms that we study in this paper make the system follow a sequence
of states that are deterministically predictable because of the underlying
lattice structure in the state space. Thus, they exhibit snap-stabilization.

Modulo: In [12], authors have presented parallel processing algorithms for
inverse, discrete roots, or a large power modulo a number that has only small
prime factors. A hardware circuit implementation for mod is presented in [2].

In this paper, we present several parallel processing algorithms which are
self-stabilizing. Some require critical preprocessing, and some do not. The
general algorithm for modulo (Algorithm 2) utilizes the power of (sequential
or parallel) modulo and multiplication operations on smaller operands.

Multiplication: In [3], authors presented three parallel implementations of
the Karatsuba algorithm for long integer multiplication on a distributed
memory architecture. Two of the implementations have time complexity of $O(n)$
on $n^{\lg 3}$ processors. The third algorithm has complexity $O(n\lg n)$ on
$n$ processors.

In this paper, we study lattice linearity of parallelized Karatsuba algorithm.

##  4 Parallelized Multiplication Operation

In this section, we study the parallelization of the classic Karatsuba’s
algorithm for multiplication [9]. The parallelization was presented in [3]; we
study the lattice linearity of this algorithm. First we discuss the idea
behind sequential Karatsuba’s algorithm, and then we elaborate on its lattice
linearity.

###  4.1 Key Idea of the Sequential Karatsuba’s Algorithm

In this part, we discuss the sequential Karatsuba’s algorithm. The input is a
pair of bitstrings $n$ and $m$. This algorithm is recursive in nature. As the
base case, if the length of $n$ and $m$ equals 1 then, the multiplication
result is trivial. When the length is greater than 1, we let $m=ab$ and
$n=cd$, where $a,b,c$ and $d$ are half the length of $m$ and $n$ respectively.
Here, $ab$, for example, represents concatenation of $a$ and $b$, which equals
$m$.

Let $z=2^{|b|}$. $n\times m$ can be computed as $ac\times z^{2}+(ad+bc)\times
z+bd$. $ad+bc$ can be computed as $(a+b)\times(c+d)-ac-bd$. Thus, to compute
$m\times n$, it suffices to compute 3 multiplications $a\times c$, $b\times d$
and $(a+b)\times(c+d)$. Hence, we can eliminate one of the multiplications. In
the following section, we analyse the lattice linearity of the parallelization
of this algorithm as described in [3].

###  4.2 The CM parallelization [3] for Karatsuba’s protocol

This algorithm induces a tree among the nodes. Every non-leaf node has three
children. This algorithm works in two phases, top-down and bottom-up.

In the top-down phase, if $|m.i|>1$ or $|n.i|>1$ a node $i$ writes (1) $a$ and
$c$ to its left child, node $3i-1$ ($m.(3i-1)=a.i$ and $n.(3i-1)=c.i$), (2)
$b$ and $d$ to its middle child, node $3i$ ($m.(3i)=b.i$ and $n.(3i)=d.i$),
and (3) $a+b$ and $c+d$ to its right child, node $3i+1$ ($m.(3i+1)=a.i+b.i$
and $n.(3i+1)=c.i+d.i$). If $|m.i|=|n.i|=1$, i.e., in the base case, the
bottom-up phase begins and node $i$ sets $ans.i=m.i\times n.i$ that can be
computed trivially since $|m.i|=|n.i|=1$.

In the bottom-up phase, node $i$ sets $ans.i=ans.(3i-1)\times
z^{2}+(ans.(3i+1)-(ans.(3i-1)+ans.(3i))\times z+ans.(3i)$. Notice that
multiplication by $z$ and $z^{2}$ corresponds to bit shifts and does not need
an actual multiplication. Consequently, the product of $m\times n$ for node
$i$ is computed by this algorithm. With some book-keeping (storing the place
values of most significant bits of $a+b$ and $c+d$), a node $i$ only needs to
write the rightmost $\frac{|m.i|}{2}$ and $\frac{|n.i|}{2}$ bits to its
children. Thus, we can safely assume that when a node writes $m$ and $n$ to
any of its children, then $m$ and $n$ of that child are of equal length and
are of length in some power of 2. (If we do not do the book-keeping, the
number of nodes required increases, this number is upper bounded by $|n|^{2\lg
3}$ as the number of levels is upper bounded by $2\lg|n|$; this observation
was not made in [3].) However, we do not show the same in the algorithm for
brevity.

This algorithm utilizes $shift$ to compute $z$. A node $i$ evalautes $shift.i$
by doubling the value of $shift$ from its children. A node $i$ evaluates that
it is forbidden because of incorrect value of $shift.i$ by evaluating the
following macro.

$\begin{array}[]{c}\textsc{Forbidden-Multiplication-Karatsuba-
Shift}(i)\equiv\\\ \begin{array}[]{l}\begin{cases}|m.i|=1\land|n.i|=1\land
shift.i\neq 0~{}~{}OR\\\ shift.(3i)=shift.(3i-1)=0\leq shift.(3i+1)\land
shift.i\neq 1~{}~{}OR\\\ 0<shift.(3i)=shift.(3i-1)\leq shift.(3i+1)\land
shift.i\neq shift.(3i)*2.\\\ \end{cases}\end{array}\end{array}$

To ensure that the data flows down correctly, we declare a node $i$ to be
forbidden as follows.

$\begin{array}[]{c}\textsc{Forbidden-Multiplication-Karatsuba-
TopDown}(i)\equiv\\\ \begin{array}[]{l}\begin{cases}i=1\land(m.i\neq m\lor
n.i\neq n)~{}~{}OR\\\ ((|m.i|>1\land n.i>1)\land\\\ (m.(3i-1)\neq
m.i\Big{[}1:\frac{|m.i|}{2}\Big{]}~{}~{}OR\\\ n.(3i-1)\neq
n.i\Big{[}1:\frac{|n.i|}{2}\Big{]}~{}~{}OR\\\ m.(3i)\neq
m.i\Big{[}\frac{|m.i|}{2}+1:|m.i|\Big{]}~{}~{}OR\\\ n.(3i)\neq
n.i\Big{[}\frac{|n.i|}{2}+1:|n.i|\Big{]}~{}~{}OR\\\ m.(3i+1)\neq
m.i\Big{[}1:\frac{|m.i|}{2}\Big{]}+m.i\Big{[}\frac{|m.i|}{2}+1:|m.i|\Big{]}~{}~{}OR\\\
n.(3i+1)\neq
n.i\Big{[}1:\frac{|n.i|}{2}\Big{]}+n.i\Big{[}\frac{|n.i|}{2}+1:|n.i|\Big{]})).\\\
\end{cases}\end{array}\end{array}$

To determine if a node $i$ has stored $ans.i$ incorrectly, it evaluates to be
forbidden as follows.

$\begin{array}[]{c}\textsc{Forbidden-Multiplication-Karatsuba-
BottomUp}(i)\equiv\\\ \begin{array}[]{l}\begin{cases}|m.i|=1\land|n.i|=1\land
ans.i\neq m.i\times n.i~{}~{}OR\\\ |m.i|>1\land|n.i|>1\land(ans.i\neq
ans.(3i-1)$ << $shift.i\\\ \quad\quad+(ans.(3i+1)-ans.(3i-1)-ans.(3i))$ <<
$shift.(3i)\\\ \quad\quad+ans.(3i+1))\end{cases}\end{array}\end{array}$

Thus, Algorithm 1 is described as follows:

###### Algorithm 1

Parallel processing version of Karatsuba’s algorithm.

$\begin{array}[]{|l|}\hline\cr\text{{Rules for node $i$}}.\\\ {}\\\
\textsc{Forbidden-Multiplication-Karatsuba-Shift}(i)\longrightarrow\\\
\begin{array}[]{l}\begin{cases}shift.i=0&\text{if $|m.i|=1\land|n.i|=1\land
shift.i\neq 0$.}\\\ shift.i=1&\text{if $shift.(3i)=shift.(3i-1)=0\leq
shift.(3i+1)\land shift.i\neq 1$}\\\ shift.i=shift.(3i)\times
2&\text{otherwise}\end{cases}\end{array}{}\\\ \textsc{Forbidden-
Multiplication-Karatsuba-TopDown}(i)\longrightarrow\\\
\begin{array}[]{l}\begin{cases}m.i=m,n.i=n&\text{if $i=1\land(m.i\neq m\lor
n.i\neq n)$.}\\\ m.(3i-1)=m.i\Big{[}1:\frac{|m.i|}{2}\Big{]}&\text{if
$m.(3i-1)\neq m.i\Big{[}1:\frac{|m.i|}{2}\Big{]}$.}\\\
n.(3i-1)=n.i\Big{[}1:\frac{|n.i|}{2}\Big{]}&\text{if $n.(3i-1)\neq
n.i\Big{[}1:\frac{|n.i|}{2}\Big{]}$.}\\\
m.(3i)=m.i\Big{[}\frac{|m.i|}{2}+1:|m.i|\Big{]}&\text{if $m.(3i)\neq
m.i\Big{[}\frac{|m.i|}{2}+1:|m.i|\Big{]}$.}\\\
n.(3i)=n.i\Big{[}\frac{|n.i|}{2}+1:|n.i|\Big{]}&\text{if $n.(3i)\neq
n.i\Big{[}\frac{|n.i|}{2}+1:|n.i|\Big{]}$.}\\\
m.(3i+1)=m.i\Big{[}1:\frac{|m.i|}{2}\Big{]}+m.i\Big{[}\frac{|m.i|}{2}+1:|m.i|\Big{]}&\text{if
$m.(3i+1)\neq m.i\Big{[}1:\frac{|m.i|}{2}\Big{]}$.}\\\
&\text{\quad\quad$+m.i\Big{[}\frac{|m.i|}{2}+1:|m.i|\Big{]}$.}\\\
n.(3i+1)=n.i\Big{[}1:\frac{|n.i|}{2}\Big{]}+n.i\Big{[}\frac{|n.i|}{2}+1:|n.i|\Big{]}&\text{otherwise}\\\
\end{cases}\end{array}{}\\\ \textsc{Forbidden-Multiplication-Karatsuba-
BottomUp}(i)\longrightarrow\\\ \begin{array}[]{l}\begin{cases}ans.i=m.i\times
n.i&\text{if $|m.i|=1\land|n.i|=1$}\\\ &\text{\quad\quad$\land ans.i\neq
m.i\times n.i$}\\\ ans.i=ans.(3i-1)$ << $shift.i+(ans.(3i+1)-ans.(3i-1)\\\
\quad\quad-ans.(3i))$ <<
$shift.(3i)+ans.(3i+1))&\text{otherwise.}\end{cases}\end{array}{}\\\
\hline\cr\end{array}$

From the above description, we see that Algorithm 1 satisfies the constraints
of lattice linearity. This algorithm executes in $O(|n|)$ time. Its work
complexity is $O(n^{\lg 3})$, which is same as the time complexity of the
sequential Karatsuba’s algorithm. The following example illustrates the
working of Algorithm 1.

###### Example 1

Figure 1 evaluates $100\times 100$ following Algorithm 1. ∎

$m=0100$
---
$n=0100$
node 1 $m=01$
---
$n=01$
node 2 $m=00$
---
$n=00$
node 3 $m=01$
---
$n=01$
node 4 $m=0$
---
$n=0$
node 5 $m=1$
---
$n=1$
node 6 $m=1$
---
$n=1$
node 7 $m=0$
---
$n=0$
node 8 $m=0$
---
$n=0$
node 9 $m=0$
---
$n=0$
node 10 $m=0$
---
$n=0$
node 11 $m=1$
---
$n=1$
node 12 $m=1$
---
$n=1$
node 13

(a) $m=0100$
---
$n=0100$
$mn=10000$
node 1 $m=01$
---
$n=01$
$mn=1$
node 2 $m=00$
---
$n=00$
$mn=0$
node 3 $m=01$
---
$n=01$
$mn=1$
node 4 $m=0$
---
$n=0$
$mn=0$
node 5 $m=1$
---
$n=1$
$mn=1$
node 6 $m=1$
---
$n=1$
$mn=1$
node 7 $m=0$
---
$n=0$
$mn=0$
node 8 $m=0$
---
$n=0$
$mn=0$
node 9 $m=0$
---
$n=0$
$mn=0$
node 10 $m=0$
---
$n=0$
$mn=0$
node 11 $m=1$
---
$n=1$
$mn=1$
node 12 $m=1$
---
$n=1$
$mn=1$
node 13 (b)

Figure 1: Demonstration of multiplication of 100 and 100 in base 2: (a) top
down (b) bottom up.

###### Theorem 4.1

Given the input bitstrings $n$ and $m$, the predicate

$\forall i\lnot(\textsc{Forbidden-Multiplication-Karatsuba-Shift}(i)\lor$

$\textsc{Forbidden-Multiplication-Karatsuba-TopDown}(i)\lor$

$\textsc{Forbidden-Multiplication-Karatsuba-BottomUp}(i))$

is lattice linear on $|n|^{\lg 3}$ threads.

###### Proof

Let us assume that node 1 does not have the correct value of $ans.1=n\times
m$. This implies that (1) node 1 does not have an updated value in $n.1$ or
$m.1$, or (2) node 1 has a non-updated value of $ans.1$, (3) node 1 has not
written the updated values to $n.2$ & $m.2$ or $n.3$ & $m.3$ or $n.4$ & $m.4$,
(4) node 1 has a non-updated value in $shift.1$, or (5) nodes 2, 3 or 4 have
incorrect values in their respective $n$, $m$, $ans$ or $shift$ variables. In
cases (1),…,(4), node 1 is forbidden.

Recursively, this can be extended to any node $i$. Let that node $i$ has
stored an incorrect value in $ans.i$ or $shift.i$. Let that $i>1$. Then (1)
node $i$ has a non-updated value in $shift.i$, $ans.i$, $n.i$ or $m.i$, or (2)
(if $m.i>1$ or $n.i>1$) node $i$ has not written updated values to $n.(3i-1)$
& $m.(3i-1)$ or $n.(3i)$ & $m.(3i)$ or $n.(3i+1)$ & $m.(3i+1)$, in which case
node $i$ is forbidden. In both these cases, node $i$ is forbidden. It is also
possible that at least one of the children of node $i$ has incorrect values in
its respective $n$, $m$, $ans$ or $shift$ variables. From these cases, we have
that given a global state $s$, where $s=\langle\langle n.1,$ $m.1$,
$ans.1\rangle$, $\langle n.2$, $m.2$, $ans.2\rangle$, $...$, $\langle
n.(|n|^{\lg 3})$, $m.(|n|^{\lg 3})$, $ans.(|n|^{\lg 3})\rangle\rangle$, if $s$
is forbidden, there is at least one node which is forbidden. This shows that
if the global state is forbidden, then there exists some node $i$ which is
forbidden.

Next, we show that if some node is forbidden, then node 1 will not store the
correct answer. node 1 is forbidden if it has not read the correct value $m.1$
and $n.1$. Additionally, $\forall i:i\in[1:n^{2\lg 3}]$ node $i$ is forbidden
if (1) it has a non-updated value in $ans.i$ or $shift.i$, (2) $i$ has not
written the correct values to $ans.(3i-1)$ or $ans.(3i)$ or $ans.(3i+1)$. This
implies that the parent of node $i$ will also store incorrect value in its
$ans$ or $shift$ variable. Recursively, we have that node 1 stores an
incorrect value in $ans.1$. Thus, the global state is forbidden.

###### Corollary 1

Algorithm 1 computes multiplication of two numbers with $n^{\lg 3}$ threads
without using synchronization.

In Algorithm 1, we have allowed a parent update its children. It can be
transformed such that all nodes update themselves only, which will stay
lattice linear. This can be done by making all nodes to copy $n$ and $m$
values from their parents update their $n$ and $m$ values. With this change,
the definition of Forbidden-Multiplication-Karatsuba-TopDown$(i)$ will change
accordingly.

##  5 Modulo: some classic sequential models

In this section, we discuss some sequential algorithms for computing $n\mod
m$. We consider two instances, one where both $n$ and $m$ are inputs and
another where $n$ is an input but $m$ is hardcoded. We utilize these
algorithms in Section 6.

The latter algorithm is motivated by algorithms such as RSA [10] where the
value of $n$ changes based on the message to be encrypted/decrypted, but the
value of $m$ is fixed once the keys are determined. Thus, some pre-processing
can potentially improve the performance of the modulo operation; we observe
that certain optimizations are possible. While the time and space complexities
required for preprocessing in this algorithm are high, thereby making it
impractical, it demonstrates a gap between the lower and upper bound in the
complexity.

###  5.1 Modulo by Long Division

The standard long-division algorithm for computing $n$ mod $m$ is as shown
below.

Division-Modulo($n,m$):
---
$ans=n[1:|m|-1]$, $i=|m|$.
while $(i<|n|)$, {$ans=append(ans,n[i])$; if $ans>m$, then $ans=ans-m$; i++}.
Return $ans$.

Clearly, the number of iterations in this algorithm is bounded by $|n|$. In
each of these iterations, the worst case complexity is $O(|m|)$ to perform the
subtraction operation. Thus, the time complexity of standard long division is
$O(|n|\times|m|)$ when $m$ and $n$ are both inputs to the algorithm.

###  5.2 Modulo by constructing DFA

If the value of $m$ is hardcoded in the algorithm, we use it to reduce the
cost of the modulo operation by creating a DFA $M=(Q$, $\Sigma$, $\delta$,
$q_{0})$ where $Q=\\{q_{0}..q_{m-1}\\}$. If $M$ has read the first $k$ digits
of $n$ then its state would be $n[1:k]\mod m$. When $M$ reads the next digit
of $n$, $n[1:(k+1)]\mod m$ is evaluated depending on the next digit. If the
next digit is 0, the next state would be $(2\times(n[1:k]\mod m))\mod m$,
otherwise, it would be $(2\times(n[1:k]\mod m)+1)\mod m$. Since the value of
$m$ is hardcoded, this DFA can be pre-computed. As an example, for $m=3$, the
corresponding DFA is in the following example.

###### Example 2

A finite automaton $M_{3}$ computing $n\mod 3$ ($m=3$ is fixed) for any
$n\in\mathbb{N}$ is shown in Figure 2. ∎

$q_{0}$$q_{1}$$q_{2}$110001 Figure 2: A finite automaton $M_{3}$ computing
$n\mod 3$ for any $n\in\mathbb{N}$.

With this DFA, the cost of computing the modulo operation corresponds to one
DFA transition for each digit of $n$. Hence, the complexity of the
corresponding operation is $O(|n|)$. $M$ does not have any accepting states,
unlike any usual finite automaton. The final state of $M$ only tells us the
value of the remainder which we would obtain after the computation of the
expression $n\mod m$. We define the construction of the transition table
$\delta$ of $M$ in Figure 3.

Elaborated definition of $\delta$.
---
1\. $v=0$. $m_{1/2}=m$ >> $1$.
$b=$ last bit of $m$.
2\. for each $i:0\leq i\leq m_{1/2}$,
3\. $\delta(q_{i},0)=q_{v}$. $v=v+1$.
4\. $\delta(q_{i},1)=q_{v}$. $v=v+1$.
5\. if $b=1$, then
//$m-1$ is even, so
//$m-1$ is $2+2\times(m_{1/2}-1)$
//assuming that $i=m_{1/2}+1$ now
6\. $\delta(q_{i},0)=q_{v}$. $v=v+1$.
7\. $v=0$.
8\. $\delta(q_{i},1)=q_{v}$. $v=v+1$.
9\. $i=i+1$.
10\. else, then
11\. $v=0$.
12\. for each $m_{1/2}+1\leq i:i\leq m$,
13\. $\delta(q_{i},0)=q_{v}$. $v=v+1$.
14\. $\delta(q_{i},1)=q_{v}$. $v=v+1$.
Figure 3: Definition of the transition function $\delta$.

While this approach takes $|n|$ steps, each step taking a constant amount of
time, the time complexity (as well as the space complexity) of the required
preprocessing is $O(m\times|m|)$, which is very high. Therefore, this approach
is not practical when $m$ is large. However, we consider it to observe that we
obtain a time complexity of $O(|n|)$ to run this automaton to evaluate for the
modulo operation if $m$ can be hardcoded.

##  6 Parallel processing modulo operation

In this section, we present a parallel processing algorithm to compute $n\mod
m$ using $4\times|n|/|m|-1$ processors. This method induces a binary tree
among the computation nodes based on their ids; there are $2\times|n|/|m|$
processors in the lowest level (level 1). This algorithm uses three variables
to represent the state of each node $i$: $shift.i$, $pow.i$ and $ans.i$.

Computation of shift.i:  The variable $shift$ stores the required power of 2.
The value of $shift$ at any node at level 1 is $0$. At level 2, the value of
$shift$ at any node is $|m|/2$. At any higher level, the value of $shilft$ is
twice the value of shift of its children. Forbidden-Log-Modulo-Shift, in this
context, is defined below.

$\begin{array}[]{c}\textsc{Forbidden-Log-Modulo-Shift}(i)\equiv\\\
\begin{cases}shift.i\neq 0&\text{if $i\geq 2\times|n|/|m|$}\\\
shift.i\neq|m|/2&\text{if $shift.(2i)=shift.(2i+1)=0$}\\\ shift.i\neq 2\times
shift.(2i)&\text{if $shift.(2i)=shift.(2i+1)\geq|m|$}\end{cases}\end{array}$

Computation of pow.i:  The goal of this computation is to set $pow.i$ to be
$2^{shift.i}\mod m$, whenever $level.i\geq 2$. This can be implemented using
the following definition for Forbidden-Log-Modulo-Pow.

$\begin{array}[]{l}\textsc{Forbidden-Log-Modulo-Pow}(i)\equiv\\\
\begin{cases}pow.i\neq 1&\text{if $shift.i=0$}\\\ pow.i\neq
2^{\frac{|m|}{2}}&\text{if $shift.i=|m|/2$}\\\ pow.i\neq(pow.(2i))^{2}\mod
m&\text{otherwise}\end{cases}\end{array}$

By definition, $pow.i$ is less than $m$. Also, computation of $pow$ requires
multiplication of two numbers that are bounded by $|m|$. Hence, this
calculation can benefit from parallelization of Algorithm 1. However, as we
will see later that the complexity of this algorithm (for modulo) is dominated
by the modulo operation happening in individual nodes which is $O(|m|^{2})$,
we can use the sequential version of Karatsuba’s algorithm for multiplication
as well, without affecting the order of the time complexity of this algorithm.

Computation of ans.i:  We split $n$ into strings of size $\frac{|m|}{2}$.
Clearly, the number representing this substring is less than $m$. At the
lowest level (level 1), $ans.i$ is set to the corresponding substring. At
higher levels, $ans.i$ is set to $(pow.i\times ans.(2i)+ans.(2i+1))\mod m$.
This computation also involves multiplication of two numbers whose size is
upper bounded by $|m|$. A forbidden node $i$ from a non-updated $ans.i$ can be
evaluated using Forbidden-Log-Modulo-Ans$(i)$.

$\begin{array}[]{c}\textsc{Forbidden-Log-Modulo-Ans}(i)\equiv\\\
\begin{cases}ans.i\neq
n[(i-2\times\dfrac{|n|}{|m|})\times\dfrac{|m|}{2}+1:(i-2\times\dfrac{|n|}{|m|}+1)\times\dfrac{|m|}{2}]&\text{if
$shift.i=0$}\\\
ans.i\neq\textsc{Mod}(\textsc{Sum}(\textsc{Mul}(ans.(2i),pow.i),ans.(2i+1)),m)&\text{otherwise}\end{cases}\end{array}$

We describe the algorithm as Algorithm 2.

###### Algorithm 2

Modulo computation by inducing a tree among the nodes.

$\begin{array}[]{|l|}\hline\cr\text{{Rules for node} $i$.}\\\ {}\\\
\textsc{Forbidden-Log-Modulo-Shift}(i)\longrightarrow\\\
\begin{array}[]{ll}\begin{cases}shift.i=0&\text{if $i\geq 2\times|n|/|m|$}\\\
shift.i=\dfrac{|m|}{2}&\text{if $shift.(2i)=shift.(2i+1)=0$}\\\
shift.i=2\times shift.(2i)&\text{if
$shift.(2i)=shift.(2i+1)\geq|m|$}\end{cases}\end{array}{}\\\
\textsc{Forbidden-Log-Modulo-Pow}(i)\longrightarrow\\\
\begin{array}[]{ll}\begin{cases}pow.i=1&\text{if $shift.i=0$}\\\
pow.i=2^{\frac{|m|}{2}}&\text{if $shift.i=|m|/2$}\\\
pow.i=\textsc{Mod}(\textsc{Mul}(pow.i,pow.i),m)&\text{otherwise}\end{cases}\end{array}{}\\\
\textsc{Forbidden-Log-Modulo-Ans}(i)\longrightarrow\\\
\begin{array}[]{ll}\begin{cases}ans.i=n[(i-2\times\dfrac{|n|}{|m|})\times\dfrac{|m|}{2}+1:(i-2\times\dfrac{|n|}{|m|}+1)\times\dfrac{|m|}{2}]&\text{if
$shift.i=0$}\\\
ans.i=\textsc{Mod}(\textsc{Sum}(\textsc{Mul}(ans.(2i),pow.i),ans.(2i+1)),m)&\text{otherwise}\end{cases}\end{array}{}\\\
\hline\cr\end{array}$

The following example demonstrates the working of Algorithm 2.

###### Example 3

Figure 4 shows the computation of $11011\mod 11$ as performed by Algorithm 2.
∎

1+10 mod 11node 11 mod 11node 210+0 mod 11node 300 mod 11node 401 mod 11node
510 mod 11node 611 mod 11node 70node 80node 90node 101node 111node 120node
131node 141node 15$pow=1$$pow=10$$pow=100\mod 11~{}(=1)$$pow=10000\mod
11~{}(=1)$(=0) Figure 4: Processing $11011\mod 11$ following Algorithm 2.

###  6.1 Proof of correctness.

###### Theorem 6.1

Given the input bitstrings $n$ and $m$, the predicate

$\forall i\lnot(\textsc{Forbidden-Log-Modulo-Shift}(i)\lor$

$\textsc{Forbidden-Log-Modulo-Pow}(i)\lor$

$\textsc{Forbidden-Log-Modulo-Ans}(i))$

is lattice linear on $4|n|/|m|-1$ threads.

###### Proof

Let us assume that node 1 has incorrect value in $ans.1$. This implies that
(1) node 1 has a non-updated value in $shift.i$ or $pow.i$, or (2) node 1 does
not have an updated value in $ans.i$. in both these cases node 1 is forbidden.
It is also possible that node 2 or node 3 have an incorrect value in their
variables.

Recursively, this can be extended to any node $i$. Let that node $i$ has
stored an incorrect value in $ans.i$. If $i<2(|n|/|m|)$, then (1) node $i$ has
a non-updated value in $ans.i$ $pow.i$ or $shift.i$, in which case node $i$ is
forbidden or (2) node $2i$ or node $2i+1$ have incorrect value in their
respective $shift$, $pow$ or $ans$ variables. If $i\geq 2(|n|/|m|)$, then $i$
do not have values $shift.i=0$, $pow.i=1$ or a correct $ans.i$ value, in which
case $i$ is forbidden. From these cases, we have that given a global state
$s$, where $s=\langle\langle shift.1$, $pow.1$, $ans.1\rangle$, $\langle
shift.2$, $pow.2$, $ans.2\rangle$, $...$, $\langle shift.(4|n|/|m|-1)$,
$pow.(4|n|/|m|-1)$, $ans.(4|n|/|m|-1)\rangle\rangle$, if $s$ is forbidden,
there is at least one node which is forbidden.

This shows that if the global state is forbidden, then there exists some node
$i$ which is forbidden.

Next, we show that if some node is forbidden, then node 1 will not store the
correct answer. $\forall i:i\in[1:4|n|/|m|-1]$ node $i$ is forbidden if it has
a non-updated value in $ans.i$, $pow.i$ or $shift.i$. This implies that the
parent of node $i$ will also store incorrect value in its $ans$ variable.
Recursively, we have that node 1 stores an incorrect value in $ans.1$, and
thus the global state is forbidden.

###### Corollary 2

Algorithm 2 computes multiplication of two numbers with $n^{\lg 3}$ threads
without using synchronization.

###  6.2 Time complexity analysis

Algorithm 2 is a general algorithm that uses the Mod( Mul($\cdots$)) and Mod(
Sum($\cdots$)). These functions can be implemented in different ways. Choices
for these implementations affect the time complexity, which we discuss as
follows.

####  6.2.1 Modulo via Long Division.

First, we consider the standard approach for computing Mod(Mul($\cdots$)) and
Mod(Sum($\cdots$)). Observe that in Algorithm 2, if we compute Mod(Mul$(x,y)$)
then $x,y<m$. Hence, we can use Karatsuba’s parallelized algorithm from
Section 4 where the input numbers are less than $m$. Using the analysis from
Section 4, we have that each multiplication operation has a time complexity of
$O(|n|)$.

Subsequently, to compute the mod operation, we need to compute $xy\mod m$
where $xy$ is upto $2|m|$ digits long. Using the standard approach of long
division, we will need $m$ iterations where in each iteration, we need to do a
subtraction operation with numbers that $|m|$ digits long. Hence, the
complexity of this approach is $O(|m|^{2})$ per modulo operation. Since this
complexity is higher than the cost of multiplication, the overall time
complexity is $O(|m|^{2}\times\lg(|n|/|m|))$.

This approach uses $m$ and $n$ as inputs. Next, we consider the case where $m$
is hardcoded. As discussed in Section 5, these problems occur in practice.
Similar to Section 5.2, the pre-processing required in these algorithms makes
them impractical in practice. However, we present them to provide lower bounds
on the complexity of the modulo operation when $m$ is hardcoded, and to show
that there is a potential to reduce the complexity.

####  6.2.2 Modulo by using $M$.

We can use Algorithm 1 for multiplication; each multiplication operation will
have a time complexity of $O(|m|)$. Subsequently, to compute the mod
operation, we need to compute $xy\mod m$ where $xy$ is upto $2|m|$ digits
long. Using $M$, we will need $2|m|$ iterations; each iteration takes a
constant amount of time. Hence, the complexity of this approach is $O(|m|)$
per modulo operation. Thus, the overall time complexity is
$O(|m|\times\lg(|n|/|m|))$.

####  6.2.3 Modulo by Constructing Transition Tables.

If $m$ is fixed, we can create a table $\delta_{sum}$ of size $m\times m$
where an entry at location $(i,j)$ represents $i+j\mod m$ in $O(m^{2})$ time.
Using $\delta_{sum}$, we can create another transition table $\delta_{mul}$ of
size $m\times m$ where an entry at location $(i,j)$ represents $i\times j\mod
m$ in $O(m^{2})$ time. Using $\delta_{mul}$, the time complexity of a
Mod(Mul($\cdots$)) operation becomes $O(1)$. Hence, the overall complexity of
the modulo operation becomes $O(\lg(|n|/|m|))$.

##  7 Conclusion

Multiplication and modulo operations are among the fundamental mathematical
operations. Fast parallel processing algorithms for such operations reduce the
execution time of the applications which they are employed in. In this paper,
we showed that these problems are lattice linear. In this context, we studied
an algorithm by Cesari and Maeder [3] which is a parallelization of
Karatsuba’s algorithm for multiplication. In addition, we presented a parallel
processing algorithm for the modulo operation.

Presence of lattice linearity in problems and algorithms allow nodes to
execute asynchronously, where they can perform executions based on old data of
other nodes. This is specifically valuable in parallel algorithms where the
overhead of synchronization can be very high.

In these algorithms, even if the initial states of the nodes are arbitrary,
the system follows the trace strictly according to its specification, and
within the same time complexity bounds. This property is the same as that
envisioned in snap-stabilization [1]. When viewed with the fact that the
algorithms are also lattice linear, this implies that the supremum state in
the lattice is also the optimal state.

## References

  * [1] A. Bui, A.K. Datta, F. Petit, and V. Villain.  State-optimal snap-stabilizing pif in tree networks.  In Proceedings 19th IEEE International Conference on Distributed Computing Systems, pages 78–85, 1999.
  * [2] J.T. Butler and T. Sasao.  Fast hardware computation of x mod z.  In 2011 IEEE International Symposium on Parallel and Distributed Processing Workshops and Phd Forum, pages 294–297, 2011.
  * [3] Giovanni Cesari and Roman Maeder.  Performance analysis of the parallel karatsuba multiplication algorithm for distributed memory architectures.  Journal of Symbolic Computation, 21(4):467–473, 1996.
  * [4] Vijay Garg.  A lattice linear predicate parallel algorithm for the dynamic programming problems.  In 23rd International Conference on Distributed Computing and Networking, ICDCN 2022, page 72–76, New York, NY, USA, 2022. Association for Computing Machinery.
  * [5] Vijay K. Garg.  Predicate detection to solve combinatorial optimization problems.  In Christian Scheideler and Michael Spear, editors, SPAA ’20: 32nd ACM Symposium on Parallelism in Algorithms and Architectures, Virtual Event, USA, July 15-17, 2020, pages 235–245. ACM, 2020.
  * [6] Vijay K. Garg.  A lattice linear predicate parallel algorithm for the housing market problem.  In Colette Johnen, Elad Michael Schiller, and Stefan Schmid, editors, Stabilization, Safety, and Security of Distributed Systems, pages 108–122, Cham, 2021. Springer International Publishing.
  * [7] Arya Tanmay Gupta and Sandeep S. Kulkarni.  Extending lattice linearity for self-stabilizing algorithms.  In Colette Johnen, Elad Michael Schiller, and Stefan Schmid, editors, Stabilization, Safety, and Security of Distributed Systems, pages 365–379, Cham, 2021. Springer International Publishing.
  * [8] Arya Tanmay Gupta and Sandeep S. Kulkarni.  Brief announcement: Fully lattice linear algorithms.  In Stabilization, Safety, and Security of Distributed Systems: 24th International Symposium, SSS 2022, Clermont-Ferrand, France, November 15–17, 2022, Proceedings, page 341–345, Berlin, Heidelberg, 2022. Springer-Verlag.
  * [9] A. Karatsuba and Y. Ofman.  Multiplication of many-digital numbers by automatic computers.  Doklady Ak- ademii Nauk SSSR, 14(145):293–294, 1962.
  * [10] R. L. Rivest, A. Shamir, and L. Adleman.  A method for obtaining digital signatures and public-key cryptosystems.  Commun. ACM, 21(2):120–126, feb 1978.
  * [11] Michael Sipser.  Introduction to the Theory of Computation.  Cengage Learning, 3 edition, 2013.
  * [12] Thomas Zeugmann.  Highly parallel computations modulo a number having only small prime factors.  Information and Computation, 96(1):95–114, 1992.