11institutetext: 1 Hankyong University 2 Universität Trier 3 Aston University
4 KAIST 5 Accelia Inc.

# On Formal Verification in Imperative Multivalued Programming over Continuous
Data Types††thanks: Based on preprint arXiv:1608.05787v1, we thank Franz
Brauße, Cyril Cohen, Pieter Collins, Jeehoon Kang, Johannes Kanig, Sunyoung
Kim, and Michal Konečný for seminal discussions — and the referees for helpful
feedback. This work was supported by the International Research & Development
Program of the National Research Foundation of Korea (NRF) funded by the
Ministry of Science and ICT (grant NRF-2016K1A3A7A03950702) and by the
European Union’s Horizon 2020 MSCA IRSES project 731143

Gyesik Lee1 Norbert Müller2 Eike Neumann3 Sewon Park4 Norbert Preining5
Martin Ziegler4

###### Abstract

Inspired and guided by the iRRAM C++ library (Müller 2001), we formally
specify a programming language for the paradigm of _Exact Real Computation_
(ERC): reliably operating on encapsulated _continuous_ data types such as (not
necessarily algebraic) real numbers — imperatively and exactly (no rounding
errors) with primitives _computable_ in the sense of Recursive Analysis
including a necessarily modified multivalued (=non-functional) semantics of
tests. Three simple numerical problems demonstrate the elegance and
convenience of writing programs handling _real_ real numbers: integer
rounding, solving systems of linear equations, and continuous root finding. We
establish Turing-_completeness_ over the reals: a partial function
$f:\subseteq\mathbb{R}^{d}\times\mathbb{Z}\to\mathbb{R}$ is computable (in the
sense of Recursive Analysis) iff it can be expressed in ERC. For rigorously
specifying and arguing about such computations in Mathematical Logic, we then
propose a _decidable_ first-order theory over two sorts, integers and real
numbers. We extend the rules of Hoare Logic to support the formal derivation
of correctness proofs in ERC; and we have them, including their real
quantification, verified in the Coq Proof Assistant.

##  1 Motivation, Introduction, Overview

Correctness is a core of Computer Science: reliably building from single
silicon gates via circuits, microprocessor hardware, operating system kernel,
libraries, compilers/interpreters, to high-level user programs. Based on
mathematical logic, the Theory of Computation collects fundamental concepts
and methods devised to help assert correctness on each level and allow for
modern modular software engineering. These include specification, algorithm
design and analysis, models of computation, semantics of operational
primitives, measures of cost, and formal verification — for problems over
discrete structures: common continuous realms like real numbers and functions,
regularly attributed to numerics, arguably lack behind regarding a rigorous
treatment [22, p.412] of its heuristics and ‘recipes’ [32]. Although
empirically often successful, numerical computations may be errorenous — with
dramatical consequences [34, 19, 18], yet facilitated by complex [26] or vague
specifications [30, e04bbc].

_Recursive Analysis_ offers a sound foundation of numerical computation over
real numbers, functions, and compact Euclidean subsets [36]:

> A real number $x$ is computable if some Turing machine can, given
> $p\in\mathbb{Z}$, produce the numerator $a_{p}\in\mathbb{Z}$ of a dyadic
> rational $a_{p}\cdot 2^{p}$ approximating $x$ up to absolute error $2^{p}$.
> Computing a partial function
> $f:\subseteq\mathbb{R}^{d}\times\mathbb{Z}\to\mathbb{R}$ means to convert
> any integer sequence $(\vec{a}_{p})$ and $j\in\mathbb{Z}$ with
> $\|\vec{x}-\vec{a}_{p}\cdot 2^{p}\|\leq 2^{p}$ for
> $(\vec{x},j)\in\operatorname{dom}(f)$ into a sequence
> $(b_{q})\subseteq\mathbb{Z}$ with $|f(\vec{x},j)-b_{q}\cdot 2^{q}|\leq
> 2^{q}$.

This notion has pleasant properties (such as closure under composition) and
renders arithmetic as well as many transcendental functions computable.
Moreover it reflects the cost of contemporary approaches determining $\pi$ up
to billions of digits, that is, up to absolute error $\varepsilon=2^{-n}$ for
$n\approx 10^{10}$ within time polynomial in the _binary_ †††A runtime
polynomial in the _unary_ precision parameter $1/\varepsilon$ would be
meaningless here. Theorem 4.1 gives another reason for preferring binary over
unary precision. precision parameter $n$. However the underlying Turing
machine model is inconvenient to code in practice [3].

The algebraic model [35, 5] on the other hand is intuitive and implicitly
prevalent but neglects the cost of internal precision, and its test for
equality exhibits superrecursive power [6]. Reconciling and combining the best
of both worlds, and inspired by [8], the free C++ library iRRAM [28] provides
via object-oriented overloading a data type REAL for (not necessarily
algebraic) real numbers: with simultaneously exact and computable primitives
by introducing a modified — namely either partial or multivalued (aka non-
functional) — semantics of comparisons. Indeed multivaluedness is well-known
inherent to reliable real computation [23, 31].

The present work establishes a theoretical foundation to the above
implementation. Section 2 below captures the _Exact Real Computation_ (ERC)
paradigm in a formal language for imperative programming over real numbers as
abstract data type. The convenience of programming in ERC with modified
semantics of comparisons and dealing with multivaluedness/non-functionality in
loops and conditional statements is demonstrated in Section 3 with rigorous
algorithms for three numerical example problems: (I) multivalued integer
rounding, (II) Gaussian Elimination with matrix rank given, and (III) simple
single 1D root finding. Theorem 2.1 asserts that a real function is computable
in the sense of Recursive Analysis iff it can be expressed in ERC: Turing-
completeness over the reals. We then propose (Section 4) a two-sorted logical
structure for rigorously specifying and arguing about the behaviour of such
programs; and show its first-order theory decidable (Theorem 4.1). For
verifying program correctness formally, Section 5 adapts and extends the
classical Floyd-Hoare Logic to this structure. As proof (pun) of concept,
Section 6 applies said extension to the aforementioned root-finding Algorithm
(III); and we have this correctness proof, based on the inherently real
Intermediate Value Theorem, formally verified in Section 7 using the Coq Proof
Assistant.

We are convinced that the elegance and verifiability of algorithms in this
imperative programming language for real numbers (as opposed to floating
point, rational, algebraic numbers, intervals, or sequences of approximations)
will simplify rigorous numerical software development [29] and enable rigorous
software engineering in Reliable Numerics.

###  1.1 Related Work

IEEE754 is the common standard for hardware-supported (and thus extremely
fast) floating-point data types. Their fixed precision and truncation errors
cause violations of (or at least heavily complicate) common mathematical
properties such as distributivity or the Intermediate Value Theorem. The many
rounding modes and other cases (underflow, overflow, NaN) make algorithm
verification a real (pun) challenge [7]. ERC introduces a new layer of
abstraction: building and relying on correct implementation of the operational
semantics on real numbers axiomatized in Definition 1 below — similarly to
regular high-level programs building and relying on correct compilers,
libraries, operating systems, and hardware. It roughly corresponds to Level
III of the core library where, however, the classical/naïve semantics of tests
spoils closure under composition [39, p.325].

Indeed, equality of real (as opposed to, say algebraic) numbers is a well-
known equivalent to the complement of the Halting Problem [36, Exercise
4.2.9]; whereas algebraic numbers, as supported for instance by the LEDA
library [25], are insufficient for important transcendental functions like
$\exp$ or $\ln$. Multiprecision calculations, interval arithmetic, Domain
Theory, and operating on streams of approximations, are common approaches to
transcendental computation — yet arguably inconvenient due to the necessity to
keep track of the error accumulation and providing only a-posteriori bounds
[9]. Functional programming can avoid these disadvantages [17, 14], but still
lacks broad practial dissemination.

##  2 Programming Language for Exact Real Computation

“Don’t test for equality!” may be the first lesson in _Numerical Programming
101_. Indeed, inequality of real (and not just, say, algebraic) numbers is
well-known equivalent to the Halting Problem [36, Exercise 4.2.9]. But which
comparisons are permitted, then? Strict inequality “$x>0$” would allow to
express equality via the Boolean combination “$\neg(x>0)\wedge\neg(-x>0)$”.
This paradox can be avoided, while maintaining exact arithmetic on real
numbers, by replacing the naïve semantics with either a partial or a
multivalued/non-functional one.

Recall that a partial multivalued‡‡‡as opposed to the _nondeterminism_
defining, e.g., the complexity class $\mathcal{NP}$ [40, §5]. mapping
$f:\subseteq X\rightrightarrows Y$ (aka search problem) is simply a relation
$f\subseteq X\times Y$, considered as total function $f:X\ni x\mapsto\\{y\in
Y:(x,y)\in f\\}$; with the understanding that an algorithm computing $f$ may,
given $x\in\operatorname{dom}(f)=\\{x:f(x)\neq\emptyset\\}$, return any $y\in
f(x)$. Single-valuedness/functionality thus means singleton values.

######  Definition 1 (Exact Real Computation)

Consider WHILE programs [15, §1.2.3] with two separate data types: integers
$\mathbb{Z}$ for counters/Booleans (0=false, 1=true) and reals $\mathbb{R}$,
connected by the _binary precision_ embedding $\imath:\mathbb{Z}\ni p\mapsto
2^{p}\in\mathbb{R}$. Perusing C++ notation for illustration purposes, a
(possibly partial, but functional) real-valued function
$f:\subseteq\mathbb{R}^{d}\times\mathbb{Z}\to\mathbb{R}$ is declared as

$\textup{{REAL}}\>f(\textup{{INTEGER}}\>p,\;\textup{{REAL}}\>x_{1},\ldots,x_{d},\;\textup{{INTEGER}}\>j);$

with one additional integer precision parameter $p$. It receives arguments
$j\in\mathbb{Z}$ and $x_{1},\ldots,x_{d}\in\mathbb{R}$ exactly; and has to
return a real approximation $z$ to $y:=f(x_{1},\ldots,x_{d},j)$ up to absolute
error $\leq 2^{p}$: provided $(x_{1},\ldots,x_{d},j)\in\operatorname{dom}(f)$;
otherwise it may behave arbitrarily. Let +, -, *, / denote the usual
arithmetic on real numbers: exactly! Moreover,
$\mathrel{\textup{{\large{>}}}_{\\!}}$ and
$\operatorname{\textup{{choose}}}$§§§We deliberately deviate from iRRAM where
$\operatorname{\textup{{choose}}}$ has codomain $\\{1,2,\uparrow\\}$ denote
the following singlevalued _partial_ comparison and _multivalued_ partial
predicate:

$(x\mathrel{\textup{{\large{>}}}_{\\!}}y)\;=\;\left\\{\begin{array}[]{cl}1&:\;x>y,\\\
0&:\;x<y,\\\
\uparrow&:\;x=y\end{array}\right.\qquad\qquad\operatorname{\textup{{choose}}}(a,b)\;=\;\left\\{\begin{array}[]{cl}0&:\;a=\textup{{true}},\\\
1&:\;b=\textup{{true}},\\\ \uparrow&:\;\text{\rm otherwise}\end{array}\right.$
(1)

for lazy/_un_ evaluated logical expressions $a,b$, some of which may be
undefined (“$\uparrow$”) results of other partial tests/predicates.

To emphasize, the partial test “$x\mathrel{\textup{{\large{>}}}_{\\!}}y$” does
not return anything/at all in case $x=y$!
A total but multivalued so-called _soft_ test similar to [38, §6] is obtained
as

$(x\mathrel{\textup{{\large{>}}}_{\\!p}}0)\quad:=\quad\left\\{\begin{array}[]{cl}1&:\;x>-2^{p}\\\
0&:\;x<+2^{p}\end{array}\right\\}\quad\equiv\quad\operatorname{\textup{{choose}}}\big{(}\>2^{p}\mathrel{\textup{{\large{>}}}_{\\!}}x\>,\>x\mathrel{\textup{{\large{>}}}_{\\!}}-2^{p}\>\big{)}\enspace.$

The non-functional binary predicate $\operatorname{\textup{{choose}}}$ is
known in Logic as _parallel OR_ [16]. Higher, yet finite, arities can be
expressed for instance in the ternary case as follows:

$\operatorname{\textup{{choose}}}(a,b,c)\;\;:\equiv\;\;\operatorname{\textup{{choose}}}\big{(}\operatorname{\textup{{choose}}}(a,b)\geq
0,c\big{)}\;\;\texttt{?}\;\;2\;\;\texttt{:}\;\;\operatorname{\textup{{choose}}}(a,b)$

Here and in the sequel we presume this conditional operator and instructions
if then else expressed via while as usual. Note that no consistency is
required among the first and second call
$\operatorname{\textup{{choose}}}(a,b)$ in case $a\neq\uparrow\neq b$. Indeed,
$a\Leftrightarrow a^{\prime}\neq\uparrow$ and $b\Leftrightarrow
b^{\prime}\neq\uparrow$ does not imply
$\operatorname{\textup{{choose}}}(a,b)=\operatorname{\textup{{choose}}}(a^{\prime},b^{\prime})$.
Assignments are _not_ lazy: storing a logical expression in an integer
variable requires the former to be defined. Also ‘multivalued assignment’ like
$x:=\operatorname{\textup{{choose}}}(a,b,c)$ is understood as abbreviation for
the following:

$\textup{{\boldmath if
}}\operatorname{\textup{{choose}}}(a,b,c)=2\textup{{\boldmath then
}}x:=2\textup{{\boldmath else }}\textup{{\boldmath if
}}\operatorname{\textup{{choose}}}(a,b)\textup{{\boldmath then
}}x:=1\textup{{\boldmath else }}x:=0\textup{{\boldmath end if
}}\textup{{\boldmath end if }}$

Definition 1 captures numerical practice and iterative methods like Newton,
Gauss-Seidel, or Runge-Kutta. There _floating point_ inputs are (often
silently) considered exact, and the output produced by approximation — and
hope [22, p.412], for instance that cancellation is rare; whereas here the
programmer can _rely_ on her code and rigorous correctness proofs (Section 5).

ERC by definition converts the convergent sequence
$\big{(}f(p,\vec{x},j)\big{)}_{{}_{p}}$ of output approximations implicitly
and automatically¶¶¶iRRAM requires the user to explicitly select and apply one
of its limit operators. into a single exact object of type REAL, enabling
calls/assignments like $t:=f(x,y,j)$. This makes ERC Turing-complete over the
Reals:

###### Theorem 2.1

Any algorithm according to Definition 1 describes a partial real function
computable in the sense of Recursive Analysis — and vice versa: Every partial
function $f:\subseteq\mathbb{R}^{d}\times\mathbb{Z}\to\mathbb{R}$ computable
in the sense of Recursive Analysis can be implemented in ERC.

###### Proof

Regarding the first claim refer to [37]. ERC can implement integer
multiplication by repeated addition, producing integer and real results along;
similarly for powering and more generally Ackermann’s Function. In fact we can
simulate a Counter Machine, and thus any Turing Machine, on discrete input.
Regarding a real argument $x$, Algorithm (I) in Subsection 3.1 implements a
multivalued integer ‘function’
$\operatorname{round}:\mathbb{R}\rightrightarrows\mathbb{Z}$ satisfying
$\forall y\in\operatorname{round}(x):\big{|}x-y\big{|}\leq 1$. Composing with
the binary precision embedding,
$a_{n}:=\operatorname{round}\big{(}x\cdot\imath(n)\big{)}\in\mathbb{Z}$ yields
the numerators of a sequence $a_{n}/2^{n}$ of dyadic approximations to $x$ up
to absolute error $\leq 2^{-n}$: by [21, Definition 2.11] the way of
presenting real argument to a Turing machine computing $f(x)$. Similarly for
higher arities and mixed real/integer arguments…∎

##  3 Programming in Exact Real Computation

The computable but partial/multivalued semantics of tests in ERC might take
some getting used to by classically-trained programmers; for instance to
guarantee not only termination of loops, but also assert the loop condition to
be total. The present section demonstrates this and related issues at
algorithms for three numerical problems: (I) multivalued/non-functional
integer rounding
$\operatorname{round}:\mathbb{R}\;\ni\;x\;\mapsto\;\\{k\in\mathbb{Z}:|k-x|<1\\}$,
(II) Gaussian Elimination with full pivoting, and (III) simple root finding.
These examples demonstrate the practicality of the formal programming language
from Definition 1 with partial tests, multivalued predicates, restriction to
Presburger Arithmetic for integers, and binary precision embedding $\imath$.
We include comments for precondition, loop in/variant, and postcondition as
preparation for the formal Hoare triples in Section 5.

###  3.1 Integer Rounding

Like comparing real numbers, rounding down/up or to the nearest integer is a
common but uncomputable operation: functionally. The multivalued/non-
functional variant

$\operatorname{round}:\mathbb{R}\;\ni\;x\;\mapsto\;\\{k\in\mathbb{Z}:x-1<k<x+1\\}\;\subseteq\;\mathbb{Z}$
(2)

on the other hand is computable. A trivial implementation in ERC might, given
$x>0$, initialize $k:=0$ and, $\textup{{\boldmath while
}}\operatorname{\textup{{choose}}}(k\>\texttt{>}\>x-1\>,\>x\>\texttt{>}\>k)$,
increment $k:=k+1$. However this ‘unary’ approach would take a number of
iterations proportional to $x$, that is, exponential in its (output) binary
length. A faster idea might be to determine the binary expansion of $k$
bitwise, by comparing $x$ with appropriate powers of two — naïvely: over the
reals determining any digit of ‘the’ binary expansion (or one of the at most
two possible ones) is well-known uncomputable since at least Turing (1937).
Instead, the following algorithm in ERC determines some signed-digit expansion
[36, Definition 7.2.4] of the argument $x\in\mathbb{R}$:

Algorithm (I) $\textup{{INTEGER}}\;\operatorname{round}(\textup{{REAL}}\;x)$

1:$\textup{{INTEGER}}\ni b;\quad\textup{{INTEGER}}\ni
k:=0;\quad\textup{{INTEGER}}\ni l:=0;\quad\textup{{REAL}}\ni y:=x$

2:while
$\operatorname{\textup{{choose}}}\big{(}|y|\texttt{<}1\>,\>|y|\texttt{>}\tfrac{1}{2}\big{)}$
do invar= $\boldsymbol{\big{[}}x=y\cdot\imath(l)\boldsymbol{\big{]}}$

3: $l:=l+1;\;y:=y\texttt{/}2;$

4:end while$\boldsymbol{\big{[}}|y|<1\boldsymbol{\big{]}}$

5:while $l>0$ do invar=
$\boldsymbol{\big{[}}x=(y+k)\cdot\imath(l)\;\wedge\;|y|<1\boldsymbol{\big{]}}$

6: $y\>:=\>y\>\texttt{*}\>2$

7:
$b:=\operatorname{\textup{{choose}}}(\>y\texttt{<}0\>,\>-1\texttt{<}y\texttt{<}1\>,\>y\texttt{>}0\>)\>-\>1$
// most significant _signed_ binary digit $-1,0,+1$ of $y$

8: $y\>:=\>y\>-\>b$

9: $k\>:=\>k\>+\>k\>+\>b$ variant= $\boldsymbol{\big{[}}l\boldsymbol{\big{]}}$

10: $l\>:=\>l\>-\>1$

11:end while; return $k$ post=
$\boldsymbol{\big{[}}|x-k|<1\boldsymbol{\big{]}}$

We wrote $y\texttt{<}0$ for $0\texttt{>}y$. Loop invariants have been included
in order to convey partial correctness with respect to the specified
postcondition; and a loop _variant_ asserts the (here trivial) termination.
Due to multivalueness of the test, after the real while loop (lines 2 to 4)
has ended, the second argument
“$y\mathrel{\textup{{\large{>}}}_{\\!}}\tfrac{1}{2}$” _may_ still, whereas the
first “$|y|\texttt{<}1$” _must_ be true; but always at least one of both is
valid, thus guaranteeing total correctness. In the integer loop (lines 5 to
10) multivalueness ‘strikes’ only at line 7: Recall that this assignment is
meant as abbreviation. Similarly, since ERC deliberately prevents mixing real
and integer arithmetic (in order for the below Theorem 4.1 to hold), line 8 is
understood as splitting and treating separately the three possible
cases/values of $b$.

###  3.2 Gaussian Elimination

Consider the task of finding, given a singular matrix $A$, some non-zero
vector $\vec{x}$ in $\operatorname{kernel}(A)$: a natural, multivalued/non-
functional problem solved by Gaussian Elimination.

This classical algorithm has been formally verified over _decidable_ fields
[2], not over the reals with _un_ decidable equality. Here one lacks naïve
pivot tests, employed to determine termination when the outer/column loop
counter exhausts the matrix’ rank. The latter being uncomputable in the input
$A$’s entries, $r:=\operatorname{rank}(A)\in\mathbb{N}$ instead has to be
provided as additional argument [42, 41, 20]:

Algorithm (II) $\textup{{REAL}}[d]\;$ Gauss ($\textup{{REAL}}[d\times
d]\;A,\;\;\textup{{INTEGER}}\;r$)

1:$\textup{{INTEGER}}\ni i,j,k,pi,pj,\pi[d];\;\textup{{REAL}}\ni
t,\vec{x}[d],B[d\times d]$
$\textbf{pre}=\;\boldsymbol{\big{[}}\operatorname{rank}(A)\underline{=r}<d\boldsymbol{\big{]}}$

2:for $i:=0\textup{{\boldmath to }}d-1$ do // Initialization $B:=A$,
$\pi:=\operatorname{id}$

3: $\pi[i]:=i$ ; for $j:=0$ to $d-1$ do $B[i,j]:=A[i,j]$ end for

4:end for

5:for $k:=0\textup{{\boldmath to }}r-1$ do // Convert $B$ to reduced row
echelon form:

6: $pi\cdot d+pj:=$Pivot($(d-k)^{2},B[k\ldots d-1,k\ldots d-1]$) // Find
$pi,pj$ such that $B[pi,pj]\neq 0$.

7: for $j:=0$ to $d-1$ do $\textup{{\boldmath swap}}(B[k,j],B[pi,j])$ end
for// Exchange rows #$k$ and #$pi$.

8: for $i:=0$ to $d-1$ do $\textup{{\boldmath swap}}(B[i,k],B[i,pj])$ end
for// Exchange columns #$k$ and #$pj$.

9: $\textup{{\boldmath swap}}(\pi[k],\pi[pj])$

10: for $j:=k\textup{{\boldmath to }}d-1$ do // Scale row #$k$ by $1/B[k,k]$

11: $B[k,j]:=B[k,j]/B[k,k]$ // and subtract $B[i,k]$-fold from rows
#$i=k+1\ldots d-1$:

12: for $i:=k+1$ to $d-1$ do $B[i,j]:=B[i,j]-B[i,k]*B[k,j]$ end for

13: end for

14:end for

15:$\vec{x}\big{[}\pi[d]\big{]}:=1$ // Back-substitution for $\vec{x}$, taking
into account permutation $\pi$:

16:for $i:=d-2\textup{{\boldmath downto }}0$ do

17: $t:=0$ ; for $j:=d-1$ downto $i$ do
$t:=t+B[i,j]*\vec{x}\big{[}\pi[j]\big{]}$ end for

18: $\vec{x}\big{[}\pi[i]\big{]}:=-t$

19:end for; return $\vec{x}$

Algorithm (II) employs full pivoting. We emphasize that the index operations,
such as identifying two-dimensional array positions
$(i,j)\in\\{0,\ldots,d-1\\}^{2}$ with one-dimensional $i\cdot
d+j\in\\{0,\ldots,d^{2}-1\\}$, can be expressed within Presburger Arithmetic
for arbitrary but _fixed_ matrix dimension $d\in\mathbb{N}$. Also note that
pivot search in line 6 is guaranteed to succeed in that the $(d-k)\times(d-k)$
submatrix $M:=B[k\ldots d-1,k\ldots d-1]$ under consideration will indeed
contain at least one — but usually far from unique — non-zero element. Finding
the index of such a _real_ pivot is therefore an inherently multivalued
problem, solvable in ERC as follows:

Algorithm
$\textup{{INTEGER}}\;\text{Pivot}(\textup{{INTEGER}}\;m,\;\textup{{REAL}}[m]\;M$)

1:$\textup{{REAL}}\ni x:=\max\big{(}|M[0]|,\ldots,|M[m-1]|)$;
$\textup{{INTEGER}}\ni i:=0$
$\textbf{pre}=\;\boldsymbol{\big{[}}M\in\mathbb{R}^{m}\land
M\neq\vec{0}\boldsymbol{\big{]}}$

2:while
$\operatorname{\textup{{choose}}}\big{(}\;|M[i]|\mathrel{\textup{{\large{>}}}_{\\!}}x/2\;,\;x\mathrel{\textup{{\large{>}}}_{\\!}}|M[i]|\;\big{)}$
do $\textbf{invar}=\;\boldsymbol{\big{[}}0\leq i<m\land\forall j(0\leq
j<i\Rightarrow|M[j]|<x)\boldsymbol{\big{]}}$

3: $i:=i+1$ $\textbf{variant}=\;\boldsymbol{\big{[}}m-i\boldsymbol{\big{]}}$

4:end while; return $i$ $\textbf{post}=\;\boldsymbol{\big{[}}M[i]\neq
0\;\wedge\;0\leq i<m\boldsymbol{\big{]}}$

###  3.3 Continuous Root Finding

Let $f:[0,1]\to\mathbb{R}$ denote a given function, supposed to be continuous
(predicate $\operatorname{cont}(f)$ in the below formal prerequisite) with a
unique and simple root:

$\operatorname{uniq}(f,a,b)\quad:\equiv\quad\big{(}\exists!x\in[a,b]:\;\;f(x)=0\;\wedge\;f(a)\cdot
f(b)<0\big{)}\enspace.$ (3)

_Bi_ section proceeds according to the sign of $f(1/2)$; but fails in ERC in
case $1/2$ already is a root! Instead, _tri_ section tests the signs of both
$f(1/3)$ and $f(2/3)$ in parallel, knowing from the hypothesis on $f$ that at
most one of both can be zero [36, Theorem 6.3.7]:

Algorithm (III) REAL
Trisection($\textup{{INTEGER}}\;p,\;\;\textup{{REAL}}\to\textup{{REAL}}\;f$)

1:$\textup{{REAL}}\ni x:=0$; $\textup{{REAL}}\ni y:=1$
$\textbf{pre}=\;\boldsymbol{\big{[}}\operatorname{uniq}(f,0,1)\;\wedge\;\operatorname{cont}(f)\boldsymbol{\big{]}}$

2:while
$\operatorname{\textup{{choose}}}\big{(}\;\imath(p)\mathrel{\textup{{\large{>}}}_{\\!}}y-x\;,\;y-x\mathrel{\textup{{\large{>}}}_{\\!}}\imath(p-1)\;\big{)}$
do $\textbf{invar}=\;\boldsymbol{\big{[}}0\leq x<y\leq
1\land\operatorname{uniq}(f,x,y)\boldsymbol{\big{]}}$

3: if
$\operatorname{\textup{{choose}}}\big{(}\;0\mathrel{\textup{{\large{>}}}_{\\!}}f((2x+y)/3)*f(y)\;,\;0\mathrel{\textup{{\large{>}}}_{\\!}}f(x)*f((x+2y)/3)\big{)}$
$\textbf{variant}=\;\boldsymbol{\big{[}}y-x-2^{p-1}\boldsymbol{\big{]}}$

4: then  $y:=(x+2*y)/3$ else $x:=(2*x+y)/3$ end if

5:end while; return $x$
$\textbf{post}=\;\boldsymbol{\big{[}}\operatorname{uniq}(f,x,y)\;\wedge\;|y-x|\leq
2^{p}\boldsymbol{\big{]}}$

We have introduced the precondition
$P\;:\equiv\;\;\operatorname{cont}(f)\;\wedge\;\operatorname{uniq}(f,0,1)$ and
postcondition $Q\;:\equiv\;\;\operatorname{uniq}(f,x,y)\;\wedge\;|y-x|\leq
2^{p}$ for use in Section 6.

##  4 Logic of Exact Real Computation

Reflecting the two data types supported by ERC, and in order to formally
specify and describe algorithms operating on them, consider the two-sorted
structure consisting of Presburger Arithmetic $(\mathbb{Z},0,1,+,>)$ and real-
closed field $(\mathbb{R},0,1,+,\times,>)$ of characteristic 0 together with
the ‘binary precision’ embedding $\imath:\mathbb{Z}\ni p\mapsto
2^{p}\in\mathbb{R}$ in order to express the output error specification.

###### Theorem 4.1

The first-order theory of the above two-sorted structure is decidable; but not
when replacing $\imath$ with the ‘unary’ embedding $\mathbb{N}_{+}\ni n\mapsto
1/n\in\mathbb{R}$.

The basically means that every valid formula, such as loop invariants of an
algorithm, can be formally derived from a fixed recursively enumerable and
sound set of axioms: thus justifying and promising the fundamental feasibility
of the rest of this work — while WHILE programs over integers classically
suffer from Gödel’s Undecidability [10, §6].

######  Proof (Theorem 4.1)

It is well-known that, as opposed to full integer arithmetic
$(\mathbb{Z},0,1,+,\times)$, Presburger $(\mathbb{Z},0,1,+,>)$ is decidable.
By elementary equivalence, the same holds for its multiplicative version
$\big{(}2^{\mathbb{Z}},1,2,\times\big{)}$ called _Skolem Arithmetic_. By
Tarski, the first-order theory of real-closed fields of characteristic zero
$(\mathbb{R},0,1,+,\times,>)$ is also decidable. Of course this becomes
blatantly wrong when adding an integrality test, that is, the unary predicate
$\mathbb{Z}$. But a celebrated result of Lou van den Dries [13] asserts
decidability of real-closed fields expanded with a predicate for integral
powers of two $\big{(}\mathbb{R},0,1,+,\times,>,2^{\mathbb{Z}}\big{)}$; cmp.
[4]. Finally observe that any formula in the two-sorted structure
$\big{(}2^{\mathbb{Z}},1,2,\times\big{)}$ and $(\mathbb{R},0,1,+,\times,>)$
with the identity embedding $\operatorname{id}:2^{\mathbb{Z}}\to\mathbb{R}$
can be interpreted in $\big{(}\mathbb{R},0,1,+,\times,>,2^{\mathbb{Z}}\big{)}$
by restricting $\operatorname{id}:\mathbb{R}\to\mathbb{R}$ to the thus
definable domain $2^{\mathbb{Z}}$: Decidability therefore carries over. ∎

By Theorem 2.1, many transcendental functions like $\exp:x\mapsto e^{x}$ can
be realized in ERC. Still, we evade for instance Tarski’s Exponential Function
Problem: by considering specifications only in the language of the above
logic, verified algorithms cannot rely on any non-algebraic relations among
inputs.

###### Example 4.2

The exponential function is commonly computed via its Taylor expansion
$\exp(x)=\sum_{n}x^{n}/n!$ and some tail bound such as

$\forall
N:\;\Big{|}\sum\limits_{j=-p+5}^{N}x^{j}/j!\Big{|}\;\leq\;\sum\limits_{j\geq-p+5}\underbrace{\tfrac{2}{1}\cdot\tfrac{2}{2}\cdot\tfrac{2}{3}\tfrac{2}{4}}_{=2/3}\cdot\underbrace{\tfrac{2}{5}}_{\leq
1/2}\cdots\underbrace{\tfrac{2}{j}}_{\leq
1/2}\;\leq\;\sum\limits_{j\geq-p+5}2^{-j+4}\;=\;2^{p}$ (4)

whenever $x\leq 2$ and $p\leq 5$. Thus calculating and returning
$z:=\sum_{j=0}^{4-p}x^{j}/j!$ as approximation to $y:=\exp(x)$ is easily
implemented in ERC; however for the purpose of formally verifying its
correctness, due to the lack of unbounded summation, we could not express the
Cauchy Property (4) within the above logic. Instead, recall that

$\exp(x)\;\overset{n\to\infty}{\longleftarrow}\;\big{(}1+\tfrac{x}{n}\big{)}^{n}\;\leq\;\exp(x)\;\leq\;\big{(}1+\tfrac{x}{n}\big{)}^{n+1}\;\overset{n\to\infty}{\longrightarrow}\;\exp(x)$
(5)

holds for every $x\geq 0$; in which case it suggests the following
_verifiable_ algorithmic alternative:

$\displaystyle\textup{{INTEGER}}\ni n:=1;\;\textup{{REAL}}\ni
c:=1+x;\;\textup{{REAL}}\ni a:=c;\;\textup{{REAL}}\ni b:=a*c;$
$\displaystyle\textup{{\boldmath while
}}\operatorname{\textup{{choose}}}\big{(}\imath(p)\mathrel{\textup{{\large{>}}}_{\\!}}b-a\;,\;b-a\mathrel{\textup{{\large{>}}}_{\\!}}\imath(p-1)\;\big{)}$
$\displaystyle n:=n+n;\;c:=1+x/n;\;a:=c^{n};\;b:=a*c;\;\textup{{\boldmath end
while }}$

It extends to the (deliberately overlapping) case $x\leq 1$ via
$\exp(x)=\exp(1)/\exp(1-x)$. ∎

Also Theorem 4.1 applies only to algorithmic properties _expressible_ in our
logic. It obviously excludes any undecidable, such as the Totality aka
Universal Halting, problem:

###### Example 4.3

ERC can implement integer multiplication by repeated addition. Therefore, to
every $d$-variate integer polynomial $p\in\mathbb{Z}[\vec{X}]$, there exists
an algorithm $\mathcal{A}_{p}$ that, on input $\vec{x}\in\mathbb{Z}^{d}$,
evaluates $p(\vec{x})$ and stops if the result is non-zero while deliberately
looping indefinitely in case $p(\vec{x})=0$. By Davis, Robinson, and
Matiyasevich (1970) the set $\\{p:\mathcal{A}_{p}\text{
total}\\}\subseteq\mathbb{Z}^{*}$ is _not_ recursively enumerable. Therefore
there exists a $p\in\mathbb{Z}[\vec{X}]$ such that
$\forall\vec{x}:\;p(\vec{x})\neq 0$ — equivalently: totality of
$\mathcal{A}_{p}$ — is true but not provable (more precisely: not
algorithmically derivable from any fixed recursively enumerable sound axiom
system). ∎

In practice the expressive power of Theorem 4.1 seems sufficient for practical
applications, as demonstrated in the sequel.

##  5 Hoare-Style Deduction for Exact Real Computation

Having a computation or computing system fully correct is crucial in many
cases, including mission-critical or life-critical applications as well as
computer-assisted proofs [33]. Formal verification splits into several layers,
ranging from hardware to user programme [27, p.169 level 4]. Algorithms based
on ERC benefit from its mathematical elegance (compared to, say, IEEE 754). So
let us put Theorem 4.1 into practice!

(Floyd-)Hoare Logic is a well-known formal system for reasoning about partial
and total correctness of imperative programs [1]. It considers each individual
command $C$ enriched with two comments:

###### Definition 2

A _Hoare Triple_ $\\{P\\}\;C\;\\{Q\\}$ consists of a command $C$ of the
programming language under consideration, along with two formulae $P$ and $Q$
in predicate logic over the program’s state/variables such that the
_postcondition_ $Q$ holds after executing $C$ whenever the _precondition_ $P$
was met before $C$.

Partial correctness of an entire program thus decomposes into a sequence of
Hoare Triples whose conditions can be derived formally, step by step according
to purely syntactic rules. The following, for instance, capture the logics of
assignment “$x:=e$”, of “$\textup{{\boldmath if }}b\textup{{\boldmath then
}}\ldots\textup{{\boldmath else }}\ldots$”, and of “$\textup{{\boldmath while
}}b\textup{{\boldmath do }}\ldots$”, where $e,b$ denotes some expressions and
$P[e/x]$ means substituting all occurrences of the variable $x$ in the
proposition $P$ with $e$.

$\displaystyle\boldsymbol{\big{[}}Q[e/x]\boldsymbol{\big{]}}\;\texttt{x:=e}\;\boldsymbol{\big{[}}Q\boldsymbol{\big{]}}$
(6)
$\displaystyle\boldsymbol{\big{[}}P\boldsymbol{\big{]}}\;\textup{{\boldmath if
}}b\textup{{\boldmath then }}C\textup{{\boldmath else
}}A\;\boldsymbol{\big{[}}R\boldsymbol{\big{]}}\boldsymbol{\big{[}}P\land
b\boldsymbol{\big{]}}\;C\;\boldsymbol{\big{[}}R\boldsymbol{\big{]}}\qquad\boldsymbol{\big{[}}P\land\neg
b\boldsymbol{\big{]}}\;A\;\boldsymbol{\big{[}}R\boldsymbol{\big{]}}$ (7)
$\displaystyle\boldsymbol{\big{[}}I\boldsymbol{\big{]}}\;\textup{{\boldmath
while }}b\textup{{\boldmath do }}C\;\boldsymbol{\big{[}}I\land\neg
b\boldsymbol{\big{]}}\boldsymbol{\big{[}}I\land
b\land(V=N)\boldsymbol{\big{]}}\;C\;\boldsymbol{\big{[}}I\land(V<N)\boldsymbol{\big{]}},\quad
I\wedge(V\leq 0)\Rightarrow\neg b$ (8)
$\displaystyle\boldsymbol{\big{[}}P\boldsymbol{\big{]}}\;C\;\boldsymbol{\big{[}}Q\boldsymbol{\big{]}}P\Rightarrow
P^{\prime}\quad\boldsymbol{\big{[}}P^{\prime}\boldsymbol{\big{]}}\;C\;\boldsymbol{\big{[}}Q^{\prime}\boldsymbol{\big{]}}\quad
Q^{\prime}\Rightarrow Q$ (9)

We follow the convention of writing
$\boldsymbol{\big{[}}\cdot\boldsymbol{\big{]}}$ for an assertion of total
correctness. The predicate $I$ in (8) is called _loop _in_ variant_, whereas
integer expression $V$ is the _loop _var_ iant_: having some initial value
$N\in\mathbb{Z}$, required to strictly decrease in each iteration, and to
satisfy the implication $I\wedge(V\leq 0)\Rightarrow\neg b$ for asserting
termination. Inference (9) is the rule of precondition strengthening and
postcondition weakening, used to connect between Hoare triples and usually
applied implicitly.

Over integers with arithmetic, Gödel prevents (9) to admit a sound, complete,
recursively enumerable axiomatization; but for ERC, we have Theorem 4.1! On
the other hand the real realm requires adapting Hoare Logic for three reasons:

  1. i)

The expression $b$ might be undefined; in which case merely evaluating it in
the while condition of an otherwise terminating loop would spoil total
correctness.

  2. ii)

