#  Solution to the Riemann Hypothesis from geometric analysis of component
series functions in the functional equation of zeta

Jeet Kumar Gaur
PhD in Mechanical Eng. Dept.
Indian Institute of Science, Bangalore, India
email id<EMAIL_ADDRESS>
###### Abstract

This paper presents a new approach towards the Riemann Hypothesis. On
iterative expansion of integration term in functional equation of the Riemann
zeta function we get sum of two series function. At the ‘non- trivial’ zeros
of zeta function, value of the series is zero. Thus, Riemann hypothesis is
false if that happens for an ‘s’ off the line $\Re(s)=1/2$ ( the critical
line). This series has two components $f(s)$ and $f(1-s)$. For the hypothesis
to be false one component is additive inverse of the other. From geometric
analysis of spiral geometry representing the component series functions $f(s)$
and $f(1-s)$ on complex plane we find by contradiction that they cannot be
each other’s additive inverse for any $s$, off the critical line. Thus,
proving truth of the hypothesis.

##  1 Introduction

The Riemann Hypothesis has become one of the most centralized problems in
mathematics today. It has to do with the position of ‘non-trivial’ zeros of
the zeta function that encrypt information about prime numbers. The hypothesis
has its roots predominantly in number theory. However, over time the ways to
approach the problem have diversified with new evidence suggesting connections
with seemingly unrelated fields like random matrices, chaos theory and quantum
physics [1][2]. In over 160 years of period, from the time it was first
formulated in 1859, Riemann hypothesis has been attempted by several great
mathematical minds. But, yet it has remained unsolved. Solving Riemann
hypothesis would also verify hundreds of theorems and solutions which are
based on the assumption of it’s truth.

Owing to the works of Turing, Hardy, Siegel, Chebyshev and others, better
algorithms to numerically evaluate the zeta function have been developed [3].
As of now the first 10 billion zeros have been numerically verified to follow
the hypothesis without a single exception [3]. Although computational pieces
of evidence speak in favour of the hypothesis, an analytic approach to the
hypothesis is necessary to prove it for each and every zero.

The Riemann zeta function, $\zeta(s)$ is defined by a Dirichlet series as

$\zeta(s)=\sum\limits_{n=1}^{\infty}\frac{1}{n^{s}}=1+\frac{1}{2^{s}}+\frac{1}{3^{s}}+\frac{1}{4^{s}}+...$

which is only convergent in the region $\Re(s)>1$, [4].

The zeta function was first known to the Swiss Mathematician Leonhard Euler.
In 1737, Euler gave the famous form of the zeta function as a product of
functions of prime numbers [5][6].

$\zeta(s)=\prod\limits_{p}\left(1-\frac{1}{p^{s}}\right)^{-1}.$

The above equation also represents the secret association of prime numbers
with the zeta function. However, it was Bernhard Riemann who first studied the
complex form of $\zeta(s)$ and also found a functional equation that
analytically continues it as a meromorphic function. He found that the
concrete prediction for the occurrence of prime numbers depends on the
placement of ‘non-trivial’( complex) zeros of this function. If the Riemann
hypothesis is true then all these zeros lie collinearly. The values of the
function on negative integers is associated to the corresponding Bernoulli
number. The Bernoulli numbers can be obtained from Euler-Maclauren
series(1735) [7][8]. The $n$th Bernoulli number $B^{+}_{n}$ can be defined as

$\displaystyle{B_{n}^{+}=-n*\zeta(1-n)},$

for all $n\geq 1$. Every odd Bernoulli number $B^{+}_{2n+1}$ is zero so, the
value of zeta function at every negative even integer. These are the ‘trivial’
(real) zeros of zeta function

In 1850’s Bernhard Riemann, a theology student in Germany turned back to his
first love mathematics and became a student of Friedrich Gauss. In his
Legendary 1859 paper titled "Ueber die Anzahl der Primzahlem unter eine
gegebener Gr${\ddot{o}}$sse" that is "On the Number of Primes Less Than a
Given Magnitude", Riemann conjectured his hypothesis. He presented the
unnatural connection between the position of zeros of the zeta function with
order of the error term in prime counting functions such as conjectured by
Gauss. The Prime Number Theorem( PNT) states that the number of primes less
than or equal to $n$ approaches $n/\log{n}$ as $n$ approaches infinity [9].

$\lim\limits_{n\to\infty}\frac{\pi(n)}{n/\log(n)}=\;1.$

Gauss shared his prime counting formula with Dirichlet who came up with a more
accurate counting function, the $Li$ function.

$Li(x):=\int\limits_{2}^{x}\frac{dx}{\log(x)}.$

Although $Li(x)$ and $\pi(x)$ asymptotically become more and more accurate as
$n$ increases but none of them is a true prime counting function. The Riemann
hypothesis is connected to the order of the error term. If the Riemann
hypothesis is true, there is an error of order $\sqrt{x}$ in counts of primes
as given by the $Li$ function [10]. That is for arbitrary $\epsilon>0$,

$\pi_{true}(x)=\int\limits_{2}^{x}\frac{dx}{\log(x)}+O(x^{\frac{1}{2}+\epsilon}).$

There are several other equivalent statements for RH made using it’s
connection with functions from number theory. For example, the
$Liouvillefunction$, $\lambda(n)$ which is defined by

$\lambda(n)\;=\;(-1)^{\omega(n)},\\\ $

where $\omega(n)$ is the number of prime factors, counted with multiplicity.
Thus, $\omega(6)=2$ and $\lambda(6)=\;1$,
$\omega(9)=2$ and $\lambda(9)=1$ and
$\omega(p)=1$ and $\lambda(p)=-1$ for all primes $p$.
In his doctoral thesis, Laudau (1899) showed that RH is equivalent to the
statement

$\lim\limits_{n\to\infty}\frac{\lambda(1)\;+\;\lambda(2)\;+\;\lambda(3)\;+\;.\;.\;.+\;\lambda(n)}{n^{\frac{1}{2}+\epsilon}}\;=\;0,$

for all $\epsilon>0$ [11]. The theory of RH is vast and connections to
seemingly unrelated fields have been found. One such example is the connection
to Random Matrix Theory as explored by Montgomery. In 1972, Hugh Montgomery
investigated the spacing between the zeros of the zeta function. This was just
after the spectral interpretation of RH provided by Hilbert-Pólya conjecture
[1]. Assuming RH to be true, Montgomery conjectured that the correlation for
the zeros of the zeta function is

$1-\frac{\sin^{2}(\pi x)}{(\pi x)^{2}}.$

This is the pair correlation Conjecture given by Montgomery. Later Freeman
Dyson identified that the above correlation function is the same as that for
pairs of eigenvalues of random Hermitian matrices with Gaussian measure used
by physicists. This gave birth to the Gaussian Unitary Ensemble( GUE)
conjecture[12].
Till now several equivalent statements and extensions of RH have been
formulated. New findings on association of zeta function with other concepts
keep redefining the hypothesis and encourage newer approach to the problem.
Despite of having acquired several pieces of information about the zeta
function, the problem has remained unsolved till date. However, famous
attempts to solve RH are worth mentioning.

In 1885, T. J Stieltjes claimed to have solved the Mertens’ conjecture which
implied the Riemann hypothesis. The results were proven wrong after his death
as Mertens’ conjecture was proven false by Andrew Odlyzko and Herman J.J. te
Riele in 1985 [13].

In the meanwhile, German-born American mathematician Hans Rademacher came up
with a proof of the falsehood of Riemann Hypothesis. He reported his proof in
Times magazine in 1945 but it was found erroneous. Legend has it that
Rademacher was terribly embarrassed by the ordeal and it was well known that
no one was to mention the words, Riemann Hypothesis in his presence[14].

In 1959 the famous mathematician John Nash presented a proof for the Riemann
hypothesis in a lecture at Columbia University. The attendees had high
expectations but were disappointed as the lecture didn’t make any sense. This
was later found to be due to his struggle with schizophrenia [15].

In 2004, Louis de Bourcia claimed to have the proof of Riemann hypothesis
which remained under skepticism for a long time. Later a counterexample from a
1998 paper by one of his student mathematicians Xian-Jin Li and Brian Conrey
showed his proof was false. Li, in July 2008 reported the proof of RH but it
was soon retracted after mathematicians pointed out a crucial flaw[16].

The most recent famous claim came in September, 2018 from British-Lebenese
mathematician Sir Micheal Atiyah. He used Todd functions and mentioned the
Fine Structure constant a fundamental physical constant in his proof by
contradiction [17]. The Fine Structure constant also called the Sommerfeld’s
constant is a dimensionless constant that equals $1/137$. A concept entirely
from physics. His proof was based on the works of leading 20th century
mathematicians, John von Neumann and Friedrich Hirzebruch. His proof remained
under skepticism and to be rigorously verified. Although few mathematicians
have declared it flawed [18].

##  2 Expansion of functional equation of zeta function into a series

The Dirichlet series form of the zeta function is only valid for $\Re(s)>1$.
Riemann gave the functional equation of Zeta function that continued it as a
meromorphic function with a simple pole at 1 [1].

$\zeta(s):=\frac{\pi^{s/2}}{\Gamma(s/2)}\Bigg{\\{}\frac{1}{s(s-1)}+\int\limits_{1}^{\infty}(x^{s/2-1}+x^{-s/2-1/2}).\left(\frac{\vartheta(x)-1}{2}\right)dx\Bigg{\\}}$

defined for Re(s) > 0, where

$\vartheta(x):=\sum\limits_{n=-\infty}^{\infty}e^{-n^{2}{\pi}x}=1+2\sum\limits_{n=1}^{\infty}e^{-n^{2}{\pi}x}$

is the Jacobi-theta function.

From analytic continuation of zeta function, one also gets the Riemann
reflection formula

$\xi(s)=\xi(1-s).$

Where $\xi(s)$ is an entire function and is defined as

$\xi(s)=\frac{s}{2}(s-1)\pi^{-\frac{s}{2}}\Gamma(\frac{s}{2})\zeta(s).$

In his extensive work on zeta function in 1859, Riemann also gave his famous
hypothesis about the position of non-trivial zeroes of the zeta function. The
Riemann hypothesis is essentially that all the values of $s$ in the critical
strip for which $\zeta(s)=0$, lie on the critical line ($\Re(s)=\frac{1}{2}$)
[1].
Let us define:

${\textstyle\Omega(s)=\frac{1}{s(s-1)}+\int\limits_{1}^{\infty}(x^{\frac{s}{2}-1}+x^{-\frac{s}{2}-\frac{1}{2}}).\left(\frac{\vartheta\left(x\right)-1}{2}\right)dx}$

defined in $(0<\Re(s)<1)$.
On expanding the integral term in $\Omega(s)$, by ‘integration by parts’
taking $x^{\frac{s}{2}-1}+x^{-\frac{s}{2}-\frac{1}{2}}$ as the first function
and the Jacobi-theta function term as the second,

$\displaystyle\textstyle\Omega(s)$
$\displaystyle=\frac{1}{s(s-1)}+\sum\limits_{n=1}^{\infty}\left(0-\frac{2}{s}-\frac{2}{1-s}\right)e^{-n^{2}\pi}$
$\displaystyle-\int\limits_{1}^{\infty}\left(\frac{2x^{\frac{s}{2}}}{s}+\frac{2x^{-\frac{s}{2}+\frac{1}{2}}}{1-s}\right)\sum\limits_{n=1}^{\infty}e^{-n^{2}\pi
x}(-n^{2}\pi)dx$

Repeating integration ‘by parts’ again and again, taking the Jacobi-theta
function part as the second function, one gets a series function:

$\displaystyle\Omega(s)=$
$\displaystyle\frac{1}{s(s-1)}+\sum\limits_{n=1}^{\infty}\left(-\frac{2}{s}-\frac{2}{1-s}\right)e^{-n^{2}\pi}$
$\displaystyle+\sum\limits_{n=1}^{\infty}\left(-\frac{1}{\frac{s}{2}(\frac{s}{2}+1)}-\frac{1}{(\frac{1-s}{2})(\frac{1-s}{2}+1)}\right)e^{-n^{2}\pi}.n^{2}\pi$
$\displaystyle+\sum\limits_{n=1}^{\infty}\left(-\frac{1}{\frac{s}{2}(\frac{s}{2}+1)}-\frac{1}{(\frac{1-s}{2})(\frac{1-s}{2}+1)}\right)e^{-n^{2}\pi}.(n^{2}\pi)^{2}+...$

We define:

$f(s)=-\frac{1}{s}+\sum\limits_{n=1}^{\infty}e^{-n^{2}\pi}\left(-\frac{1}{\left(\frac{s}{2}\right)^{\overline{1}}}-\frac{n^{2}\pi}{\left(\frac{s}{2}\right)^{\overline{2}}}-\frac{(n^{2}\pi)^{2}}{\left(\frac{s}{2}\right)^{\overline{3}}}-\frac{(n^{2}\pi)^{3}}{\left(\frac{s}{2}\right)^{\overline{4}}}-...\right)$

where
$(\frac{s}{2})^{\overline{n}}=\frac{s}{2}(\frac{s}{2}+1)(\frac{s}{2}+2)...(\frac{s}{2}+n-1)$
is the rising factorial power function.