A real loop variant $V$ can strictly decrease infinitely often without ever
reaching zero.

  3. iii)

Even when defined, a multivalued predicate may violate
$\operatorname{\textup{{choose}}}(a,c)\lor\neg\operatorname{\textup{{choose}}}(a,c)=\textup{{true}}$
and thus render vacuous both hypotheses of Rule (7).

To mend (i), add the hypothesis that $b$ be defined. To fix (ii), require that
$V$ decreases in each round by at least some fixed (within the loop, not
necessarily globally) $\varepsilon>0$. And also (iii) is taken care of in the
following real variant of the latter two above deductive rules:

$\displaystyle\boldsymbol{\big{[}}P\wedge(a\vee
c)\boldsymbol{\big{]}}\;\textup{{\boldmath if
}}\operatorname{\textup{{choose}}}(a,c)\textup{{\boldmath then
}}C\textup{{\boldmath else
}}A\;\boldsymbol{\big{[}}R\boldsymbol{\big{]}}\boldsymbol{\big{[}}P\land
c\boldsymbol{\big{]}}\;C\;\boldsymbol{\big{[}}R\boldsymbol{\big{]}}\qquad\boldsymbol{\big{[}}P\land
a\boldsymbol{\big{]}}\;A\;\boldsymbol{\big{[}}R\boldsymbol{\big{]}}$ (10)
$\displaystyle\boldsymbol{\big{[}}I\land(a\vee
c)\boldsymbol{\big{]}}\;\textup{{\boldmath while
}}\operatorname{\textup{{choose}}}(a,c)\textup{{\boldmath do
}}C\;\boldsymbol{\big{[}}I\land
a\boldsymbol{\big{]}}\boldsymbol{\big{[}}I\land
c\land(V=R)\boldsymbol{\big{]}}\;C\;\boldsymbol{\big{[}}(a\vee c)\land
I\land(V\leq R-\varepsilon)\boldsymbol{\big{]}},\quad I\wedge(V\leq
0)\Rightarrow\neg c$ (11)

In terms of _predicate transformer semantics_ [12] and with respect to the
statement “$b:=\operatorname{\textup{{choose}}}(a,c)$”,
$\boldsymbol{\big{[}}c=1\wedge a\neq 1\boldsymbol{\big{]}}$ is the _weakest
precondition_ of $\boldsymbol{\big{[}}b=1\boldsymbol{\big{]}}$; see Appendix
0.A.

##  6 Formal Correctness Proof of Trisection

The present section picks up from Subsection 3.3 to illustrate formal
verification in ERC. To emphasize, our purpose here is not to actually
establish correctness of the long-known trisection method, but to demonstrate
the extended Hoare rules from Section 5 using a toy example. Since trisection
relies on the Intermediate Value Theorem, any correctness proof must make full
use of real (as opposed to, say, floating point, rational, or algebraic)
numbers: roots of computable functions can be very general.

Let us define some abbreviations such that Algorithm (III) becomes of the form
$C_{1};C_{2}$.

$\displaystyle\tilde{a}$ $\displaystyle:=\;$
$\displaystyle\imath(p)\mathrel{\textup{{\large{>}}}_{\\!}}y-x\;,\;\;\tilde{c}:=\;y-x\mathrel{\textup{{\large{>}}}_{\\!}}\imath(p-1)$
$\displaystyle a$ $\displaystyle:=\;$ $\displaystyle
f\big{(}\tfrac{2x+y}{3}\big{)}\cdot
f\big{(}b\big{)}<0\;,\;\;c:=\;f\big{(}x\big{)}\cdot
f\big{(}\tfrac{x+2y}{3}\big{)}<0$ $\displaystyle C_{1}$
$\displaystyle:\equiv\;$ x:=0; y:=1 $\displaystyle C_{2}$
$\displaystyle:\equiv\;$ $\displaystyle\textup{{\boldmath while
}}\operatorname{\textup{{choose}}}(\tilde{a},\tilde{c})\;\textup{{\boldmath do
}}C_{3}\;\textup{{\boldmath end while }}$ $\displaystyle C_{3}$
$\displaystyle:\equiv\;$ $\displaystyle\textup{{\boldmath if
}}\operatorname{\textup{{choose}}}(a,c)\textup{{\boldmath then
}}C_{5}\textup{{\boldmath else }}C_{4}\textup{{\boldmath end if }}$
$\displaystyle C_{4}$ $\displaystyle:\equiv\;$ x:=(2*x+y)/3 $\displaystyle
C_{5}$ $\displaystyle:\equiv\;$ y:=(x+2*y)/3