So,

$\Omega(s)=f(s)+f(1-s),$

and

$\textstyle f(s)=-\left(\frac{1}{s}+J(s)\right).$

Where,

$J(s):=\sum\limits_{m=1}^{\infty}\frac{\sum\limits_{n=1}^{\infty}e^{-n^{2}\pi}\left(n^{2}\pi\right)^{m-1}}{\left(\frac{s}{2}\right)^{\overline{m}}}.$

##  3 Proof of the truth of RH

The zeros of $\Omega(s)$ are the non-trivial zeros of the zeta function.
$\Omega(s)$ can be equal to zero in either of the two cases:

  * •

if $f(s)$ and $f(1-s)$ are individually zero.

Or

  * •

if $f(s)=-f(1-s)$ and none equals zero.

Therefore, proving that $f(s)$ is not the additive inverse of $f(1-s)$
anywhere off the critical line is sufficient to proof the Riemann Hypothesis.

In the next section deductions about the allowed values of the series function
$J(s)$ are made by geometric analysis of simpler series functions and $J(s)$.
From the allowed values of $J(s)$, the existence of the first condition for
$\Omega(s)=0$ is obtained.

###  3.1  Geometry of series $J(s)$ on the complex plane

Plotting $J(s)$ term by term on the complex plane for a given s, leads us to
value of the series. In the present section we will find that $J(s)$ converges
as an anticlockwise inward spiral starting from origin, for values of s in the
lower half region of the critical strip.
Let’s define $\alpha\in(0,\frac{1}{2})$. Thus, $s=\frac{1}{2}\pm\alpha-i|t|$
represents any point in the lower region of the critical strip to the real
axis (of course, excluding the critical line $\Re(s)=\frac{1}{2}$). This is
our region of interest because if Riemann hypothesis is false then $\Omega(s)$
is zero in this region.
Now, for $s=\frac{1}{2}\pm\alpha-i|t|$, the argument of the $m^{th}$ term of
$J(s)$ is given by:

$\displaystyle\theta_{m,J}$
$\displaystyle=arg\left(\frac{1}{\frac{s}{2}}\right)^{\overline{m}}=-arg\left(\frac{s}{2}\right)^{\overline{m}}$
$\displaystyle=\arctan\left(\frac{\frac{|t|}{2}}{\frac{\Re(s)}{2}}\right)+\arctan\left(\frac{\frac{|t|}{2}}{\frac{\Re(s)}{2}+1}\right)+...+\arctan\left(\frac{\frac{|t|}{2}}{\frac{\Re(s)}{2}+m-1}\right)$

So, $\theta_{m,J}<m\theta$, where $\theta$ is $\theta_{1,J}$ (argument of the
first term OA).
$0<\theta_{m+1,J}-\theta_{m,J}=\arctan\left(\frac{\frac{|t|}{2}}{\frac{\Re(s)}{2}+m}\right)<\frac{\pi}{2}$
i.e., is positive and acute $\forall\;\;m$ . Thus, as $m$ increases, the
difference between the arguments of consecutive terms of $J(s)$ decreases.

Since, $J(s)$ is a component function in the functional equation of zeta
function, it is also analytic throughout the critical strip. A proof of its
convergence is provided in the appendix. Form its convergence we have,

$\frac{|T_{m+1,J}|}{|T_{m,J}|}=\frac{\sum\limits_{n=1}^{\infty}e^{-n^{2}\pi}\left(n^{2}\pi\right)^{m}}{|\frac{s}{2}+m|\sum\limits_{n=1}^{\infty}e^{-n^{2}\pi}\left(n^{2}\pi\right)^{m-1}}<1\;\;\;\forall\;\;m$

(a)
$\frac{|T_{m+1}|}{|T_{m}|}=1,\;\theta_{m}=m\theta\;and\;\frac{2\pi}{\theta}\in\;\mathbb{N}^{+}$

(b)
$\frac{|T_{m+1}|}{|T_{m}|}=1,\;\theta_{m}=m\theta\;and\;\frac{2\pi}{\theta}\notin\;\mathbb{N}^{+}$

Figure 1: Different cases of series $U(s)$

Now, let’s try to plot $J(s)$, by plotting some simpler series.
Take a series,

$U(s)=\sum\limits_{m=1}^{\infty}T_{m}(s)=\sum\limits_{m=1}^{\infty}{|T_{m}(s)|}e^{i\theta_{m}}$

such that

$|T_{m+1}(s)|=|T_{m}(s)|=T$

Therefore, we have

$OA=T_{1}(s)=|T_{1}(s)|e^{i\theta}=Te^{i\theta}$
$AB=|T_{2}(s)|e^{2i\theta}=Te^{2i\theta}$
$BC=T_{3}(s)=|T_{3}(s)|e^{3i\theta}=Te^{3i\theta}$

Fig. 1 shows U(s) series with $\frac{|T_{m+1}|}{|T_{m}|}=1$,
$\theta_{m}=m\theta$, i.e. $\frac{T_{m+1}(s)}{T_{m}(s)}=e^{i\theta}$. Here,
either of the following two cases will occur:

  * •

A uniform polygon of sides $\frac{2\pi}{\theta}$ : the whole series plots as a
single polygon of $\frac{2\pi}{\theta}$ sides (fig. a).

  * •

if $\theta$ does not divide $2\pi$ i.e., the circle into $n$ equal sectors and
we have a geometry as shown (fig. b) with consecutive cycles going out of
phase and there exists $\beta:0<\beta<\theta$, and $\beta+(n-1)\theta=2\pi$,
which is the angle formed by the last sector at the center of the polygon on
completion of $2\pi$.

Figure 2: $\mu(s):\frac{|T_{m+1}|}{|T_{m}|}<1$, $\theta_{m}=m\theta$

But, in none of the cases above, does the series converge. Let’s define a new
series function $\mu(s)=\sum\limits_{m=1}^{\infty}T_{m,\mu}(s)$, such that
$\frac{|T_{m+1,\mu}|}{|T_{m,\mu}|}<1$ and,
$\theta_{m,\mu}=\theta_{m}=m\theta$, for $\mu(s)$, where $\theta$ is acute and
positive.

As shown in fig. 2, we get an inwardly spiraling geometry going anticlockwise
from the origin. The value of $\mu(s)$ always lie inside the circle passing
through points O, A, B, C, D, …. Fig. 2 shows that any series with the
attributes of $\mu(s)$ necessarily lies in the circle with centre $P_{1}$ and
radius $AP_{1}=BP_{1}$. Take the first term $OA$ as in the $U(s)$ series.
Since,

$|T_{2,\mu}(s)|<|T_{1,\mu}(s)|$

thus B’ lies between points A and B on chord $AB$. Thus, the point B’ lies
inside circle with centre $P_{1}$. Now, draw a circle with AB’ as its chord
such that it projects the same angle at it’s center $(P_{2})$ as chord AB
projects at the center ($P_{1}$) of outer circle. This new circle with centre
$P_{2}$ and radius $B^{\prime}P_{2}=AP_{2}\;<AP_{1}$ is drawn by moving center
$P_{1}$ along $AP_{1}$ to $P_{2}$ such that it’s tangential at point A to the
outer circle and $B^{\prime}P_{2}$ is parallel to $BP_{1}$. Using the
similarity of $\bigtriangleup AP_{1}B$ with $\bigtriangleup AP_{2}B^{\prime}$,
the length $AP_{2}$ follows
$\frac{AP_{2}}{AP_{1}}=\frac{AB^{\prime}}{AB}\;<\;1$. Since, circle with
centre $P_{2}$ touches the outer circle only tangentially at point A, it’s
area always lie inside that of the outer circle. Same treatment can be done
for the third term $T_{3,\mu}(s)$ represented by $B^{\prime}C^{\prime}$ which
is less than $BC$ by generating new circle tangential to the previous one at
point B’. Thus circle with center $P_{3}$ and radius
$C^{\prime}P_{3}=B^{\prime}P_{3}$ ($<\;B^{\prime}P_{2}$ since
$|T_{3,\mu}(s)|<|T_{2,\mu}(s)|$) is generated by moving poimt $P_{2}$ along
$B^{\prime}P_{2}$ to point $P_{3}$. Again, area of the circle with center
$P_{3}$ and radius $C^{\prime}P_{3}$ lies inside the circle with center
$P_{2}$. If $C_{n}$ represents the area under the $nth$ circle with center
$P_{n}$ . We have,

$C_{n}\subset\;C_{n-1}\subset\;...C_{2}\subset\;C_{1}$

. The value of the series $\mu(s)$ lies in the $nth$ circle as $n$ approaches
infinity. The area of circles with center $P_{2},P_{3},..,P_{n}$ lie inside
area $C_{1}$. This shows that the value of all series with attributes of
$\mu(s)$ must lie inside the circle with center $P_{1}$ circumscribing the
polygon formed for the case $T_{m+1}(s)=T_{m}(s)$ as represented by $U(s)$
above. The red colored spiral is formed by the arcs
$\stackrel{{\scriptstyle\mbox{\large$\frown$}}}{{OA}}$,
$\stackrel{{\scriptstyle\mbox{\large$\frown$}}}{{AB^{\prime}}}$,
$\stackrel{{\scriptstyle\mbox{\large$\frown$}}}{{B^{\prime}C^{\prime}}}$,
$\stackrel{{\scriptstyle\mbox{\large$\frown$}}}{{C^{\prime}D^{\prime}}}$ … and
converges at the same point as the series does.
The rate of convergence depends on the ratio of the consecutive terms. The
lower the ratio of consecutive terms $\frac{|T_{m+1}|}{|T_{m}|}$, the higher
the rate of convergence of the series $\mu(s)$. Thus, on choosing the first
term in $\mu(s)$ same as for $U(s)$, it converges at a point inside the circle
circumscribing the geometry formed for the $\frac{|T_{m+1}|}{|T_{m}|}=1$ cases
above. The grey to black circular regions represent necessary region in which
the value converges after first few term.

Now, coming back to our original series

$J(s)=\sum\limits_{m=1}^{\infty}T_{m}(s)=\sum\limits_{m=1}^{\infty}|T_{m}(s)|e^{i\theta_{m,J}}$

.
where,
$\theta_{m,J}=arg\left(\frac{1}{\frac{s}{2}}\right)^{\overline{m}}=-arg\left(\frac{s}{2}\right)^{\overline{m}}$

Figure 3: $J(s):\frac{|T_{m+1}|}{|T_{m}|}<1$, $\theta_{m,J}<m\theta$ and
$\mu(s)$

We exploit the fact that for both $\mu(s)$ and $J(s)$ we have,
$\frac{|T_{m+1}|}{|T_{m}|}<1$ and draw $\mu(s)$ with terms having magnitude as
that for $J(s)$ but $\theta_{m}=m\theta$. This corresponds to the geometry
OAB’C’D’E’… in fig. $(3)$ below. Now, we just have to rotate each term
(starting from the second term) relative to the preceding one, clockwise. This
rotation or unfastening is fixed for each term. So, from the last case we fix
OA and rotate AB’ by ($\theta_{2,J}-\theta$) clockwise about A, to new
position AB”, so that the acute angle between AB” and line colinear with OA is
now $\theta_{2,J}$. Consecutively, rotate C’B” clockwise about B” by
$\theta_{3,J}-\theta$ to new position B”C” and D’C” about C” by
$\theta_{4,J}-\theta$ and so on. As shown in fig. 3, $J(s)$ is represented by
OAB”C”D”E”… which is unfastened form of the spiral OAB’C’D’E’…. Thus the rate
of convergence is reduced. The unfastening is not done uniformly to all line
segments representing series terms. Rather, the unfastening monotonously
decreases over terms, from $arctan\left(\frac{|t|/2}{\Re(s)/2+1}\right)$ at m=
2 (since OA is fixed) to 0 as $m\rightarrow\infty$.

In fig. 3, the spiral formed in fig. 2 is unfastened so that, argument of
$mth$ term is $\theta_{m,J}(<m\theta)$ and depending on the difference
$\theta_{m,J}-m\theta,$ we get the extent of unfastening( clockwise rotation)
of each term of the spiral. Nonetheless this unfastened spiral also has a
monotonous change in curvature from its centre to it’s farthest point at the
origin, $O$. This is true for all series with $|T_{m+1}|<|T_{m}|$ and
$\theta_{m+1}>\theta_{m}$. The same applies to $J(s)$.

Using the above deductions about $f(s)$ and $J(s)$, we look at the conditions
for $\Omega(s)$ to be zero for any $\alpha\in\left(0,\frac{1}{2}\right)$.

###  3.2 Conditions for the Falsification of the Riemann Hypothesis

####  3.2.1 Case-1: RH is false if $f(s)=f(1-s)=0$, for
$s=\frac{1}{2}\pm\alpha-i|t|$

The condition $f(s)=0$ implies $J(s)=-\frac{1}{s}$. For any monotonously
converging spiral, the distance from the centre of the spiral to a point on
the spiral changes monotonously along the spiral.

Let point Q represent the value $-\frac{1}{s}$ on the complex plane.

$OA=|\;T_{1,J}(s)\;|=\;\frac{2\sum\limits_{n=1}^{\infty}e^{-n^{2}\pi}}{|s|}=\frac{\frac{\pi^{\frac{1}{4}}}{\Gamma(\frac{3}{4})}-1}{|s|}=\frac{0.0864348}{|s|}$

Figure 4: Value of $J(s)$ lies col-linearly with chord OA if J(s)
=-$\frac{1}{s}$

If (see fig. 6) $J(s)=-\frac{1}{s}$ then, the the centre of spiral
$OAB"C"D"E"...$ coincides point Q. Since, $J(s)$ is a monotonously converging
spiral starting anticlockwise from first term $OA$ as discussed in previous
section. Thus,

$...<QD"<QC"<QB"<QA<QO$

as distance from Q to points on the spiral change monotonously as you move
along the spiral. But, one finds that distance $QO=\frac{1}{|s|}$ and
$QA=OA+OQ=\frac{1.0864348}{|s|}$. Thus,

$QO<QA.$

Thus by contradiction $\frac{1}{s}+J(s)\neq 0$,

$\Rightarrow f(s)\neq 0.$

Also, since $(f(s))^{*}=f(s^{*})$ thus, $f(s)\neq O$ for any $s$, off the
critical line, $\Re(s)=\frac{1}{2}$.

Thus, no such case exists to falsify RH.

####  3.2.2 Case-2: RH is false if $f(s)=-f(1-s)$

For this case we first show that the magnitude of our complex series function,
$|f(s)|$ always decreases as the real part of $s$ changes from zero to one
keeping the imaginary part of $s$ constant. Since, showing that for $\Re(s)$
less than $\frac{1}{2}$, the magnitude $|f(1-s*)|$ is always less than
$|f(s)|$, is sufficient to show that $f(s)\neq-f(1-s)$ anywhere in the
critical strip off the critical line.

We know that,

$f(s)=-\left(\frac{1}{s}+J(s)\right)=\sum\limits_{m=1}^{\infty}T_{m}(s)=\sum\limits_{m=1}^{\infty}\frac{C_{m}}{\left(\frac{s}{2}\right)^{\overline{m}}}$

.

Figure 5: Magnitude and argument (not necessarily less than $\pi/2$ as
depicted) of the $nth$ series term $T_{n}(s_{2})$ of $-f(s_{2})$ is always
less than the corresponding term $T_{n}(s_{1})$ of $-f(s_{1})$ if
$\Re(s_{1})<\Re(s_{2})$ and $\Im(s_{1})=\Im(s_{2})$.

Note that $-f(s)$ can be plotted by plotting $J(s)$ but with the first term’s
coefficient $C_{1}$ taken as
$\left(\frac{1}{2}+\sum\limits_{n=1}^{\infty}e^{-n^{2}\pi}\right)$ instead of
$\sum\limits_{n=1}^{\infty}e^{-n^{2}\pi}$. We, know already that from the
center of the spiral the radial distance monotonously increases with each term
such that it’s maximum after counting all the terms i.e. from the spiral
centre to the origin. The point Q is the center of the spiral. It has value
$-f(s_{1})$ which is reached by laying the terms $OA=T_{1}(s_{1})$,
$AB=T_{2}(s_{1})$, $BC=T_{3}(s_{1})$ and so on in the complex plane. We also
draw concentric circles centered at the center of the spiral and having radius
equal to the radial distance from the center after the addition of each term
(see fig. 6). Further, we decorate the corresponding terms
$OA^{\prime}=T_{1}(s_{2})$, $AB^{\prime}=T_{2}(s_{2})$,
$BC^{\prime}=T_{3}(s_{2})$ and so on of $-f(s_{2})$ at points $O,A,B,C....$
respectively. As a term with relatively smaller magnitude and argument
replaces a term in the spiral, shorter radial distance covered by the spiral
in that term. From the figure one can see that to draw the spiral geometry
$-f(s_{2})$ about point Q one has to move it’s terms across the light red
annular regions to move close to Q. So, the term OA’ is moved to coincide
point A’ and point A across the last light red circular band. Similarly,
second term AB’ has to be moved across the second annular region (in light
red) to coincide points B and B’. Repeating the same for rest of the terms to
form a continuous spiral geometry $-f(s_{2})$. Thus, to draw the spiral
geometry for $-f(s_{2})$ about Q all the light red coloured annular regions
are eliminated. This implies that the radial distance covered by each term of
$-f(s_{2})$ (in yellow) about Q is less than that for $-f(S_{1})$. Thus, we
conclude
$|\sum\limits_{i=n_{1}}^{n_{2}}T_{n}(s_{1})|>|\sum\limits_{i=n_{1}}^{n_{2}}T_{n}(s_{2})|$
and therefore, $|f(s_{1})|>|f(s_{2})|$.

Figure 6: Spiral geometry O-A-B-C-..-Q is term by term plot of $-f(s_{1})$ and
the corresponding terms OA’, AB’, BC’,… of $-f(s_{2})$ (in yellow) are placed
from first term till the last term. The light red annular regions correspond
(not equal to) the term wise difference in radial distance about Q, overall
resulting in $|f(s_{1})|>|f(s_{2})|$.Given, $\Re(s_{1})<\Re(s_{2})$ and
$\Im(s_{1})=\Im(s_{2})$.

If for two variables $s_{1}$ and $s_{2}$, where
$s_{1}=\frac{1}{2}-\alpha-i|t|$ and
$s_{2}=1-s_{1}^{*}=\frac{1}{2}+\alpha-i|t|$ where $\alpha\in(0,\frac{1}{2})$
and $t\in\mathbb{R}$ then as shown in the fig. 5,

$\bigg{|}T_{n}(s_{1})\bigg{|}>\bigg{|}T_{n}(s_{2})\bigg{|}$

.

$\arg\left(T_{n}(s_{1})\right)>\arg\left(T_{n}(s_{2})\right)$

.

Thus,

$|f(\frac{1}{2}-\alpha-i|t|)|>|f(\frac{1}{2}+\alpha-i|t|)|=|f(\frac{1}{2}+\alpha+i|t|)|$

.

This means nowhere in the critical strip can $f(s)$ be equal to $-f(1-s)$ off
the critical line. Here we can see that the non-trivial zeros correspond to
purely imaginary values of $f(\frac{1}{2}+it)$ such that
$f(s)+f(1-s)=f(s)+f(s^{*})=0$.

Since, both the cases for falsification of the hypothesis are non-existent.
Thus, the Riemann hypothesis is true.

## APPENDIX: On convergence of $J(s)$

In this section it’s shown that the series function $J(s)$ converges in the
critical strip region for $|\frac{s}{2}|\;\geq\;1$. This is essentially the
region of interest as in the region under the disk $|\frac{s}{2}|\;<\ 1$, it’s
known that no zero exists.
We use the Taylor series expansion of a new series $\Lambda(s)$. Using the
integration test for convergence of $\Lambda(s)$ in the critical region and
comparing the magnitude of coefficients of corresponding terms in $J(s)$, it’s
shown that $|J(s)|<|\Lambda(\frac{s}{2})|$. Thus, $\;J(s)\;$ converges.

Let $\Lambda(x)=\sum\limits_{n=1}^{\infty}e^{-n^{2}\pi x}$
Then, Taylor series expansion of $\Lambda(x)$ about point $x=1$ is

$\Lambda(x)=\sum\limits_{n=1}^{\infty}e^{-n^{2}\pi}+\left(1-x\right)\sum\limits_{n=1}^{\infty}e^{-n^{2}\pi}\left(n^{2}\pi\right)+\frac{\left(1-x\right)^{2}}{2}\sum\limits_{n=1}^{\infty}e^{-n^{2}\pi}\left(n^{2}\pi\right)^{2}+...$

Also notice that all terms are positive for $0<x<1$.

Now, using the integration test for convergence of $\Lambda(x)$ [19],

$a_{1}\geq\lim_{n\to\infty}\left(S_{n}\left(x\right)-I_{n}\left(x\right)\right)\geq
0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;$ (1)
here, $S_{n}\left(x\right)=\sum\limits_{j=1}^{n}e^{-j^{2}\pi x}$
and $a_{1}=e^{-\pi x}$

$\lim_{n\to\infty}I_{n}\left(x\right)=\int\limits_{1}^{\infty}e^{-k^{2}\pi
x}dk=\frac{1}{2\sqrt{x}}erfc(1)$

Where, erfc(x) is the complimentary error function defined as,
$erfc(x)=\frac{2}{\sqrt{\pi}}\int\limits_{x}^{\infty}e^{-t^{2}}dt$
thus, from (1) we get

$e^{-\pi
x}+\frac{1}{2\sqrt{x}}erfc(1)\geq\Lambda(x)\geq\frac{1}{2\sqrt{x}}erfc(1)$

Thus, $\Lambda(s)$ is convergent only for $x>0$

on expanding $\Lambda(x)$ from a real function to the complex plane

$\Lambda(s)=\sum\limits_{n=1}^{\infty}e^{-n^{2}\pi}+\left(1-s\right)\sum\limits_{n=1}^{\infty}e^{-n^{2}\pi}\left(n^{2}\pi\right)+\frac{\left(1-s\right)^{2}}{2}\sum\limits_{n=1}^{\infty}e^{-n^{2}\pi}\left(n^{2}\pi\right)^{2}+...$

Where $s=\sigma+it\;\;\;\;\;$, $\sigma\;\in(0,1)$ and $t\in\mathbb{R}$

but,

$|\sum\limits_{n=1}^{\infty}e^{-n^{2}\pi(\sigma+it)}|\leq|\sum\limits_{n=1}^{\infty}e^{-n^{2}\pi(\sigma+it)}|=\sum\limits_{n=1}^{\infty}e^{-n^{2}\pi\sigma}\leq
e^{-\pi\sigma}+\frac{1}{2\sqrt{\sigma}}erfc(1)$
thus,
$\Lambda(s)=\sum\limits_{n=1}^{\infty}e^{-n^{2}\pi s}$
is convergent for $\sigma>0$

Coming back to $J(s)$,

$J(s):=\sum\limits_{m=1}^{\infty}\frac{\sum\limits_{n=1}^{\infty}e^{-n^{2}\pi}\left(n^{2}\pi\right)^{m-1}}{\left(\frac{s}{2}\right)^{\overline{m}}}$

Since, $\;\forall\;s\;$ in the critical strip :$|s|\geq 1$ and
$k\in\mathbb{N}^{+}$, one can say that

$\;\;\;\;\;|s||s+k|>k$

$\Rightarrow\frac{|s|}{k}>\frac{1}{|s+k|}$

$\Rightarrow\frac{|\left(\frac{s}{2}\right)^{k}|}{k!}\;=\;\frac{|\frac{s}{2}|^{k}}{k!}\;>\;\frac{1}{\big{|}\left(\frac{s}{2}\right)\big{|}\big{|}\left(\frac{s}{2}+1\right)\big{|}...\big{|}\left(\frac{s}{2}+k-1\right)\big{|}\big{|}\left(\frac{s}{2}+k\right)\big{|}}\;=\;\frac{1}{\big{|}\left(\frac{s}{2}\right)\left(\frac{s}{2}+1\right)...\left(\frac{s}{2}+k-1\right)\left(\frac{s}{2}+k\right)\big{|}}$

$\Rightarrow|\Lambda\left(\frac{s}{2}\right)|>|J(s)|\;\;\;\forall\;s$ in the
critical strip outside the disk $|\frac{s}{2}|\leq 1$ centered at origin.

Thus,

$J(s):=\sum\limits_{m=1}^{\infty}\frac{\sum\limits_{n=1}^{\infty}e^{-n^{2}\pi}\left(n^{2}\pi\right)^{m-1}}{\left(\frac{s}{2}\right)^{\overline{m}}}\;\;$
converges in our region of interest.

Hence, proved!

Therefore, the ratio of it’s consecutive terms,
$\frac{|T_{m+1}|}{|T_{m}|}=\frac{\sum\limits_{n=1}^{\infty}e^{-n^{2}\pi}\left(n^{2}\pi\right)^{m}}{|\frac{s}{2}+m|\sum\limits_{n=1}^{\infty}e^{-n^{2}\pi}\left(n^{2}\pi\right)^{m-1}}<1$.

Another way to conclude the convergence of $J(s)$ in the critical strip is to
use the fact that the Functional equation of zeta function converges in the
critical strip and therefore all its component functions must.