Furthermore, with the specified loop invariant, $I:\equiv 0\leq x<y\leq
1\land\operatorname{uniq}(f,x,y)\land\text{cont}(f)$, and loop variant,
$V:=y-x-2^{p-1}$, we define predicates $\tilde{P}:\equiv
I\wedge\tilde{c}\wedge\big{(}V=R\big{)}$ and
$\tilde{Q}:\equiv(\tilde{a}\vee\tilde{c})\wedge I\wedge(V\leq R-\varepsilon)$,
which are repectively the pre and post conditions of the premise of rule (11)
on $C_{3}$. We omit the condition $I\wedge V\leq 0\Rightarrow\neg\tilde{c}$ in
the premise, as $V\leq 0\equiv y-x-2^{p-1}\leq 0\equiv\neg\tilde{c}$ holds.

With $\varepsilon:=\tfrac{1}{3}\cdot 2^{p-1}$ (note that $\varepsilon$ is
fixed within the loop) and $R\in\mathbb{R}$ universally quantified, the proof
tree (12) constitutes the formal verification of the Algorithm (III). Note
that the inference Rule (9) is applied implicitly here, and explicitly with
formal proof in Section 7.

(6)
$\boldsymbol{\big{[}}\big{(}I\wedge(\tilde{a}\vee\tilde{c})\big{)}\big{[}0\big{/}x,1\big{/}y\big{]}\boldsymbol{\big{]}}\;C_{1}\;\boldsymbol{\big{[}}I\wedge(\tilde{a}\vee\tilde{c})\boldsymbol{\big{]}}$
(6)
$\boldsymbol{\big{[}}\tilde{Q}\big{[}\tfrac{x+2y}{3}\big{/}y\big{]}\boldsymbol{\big{]}}\;C_{5}\;\boldsymbol{\big{[}}\tilde{Q}\boldsymbol{\big{]}}$
(9) $\boldsymbol{\big{[}}\tilde{P}\wedge
c\boldsymbol{\big{]}}\;C_{5}\;\boldsymbol{\big{[}}\tilde{Q}\boldsymbol{\big{]}}$
(6)
$\boldsymbol{\big{[}}\tilde{Q}\big{[}\tfrac{2x+y}{3}\big{/}x\big{]}\boldsymbol{\big{]}}\;C_{4}\;\boldsymbol{\big{[}}\tilde{Q}\boldsymbol{\big{]}}$
(9) $\boldsymbol{\big{[}}\tilde{P}\wedge
a\boldsymbol{\big{]}}\;C_{4}\;\boldsymbol{\big{[}}\tilde{Q}\boldsymbol{\big{]}}$
(10) $\boldsymbol{\big{[}}\tilde{P}\wedge(a\vee
c)\boldsymbol{\big{]}}\;C_{3}\;\boldsymbol{\big{[}}\tilde{Q}\boldsymbol{\big{]}}$
(9)
$\boldsymbol{\big{[}}\tilde{P}\boldsymbol{\big{]}}\;C_{3}\;\boldsymbol{\big{[}}\tilde{Q}\boldsymbol{\big{]}}$
(11)
$\boldsymbol{\big{[}}I\wedge(\tilde{a}\vee\tilde{c})\boldsymbol{\big{]}}\;C_{2}\;\boldsymbol{\big{[}}I\wedge\tilde{a}\boldsymbol{\big{]}}$
(9)
$\boldsymbol{\big{[}}P\boldsymbol{\big{]}}\;C_{1};C_{2}\;\boldsymbol{\big{[}}Q\boldsymbol{\big{]}}$
(12)

In the above proof tree, each derivation is labelled with its inference rule.
Note that
$P\;:\equiv\;\;\operatorname{cont}(f)\;\wedge\;\operatorname{uniq}(f,0,1)$ and
$Q\;:\equiv\;\;\operatorname{uniq}(f,x,y)\;\wedge\;|y-x|\leq 2^{p}$.

Readers can easily verify the correctness of the proof tree by syntactic
application of the rules (6), (10), and (11). However, for the applications of
rule (9), the following implications remain to be proven:

$\displaystyle\tilde{P}\wedge a$
$\displaystyle\Rightarrow\tilde{Q}\big{[}\tfrac{2x+y}{3}\big{/}x\big{]}$ (13)
$\displaystyle\tilde{P}\wedge c$
$\displaystyle\Rightarrow\tilde{Q}\big{[}\tfrac{x+2y}{3}\big{/}y\big{]}$ (14)
$\displaystyle\tilde{P}$ $\displaystyle\Rightarrow\tilde{P}\wedge(a\vee c)$
(15) $\displaystyle P$
$\displaystyle\Rightarrow\big{(}I\wedge(\tilde{a}\vee\tilde{c})\big{)}\big{[}0\big{/}x,1\big{/}y\big{]}$
(16) $\displaystyle I\wedge\tilde{a}$ $\displaystyle\Rightarrow Q$ (17)

Their formal verification in the Coq proof assistant is described in the
following section.

##  7 Formal Verification of the Correctness Proof in Coq

We present here a formal verification of the five applications of Hoare Rule
(9) in Section 6, that is, of the implications (13), (14), (15), (16), (17).
As before, the purpose is not the assertion itself, but to complete our proof-
of-concept and example demonstrating fully formalized correctness of
transcendental computing in ERA.

To this end we employ the Coq Proof Assistant [24]. Note that Implications
(13), (14), and (15) involve a function symbol $f$ with (obviously incomplete)
first-order axioms expressing continuity $\operatorname{cont}(f)$, existence
and uniqueness of a furthermore simple root $\operatorname{uniq}(f,0,1)$.
Fortunately the Coq standard library on real numbers already provides all the
necessary materials and instruments, including the Intermediate Value Theorem.

For execution by the reader’s own Coq installation, the full source code of
our formal proof is available at http://complexity.kaist.edu/REALVERIFICATION.
The rest of this section provides some high-level explanations of our
formalization with some inserts of Coq pseudo code.

###  7.1 Binary Precision Embedding

Our formalization starts with the definition of the _binary precision_
embedding $\imath:\mathbb{Z}\ni p\mapsto 2^{p}\in\mathbb{R}$, defined in
Definition 1. It is called iota : $\mathbb{Z}\rightarrow\mathbb{R}$, a trivial
extension of the exponentiation with base 2 from $\mathbb{N}$ to $\mathbb{R}$
which is in turn defined by a primitive recursion. The extention is done based
on a case distinction on integers: zero, positive, or negative.

Fixpoint iota  nat (n : $\mathbb{N}$) : $\mathbb{R}$ :=

match n with

$|$ O $\Rightarrow$ 1

$|$ S n $\Rightarrow$ 2 $\times$ iota  nat n

end.

Definition iota (z : $\mathbb{Z}$) : $\mathbb{R}$ :=

match z with

$|$ 0%Z $\Rightarrow$ 1

$|$ Zpos n $\Rightarrow$ iota  nat (Pos.to  nat n)

$|$ Zneg n $\Rightarrow$ 1 / (iota  nat (Pos.to  nat n))

end.

Then we get the homomorphic properties such as the following two lemmas.

Lemma iota  mult : $\forall$ z1 z2 : $\mathbb{Z}$, iota (z1 + z2) = iota z1
$\times$ iota z2.

Lemma iota  div : $\forall$ z1 z2 : $\mathbb{Z}$, iota (z1 - z2) = iota z1 /
iota z2.

###  7.2 Generalized Intermediate Value Theorem