## References

  * [1] B Rooney P Borwein.  The Riemann Hypothesis: A Resource for the Afficionado and Virtuoso Alike, chapter Empirical Evidence, pages 40–43.  Springer, june 2007.
  * [2] John Derbyshire  Prime obsession : Bernhard Riemann and the greatest unsolved problem in mathematics, chapter Number Theory Meets Quantum Mechanics, pages 280–295.  Joseph Henry Press, Washington, D.C., august 2003.
  * [3] B Rooney P Borwein.  The Riemann Hypothesis: A Resource for the Afficionado and Virtuoso Alike, chapter Algorithms for calculating $\zeta(s)$, pages 28–36.  Springer, june 2007.
  * [4] B Rooney P Borwein.  The Riemann Hypothesis: A Resource for the Afficionado and Virtuoso Alike, chapter Analytic Preliminaries, pages 9–11.  Springer, june 2007.
  * [5] John Derbyshire.  Prime obsession : Bernhard Riemann and the greatest unsolved problem in mathematics, chapter The Golden Key, and an Improved Prime Number Theorem, pages 102–104.  Joseph Henry Press Washington, D.C., august 2003.
  * [6] B Rooney P Borwein.  The Riemann Hypothesis: A Resource for the Afficionado and Virtuoso Alike, chapter 1, pages 7–8.  Springer, june 2007.
  * [7] B Rooney P Borwein.  The Riemann Hypothesis: A Resource for the Afficionado and Virtuoso Alike, chapter Algorithms for calculating $\zeta(s)$ pages 29–30.  Springer, june 2007.
  * [8] P Sebah, X Gourdon  "Introduction to Bernaulli’s Numbers",  https://math.ucr.edu/~res/math153/s12/bernoulli-numbers.pdf.  june 2002.
  * [9] John Derbyshire.  "Prime obsession : Bernhard Riemann and the greatest unsolved problem in mathematics", chapter The Prime Number Theorem, pages 32–47.  Joseph Henry Press, Washington, D.C., august 2003.
  * [10] B Rooney P Borwein.  The Riemann Hypothesis: A Resource for the Afficionado and Virtuoso Alike, chapter 1, pages 17–18.  Springer, june 2007.
  * [11] B Rooney P Borwein.  The Riemann Hypothesis: A Resource for the Afficionado and Virtuoso Alike, chapter 1, pages 5–6.  Springer, june 2007.
  * [12] J. Brian Conrey.  The Riemann Hypothesis.  Notices of the AMS, 50(3), 2003.
  * [13] Juan Miguel Marín.  Marvelous cancellations: T.J. Stieltjes’ letters concerning the zeta function.  Arxiv e-print, 01 2020.
  * [14] Daniel Nahum Rockmore.  Stacking The Riemann Hypothesis.  Random House, 2011.
  * [15] Karl Sabbagh.  The Riemann Hypothesis : The Greatest Unsolved Problem in Mathematics, page 105.  Random House, 2011.
  * [16] Columbia University Dept. of mathematics.  Proof of the riemann hypothesis?  https://www.math.columbia.edu/~woit/wordpress/?p=707.  2008\.
  * [17] Michael Atiyah.  The Riemann Hypothesis.  2018\.
  * [18] ASM blogs.  On michael atiyah and the riemann hypothesis.  https://blogs.ams.org/blogonmathblogs/2018/10/01/on-michael-atiyah-and-the-riemann-hypothesis/.  Last accessed: 25 Feb 2021.
  * [19] Thomas John I’Anson Bromwich.  An introduction to the theory of infinite series.  London Macmillan, 1908. pages 29–30. ISBN : 978-935128-594-6.


# A blueprint for precise and fault-tolerant analog neural networks

Cansu Demirkiran<EMAIL_ADDRESS>Boston University  Lakshmi Nair  Lightmatter
Darius Bunandar  Lightmatter  Ajay Joshi  Boston University

Abstract

Analog computing has reemerged as a promising avenue for accelerating deep
neural networks (DNNs) due to its potential to overcome the energy efficiency
and scalability challenges posed by traditional digital architectures.
However, achieving high precision and DNN accuracy using such technologies is
challenging, as high-precision data converters are costly and impractical. In
this paper, we address this challenge by using the residue number system
(RNS). RNS allows composing high-precision operations from multiple low-
precision operations, thereby eliminating the information loss caused by the
limited precision of the data converters. Our study demonstrates that analog
accelerators utilizing the RNS-based approach can achieve ${\geq}99\%$ of FP32
accuracy for state-of-the-art DNN inference using data converters with only
$6$-bit precision whereas a conventional analog core requires more than
$8$-bit precision to achieve the same accuracy in the same DNNs. The reduced
precision requirements imply that using RNS can reduce the energy consumption
of analog accelerators by several orders of magnitude while maintaining the
same throughput and precision. Our study extends this approach to DNN
training, where we can efficiently train DNNs using $7$-bit integer arithmetic
while achieving accuracy comparable to FP32 precision. Lastly, we present a
fault-tolerant dataflow using redundant RNS error-correcting codes to protect
the computation against noise and errors inherent within an analog
accelerator.

## Introduction

Deep Neural Networks (DNNs) are widely employed across various applications
today. Unfortunately, their compute, memory, and communication demands are
continuously on the rise. The slow-down in CMOS technology scaling, along with
these increasing demands has led analog DNN accelerators to gain significant
research interest. Recent research has been focused on using various analog
technologies such as photonic cores 1, 2, 7, 3, 4, 5, 6, resistive arrays 8,
9, 10, 11, 12, switched capacitor arrays 13, 14, Phase Change Materials (PCM)
15, Spin-Transfer Torque (STT)-RAM 16, 17, etc., to enable highly parallel,
fast, and efficient matrix-vector multiplications (MVMs) in the analog domain.
These MVMs are fundamental components used to build larger general matrix-
matrix multiplication (GEMM) operations, which make up more than $90\%$ of the
operations in DNN inference and training 18.

The success of this approach, however, is constrained by the limited precision
of the digital-to-analog and analog-to-digital data converters (i.e., DACs and
ADCs). In an analog accelerator, the data is converted between analog and
digital domains using data converters before and after every analog operation.
Typically, a complete GEMM operation cannot be performed at once in the analog
domain due to the fixed size of the analog core. Instead, the GEMM operation
is tiled into smaller MVM operations. As a result, each MVM operation produces
a partial output that must be accumulated with other partial outputs to obtain
the final GEMM result. Concretely, an MVM operation consists of parallel dot
products between $b_{w}$-bit signed weight vectors and $b_{in}$-bit signed
input vectors—each with $h$ elements—resulting in a partial output containing
$b_{\operatorname{\text{out}}}$ bits of information, where
$b_{\operatorname{\text{out}}}=b_{\operatorname{\text{in}}}+b_{w}+\log_{2}(h)-1$.
An ADC with a precision greater than $b_{\operatorname{\text{out}}}$ (i.e.,
$b_{\operatorname{\text{ADC}}}\geq b_{\operatorname{\text{out}}}$) is required
to ensure no loss of information when capturing these partial outputs.
Unfortunately, the energy consumption of ADCs increases exponentially with bit
precision (often referred to as effective number of bits (ENOB)). This
increase is roughly 4$\times$ for each additional output bit 19.

As a result, energy-efficient analog accelerator designs typically employ ADCs
with lower precision than $b_{\operatorname{\text{out}}}$ and only capture the
$b_{\operatorname{\text{ADC}}}$ most significant bits (MSBs) from the
$b_{\operatorname{\text{out}}}$ bits of each partial output 20. Reading only
MSBs causes information loss in each partial output leading to accuracy
degradation in DNNs, as pointed out by Rekhi et al. 20. This degradation is
most pronounced in large DNNs and large datasets. Fig. 1 shows the impact of
this approach on DNN accuracy in two tasks: (1) a two-layer convolutional
neural network (CNN) for classifying the MNIST dataset 21: a simple task with
only 10 classes, and (2) the ResNet50 CNN 22 for classifying the ImageNet
dataset 23: a more challenging task with 1000 classes. As the vector size $h$
increases, higher precision is needed at the output to maintain the accuracy
in both DNNs. Moreover, ResNet50 experiences accuracy degradation at smaller
values of $h$ compared to the two-layer CNN. While using a higher precision
ADC can help recover from this accuracy degradation, it significantly reduces
the energy efficiency of the analog hardware. Essentially, to efficiently
execute large DNNs using analog accelerators, it is crucial to find a better
way to achieve high accuracy than simply increasing the bit precision of the
data converters.

In this work, we present a universal residue number system (RNS)-based
framework to overcome the abovementioned challenge in analog DNN inference as
well as DNN training. RNS represents high-precision values using multiple low-
precision integer residues for a selected set of moduli. As such, RNS enables
high-precision arithmetic without any information loss on the partial
products, even when using low-precision DACs and ADCs. Utilization of RNS
leads to a significant reduction in the data converter energy consumption,
which is the primary contributor to energy usage in analog accelerators. This
reduction can reach up to six orders of magnitude compared to a conventional
fixed-point analog core with the same output bit precision.

Figure 1: Inference accuracy versus vector size for varying data bit-width in
a conventional analog core. a Inference accuracy for a two-layer CNN
classifying handwritten digits from the MNIST dataset. b Inference accuracy
for ResNet50 classifying images from the ImageNet dataset evaluated in an
analog core with varying precision $b$ and vector sizes $h$. For both a and b,
$b$-bit precision means
$b=b_{\operatorname{\text{DAC}}}=b_{\operatorname{\text{ADC}}}<b_{\text{out}}$
where $b$ varies between $2$ and $8$.

Our study shows that the RNS-based approach enables $\geq 99\%$ FP-32
inference accuracy by using only 6-bit data converters for state-of-the-art
MLPerf (Inference: Datacenters) benchmarks 24 and Large Language Models
(LLMs). We also demonstrate the applicability of this approach in training and
fine-tuning state-of-the-art DNNs using low-precision analog hardware. The RNS
approach, however, is susceptible to noise as small errors in the residues
scale up during output reconstruction, leading to larger errors in the
standard representation. To address this issue, we incorporate the Redundant
RNS (RRNS) error-correcting code 25, 26, 27 to introduce fault-tolerance
capabilities into the dataflow.

As RNS is closed under multiplication and addition, no significant changes are
required in the design of the analog core or in how GEMM operations are
performed. Unlike a conventional analog core design, performing RNS operations
necessitates an analog modulo operation. This operation can be implemented by
using ring oscillators 28 in an analog electrical core or by using optical
phase shifters in an analog optical core. Our proposed framework, however,
remains agnostic to the underlying technology. Importantly, arbitrary fixed-
point precision can be achieved by combining the positional number system
(PNS) and RNS in analog hardware. Overall, our presented RNS-based methodology
offers a solution combining high accuracy, high energy efficiency, and fault
tolerance in analog DNN inference and training.

## Results

### DNN Inference and Training Using RNS

The RNS represents an integer as a set of smaller (integer) residues. These
residues are calculated by performing a modulo operation on the said integer
using a selected set of $n$ _co-prime_ moduli. Let $A$ be an integer. $A$ can
be represented in the RNS with $n$ residues as $\\{a_{1},\dots,a_{n}\\}$ for a
set of co-prime moduli $\mathcal{M}=\\{m_{1},\dots,m_{n}\\}$ where
$a_{i}=|A|_{m_{i}}\equiv A\mod m_{i}$ for $i\in\\{1\dots n\\}$. $A$ can be
uniquely reconstructed using the Chinese Remainder Theorem (CRT):

$A=\bigg{|}\sum_{i=1}^{n}a_{i}M_{i}T_{i}\bigg{|}_{M},$ (1)

if $A$ is within the range $[0,M)$ where $M=\prod_{i}m_{i}$. Here,
$M_{i}=M/m_{i}$ and $T_{i}$ is the multiplicative inverse of $M_{i}$, i.e.,
$\left|M_{i}T_{i}\right|_{m_{i}}\equiv 1$. Hereinafter, we refer to the
integer $A$ as the _standard representation_ , while we refer to the set of
integers $\\{a_{1},\dots,a_{n}\\}$ simply as the residues.

A DNN consists of a sequence of $L$ layers. During inference, where the DNN is
previously trained and its parameters are fixed, only a forward pass is
performed. Generically, the input $X$ to $({\ell}+1)$-th layer of a DNN during
the forward pass is the output generated by the previous $\ell$-th layer:

$X^{(\ell+1)}=f^{(\ell)}\big{(}W^{(\ell)}X^{(\ell)}\big{)},$ (2)

where $W^{(\ell)}X^{(\ell)}=O^{(\ell)}$ is a GEMM operation and $f(\cdot)$ is
an element-wise nonlinear function.

DNN training requires both forward and backward passes as well as weight
updates. The forward pass in the training is performed the same way as in Eq.
(2). After the forward pass, a loss value $\mathcal{L}$ is calculated using
the output of the last layer and the ground truth. The gradients of the DNN
activations and parameters with respect to $\mathcal{L}$ for each layer are
calculated by performing a backward pass after each forward pass:

$\frac{\partial{\mathcal{L}}}{\partial{X^{(\ell)}}}={{W^{(\ell)}}^{T}}\frac{\partial{\mathcal{L}}}{\partial{O^{(\ell)}}},$
(3)
$\frac{\partial{\mathcal{L}}}{\partial{W^{(\ell)}}}=\frac{\partial{\mathcal{L}}}{\partial{O^{(\ell)}}}{X^{(\ell)}}^{T}.$
(4)

Using these gradients $\Delta
W^{(\ell)}={\frac{\partial{\mathcal{L}}}{\partial{W^{(\ell)}}}}$, the DNN
parameters are updated in each iteration $i$:

$W^{(\ell)}_{i+1}=W^{(\ell)}_{i}-\eta\Delta W^{(\ell)}_{i}$ (5)

with a step size $\eta$ for a simple stochastic gradient descent (SGD)
optimization algorithm.

Essentially, for each layer, one GEMM operation is performed in the forward
pass and two GEMM operations are performed in the backward pass. Because RNS
is closed under addition and multiplication operations, GEMM operations can be
performed in the RNS space. Using the RNS, Eq. (2) can be rewritten as:

$X^{(\ell+1)}=f^{(\ell)}\Bigg{(}\text{CRT}\bigg{(}\Big{|}\big{|}W^{(\ell)}\big{|}_{\mathcal{M}}\big{|}{X}^{(\ell)}\big{|}_{\mathcal{M}}\Big{|}_{\mathcal{M}}\bigg{)}\Bigg{)}.$
(6)

The same approach applies for Eqs. (3) and (4) in the backward pass.

The moduli set $\mathcal{M}$ must be chosen to ensure that the outputs of the
RNS operations are smaller than $M$, which means that

$\log_{2}M\geq
b_{\operatorname{\text{out}}}=b_{\operatorname{\text{in}}}+b_{w}+\log_{2}(h)-1,$
(7)

should be guaranteed for a dot product between
$b_{\operatorname{\text{in}}}$-bit input and $b_{w}$-bit weight vectors with
$h$-elements. This constraint prevents overflow in the computation.

### Precision and Energy Efficiency in the RNS-based Analog Core

The selection of moduli set $\mathcal{M}$, which is constrained by Eq. (7),
impacts the achievable precision at the output as well as the energy
efficiency of the RNS-based analog core. Table 1 compares RNS-based analog
GEMM cores with example moduli sets and regular fixed-point analog GEMM cores.
Here, we show two cases for the regular fixed-point representation: (1) the
low-precision (LP) case where
$b_{\operatorname{\text{out}}}>b_{\operatorname{\text{ADC}}}=b_{\operatorname{\text{DAC}}}$,
and (2) the high-precision (HP) case where
$b_{\operatorname{\text{out}}}=b_{\operatorname{\text{ADC}}}>b_{\operatorname{\text{DAC}}}$.
It should be noted that all three analog cores represent data as fixed-point
numbers. We use the term ‘regular fixed-point core’ to refer to a typical
analog core that performs computations in the standard representation (without
RNS). ‘RNS-based core’ refers to an analog core that performs computations on
the fixed-point residues.

While the LP approach introduces
$b_{\operatorname{\text{out}}}-b_{\operatorname{\text{ADC}}}$ bits of
information loss in every dot product, the HP approach uses high-precision
ADCs to prevent this loss. For the RNS-based core, we picked
$b_{\operatorname{\text{in}}}=b_{w}=b_{\operatorname{\text{ADC}}}=b_{\operatorname{\text{DAC}}}=\lceil\log_{2}m_{i}\rceil\equiv
b$ for ease of comparison against the fixed-point cores. Table 1 shows example
moduli sets that are chosen to guarantee Eq. (7) for $h=128$ while keeping the
moduli under the chosen bit-width $b$. In this case, for $n$ moduli with bit-
width of $b$, $M$ covers $\approx n\cdot b$ bits of range at the output. $h$
is chosen to be 128 as an example considering the common layer sizes in the
evaluated MLPerf (Inference: Datacenter) benchmarks. The chosen $h$ provides
high throughput with high utilization of the GEMM core.

Fig. 2a compares the error (with respect to FP32 results) observed when
performing dot products with the RNS-based core and the LP fixed-point core
with the same bit precision. Both cores use the configurations described in
Table 1 for the example vector size $h=128$. The larger absolute error
observed in the LP fixed-point case illustrates the effect of the information
loss mentioned above. HP fixed-point case is not shown as it is equivalent to
the RNS case.

Fig. 2b shows the energy consumption of DACs and ADCs per dot product for the
three aforementioned analog hardware configurations. To achieve the same MVM
throughput as the (LP/HP) fixed-point cores, the RNS-based core with $n$
moduli must use $n$ distinct MVM units and $n$ sets of DACs and ADCs. This
makes the energy consumption of the RNS-based core $n\times$ larger compared
to the LP fixed-point approach. However, the LP fixed-point approach with low-
precision ADCs experiences information loss in the partial outputs and hence
has lower accuracy.

The RNS-based approach and the HP fixed-point approach provide the same bit
precision (i.e., the same DNN accuracy). Yet, using the RNS-based approach is
orders of magnitude more energy-efficient than the HP fixed-point approach.
This is mainly because of the high cost of high-precision ADCs required to
capture the full output in the HP fixed-point approach. ADCs dominate the
energy consumption with approximately three orders of magnitude higher energy
usage than DACs with the same bit precision. In addition, energy consumption
in ADCs increases exponentially with increasing bit precision 19. This favors
using multiple DACs and ADCs with lower precision in the RNS-based approach
over using a single high-precision ADC. Briefly, the RNS-based approach
provides a sweet spot between LP and HP fixed-point approaches without
compromising on both high accuracy and high energy efficiency.

Figure 2: Comparison of the RNS-based and regular fixed-point analog approaches. a The distribution of average error observed at the output of a dot product performed with the RNS-based analog approach (pink) and the LP regular fixed-point analog approach (cyan). Error is defined as the distance from the result calculated in FP32. The experiments are repeated for 10,000 randomly generated vector pairs with vector size $h=128$. b Energy consumption of data converters (i.e., DACs and ADCs) per dot product for the RNS-based analog approach (pink) and the LP (cyan) and HP (dark blue) regular fixed-point analog approaches. See Methods for the energy estimation methodology.  Table 1: Data and data converter precision in RNS-based, LP fixed-point, and HP fixed-point analog cores.  | RNS-based Core (This work) | LP Fixed-Point Core | HP Fixed-Point Core
---|---|---|---
$b_{\operatorname{\text{in}}}$, $b_{w}$ | $b_{\operatorname{\text{DAC}}}$ | ${\log}_{2}\mathcal{M}$ | $b_{\operatorname{\text{ADC}}}$ | Moduli Set $(\mathcal{M})$ | RNS Range ($M$) | $b_{\operatorname{\text{DAC}}}$ | $b_{\operatorname{\text{out}}}$ | $b_{\operatorname{\text{ADC}}}$ | Lost Bits | $b_{\operatorname{\text{DAC}}}$ | $b_{\operatorname{\text{out}}}$ | $b_{\operatorname{\text{ADC}}}$
4 | 4 | 4 | 4 | $\\{15,14,13,11\\}$ | $\simeq 2^{15}-1$ | 4 | 14 | 4 | 10 | 4 | 14 | 14
5 | 5 | 5 | 5 | $\\{31,29,28,27\\}$ | $\simeq 2^{19}-1$ | 5 | 16 | 5 | 11 | 5 | 16 | 16
6 | 6 | 6 | 6 | $\\{63,62,61,59\\}$ | $\simeq 2^{24}-1$ | 6 | 18 | 6 | 12 | 6 | 18 | 18
7 | 7 | 7 | 7 | $\\{127,126,125\\}$ | $\simeq 2^{21}-1$ | 7 | 20 | 7 | 13 | 7 | 20 | 20
8 | 8 | 8 | 8 | $\\{255,254,253\\}$ | $\simeq 2^{24}-1$ | 8 | 22 | 8 | 14 | 8 | 22 | 22
Figure 3: Accuracy performance of the RNS-based analog core. a Inference
accuracy of regular fixed-point (LP) and RNS-based cores (See Table 1) on
MLPerf (Inference: Datacenters) benchmarks. The accuracy numbers are
normalized to the FP32 accuracy. b-d Loss during training for FP32 and RNS-
based approaches with varying moduli bit-width. ResNet50 (a) is trained from
scratch for 90 epochs using SGD optimizer with a momentum. BERT-Large (b) and
OPT-125M (c) are fine-tuned from pre-trained models. Both models are fine-
tuned using the Adam optimizer with a linear learning rate scheduler for 2 and
3 epochs for BERT-Large and OPT-125M, respectively. All inference and training
experiments use FP32 for all non-GEMM operations.  Figure 4: RNS-based analog
GEMM dataflow. The operation is shown for a moduli set
$\mathcal{M}=\\{m_{1},\dots,m_{n}\\}$. The $n$ $h\times h$ analog MVM units
are represented as generic blocks. The dataflow is agnostic of the technology.

### Accuracy in the RNS-based Analog Core

Fig. 3a compares the inference accuracy of MLPerf (Inference: Datacenters)
benchmarks 24 and OPT 29 (a transformer-based LLM) when run on an RNS-based
analog core and a fixed-point (LP) analog core. The HP fixed-point analog core
results are not shown as they are equivalent to the RNS-based results. The
evaluated DNNs, their corresponding tasks, and the datasets are listed in
Table 2. Fig. 3a shows that the RNS-based approach significantly ameliorates
the accuracy drop caused by the low-precision ADCs used in the LP fixed-point
approach for all the networks. By using the RNS-based approach, it is possible
to achieve ${\geq}99\%$ of FP32 accuracy (this cut-off is defined in the
MLPerf benchmarks 24) for all evaluated benchmarks when using residues with as
low as $6$ bits. This number can be lowered to $5$ bits for BERT-large and
RNN-T and to $4$ bits for DLRM.

Besides inference, the RNS-based approach opens the door for analog computing
to be used in tasks that require higher precision than inference such as DNN
training. Figure 3b shows the loss during DNN training/fine-tuning. Table
LABEL:table:training reports the validation accuracies after FP32 and RNS-
based low-precision training. Here, the GEMM operations during forward and
backward passes of training follow the same methodology as inference, with
weight updates carried out in FP32. Our experiments show that ${\geq}99\%$
FP32 validation accuracy is achievable after training ResNet50 from scratch
using the RNS-based approach with only $6$-bit moduli. Similarly, fine-tuning
BERT-large and OPT-125M by using $5$-bit and $7$-bit moduli, respectively, can
reach ${\geq}99\%$ FP32 validation accuracy. The results are noticeably
promising as the previous efforts on analog DNN hardware that adopted the LP
fixed-point approach had never successfully demonstrated the training of
state-of-the-art DNNs due to the limited precision of this approach.

Fig. 4 illustrates the dataflow of the RNS-based analog core when performing
MVM as part of the DNN inference/training. An input vector $X$ and a weight
matrix $W$ to be multiplied in the MVM unit are first mapped to signed
integers. To mitigate the quantization effects, $X$ and each row in $W$ are
scaled by an FP32 scaling factor that is unique to the vector (See Methods).
The signed integers are then converted into RNS residues through modulo
operation (i.e., forward conversion). By construction, each residue is within
the range of $[0,m_{i})$. To achieve the same throughput as a fixed-point
analog core, the RNS-based analog core with $n$ moduli requires using $n$
analog MVM units—one for each modulus—and running them in parallel. Each
analog MVM unit requires a set of DACs for converting the associated input and
weight residues into the analog domain. The MVM operations are followed by an
analog modulo operation on each output residue vector. Thanks to the modulo
operation, the output residues—to be captured by ADCs—are reduced back to the
$[0,m_{i})$ range. Therefore, a bit precision of $\lceil\log_{2}{m_{i}}\rceil$
is adequate for both DACs and ADCs to perform input and output conversions
without any information loss. The output residues are then converted back to
the standard representation in the digital domain using Eq. (1) to generate
the signed-integer output vector, which is then mapped back to an FP32 final
output $Y$. The non-linear function $f$ (e.g., ReLU, sigmoid, etc.) is then
performed digitally in FP32.

Table 2: MLPerf (Inference: Datacenters) benchmarks.  DNN |  Task |  Dataset
---|---|---
ResNet50 |  Image classification |  ImageNet 23
SSD-ResNet34 |  Object detection |  MS COCO 35
BERT-Large |  Question answering |  SQuADv1.1 36
RNN-T |  Speech recognition |  Librispeech 37
DLRM |  Recommendation |  1TB Click Logs 38
OPT-125M |  Language Modeling |  Wikitext 39
OPT-350M |  Language Modeling |  Wikitext
Table 3: Validation accuracy results after training/fine-tuning.  |  ResNet50 |  BERT-Large |  OPT-125M
---|---|---|---
Precision |  Acc.($\%$) |  F1 Score ($\%$) |  Acc.($\%$)/PPL
FP32 |  75.80 |  91.03 |  43.95/19.72
8-bit |  75.77 |  90.98 |  43.86/20.00
7-bit |  75.68 |  90.97 |  43.59/20.71
6-bit |  75.13 |  90.85 |  42.79/22.62
5-bit |  59.72 |  90.81 |  41.45/26.17
4-bit |  42.15 |  89.66 |  38.64/35.65

### Redundant RNS for Fault Tolerance

Analog compute cores are sensitive to noise. In the case of RNS, even small
errors in the residues can result in a large error in the corresponding
integer they represent. The Redundant Residue Number System (RRNS) 26, 25, 27
can detect and correct errors—making the RNS-based analog core fault tolerant.
RRNS uses a total of $n+k$ moduli: $n$ non-redundant and $k$ redundant. An
RRNS($n+k,n$) code can detect up to $k$ errors and can correct up to
$\lfloor\frac{k}{2}\rfloor$ errors. In particular, the error in the codeword
(i.e., the $n+k$ residues representing an integer in the RRNS space) can be
one of the following cases:

  * •

Case 1: Fewer than $\lfloor\frac{k}{2}\rfloor$ residues have errors—thereby
they are correctable,

  * •

Case 2: Between $\lfloor\frac{k}{2}\rfloor$ and $k$ residues have errors or
the codeword with more than $k$ errors does not overlap with another codeword
in the RRNS space—thereby the error is detectable,

  * •

Case 3: More than $k$ residues have errors and the erroneous codeword overlaps
with another codeword in the RRNS space—thereby the error goes undetected.

Errors are detected by using majority logic decoding wherein we divide the
total $n+k$ output residues into ${n+k}\choose n$ groups with $n$ residues per
group. One simple way of majority logic decoding in this context is to convert
the residues in each group back to the standard representation via CRT to
generate an output value for each group and compare the results of the
${n+k}\choose n$ groups. If more than $50\%$ of the groups have the same
result in the standard representation, then the generated codeword is correct.
This corresponds to Case 1. In contrast, not having a majority indicates that
the generated codeword is erroneous and cannot be corrected. This corresponds
to Case 2. In this case, the detected errors can be eliminated by repeating
the entire calculation. In Case 3, the erroneous codeword generated by the
majority of the groups overlaps with another codeword. As a result, more than
$50\%$ of the groups have the same incorrect result and the error goes
undetected. To optimize the hardware performance of this process, more
efficient base-extension-based algorithms 30 instead of CRT can be used for
error detection.

The final error probability in an RRNS code depends on the percentage of the
_non-correctable_ errors observed in the residues. This probability is
influenced by the chosen moduli set and the number of error correction
iterations (See Methods). Let $p_{c},p_{d},$ and $p_{u}$ be the probabilities
that Cases 1, 2, and 3 occur respectively when computing a single output.
Overall, $p_{c}+p_{d}+p_{u}=1$. For a single attempt (i.e., $R=1$), the
probability of producing the incorrect output integer is
$p_{\text{err}}(R=1)=1-p_{c}=p_{u}+p_{d}$. Generally, it is possible to repeat
the calculations $R$-times until no detectable error is found at the expense
of increasing compute latency. In this case, the probability of having an
incorrect output after $R$ attempts of error correction is

$p_{\text{err}}(R)=1-p_{c}\sum_{r=0}^{R-1}(p_{d})^{r}.$ (8)

As the number of attempts increases, the output error probability decreases
and converges to $\lim_{R\to\infty}p_{\text{err}}(R)=p_{u}/(p_{u}+p_{c})$.

Fig. 5 shows $p_{\text{err}}$ for different numbers of redundant moduli $(k)$,
attempts $(R)$, and moduli sets with different bit-widths. Broadly, as the
probability of a single residue error $p$ increases, the output error
probability tends to $1$. For a given number of attempts, increasing bit
precision and the number of redundant moduli decreases $p_{\text{err}}$. For a
fixed number of redundant moduli and a fixed number of bits per moduli,
$p_{\text{err}}$ decreases as the number of attempts increases.

Fig. 6 investigates the impact of noise on the accuracy of two large and
important MLPerf benchmarks—ResNet50 and BERT-Large—when using RRNS. The two
networks show similar behavior: adding extra moduli and increasing the number
of attempts decrease $p_{\text{err}}$ at the same value of $p$. ResNet50
requires ${\sim}3.9$ GigaMAC operations (GOp) for one inference on a single
input image. For a $128\times 128$ MVM unit, inferring an ImageNet image
through the entire network involves computing ${\sim}29.4$M partial output
elements. Therefore, we expect the transition point from an accurate network
to an inaccurate network to occur at $p_{\text{err}}$ to be $\leq 1/29.4$M =
$3.4\times 10^{-8}$. This $p_{\text{err}}$ transition point is $\leq 1/358.6$M
= $2.8\times 10^{-9}$ for BERT-Large. Fig. 6, however, shows that the
evaluated DNNs are more resilient to noise than expected: it is able to
tolerate higher $p_{\text{err}}$ while maintaining good accuracy. The accuracy
of ResNet50 only starts degrading (below $99\%$ FP32) when
$p_{\text{err}}\approx 4.5\times 10^{-5}$ ($1000\times$ higher than the
estimated value) on average amongst the experiments shown in Figure 6. This
transition probability is $p_{\text{err}}\;{\approx}\;4\times 10^{-4}$ for
BERT-Large (on average $100,000\times$ higher than the estimated value).

Figure 5: Calculated output error probability ($\mathbf{p_{\text{err}}}$)
versus single residue error probability ($\mathbf{p}$). a-c $p_{\text{err}}$
for one (a), two (b), and infinite (c) error correction attempts and a varying
number of redundant moduli $(k)$. Figure 6: Inference accuracy versus single
residue error probability ($\mathbf{p}$). a-f The plots show ResNet-50 (a-c)
and BERT-Large (d-f) inference accuracy results under varying $p$ for one (a
and d), two (b and e), and infinite (c and f) error correction attempts and a
varying number of redundant moduli ($k$).

## Discussion

The RNS (and the fault-tolerant RRNS) framework are agnostic to the analog
technology employed. Generally, the RNS GEMM operations can be performed as a
regular GEMM operation followed by a modulo operation in the analog domain.
Analog GEMM is well-explored in the literature. Previous works leveraged
photonics 1, 2, 7, 3, 4, 5, 6, crossbar arrays consisting of resistive RAM 8,
9, 10, 11, 12, switched capacitors 13, 14, PCM cells 15, STT-RAM 16, 17, etc.

The analog modulo operation can be performed electrically or optically. In the
electronic domain, one can use ring oscillators: a circuit that generates a
continuous waveform by cycling through a series of inverters 28 to perform
modulo operations. By carefully designing the parameters of the ring
oscillator, it is possible to create an output frequency that corresponds to
the desired modulus value. Alternatively, the phase of a signal can be used
for performing modulo due to the periodicity of phases in optical systems.
Optical phase is inherently modular against $2\pi$. By modulating the phase of
an optical signal, one can achieve modulo operations in the analog domain.
Using RNS requires forward and reverse conversion circuits to switch between
the RNS and the standard number system. The forward conversion is a modulo
operation while the reverse conversion can be done using the CRT, mixed-radix
conversion, or look-up tables. The (digital) hardware costs of these circuits
can be reduced by choosing special moduli sets 31, 32.

The RNS framework can be extended with the PNS to work with arbitrary
precision, despite having DACs and ADCs with limited precision. For
applications requiring higher-precision arithmetic than the example cases in
this study (e.g., some high-performance computing applications, homomorphic
encryption, etc.), a higher $M$ value and therefore moduli with higher bit-
width might be necessary, which will be bound by the same limitations
discussed in this paper. Instead, one can represent an integer value as $D$
separate digits where each digit is represented as a set of residues in the
RNS domain and has an RNS range of $M$. This hybrid scheme can achieve
$D\log_{2}M$ bit precision where $D$ can liberally be increased without
increasing the bit precision of the data converters. Different from the RNS-
only scheme, the hybrid scheme requires overflow detection and carry
propagation from lower digits to higher digits. The overflow detection can be
achieved using two sets of residues: primary and secondary. While the
operations are performed with both sets of residues, base extension between
two sets helps detect any overflow and propagate the carry to the higher
digits if required (See Methods).

In conclusion, our work provides a methodology for precise, energy-efficient,
and fault-tolerant analog DNN acceleration. Overall, we believe that RNS is a
crucial numeral system for the development of next-generation analog hardware
capable of both inference and training state-of-the-art neural networks for
advanced applications, such as generative artificial intelligence.

## Methods

### Handling Negative Numbers with RNS

An RNS with dynamic range $M$ allows representing values within the range of
$[0,M)$. This range can be shifted to $[-\psi,\psi]$, where
$\psi=\lfloor(M-1)/2\rfloor$, to represent negative values. This is achieved
by reassigning the values in between $(0,\psi]$ to be positive, 0 to be zero,
and the numbers in between $(\psi,2\psi]$ to be negative (i.e. $[-\psi,0)$).
Then, the values can be recovered uniquely by using CRT with a slight
modification:

$A=\begin{dcases}\sum_{i=1}^{n}|a_{i}M_{i}T_{i}|_{M},&\text{if
}\sum_{i=1}^{n}|a_{i}M_{i}T_{i}|_{M}\leq\psi\\\
\sum_{i=1}^{n}|a_{i}M_{i}T_{i}|_{M}-M,&\text{otherwise}.\end{dcases}$ (9)

### Data Converter Energy Estimation

The DAC and ADC energy numbers in Fig. 2 (b) are estimated by using equations
formulated by Murmann 19, 33. The energy consumption of a DAC per conversion
is

$E_{\operatorname{\text{DAC}}}=\text{ENOB}^{2}C_{u}V_{\text{DD}}^{2},$ (10)

where $C_{u}=0.5$ fF is a typical unit capacitance and $V_{\text{DD}}$ = 1V is
the supply voltage 19. The energy consumption of an ADC per conversion is
estimated as

$E_{\operatorname{\text{ADC}}}=k_{1}\text{ENOB}+k_{2}4^{\text{ENOB}},$ (11)

where $k_{1}{\approx}100$ fJ and $k_{2}{\approx}1$ aJ.
$E_{\operatorname{\text{ADC}}}$ is dominated by the exponential term (i.e.,
$k_{2}4^{\text{ENOB}}$) at large ENOB ( $\geq 10$-bits).

### Accuracy Modeling

Both RNS-based and regular fixed-point analog cores are modeled using PyTorch
for estimating inference and training accuracy. Convolution, linear, and
batched matrix multiplication (BMM) layers are performed as tiled-GEMM
operations which are computed tile-by-tile as a set of tiled-MVM operations.
Each input, weight, and output of the tiled MVM are quantized with a desired
bit precision.

Pre-quantization, the input vectors and weight tiles are first dynamically
scaled, i.e. scaled at runtime, to mitigate the quantization effects as
follows: For an $h\times h$ weight tile $\mathcal{W}_{t}$, we denote each row
vector as ${\mathcal{W}}_{rt}$ where the subscript $r$ stands for the row and
$t$ for the tile. Similarly, an input vector of length $h$ is denoted as
$\mathcal{X}_{t}$ where $t$ indicates the tile. Each weight row
$\mathcal{W}_{rt}$ shares a single FP32 scale
$s_{rt}^{w}=\max(|\mathcal{W}_{rt}|)$ and each input vector $\mathcal{X}_{t}$
shares a single FP32 scale $s_{t}^{x}=\max(|\mathcal{X}_{t}|)$. $h$ scales per
$h\times h$ weight tile and 1 scale per input vector, in total $h+1$ scales,
are stored for each tiled-MVM operation. The tiled MVM is performed between
the scaled weight and input vectors,
$\widehat{\mathcal{W}}_{rt}=\mathcal{W}_{rt}/s_{rt}^{w}$ and
$\widehat{\mathcal{X}}_{t}=\mathcal{X}_{t}/s_{t}^{x}$, to produce
$\widehat{Y}_{rt}=\widehat{\mathcal{W}}_{t}\widehat{\mathcal{X}}_{t}$. The
output $\widehat{Y}_{rt}$ is then quantized (if required) to resemble the
output ADCs and multiplied back with the appropriate scales so that the actual
output elements $Y_{rt}=\widehat{Y}_{rt}\cdot s_{rt}^{w}\cdot s_{t}^{x}$ are
obtained.

Here, the methodology is the same for RNS-based and regular fixed-point cores.
For the RNS-based case, in addition to the description above, the quantized
input and weight integers are converted into the RNS space before the tiled-
MVM operations. MVMs are performed separately for each set of residues and are
followed by a modulo operation before the quantization step. The output
residues for each tiled MVM are converted back to the standard representation
using the CRT.

The GEMM operations (i.e., convolution, linear, and BMM layers) are sandwiched
between an input operation $O_{\text{in}}$ and an output operation
$O_{\text{out}}$. This makes the operation order
$O_{\text{in}}$-GEMM-$O_{\text{out}}$ during the forward pass, and
$O_{\text{out}}$-GEMM-$O_{\text{in}}$ in the backward pass. $O_{\text{in}}$
quantizes the input and weight tensors in the forward pass and is a null
operation in the backward pass. In contrast, $O_{\text{out}}$ is a null
operation in the forward pass and quantizes the activation gradients in the
backward pass. In this way, the quantization is always performed before the
GEMM operation. The optimizer (i.e., SGD or Adam) is modified to keep a copy
of the FP32 weights to use during the weight updates. Before each forward
pass, the FP32 weights are copied and stored. After the forward pass, the
quantized model weights are replaced by the previously stored FP32 weights
before the step function so that the weight updates are performed in FP32.
After the weight update, the model parameters are quantized again for the next
forward pass. This high-precision weight update step is crucial for achieving
high accuracy in training.