In order to prove that one can refine root intervals, it is necessary to apply
the Intermediate Value Theorem(IVT), but the version of IVT in the standard
library of Coq (see https://coq.inria.fr/library/Coq.Reals.Rsqrt_def.html#IVT)
should be generalized a little bit.

First, the existence of a unique root can be defined by exactly the same way
as described in $(\ref{e:uniq})$.

Definition uniq (f : $\mathbb{R}$ $\rightarrow$ $\mathbb{R}$) (a b :
$\mathbb{R}$ ) := ($\exists$! z : $\mathbb{R}$, a $\leq$ z $\leq$ b $\land$ f
z = 0) $\land$ f a $\times$ f b < 0.

The geneneralized form the of the IVT is defined as g  IVT. Combined with the
property of uniqueness, the lemma proves the implication (15) with deriving
contradictions; if $f\big{(}\tfrac{2x+y}{3}\big{)}\cdot f\big{(}b\big{)}\geq
0\wedge f\big{(}x\big{)}\cdot f\big{(}\tfrac{x+2y}{3}\big{)}\geq 0$, the
interval no longer contains a unique root, according to the following lemma.

Lemma g  IVT : $\forall$ (f : $\mathbb{R}$ $\rightarrow$ $\mathbb{R}$) (x y :
$\mathbb{R}$),

cont f $\rightarrow$ x < y $\rightarrow$ (f x) $\times$ (f y) < 0
$\rightarrow$ {z : $\mathbb{R}$ $|$ x $\leq$ z $\leq$ y $\land$ f z = 0}.

One can use g  IVT to prove the refinement of root intervals. The following
lemma justifies the loop invariant, thus is a key part in proving implications
(13, 14); if an interval has a unique root with a sign change, any its
sub–interval with a sign change contains a unique root.

Lemma uniq  refinement : $\forall$ (f : $\mathbb{R}$ $\rightarrow$
$\mathbb{R}$) (x a b y : $\mathbb{R}$),

cont f $\rightarrow$ uniq f x y $\rightarrow$ a < b $\rightarrow$ x $\leq$ a
$\rightarrow$ b $\leq$ y $\rightarrow$ (f a) $\times$ (f b) < 0 $\rightarrow$
uniq f a b.

###  7.3 Formal Proofs of Implications in Verification

Now we can formally prove all the five implications (13 – 17). In the rest of
this section, we assume that a function f : $\mathbb{R}$ $\rightarrow$
$\mathbb{R}$, three real numbers a, b, r, and an integer p are arbitrarily
given. We also use $\varepsilon$ = iota (p-1) / 3.

Lemma implication  13 :

f ((2 $\times$ a + b) / 3) $\times$ (f b) < 0 $\land$ 0 $\leq$ a $\land$ a < b
$\land$ b $\leq$ 1 $\land$ uniq f a b $\land$ cont f

$\land$ b - a > iota (p - 1) $\land$ b - a - iota (p - 1) = r

$\rightarrow$ (b - (2 $\times$ a + b) / 3 < iota p $\lor$ b - (2 $\times$ a +
b) /3 > iota (p - 1))

$\land$ 0 $\leq$ (2 $\times$ a + b) / 3 $\land$ (2 $\times$ a + b) / 3 < b
$\land$ b $\leq$ 1 $\land$ uniq f ((2 $\times$ a + b) / 3) b

$\land$ cont f $\land$ b - ((2 $\times$ a + b) / 3) - iota (p - 1) $\leq$ r -
$\varepsilon$.

Lemma implication  14 :

f a $\times$ f ((a + 2 $\times$ b) / 3) < 0 $\land$ 0 $\leq$ a $\land$ a < b
$\land$ b $\leq$ 1 $\land$ uniq f a b $\land$ cont f

$\land$ b - a > iota (p - 1) $\land$ b - a - iota (p - 1) = r

$\rightarrow$ (((a + 2 $\times$ b) / 3) - a < iota p $\lor$ ((a + 2 $\times$
b) / 3) - a > iota (p - 1))

$\land$ 0 $\leq$ a $\land$ a < ((a + 2 $\times$ b) / 3) $\land$ ((a + 2
$\times$ b) / 3) $\leq$ 1

$\land$ uniq f a ((a + 2 $\times$ b) / 3) $\land$ cont f

$\land$ ((a + 2 $\times$ b) / 3) - a - iota (p - 1) $\leq$ r - $\varepsilon$.

Lemma implication  15 :

0 $\leq$ a $\land$ a < b $\land$ b $\leq$ 1 $\land$ uniq f a b $\land$ cont f

$\land$ b - a > iota (p - 1) $\land$ b - a - iota (p - 1) = r

$\rightarrow$ ((f ((2 $\times$ a + b) / 3)) $\times$ (f b) < 0 $\lor$ (f a)
$\times$ (f ((a + 2 $\times$ b) / 3)) < 0)

$\land$ 0 $\leq$ a $\land$ a < b $\land$ b $\leq$ 1 $\land$ uniq f a b $\land$
cont f

$\land$ b - a > iota (p - 1) $\land$ b - a - iota (p - 1) = r.

Lemma implication  16 :

0 $\leq$ 0 $\land$ 0 < 1 $\land$ 1 $\leq$ 1 $\land$ uniq f 0 1 $\land$ cont f

$\rightarrow$ 0 $\leq$ 0 $\land$ 0 < 1 $\land$ 1 $\leq$ 1 $\land$ uniq f 0 1
$\land$ cont f $\land$ (1 - 0 < iota p $\lor$ 1 - 0 > iota (p - 1)).

Lemma implication  17 :

b - a < iota p $\land$ 0 $\leq$ a $\land$ a < b $\land$ b $\leq$ 1 $\land$
uniq f a b $\land$ cont f

$\rightarrow$ 0 $\leq$ a $\land$ a < b $\land$ b $\leq$ 1 $\land$ (($\exists$!
x, a $\leq$ x $\leq$ b $\land$ f x = 0) $\land$ Rabs (b-a) $\leq$ iota p).

##  8 Perspective

Technically, integer rounding (Subsection 3.1) and Gaussian Elimination
(Subsection 3.1) both consider multivalued problems. They refer to the
following subtle semantical refinement of Definition 1:

######  Definition 3 (Sketch)

Consider a partial _multivalued_ vector ‘function’
$f:\subseteq\mathbb{R}^{\omega}\times\mathbb{Z}\to\mathbb{R}^{d}\times\mathbb{Z}$.
An algorithm $\mathcal{A}$ in Exact Real Computation _computing_ the
evaluation mapping $(\bar{x},j)\mapsto f(\bar{x},j)$ receives, in addition to
argument $j$, random access to the components $x[i]\in\mathbb{R}$, exactly, as
well as a binary precision requirement $p\in\mathbb{Z}$. $\mathcal{A}$ then
has has to return _some_ $(\vec{z},k)\in\mathbb{R}^{d}\times\mathbb{Z}$ such
that it holds

$\exists(\vec{y},k)\in f(\bar{x},j):\quad\|\vec{z}-\vec{y}\|_{\infty}\leq
2^{p}\enspace.$ (18)

Moreover $\mathcal{A}$ must satisfy the following _consistency_ condition:
When called twice independently with the same arguments $(\bar{x},j)$, once
with $p$ and once with increased precision $p-1$, whenever those calls to
$\operatorname{\textup{{choose}}}()$ that $\mathcal{A}$ makes in both cases do
return the same multivalued choices, then the values $\vec{z}$ and
$\vec{z}^{\prime}$ returned by $\mathcal{A}$ must in both cases approximate
the _same_ $\vec{y}$ up to errors $2^{p}$ and $2^{p-1}$, respectively.

Our proof of Theorem 2.1 immediately generalizes to multivalued functions
accordingly.

Note that, as opposed to the reals, for integer arguments and values if
suffices to consider the one-dimensional case, since tuples and sequences can
be considered encoded using a pairing function. We also remark that
identifying a sequence-valued mapping $f:\subseteq X\to\mathbb{R}^{\omega}$
with $g:\subseteq X\times\mathbb{N}\to\mathbb{R}$ is possible only under
functionality hypotheses.

Future work will extend Definitions 1+3 to formalize the computation of
function _al_ s (like trisection) as higher type objects receiving real
functions as arguments; and generalize Exact Real Computation from
$\mathbb{R}$ to other structures.

## References

  * [1] Apt, K.R.: Ten years of hoare’s logic: A survey—part i. ACM Transactions on Programming Languages and Systems (TOPLAS) 3(4), 431–483 (1981)
  * [2] Aransay, J., Divason, J.: Formalisation in higher-order logic and code generation to functional languages of the gauss-jordan algorithm. Journal of Functional Programming 25 (2015)
  * [3] Arnold Schönhage, Andreas F. W. Grotefeld, E.V.: Fast algorithms. B.I. Wissenschaftsverlag (1994)
  * [4] Avigad, J., Yin, Y.: Quantifier elimination for the reals with a predicate for the powers of two. Theor. Comput. Sci. 370(1-3), 48–59 (2007)
  * [5] Blum, L., Cucker, F., Shub, M., Smale, S.: Complexity and real computation. Springer (1998), http://www.worldcat.org/oclc/37004484
  * [6] Boldi, P., Vigna, S.: Equality is a jump. Theoretical computer science 219(1-2), 49–64 (1999)
  * [7] Boldo, S., Jourdan, J.H., Leroy, X., Melquiond, G.: Verified compilation of floating-point computations. Journal of Automated Reasoning 54(2), 135–163 (2015)
  * [8] Brattka, V., Hertling, P.: Feasible real random access machines. In: International Conference on Current Trends in Theory and Practice of Computer Science. pp. 335–342. Springer (1996)
  * [9] Collins, P., Niqui, M., Revol, N.: A validated real function calculus. Mathematics in Computer Science 5(4), 437–467 (2011)
  * [10] Cook, S.A.: Soundness and completeness of an axiom system for program verification. SIAM Journal on Computing 7(1), 70–90 (1978)
  * [11] Cuoq, P., Kirchner, F., Kosmatov, N., Prevosto, V., Signoles, J., Yakobowski, B.: Frama-c. In: International Conference on Software Engineering and Formal Methods. pp. 233–247. Springer (2012)
  * [12] Dijkstra, E.W.: Guarded commands, nondeterminacy and formal derivation of programs. Communications of the ACM 18(8), 453–457 (1975)
  * [13] Dries, L.v.d.: The field of reals with a predicate for the powers of two. Manuscripta mathematica 54, 187–196 (1986), http://eudml.org/doc/155108
  * [14] Duracz, J., Farjudian, A., Konecný, M., Taha, W.: Function interval arithmetic. In: Mathematical Software - ICMS 2014 - 4th International Congress, Seoul, South Korea, August 5-9, 2014. Proceedings. pp. 677–684 (2014)
  * [15] Enderton, H.B.: Computability Theory. Academic Press (2011)
  * [16] Escardó, M.H., Hofmann, M., Streicher, T.: On the non-sequential nature of the interval-domain model of real-number computation. Mathem. Structures in Computer Science 14(6), 803–814 (2004)
  * [17] Escardó, M.H., Simpson, A.: Abstract datatypes for real numbers in type theory. In: Rewriting and Typed Lambda Calculi - Joint International Conference, RTA-TLCA 2014, Held as Part of the Vienna Summer of Logic, VSL 2014, Vienna, Austria, July 14-17, 2014. Proceedings. pp. 208–223 (2014)
  * [18] Holand, I.: The loss of the sleipner condeep platform. In: Kusters, G.M., Hendriks, M.A. (eds.) Proc. 1st International DIANA Conference on Computational Mechanics. Springer, Dordrecht (1994)
  * [19] Jézéquel, J.M., Meyer, B.: Design by contract: The lessons of ariane. Computer 30(1), 129–130 (January 1997)
  * [20] Kihara, T., Pauly, A.: Dividing by zero - how bad is it, really? In: 41st International Symposium on Mathematical Foundations of Computer Science, MFCS 2016, August 22-26, 2016 - Kraków, Poland. pp. 58:1–58:14 (2016)
  * [21] Ko, K.I.: Computational complexity of real functions. In: Complexity Theory of Real Functions, pp. 40–70. Springer (1991)
  * [22] Linz, P.: A critique of numerical analysis. Bulletin of the American Mathematical Society 19(2), 407–416 (1988)
  * [23] Luckhardt, H.: A fundamental effect in computations on real numbers. Theoretical Computer Science 5(3), 321–324 (1977)
  * [24] The Coq development team: The Coq proof assistant reference manual (2016), http://coq.inria.fr, Version 8.61
  * [25] Mehlhorn, K., Schirra, S.: Exact computation with leda  real \- theory and geometric applications. In: Symbolic Algebraic Methods and Verification Methods, pp. 163–172. Springer (2001)
  * [26] Monniaux, D.: The pitfalls of verifying floating-point computations. ACM Transactions on Programming Languages and Systems 30(3) (May 2008)
  * [27] Müller, N., Uhrhan, C.: Some steps into verification of exact real arithmetic. NASA Formal Methods pp. 168–173 (2012)
  * [28] Müller, N.T.: The iRRAM: Exact arithmetic in C++. In: Blanck, J., Brattka, V., Hertling, P. (eds.) Computability and Complexity in Analysis. Lecture Notes in Computer Science, vol. 2064, pp. 222–252. Springer, Berlin (2001), 4th International Workshop, CCA 2000, Swansea, UK, September 2000
  * [29] Müller, N.T., Ziegler, M.: From calculus to algorithms without errors. In: Proc. 4th International Congress on Mathematical Software (ICMS). LNCS, vol. 8592, pp. 718–724. Springer (2014)
  * [30] (NAG), T.N.A.G.: The nag library. Oxford, United Kingdom
  * [31] Pauly, A., Ziegler, M.: Relative computability and uniform continuity of relations. Journal of Logic and Analysis 5(7), 1–39 (2013)
  * [32] Press, W.H.: Numerical recipes: The art of scientific computing. Cambridge university press (2007)
  * [33] Rump, S.M.: Computer-assisted proofs i. Bulletin of the Japan Society for Industrial and Applied Mathematics 14(3), 2–11 (2004)
  * [34] Skeel, R.: Roundoff error and the patriot missile. SIAM News 25(4), 11 (July 1992)
  * [35] Tucker, J.V., Zucker, J.I.: Computable functions and semicomputable sets on many-sorted algebras. In: Handbook of logic in computer science. pp. 397–525. Oxford University Press (2001)
  * [36] Weihrauch, K.: Computable Analysis: An Introduction. EATCS Texts, Springer (2000)
  * [37] Weihrauch, K.: The computable multi-functions on multi-represented sets are closed under programming. J. UCS 14(6), 801–844 (2008)
  * [38] Yap, C., Sagraloff, M., Sharma, V.: Analytic root clustering: A complete algorithm using soft zero tests. In: Conference on Computability in Europe. pp. 434–444. Springer (2013)
  * [39] Yap, C.K.: On guaranteed accuracy computation. In: Geometric Computation, pp. 322–373. World Scientific (2004)
  * [40] Ziegler, M.: Revising type-2 computation and degrees of discontinuity. Electronic Notes in Theoretical Computer Science 167, 255–274 (2007)
  * [41] Ziegler, M.: Real computation with least discrete advice: A complexity theory of nonuniform computability with applications to effective linear algebra. Annals of Pure and Applied Logic 163(8), 1108–1139 (2012)
  * [42] Ziegler, M., Brattka, V.: Computability in linear algebra. Theor.Comp.Sci. 326(1-3), 187–211 (2004)

##  Appendix 0.A Weakest Precondition Calculus of Trisection

Weakest–precondition calculus, suggested by Dijkstra, is an efficient strategy
of applying Hoare logic’s inference rules, implemented and widely adapted in
modern deductive verification systems; e.g., see [11]. With
$\boldsymbol{\big{[}}P\boldsymbol{\big{]}}C\boldsymbol{\big{[}}Q\boldsymbol{\big{]}}$,
formal verification becomes proving $P\Rightarrow wp(C,Q)$, where $wp(C,Q)$
can be obtained recursively: $wp(C_{1};C_{2},Q)=wp(C_{1},wp(C_{2},Q))$. In the
appendix section, we demonstrate the weakest precondition calculus on the
examples presented in section 3.

###### Definition 4

For a _if_ command with $\operatorname{\textup{{choose}}}$, the reformulation
follows directly from [12]:

$\displaystyle wp(\textup{{\boldmath if
}}\operatorname{\textup{{choose}}}(x_{2}>y_{2},x_{1}>y_{1})\textup{{\boldmath
then }}C_{1}\textup{{\boldmath else }}C_{2},R)$ $\displaystyle=x_{1}>y_{1}\lor
x_{2}>y_{2}$ $\displaystyle\land(x_{1}>y_{1}\Rightarrow wp(C_{1},R)$
$\displaystyle\land x_{2}>y_{2}\Rightarrow wp(C_{2},R))$

The following is a weakest precondition formulation for the while  command,
with $\operatorname{\textup{{choose}}}$, whose loop variant is real where
itself and its decrement has some lower bounds, $c$ and $c_{0}$:

$\displaystyle wp(\textup{{\boldmath while
}}\operatorname{\textup{{choose}}}(x_{2}>y_{2},x_{1}>y_{1})\textup{{\boldmath
do }}C,R)$ $\displaystyle=$ $\displaystyle I$ (19) $\displaystyle\land$
$\displaystyle\exists_{c}\forall_{\bm{a},z}\big{(}I\land x_{1}>y_{1}\land
V=z\Rightarrow wp(I\land V\leq z-c,C)\big{)}$ (20) $\displaystyle\land$
$\displaystyle\exists_{c_{0}}\forall_{\bm{a}}\big{(}I\land V\leq
c_{0}\Rightarrow x_{1}\leq y_{1}\big{)}$ (21) $\displaystyle\land$
$\displaystyle\forall_{\bm{a}}\big{(}(I\land x_{2}>y_{2})\Rightarrow R\big{)}$
(22) $\displaystyle\land$
$\displaystyle\forall_{\bm{a}}\big{(}I\Rightarrow(x_{1}>y_{1}\lor
x_{2}>y_{2})\big{)}$ (23)

$\forall_{\bm{a}}$ is a universal quantifier over all variables whose values
get manipulated in the loop body, $C$. Briefly explaining a meaning of each
line of the precondition, (19) states that the loop invariant $I$ should be
satisfied in the beginning of the while command; (20) is for proving that $I$
and $V$ are valid; (21) is a condition which ensures that $V$ is bounded lower
by $c_{0}$; (22) claims that when the loop terminates, the postcondition
should be satisfied; and (23) ensure that at the every entrance of the loop
body, $\operatorname{\textup{{choose}}}$ should be safely defined.

Using a convention suggested by Isabelle theorem prover, we list a
verification condition $P\Rightarrow(wp_{1}\land\cdots\land wp_{n})$, where
$wp_{1}\land\cdots\land wp_{n}=wp(C,Q)$, $P,Q$ are pre- and postconditions and
$C$ is the program we will verify, as follow:

$\displaystyle\llbracket P\rrbracket$ $\displaystyle 1.\;wp_{1}$
$\displaystyle\vdots$ $\displaystyle n.\;wp_{n}$

Recall the precise description of the trisection algorithm in the section 3.3.
With $f\in\mathbb{R}\to\mathbb{R},p\in\mathbb{Z}$ and bounded variables
$x,y,c,c_{0},z\in\mathbb{R}$, Applying the syntactic calculus yields the
following verification condition which is needed to be proven;

$\displaystyle\llbracket\text{un}$
$\displaystyle\text{iqe}(f,0,1)\land\operatorname{cont}(f)\rrbracket$
$\displaystyle 1.\;$ $\displaystyle 0\leq 0<1\leq
1\land\operatorname{uniq}(f,0,1)$ $\displaystyle 2.\;$
$\displaystyle\exists_{c}\forall_{x,y,z}\big{[}\big{(}0\leq x<y\leq
1\land\operatorname{uniq}(f,x,y)\land y-x>2^{p-1}\land
y-x=z\big{)}\Rightarrow$ $\displaystyle\big{(}f(x)f(\frac{x+2y}{3})<0\lor
f(\frac{2x+y}{3})f(y)<0\big{)}$
$\displaystyle\land\big{(}f(x)f(\frac{x+2y}{3})<0\Rightarrow 0\leq
x<\frac{x+2y}{3}\leq
1\land\operatorname{uniq}(f,x,\frac{x+2y}{3})\land\frac{x+2y}{3}-x\leq
z-c\big{)}$ $\displaystyle\land\big{(}f(\frac{2x+y}{3})f(y)<0\Rightarrow
0\leq\frac{2x+y}{3}<y\leq 1\land\operatorname{uniq}(f,\frac{2x+y}{3},b)\land
y-\frac{2x+y}{3}\leq z-c\big{)}\big{]}$ $\displaystyle 3.\;$
$\displaystyle\exists_{c_{0}}\forall_{x,y}\big{[}\big{(}0\leq x<y\leq
1\land\operatorname{uniq}(f,x,y)\land y-x\leq c_{0}\big{)}\Rightarrow y-x\leq
2^{p-1}\big{]}$ $\displaystyle 4.\;$
$\displaystyle\forall_{x,y}\big{[}\big{(}0\leq x<y\leq
1\land\operatorname{uniq}(f,x,y)\land
y-x<2^{p}\big{)}\Rightarrow\exists!x^{\prime}\in[0,1]:f(x^{\prime})=0\;\wedge\;|x^{\prime}-x|\leq
2^{p}\big{]}$ $\displaystyle 5.\;$
$\displaystyle\forall_{x,y}\big{[}\big{(}0\leq x<y\leq
1\land\operatorname{uniq}(f,x,y)\big{)}\Rightarrow(y-x>2^{p-1}\lor
y-x<2^{p})\big{]}$

As it was shown in Section 7, the formal proofs of the above implication can
be done in a Coq Proof Assistant.


# Recent $J/\psi$ results measured with PHENIX

T. Novák (for the PHENIX Collaboration)  MATE KRC, H-3200 Gyöngyös, Mátrai út
36, Hungary
BGE, H-1054 Budapest, Alkotmány út 9-11, Hungary

###### Abstract

Finalized PHENIX results are presented on $J/\psi$ modification from the RHIC
system-size scan, which include $p$+Al, $p$+Au, and 3He+Au (small systems)
collisions as well as $AA$ (large systems) collisions at forward and backward
rapidity ($1.2<|y|<2.2$). Comparisons with state-of-the-art theory
calculations and LHC data will be discussed. The forward rapidity data in
small systems are well described by model calculations for all three systems,
however, the backward rapidity $p$+Au and 3He+Au data are not reproduced by
the shadowing calculation alone. $J/\psi$ $v_{2}$ meausrements in Au+Au
collisions at forward rapidity is consistent with zero.

##  1 Introduction

Based on the time scale in which the nuclear collisions occur, the effects
that modify the $J/\psi$ production probability can be separated into two
categories. Those that occur during the expansion and hadronization of the
energy produced in the collision are generally referred to as hot matter
effects. Those that occur on a time scale small than the crossing time of the
colliding nuclei are frequently referred to as cold nuclear matter (CNM)
effects. The latter include modification of the nuclear-parton-distribution
functions (nPDFs) in a nucleus [1], [2], initial state parton energy loss [3],
coherent gluon saturation [4], [5], breakup of the forming charmonium in
collisions with target nucleons [6], [7], and transverse momentum broadening
[8]. Here is a short list with definition of CNM effects:

Gluon Shadowing/Anit-Shadowing:

Modification (suppression/enhancement) of heavy quark crosssection due to
modifications of the gluon structure function

Parton Energy Loss:

The projectile gluon experiences multiple scattering while passing through the
target before $J/\psi$ production, reducing the rapidity of the $J/\psi$

Cronin Effect:

Modification of the $J/\psi$ $p_{\mathrm{T}}$ distribution due to multiple
elastic scattering of partons

Nuclear Break-Up:

The break up of the bound $J/\psi$ (or precursor state) in collisions with
other target nucleons that pass through $J/\psi$ production point

Co-Movers Break-Up:

Final state break up of the $J/\psi$ through interactions with producedpartons

To explore the effects of energy production on $J/\psi$ yields in light
systems, the PHENIX collaboration has recently measured the modification in
three collision systems involving small projectiles, $p$+Al, $p$+Au and 3He+Au
[9]. The $p$+Au and 3He+Au data, combined with the earlier PHENIX $d$+Au data
[10], [11], provide a direct comparison of inclusive charmonium production for
projectiles with one, two and three nucleons on a Au target.

##  2 Experimental Setup

PHENIX was planned to operate for 10 years, but finally it was running between
2000 and 2016. During the last few years of operation, PHENIX collected data
in a wide range of different collision systems.

The PHENIX Muon Arms are located parallel to the beam pipe, and measure muons
and unidentified charged hadrons at forward ($1.2<\eta<2.4$) and backward
($-2.2<\eta<-1.2$) rapidity. The Muon Arms compromise four main components:
the FVTX Detector, the Muon Tracker, the Muon Identifier, and a series of
steel hadron absorbers located throughout each Muon Arm [12], [13]. The
measurements were attainable using the Forward Vertex Silicon Detector (FVTX),
installed in the PHENIX 2012 upgrade, which covers the pseudorapidity range
$1.2<|\eta|<2.2$. The FVTX detector, a precision silicon tracking detector,
provides additional space points closest to the interaction region, supplying
the mass resolution necessary to extract the $\psi$(2S) signal.

##  3 Results

There were significant changes to the PHENIX muon arm configuration and to the
PHENIX simulation framework between the time that the d+Au data set was
recorded in 2008 and when the $p$+Al, $p$+Au, 3He+Au and their reference
$p$+$p$ data were recorded in 2014 and 2015. In addition, the $d$+Au data set
used a different reference $p$+$p$ measurement. Therefore the relative
systematic uncertainties between the modifications measured for $p$+Al, $p$+Au
and 3He+Au are smaller.

###  3.1 Small Systems

Figure 1 shows the measured rapidity dependence of the modifiation for $p$+Al,
$p$+Au and 3He+Au, integrated over $p_{\mathrm{T}}$ and for 0-100$\%$
centrality. The data are compared with model calculations by Shao et al [14],
[15], [16], [17] that include only CNM effects.

The forward rapidity data are well described for all three systems, as are the
backward rapidity $p$+Al data. However, the backward rapidity $p$+Au and
3He+Au data are not reproduced by the shadowing calculation alone. At RHIC
energies a significant suppression is expected at backward rapidity due to
breakup of the forming charmonium in collisions with nucleons that have not
yet passed through the production point (commonly referred to as “nuclear
absorption”). An estimate of this effect has been folded into the shadowing
calculation and is included in Figure 1. Including the absorption estimate
brings the backward rapidity calculation into reasonable agreement with the
$p$+Au and 3He+Au data.

Figure 1: Nuclear modification factor of inclusive $J/\psi$ as a function of
rapidity for $0-100\%$ $p$+Al (a), $p$+Au (b), and 3He+Au (c) collisions. Bars
(boxes) around data points represent point-to-point uncorrelated (correlated)
uncertainties. The theory bands are discussed in the text.

Figure 2 shows the $J/\psi$ and $\psi(2S)$ nuclear modification measurements
as a function of $\langle N_{coll}\rangle$. The nuclear modification between
the two states follows a similar trend at forward rapidity, with no clear
difference in suppression in the most central collisions. The $\psi(2S)$
shadowing predictions provided by Shao et al. shown in Figure 2(b)
underpredict the suppression at forward rapidity. Also shown is a comparison
to transport-model (TM) predictions for $\psi(2S)$ $J/\psi$ provided by Du and
Rapp [18].

A difference can be seen at backward rapidity between the $J/\psi$ and
$\psi(2S)$ nuclear modification, consistent with a 2.9$\sigma$ effect. The
$\psi(2S)$ antishadowing predictions provided by Shao et al. shown in Figure
2(a) do not predict the suppression. These predictions are purely
antishadowing and do not contain any additional CNM effects, such as nuclear
absorption [19].

Figure 2: Nuclear modification factor of inclusive $J/\psi$ as a function of
$\langle N_{coll}\rangle$ for $p$+Au collisions compared with the transport
model. Bars (boxes) around data points represent point-to-point uncorrelated
(correlated) uncertainties.

Figure 2(a) also shows a comparison to TM predictions at backward rapidity,
and final-state effects and nuclear absorption are expected to be important.
The TM includes a nuclear-absorption estimate for both the $\psi(2S)$ and
$J/\psi$. The TMs underpredict the suppression, but describe the relative
modification well and indicate that the $\psi(2S)$ suppression in $p$+Au
collisions is consistent with final state effects.

###  3.2 Large Systems

We present a preliminary result for $J/\psi$ $v_{2}$ using the PHENIX Run14
Au+Au dataset at $\sqrt{s_{NN}}=200$ GeV [20]. Figure 3 shows the
$p_{\mathrm{T}}$-dependent $J/\psi$ $v_{2}$. The measurement in this analysis
for PHENIX Run 14 at forward rapidity in a centrality range of 10 \- 60$\%$ is
shown in red. The measurement made by STAR at mid-rapidity and in a centrality
range of 10-40$\%$ is shown in black. The ALICE result at forward rapidity in
a centrality range of 20-40$\%$ is shown in blue. Boxes surrounding the data
points represent systematic uncertainties.

At LHC energies, a nonzero $v_{2}$ is observed, which is in line with $J/\psi$
formed by coalescence in the QGP medium, and carrying the azimuthal anisotropy
of the system [21]. At RHIC energies, STAR has measured $v_{2}$ that is
consistent with zero, but due to limited statistics remains inconclusive [22].
With coalescence being the dominant mechanism for nonzero $J/\psi$ $v_{2}$ it
should follow that systems where fewer $c\bar{c}$ pairs are formed should have
a smaller azimuthal anisotropy.

From the figure we can see the clear nonzero $v_{2}$ measured by ALICE.
Although the ALICE measurement is at a much higher energy, we know $v_{2}$
does not scale with energy for $J/\psi$, so it makes for a good comparison
that the ALICE result which is clearly nonzero is different from our
measurement. In our measurement, we see a $v_{2}$ that is clearly consistent
with zero across all $p_{\mathrm{T}}$ bins. The systematic uncertainties were
conservatively estimated, not taking into account cancellations or
correlations of uncertainties from different sources. Additional data from Run
16 of RHIC will be included in the final results, and we expect that both
statistical and systematic uncertainties will be significantly reduced.

Figure 3: Plot of $p_{\mathrm{T}}$ dependent $J/\psi$ $v_{2}$. The PHENIX
result in light gray/red/circle is compared to STAR in black/star and ALICE
gray/blue/square.

Partially supported by agencies that support the PHENIX collaboration
(https://www.bnl.gov/rhic/funding.php), as well as by the NKFIH grant K133046
(Hungary).

## References

## References

  * [1] K. J. Eskola, P. Paakkinen, H. Paukkunen, and C. A. Salgado, Eur. Phys. C 77, 3 (2017).
  * [2] K. Kovarik, et al, Phys. Rev. D 93, 8 (2016).
  * [3] I. Vitev, Phys. Rev. C 75, 3 (2007).
  * [4] D. Kharzeev and K. Tuchin, Nucl. Phys. A 770, 40 (2006).
  * [5] H. Fujii, F. Gelis, and R. Venugopalan, Nucl. Phys. A 780, 146 (2006).
  * [6] D. C. McGlinchey, A. D. Frawley and R. Vogt, Phys. Rev. C 87, 5 (2013).
  * [7] F. Arleo, P. B. Gossiaux, T. Gousset, and J. Aichelin, Phys. Rev. C 61, 5 (2000).
  * [8] J. W. Cronin, H. J. Frisch, M. J. Shochet, J. P. Boymond, R. Mermod, P. A. Piroue, and R. L. Sumner, Phys. Rev. D 11, 3105 (1975).
  * [9] U. A. Acharya et al, Phys. Rev. C 102, 014902 (2020).
  * [10] A. Adare et al, Phys. Rev. Lett. 107, 142301 (2011).
  * [11] A. Adare et al, Phys. Rev. C 87, 3 (2013).
  * [12] H. Akikawa et al, Nucl. Instrum. Methods Phys. Res. A 499, 537 (2003).
  * [13] C. Aidala et al, Nucl. Instrum. Methods Phys. Res. A 755, 44 (2014).
  * [14] A. Kusina, J.-P. Lansberg, I. Schienbein, and H.-S. Shao, Phys. Rev. Lett. 121, 5 (2018).
  * [15] H.-S. Shao, Comput. Phys. Commun., 184, (2013).
  * [16] H.-S. Shao, Comput. Phys. Commun., 198, (2016).
  * [17] J.-P. Lansberg and H.-S. Shao, Eur. Phys. C 77, 1 (2017).
  * [18] X. Du and R. Rapp, J. High Energy Physics, 03, (2019).
  * [19] K. Smith, Acta Physica Polonica B, 16, (2023).
  * [20] L. Bichon, 10.48550/arXiv.2301.04186. (2023).
  * [21] S. Acharya et al, J. High Energy Physics, 10, (2020).
  * [22] L. Adamczyk et al, Phys. Rev. Lett. 111, 5 (2013).