In Fig. 3b, all the convolution, linear, and BMM layers in the models were
replaced by the quantized versions. We trained ResNet50 from scratch by using
SGD optimizer for 90 epochs with a momentum of 0.9 and a learning rate
starting from 0.1. The learning rate was scaled down by 10 at epochs 30, 60,
and 80. We fine-tuned BERT-large and OPT-125M from the implementations
available in the Huggingface transformers repository 34. We used the Adam
optimizer for both models with the default settings. The script uses a linear
learning rate scheduler. The learning rate starts at 3e-05 and 5e-05 and the
models are trained for 2 and 3 epochs, respectively for BERT-Large and
OPT-125M.

### Error distribution in the RRNS code space

For an RRNS$({n+k},n)$ with $n$ non-redundant moduli $(m_{1},m_{2},...,m_{n})$
and $k$ redundant moduli $(m_{n+1},m_{n+2},...,m_{n+k})$, the probability
distributions ($p_{c}$, $p_{d}$, and $p_{u}$) of different types of errors
(Case 1, Case 2, and Case 3 that were mentioned in Redundant RNS for Fault
Tolerance) are related to the Hamming distance distribution of the RRNS code
space. In an RRNS$({n+k},n)$, every integer is represented as ${n+k}$ residues
($r_{i}$ where $i\in\\{1,...,{n+k}\\}$) and this vector of ${n+k}$ residues is
considered as an RRNS codeword. A Hamming distance of
$\eta\in\\{0,1,...,{n+k}\\}$ between the original codeword and the erroneous
codeword indicates that $\eta$ out of ${n+k}$ residues are erroneous. The
erroneous codewords create a new vector space of ${n+k}$-long vectors where at
least one $r_{i}$ is replaced with $r_{i}^{\prime}\neq r_{i}$ with
$i\in\\{1,...,{n+k}\\}$ and $r_{i}^{\prime}<m_{i}$. This vector space includes
all the RRNS$({n+k},n)$ codewords as well as other possible ${n+k}$-long
vectors that do not overlap with any codeword in the RRNS code space. A vector
represents a codeword and is in the RRNS code space if and only if it can be
converted into a value within the legitimate range $[0,M)$ of the
RRNS$({n+k},n)$ by using the CRT. The number of all vectors that have a
Hamming distance $\eta$ from a codeword in RRNS$({n+k},n)$ can be expressed as

$V_{\eta}=\sum_{Q\binom{n+k}{\eta}}\prod_{i=1}^{\eta}(m_{i}-1),$ (12)

where ${Q\binom{n+k}{\eta}}$ represents one selection of $\eta$ moduli from
${n+k}$ moduli while $\sum_{Q\binom{n+k}{\eta}}$ represents the summation over
all distinct $\binom{n+k}{\eta}$ selections. The number of codewords that are
in the RNS code space with a Hamming distance of $\eta\in\\{0,1,...,{n+k}\\}$
can be expressed as

$D_{\eta}=\sum_{h=0}^{\eta-1-k}(-1)^{h}\binom{n+k-\eta+h}{n+k-\eta}\zeta(n+k,\eta-h),$
(13)

for $k+1\leq\eta\leq n+k$. For $1\leq\eta\leq k$, $D_{\eta}=0$ and $D_{0}=1$.
$\zeta(n+k,\eta)$ represents the total number of non-zero common divisors in
the legitimate range $[0,M)$ for any $n+k-\eta$ moduli out of the $n+k$ moduli
of the RRNS$(n+k,n)$ code and can be denoted as

$\zeta(n+k,\eta)=\sum_{Q\binom{n+k}{n+k-\eta}}\Bigg{\lfloor}\frac{M-1}{m_{i_{1}}m_{i_{2}}...m_{i_{n+k-\eta)}}}\Bigg{\rfloor},$
(14)

where $(m_{i_{1}},m_{i_{2}},...,m_{i_{\lambda}})$ with $1\leq\lambda\leq n+k$
is a subset of the RRNS$(n+k,n)$ moduli set.

An undetectable error occurs only if a codeword with errors overlaps with
another codeword in the same RRNS space. Given the distance distributions for
the vector space $V$ and the codespace $D$ (Eq. (12), (13), respectively), the
probability of observing an undetectable error ($p_{u}$) for RRNS$(n+k,n)$ can
be computed as

$p_{u}=\sum_{\eta=k+1}^{n+k}\frac{D_{\eta}}{V_{\eta}}p_{E}(\eta),$ (15)

where $p_{E}(\eta)$ is the probability of having $\eta$ erroneous residues in
a codeword which can be calculated as

$p_{E}(\eta)=\sum_{Q\binom{n+k}{\eta}}p^{\eta}(1-p)^{(n+k-\eta)},$ (16)

for an error probability of $p$ in a single residue.

Eq. (13) indicates that for up to $\eta=k$ erroneous residues $D_{\eta}=0$,
and so it is not possible for an erroneous codeword to overlap with another
codeword in the RRNS code space. This guarantees the successful detection of
the observed error. If the Hamming distance of the erroneous codeword is
$\eta\leq\lfloor\frac{k}{2}\rfloor$, the error can be corrected by the
majority logic decoding mechanism. In other words, the probability of
observing a correctable error is equal to observing less or equal to
$\lfloor\frac{k}{2}\rfloor$ errors in the residues and can be calculated as

$p_{c}={\sum_{\eta=0}^{\lfloor\frac{k}{2}\rfloor}}p_{E}(\eta)=\sum_{\eta=0}^{\lfloor\frac{k}{2}\rfloor}\Big{(}\sum_{Q\binom{n+k}{\eta}}p^{\eta}(1-p)^{(n+k-\eta)}\Big{)}.$
(17)

All the errors that do not fall under the undetectable or correctable
categories are referred to as detectable but not correctable errors with a
probability $p_{d}$ where $p_{d}=1-(p_{c}+p_{d})$. The equations in this
section were obtained from the work conducted by Yang 27.

To model the error in the RNS core for the analysis shown in Fig. 6, $p_{c}$,
$p_{d}$, and $p_{u}$ are computed for a given RRNS$(n+k,n)$ using Eqs. (15)
and (17). Given the number of error correction attempts, output error
probability ($p_{err}$) is calculated according to Eq. (8). Random noise is
injected at the output of every tiled-MVM operation using a Bernoulli
distribution with the probability of $p_{err}$.

### RNS Operations

The proposed analog RNS-based approach requires modular arithmetic. In this
section, we discuss two ways of performing modular arithmetic in the analog
domain, one electrical and one optical.

#### Modular Arithmetic with Ring Oscillators

In a ring oscillator, where each inverter has a propagation delay
$t_{\text{prop}}>0$, there is always only one inverter that has the same input
and output—either 1-1 or 0-0—at any given time when the ring oscillator is on.
The location of this inverter with the same input and output propagates along
with the signal every $t_{\text{prop}}$ time and rotates due to the ring
structure. This rotation forms a modular behavior in the ring when the
location of this inverter is tracked.

Let $S_{\text{RO}}(t)$ be the state of a ring oscillator with $N$ inverters.
Here, $S_{\text{RO}}(t)\in\\{0,...,N-1\\}$ and $S_{\text{RO}}(t)=k$ means that
the $k+1$-th inverter’s input and output have the same value at time $t$.
$S_{\text{RO}}(t)$ keeps rotating between $0$ to $N-1$ as long as the
oscillator is on. Figure 7a shows a simple example where $N=3$. In the first
$t_{\text{prop}}$ time interval, the input and output of the first inverter
are both 0, therefore, the state $S_{\text{RO}}(t<t_{\text{prop}})=0$.
Similarly, when $t_{\text{prop}}<t<2t_{\text{prop}}$, the input and output of
the second inverter are 1, so
$S_{\text{RO}}(t_{\text{prop}}<t<2t_{\text{prop}})=1$. Here, the time between
two states following one another (i.e., $t_{\text{prop}}$) is fixed and
$S_{\text{RO}}(t)$ rotates ($0,1,2,0,1,...$). Assume the state of the ring
oscillator is sampled periodically with a sampling period of $T_{s}=A\cdot
t_{\text{prop}}$. Then, the observed change in the state of the ring
oscillator between two samples
($S_{\text{RO}}({t=T_{s}})-S_{\text{RO}}({t=0})$) is equivalent to $|A|_{N}$
where $A$ is a positive integer value. Therefore, to perform a modulo against
a modulus value $m$, the number of inverters $N$ should be equal to $m$. The
dividend number $A$ and the sampling period can be adjusted by changing the
analog input voltage to a voltage-to-time converter (VTC).

In a modular dot product or an MVM operation, the dividend $A$ is replaced by
the output of the dot product. Analog dot products can be performed using
traditional methods with no change with any desired analog technology where
output can be represented as an analog electrical signal (e.g., current or
voltage) before the analog modulo.

#### Modular Arithmetic with Phase Shifters

Figure 7: Analog modulo implementations. a Modulo operation performed using a
ring oscillator. A ring oscillator with $N=3$ inverters is shown to perform
modulo against a modulus $m=3$. This operation is performed after every analog
dot product to perform a modular dot product. b Modular dot product performed
using phase shifters. A modular dot product operation between two 2-element
vectors $x$ and $w$, each with 3 digits, is shown by using a dual-rail set of
cascaded phase shifters. The transistor switch turns on and supplies voltage
to the phase shifter when the corresponding digit of $w$ is 1.

The amount of phase shift introduced by a single dual-rail phase shifter when
$v$ and $-v$ voltages are applied on the upper and the bottom arms,
respectively, is

$\Delta\Phi=\frac{vL}{V_{\pi\cdot\text{cm}}},$ (18)

where $V_{\pi\cdot\text{cm}}$ is the modulation efficiency of the phase
shifter and is a constant value. $\Delta\Phi$ is proportional to both the
length of the shifter $L$ and the amount of applied voltage $v$. Figure 7b
shows an example modular dot product operation between two vectors, $x$ and
$w$, using cascaded dual-rail phase shifters. $w$ is encoded digit-by-digit
using phase shifters with lengths proportional to $2^{j}$ where $j$ represents
the digit number. In the example, each element (i.e., $w_{0}$ and $w_{1}$) of
the 2-element vector $w$ consists of $3$ digits and uses $3$ phase shifters,
each with lengths $L,2L,$ and $4L$. If the $j$-th digit of the $i$-th element
of $w$, $w_{i}^{j}=1$, a voltage $v_{i}$ is applied to the phase shifter pair
(top and bottom) with the length $2^{j}L$. If the digit $w_{i}^{j}=0$, then no
voltage is applied and therefore no phase shift is introduced to the input
signal. To encode the second operand $x$, a voltage $v_{i}$ that is
proportional to $x_{i}$ is applied to all non-zero digits of $w_{i}$. To take
modulo with a modulus $m$ instead of $2\pi$, the input $x$ and therefore the
applied voltage $v$ should be multiplied by the constant $2\pi/m$. For
encoding $x_{i}$,

$v_{i}=x_{i}\cdot\frac{V_{\pi\cdot\text{cm}}}{\pi L}\cdot\frac{2\pi}{m},$ (19)

should be applied so that the total phase shift at the end of the optical path
is

$\Delta\Phi_{\text{total}}=\big{|}\frac{2\pi}{m}\sum_{i}\big{(}\sum_{j}(2^{j}w_{i}^{j})x_{i}\big{)}\big{|}_{2\pi}=\frac{2\pi}{m}\big{|}\sum_{i}(w_{i}x_{i})\big{|}_{m}.$
(20)

The resulting output values are collected at the end of the optical path and
are in the form of the phase difference between input and output. These
outputs are then re-multiplied by $m/2\pi$ to obtain the outputs of the
modular dot products for each residue.

### Extended RNS

By combining RNS and PNS, an integer value $Z$ can be represented as $D$
separate digits, $z_{d}$ where $d\in\\{0,1,...,{D-1}\\}$ and $0\leq z_{d}<M$
where $M$ is the RNS range:

$Z=\sum_{d=0}^{D-1}z_{d}M,$ (21)

and can provide up to $D\log_{2}M$ bit precision. This hybrid scheme requires
carry propagation from lower digits to higher digits, unlike the RNS-only
scheme. For this purpose, one can use two sets of moduli, primary and
secondary, where every operation is performed for both sets of residues. After
every operation, overflow is detected for each digit and carried over to the
higher-order digits.

Let us define and pick $n_{p}$ primary moduli $m_{i}$, where
$i\in\\{1,...,n_{p}\\}$ and $n_{s}$ secondary moduli $m_{j}$, where
$j\in\\{1,...,n_{s}\\}$ and $m_{i}\neq m_{j}\ \forall\ \\{i,j\\}$. Here
$M=M_{p}\cdot M_{s}=\prod_{i=1}^{n_{p}}m_{i}\cdot\prod_{j=1}^{n_{s}}m_{j}$ is
large enough to represent the largest possible output of the operations
performed in this numeral representation and $M_{p}$ and $M_{s}$ are co-prime.

To execute an operation in this hybrid number system, the operation is
performed separately for each digit of the output. These operations for each
digit are independent of one another and can be parallelized except for the
overflow detection and carry propagation. Assume $z_{d}=z_{d}|_{p;s}$ consists
of primary and secondary residues and is a calculated output digit of an
operation before overflow detection.

$z_{d}$ can be decomposed as $z_{d}|_{p}=Q_{d}|_{p}M_{p}+R_{d}|_{p}$ where
$Q_{d}|_{p}$ and $R_{d}|_{p}$ are the quotient and the remainder of the digit,
with respect to the primary RNS. To detect a potential overflow in the digit
$z_{d}$, a base extension from primary to secondary RNS is performed on
$z_{d}|_{p}$ and the base extended residues are compared with the original
secondary residues of the digit, $z_{d}|_{s}$. If the residues are the same,
this indicates that there is no overflow, i.e., $Q_{d}|_{p;s}=0$, and both
primary and secondary residues are kept without any carry moved to the next
higher digit. As against that, if the base-extended secondary residues and the
original secondary residues are not the same, it means that there exists an
overflow (i.e., $Q_{d}|_{p;s}\neq 0$). In the case of overflow, the remainder
of the secondary RNS, $R_{d}|_{s}$, is calculated through a base extension
from primary to secondary RNS on $R_{d}|_{p}$ where $R_{d}|_{p}=z_{d}|_{p}$.
$Q_{d}|_{s}$ can then be computed as
$Q_{d}|_{s}=(z_{d}|_{s}-R_{d}|_{s})M_{p}^{-1}$ where $|M_{p}\cdot
M_{p}^{-1}|_{M_{s}}\equiv 1$. $Q_{d}|_{p}$ is calculated through base
extension from the secondary to primary RNS on the computed $Q_{d}|_{s}$. The
full quotient $Q_{d}|_{p;s}$ is then propagated to the higher-order digit.
Algorithm LABEL:alg:pseudocode shows the pseudo-code for handling an operation
$\square$ using the extended RNS representation. The operation can be replaced
by any operation that is closed under RNS.

It should be noted that $z_{d}|_{p;s}$ is not always computed as
$x_{d}|_{p;s}\square y_{d}|_{p;s}$. For operations such as addition, each
digit before carry propagation is computed by simply adding the same digits of
the operands, i.e., $z_{d}|_{p;s}=x_{d}|_{p;s}+y_{d}|_{p;s}$. However, for
multiplication, each digit of $z_{d}|_{p;s}$ should be constructed as in long
multiplication. The multiplication of two numbers in the hybrid number system
with $D_{x}$ and $D_{y}$ digits requires $D_{x}D_{y}$ digit-wise
multiplications and the output will result in $D_{z}=D_{x}+D_{y}$ digits in
total. Similarly, a dot product is a combination of multiply and add
operations. If two vectors with $h$ elements where each element has $D_{x}$
and $D_{y}$ digits, the output will require in $D_{z}=D_{x}+D_{y}+\log_{2}h$
digits.

Algorithm 1 Pseudocode for performing the operation $\square$ using the hybrid
number system. Here, $x$ and $y$ are the inputs for operation $\square$ and
$z$ is the output with $D$ digits. $z_{d}$ represents the digits of the output
where $z_{d}|_{p}$ are the primary residues and $z_{d}|_{s}$ are the secondary
residues. Primary and secondary residues together are referred to as
$z^{\prime}_{d}|_{p;s}$. $Q$ is the quotient and $R$ is the remainder where
$z_{d}=Q_{d}M_{p}+R_{d}$. $\text{p2s}()$ and $\text{s2p}()$ refer to base
extension algorithms from primary to secondary residues and from secondary to
primary residues, respectively.

1:$Q_{-1}|_{p;s}=0$

2:for $d\leftarrow 0$ to $D_{z}$ do

3: $z^{\prime}_{d}|_{p;s}={(x|_{p;s}}\ \square\ {y|_{p;s})_{d}}$

4:end for

5:for $d\leftarrow 0$ to $D_{z}$ do

6: $z_{d}|_{p;s}=z^{\prime}_{d}|_{p;s}+Q_{d-1}|_{p;s}$

7: $R_{d}|_{p}=z_{d}|_{p}$

8: $R_{d}|_{s}=\text{p2s}(R_{d}|_{p})$

9: if $R_{d}|_{s}=z^{\prime}_{d}|_{s}$ then

10: $Q_{d}|_{p;s}=0$

11: else

12: $Q_{d}|_{s}=(z^{\prime}_{d}|_{s}-R_{d}|_{s})M_{p}^{-1}$

13: $Q_{d}|_{p}=\text{s2p}(Q_{d}|_{s})$

14: end if

15:end for

## Acknowledgements

We thank Dr. Rashmi Agrawal and Prof. Vijay Janapa Reddi for their insightful
discussions.

## Author Contributions

D.B. conceived the project idea. C.D. and D.B. developed the theory. C.D.
implemented the accuracy modeling and the analytical error models with
feedback from D.B. and A.J. C.D. and L.N. conducted the experiments. D.B. and
A.J. supervised the project. C.D. wrote the manuscript with input from all
authors.

## Competing Interests

The authors declare the following patent application: U.S. Patent Application
No.: 17 / 543,676. L.N. and D.B. declare individual ownership of shares in
Lightmatter, a startup company developing photonic hardware for AI.

## References

  * 1 Shen, Y. et al. Deep learning with coherent nanophotonic circuits. Nature Photonics. 11, 441-446 (2017).
  * 2 Xu, X. et al. 11 TOPS photonic convolutional accelerator for optical neural networks. Nature. 589, 44-51 (2021).
  * 3 Peng, J. et al. A Deep Neural Network Accelerator using Residue Arithmetic in a Hybrid Optoelectronic System. ACM Journal On Emerging Technologies In Computing Systems (JETC). 18, 1-26 (2022).
  * 4 Shiflett, K., Wright, D., Karanth, A. & Louri, A. Pixel: Photonic neural network accelerator. 2020 IEEE International Symposium On High Performance Computer Architecture (HPCA). pp. 474-487 (2020).
  * 5 Shiflett, K., Karanth, A., Bunescu, R. & Louri, A. Albireo: Energy-efficient acceleration of convolutional neural networks via silicon photonics. 2021 ACM/IEEE 48th Annual International Symposium On Computer Architecture (ISCA). pp. 860-873 (2021).
  * 6 Demirkiran, C. et al. An Electro-Photonic System for Accelerating Deep Neural Networks. ACM Journal On Emerging Technologies In Computing Systems (JETC). 19 (2023).
  * 7 Tait, A. et al. Neuromorphic photonic networks using silicon photonic weight banks. Scientific Reports. 7, 7430 (2017).
  * 8 Yao, P. et al. Fully hardware-implemented memristor convolutional neural network. Nature. 577, 641-646 (2020).
  * 9 Chi, P. et al. Prime: A novel processing-in-memory architecture for neural network computation in reram-based main memory. ACM SIGARCH Computer Architecture News. 44, 27-39 (2016).
  * 10 Shafiee, A. et al. ISAAC: A convolutional neural network accelerator with in-situ analog arithmetic in crossbars. ACM SIGARCH Computer Architecture News. 44, 14-26 (2016).
  * 11 Hu, M. et al. Dot-product engine for neuromorphic computing: Programming 1T1M crossbar to accelerate matrix-vector multiplication. Proceedings Of The 53rd Annual Design Automation Conference. pp. 1-6 (2016).
  * 12 Tang, T., Xia, L., Li, B., Wang, Y. & Yang, H. Binary convolutional neural network on RRAM. 2017 22nd Asia And South Pacific Design Automation Conference (ASP-DAC). pp. 782-787 (2017).
  * 13 Bankman, D. & Murmann, B. Passive charge redistribution digital-to-analogue multiplier. Electronics Letters. 51, 386-388 (2015).
  * 14 Bankman, D. & Murmann, B. An 8-bit, 16 input, 3.2 pJ/op switched-capacitor dot product circuit in 28-nm FDSOI CMOS. 2016 IEEE Asian Solid-state Circuits Conference (a-sscc). pp. 21-24 (2016).
  * 15 Feldmann, J. et al. Parallel convolutional processing using an integrated photonic tensor core. Nature. 589, 52-58 (2021).
  * 16 Jain, S., Ranjan, A., Roy, K. & Raghunathan, A. Computing in memory with spin-transfer torque magnetic RAM. IEEE Transactions On Very Large Scale Integration (VLSI) Systems. 26, 470-483 (2017).
  * 17 Shi, Y. et al. Performance prospects of deeply scaled spin-transfer torque magnetic random-access memory for in-memory computing. IEEE Electron Device Letters. 41, 1126-1129 (2020).
  * 18 Chen, Y., Xie, Y., Song, L., Chen, F. & Tang, T. A survey of accelerator architectures for deep neural networks. Engineering. 6, 264-274 (2020).
  * 19 Murmann, B. Mixed-signal computing for deep neural network inference. IEEE Transactions On Very Large Scale Integration (VLSI) Systems. 29, 3-13 (2020).
  * 20 Rekhi, A. et al. Analog/mixed-signal hardware error modeling for deep learning inference. Proceedings Of The 56th Annual Design Automation Conference 2019. pp. 1-6 (2019).
  * 21 Deng, L. The mnist database of handwritten digit images for machine learning research. IEEE Signal Processing Magazine. 29, 141-142 (2012).
  * 22 He, K., Zhang, X., Ren, S. & Sun, J. Deep residual learning for image recognition. Proceedings Of The IEEE Conference On Computer Vision And Pattern Recognition. pp. 770-778 (2016).
  * 23 Deng, J. et al. Imagenet: A large-scale hierarchical image database. 2009 IEEE Conference On Computer Vision And Pattern Recognition. pp. 248-255 (2009).
  * 24 Reddi, V. et al. Mlperf inference benchmark. 2020 ACM/IEEE 47th Annual International Symposium On Computer Architecture (ISCA). pp. 446-459 (2020).
  * 25 James, J. & Pe, A. Error correction based on redundant Residue Number System. 2015 IEEE International Conference On Electronics, Computing And Communication Technologies (CONECCT). pp. 1-5 (2015).
  * 26 Yang, L. & Hanzo, L. Coding theory and performance of redundant residue number system codes. IEEE Transactions on Information Theory. (1999).
  * 27 Yang, L. & Hanzo, L. Redundant residue number system based error correction codes. IEEE 54th Vehicular Technology Conference. VTC Fall 2001. Proceedings. 3 pp. 1472-1476 (2001).
  * 28 Ordentlich, O., Tabak, G., Hanumolu, P., Singer, A. & Wornell, G. A modulo-based architecture for analog-to-digital conversion. IEEE Journal Of Selected Topics In Signal Processing. 12, 825-840 (2018).
  * 29 Zhang, S. et al. Opt: Open pre-trained transformer language models. Preprint at https://arxiv.org/abs/2205.01068 (2022).
  * 30 Babenko, M. et al. RRNS base extension error-correcting code for performance optimization of scalable reliable distributed cloud data storage. 2021 IEEE International Parallel And Distributed Processing Symposium Workshops (IPDPSW). pp. 548-553 (2021).
  * 31 Hiasat, A. & Abdel-Aty-Zohdy, S. Residue-to-binary arithmetic converter for the moduli set (2/sup k/, 2/sup k/-1, 2/sup k-1/-1). IEEE Transactions On Circuits And Systems II: Analog And Digital Signal Processing. 45, 204-209 (1998).
  * 32 Gallaher, D., Petry, F. & Srinivasan, P. The digit parallel method for fast RNS to weighted number system conversion for specific moduli (2/sup k/-1, 2/sup k/, 2/sup k/+ 1). IEEE Transactions On Circuits And Systems II: Analog And Digital Signal Processing. 44, 53-57 (1997).
  * 33 Murmann, B. ADC performance survey 1997-2011. http://web.stanford.edu/${\sim}$murmann/adcsurvey.html. (2011).
  * 34 Wolf, T. et al. Transformers: State-of-the-art natural language processing. Proceedings Of The 2020 Conference On Empirical Methods In Natural Language Processing: System Demonstrations. pp. 38-45 (2020).
  * 35 Lin, T. et al. Microsoft coco: Common objects in context. Computer Vision–ECCV 2014: 13th European Conference Proceedings, Part V 13. pp. 740-755 (2014).
  * 36 Rajpurkar, P., Zhang, J., Lopyrev, K. & Liang, P. Squad: 100,000+ questions for machine comprehension of text. Preprint at https://arxiv.org/abs/1606.05250 (2016).
  * 37 Panayotov, V., Chen, G., Povey, D. & Khudanpur, S. Librispeech: an asr corpus based on public domain audio books. 2015 IEEE International Conference On Acoustics, Speech And Signal Processing (ICASSP). pp. 5206-5210 (2015).
  * 38 Zhao, P., Xiao, K., Zhang, Y., Bian, K. & Yan, W. AMEIR: Automatic behavior modeling, interaction exploration and MLP investigation in the recommender system. Preprint at https://arxiv.org/abs/2006.05933 (2020).
  * 39 Merity, S., Xiong, C., Bradbury, J. & Socher, R. Pointer sentinel mixture models. Preprint at https://arxiv.org/abs/1609.07843 (2016)