# A framework to test interval arithmetic libraries and their IEEE 1788-2015
compliance

Luis Benet1, Luca Ferranti2, Nathalie Revol3
$[1]$: Instituto de Ciencias Físicas, Universidad Nacional Autónoma de México,
Cuernavaca, Mexico
$[2]$: University of Vaasa, Vaasa, Finland
$[3]$: INRIA - LIP (UMR 5668, ENS Lyon, University Lyon 1, Inria, CNRS), Lyon,
France

(Corresponding author: Nathalie Revol, INRIA - LIP, ENS de Lyon, 69364 Lyon
Cedex 07, France<EMAIL_ADDRESS>)

###### Abstract

As developers of libraries implementing interval arithmetic, we faced the same
difficulties when it comes to testing our libraries. What must be tested? How
can we devise relevant test cases for unit testing? How can we ensure a high
(and possibly 100%) test coverage? Before considering these questions, we
briefly recall the main features of interval arithmetic and of the IEEE
1788-2015 standard for interval arithmetic. After listing the different
aspects that, in our opinion, must be tested, we contribute a first step
towards offering a test suite for an interval arithmetic library. First we
define a format that enables the exchange of test cases, so that they can be
read and tried easily. Then we offer a first set of test cases, for a selected
set of mathematical functions. Next, we examine how the Julia interval
arithmetic library, IntervalArithmetic.jl, actually performs to these tests.
As this is an ongoing work, we list extra tests that we deem important to
perform.

Keywords: Unit tests for interval arithmetic libraries ; Test cases for
interval arithmetic ; Testing IEEE 1788-2015 compliance

##  1 Introduction

###  1.1 Context

Interval arithmetic is an arithmetic that operates on intervals, that is, on
sets of the form $[\underline{x},\bar{x}]$ rather than on numbers. The set
$[\underline{x},\bar{x}]$ is an interval if $\underline{x}\in\mathbb{R}$,
$\bar{x}\in\mathbb{R}$ and $\underline{x}\leq\bar{x}$. The result of the
arithmetic operation $\Diamond$ on the intervals $[\underline{x},\bar{x}]$ and
$[\underline{y},\bar{y}]$ (if $\Diamond$ takes two arguments), denoted by
$[\underline{x},\bar{x}]\Diamond[\underline{y},\bar{y}]$, is the smallest (for
inclusion) interval that contains the set $\\{x\Diamond
y,\>x\in[\underline{x},\bar{x}],\>y\in[\underline{y},\bar{y}]\\}$. For usual
operations, there are explicit formulas that express this definition in a more
amenable way, for instance:

$[\underline{x},\bar{x}]+[\underline{y},\bar{y}]=[\underline{x}+\underline{y},\bar{x}+\bar{y}],$
$[\underline{x},\bar{x}]-[\underline{y},\bar{y}]=[\underline{x}-\bar{y},\bar{x}-\underline{y}],$

and

$[\underline{x},\bar{x}]\times[\underline{y},\bar{y}]=[\min(\underline{x}\times\underline{y},\underline{x}\times\bar{y},\bar{x}\times\underline{y},\bar{x}\times\bar{y}),\max(\underline{x}\times\underline{y},\underline{x}\times\bar{y},\bar{x}\times\underline{y},\bar{x}\times\bar{y})].$

The fundamental principle of interval arithmetic is that, even if the input
arguments are known with uncertainty, such as $\pi\in[3.1415,3.1416]$, or if
the input arguments are only known to belong to a given set, the computed
result encloses the exact result (such as
$\sqrt{2\pi}\in[2.506628274631000,2.506628274631001]$) or it encloses the
range of the operation over the whole input sets (such as
$[0.5,1.5]\times\sin(\sqrt{[2,4]})\subseteq[0.4546,1.5]$). This fundamental
principle is sometimes expressed as ”Thou shalt not lie: the sought result
cannot be outside the computed result.” Indeed, interval arithmetic offers a
guarantee on the computed results, as they enclose the sought results.

However, usually interval arithmetic cannot be employed naively by replacing,
in a given mathematical expression or program, all numbers by intervals and
all operations by their interval counterparts. As a matter of fact, interval
computations are prone to overestimation, which can become so prominent as to
hinder the relevance of the result. A (stupid but archetypal) example is the
evaluation of the expression $x-x+x-x+x-x+x$, which is mathematically
equivalent to $x$ but, when evaluated using interval arithmetic with an
interval of width $w$, evaluates to an interval (containing the input
interval) of width $7w$. Interval analysis is thus the (fine) art of devising
expressions and programs that benefit from the guarantees provided by interval
arithmetic and that determine tight enclosures of the results. For a given
function $f$ and an interval $[\underline{x},\bar{x}]$, its goal is to
determine an enclosure of the range
$f([\underline{x},\bar{x}])=\\{f(x)\>:\>x\in[\underline{x},\bar{x}]\\}$
without a too large overestimation.

In this work, we do not address interval analysis but we focus on libraries
that implement interval arithmetic. First, the definition of an admissible
interval, for a library, must be clearly stated. Depending on the underlying
mathematical theory, unbounded intervals, the empty interval or intervals
$[\bar{x},\underline{x}]$ with $\bar{x}\geq\underline{x}$ are valid intervals
or not. In an attempt to clarify and to unify the definitions among the
interval community, a standard for interval arithmetic has been elaborated,
namely the IEEE 1788-2015 standard [12]. This standard first defines so-called
”common intervals” which are intervals of the form $[\underline{x},\bar{x}]$
with $\underline{x}\in\mathbb{R}$, $\bar{x}\in\mathbb{R}$ and
$\underline{x}\leq\bar{x}$, that is, connected, closed and bounded sets of
$\mathbb{R}$. The standard provides hooks to incorporate different
mathematical theories, called ”flavors”. For the time being, the only flavor
defined by the standard is the set-based flavor, based on the set theory; in
this flavor, $\emptyset$ and unbounded intervals are valid intervals.
Furthermore, for the set-based flavor, the standard defines a set of
decorations that provide additional information on the computation that has
led to the interval under consideration: there is a decoration attached to
each interval. Let us consider the example of the interval $[0,1]$ that has
been obtained as the result of $\sqrt{[-2,1]}$: $\sqrt{[-2,1]}$ is defined in
the set-based flavor as
$\sqrt{([-2,1]\cap\mathrm{Dom}_{\sqrt{~{}}})}=\sqrt{[0,1]}=[0,1]$, where
$\mathrm{Dom}_{\sqrt{~{}}}$ denotes the domain of $\sqrt{~{}}$, namely
$[0,+\infty)$. It is useful, in particular to avoid an incorrect application
of Brouwer theorem (see explanations and more details in the standard [12,
page 16]), to check whether such a intersection with the domain of the
operation (of $[-2,1]$ with $[0,+\infty)$, resulting in $[0,1]$ in our
example) has occurred. The role of the decoration, attached to the result
$[0,1]$, is to convey such information. Libraries implementing interval
arithmetic must clearly define upon which mathematical theory they are based;
for instance they can claim compliance with the IEEE 1788-2015 standard.

These libraries are based on an arithmetic on numbers, to be able to compute
for instance $\underline{x}+\underline{y}$ and $\bar{x}+\bar{y}$ when they
perform the interval addition
$[\underline{x},\bar{x}]+[\underline{y},\bar{y}]$. Often, the underlying
arithmetic is the floating-point arithmetic as defined by the IEEE 754-2008
standard [11], but not necessarily: arbitrary-precision floating-point
arithmetic or exact rational arithmetic can also be employed. Our focus is on
libraries implementing interval arithmetic, and more specifically on testing
these libraries.

Building tests along the development of a library is an important step towards
gaining confidence in the quality of this library. In this work, we focus on
unit tests for libraries implementing interval arithmetic. By a unit test, it
is meant that each operation, or function, of the library, is tested
individually; by contrast, benchmarks can be used to test a whole program that
utilizes the library, and will not be addressed here. We refer to [27] for an
example of such a benchmark. In our preliminary work [23], we presented our
vision for a set of such unit tests, that can be used as a basis to check
interval arithmetic libraries, or that can be incorporated as internal check
procedures.

Let us first recall briefly the different categories we identified in [23],
then we will sketch the content of this article, with an emphasis on our
contributions.

###  1.2 Which tests?

####  1.2.1 Definition of unit tests

We recall here our definition and notations for unit tests, following [23].
Let us denote by $\cal{L}$ the tested library, and interval quantities using
boldface: $\mathbf{x}$, $\mathbf{y}$. Here $\mathbf{x}$ can be an interval or
an interval vector, that is, a vector whose components are intervals, or any
other interval quantity such as a matrix, even if this latter case does not
occur yet in the presented work. We focus here on unit tests, for a function,
denoted by $f$. A unit test case is a pair composed of the input argument
$\mathbf{x}$ and the expected output $\mathbf{y}$.

First, the output $\mathbf{y}$ must be the tightest representable interval
enclosing $f(\mathbf{x})$; otherwise a very accurate library could compute
$\mathbf{z}$ such that $\mathbf{z}\subsetneq\mathbf{y}$ and still
$\mathbf{z}\supseteq f(\mathbf{x})$. To ensure the tightness property for
$\mathbf{y}$, typically, one computes the endpoints $\underline{y}$ and
$\bar{y}$ of $\mathbf{y}=[\underline{y},\bar{y}]$ using a precision (that is,
a number of binary – or decimal in case the library uses decimal arithmetic –
digits used for the computations) higher than the computing precision of
$\cal{L}$. Let us illustrate this with $\cal{L}$ using binary64 floating-point
numbers, on the (simple) example of the exponential function. We assume that
the given floating-point implementation of $\exp$ does not provide correct
rounding, but that, for any precision $q$, it satisfies
$RD_{q}(\exp(x))\leq\exp(x)\leq RU_{q}(\exp(x))$, for $RD_{q}$ rounding
downwards, and $RU_{q}$ rounding upwards, both in precision $q$. Given
$\mathbf{x}=[\underline{x},\bar{x}]$, to compute the infimum $\underline{y}$
of $\mathbf{y}$, one gets the approximations $RD_{q}(\exp(\underline{x}))$ and
$RU_{q}(\exp(\underline{x}))$ of $\exp(\underline{x})$ in high precision $q$,
and finally round them downwards in the target precision, $53$ for binary64.
If $RD_{53}(RD_{q}(\exp(\underline{x})))$ and
$RD_{53}(RU_{q}(\exp(\underline{x})))$ are equal, then they are the sought
value $\underline{y}$ for the infimum of $\exp(\mathbf{x})$.

Once the test cases are devised, what conclusions can be drawn from the
comparison between $\mathbf{y}$ and $\mathbf{z}$ the result computed by
$\cal{L}$? Requiring equality may be too demanding. Inclusion is required, but
we want to dismiss, for instance, an implementation of $\sin$ that returns
$[-1,1]$ for any argument. If the accuracy of the function is given in the
specification of the library, how can it be used and checked?

A last general recommendation is to incorporate the following procedure: if
$(\mathbf{x},\mathbf{y})$ is given as a pair of input and output, pick at
random (a reasonably large number of) values $x\in\mathbf{x}$ and check
whether $f(x)\in\mathbf{y}$ with $f(x)$ computed by the underlying arithmetic,
and whether $f([x,x])\subseteq\mathbf{y}$. Those tests are disconnected from
the knowledge of the implementation of $f$ and may hit a zone not considered
(forgotten) in the development of $f$, such as an overlooked quadrant for a
trigonometric function.

We now detail the different categories of unit tests we recommend. First,
there are unit tests that can be used to check any interval arithmetic
library. Then some tests are designed to check specifically IEEE 1788-2015
compliance. Finally, each library has its own unique features that must also
be tested. Let us detail these three categories.

####  1.2.2 Tests common to all libraries

As every interval arithmetic library implements usual arithmetic, mathematic
and interval-specific functions, unit tests must check these functions. The
set of unit tests must encompass:

  * •

easy test cases: these are useful especially during the early development
phase, to detect and correct “obvious” bugs;

  * •

special and exceptional values: the handling of such values is usually not
difficult, but it is easy to omit one of them; test cases must then be
exhaustive, to test every possible configuration;

  * •

cornercases: for instance when a function is not monotonic, such as the $\sin$
function, test cases must give examples where the interval includes, or not,
the optimum; cases where one of the endpoints is close to this optimum are
difficult ones and must appear in the set of test cases. Another class of
cornercases is related to the specificities of the underlying floating-point
arithmetic, if it is used to implement the interval arithmetic: it includes
subnormal values and hard-to-round cases. The test cases in this category
heavily depend on the underlying arithmetic, such as floating-point arithmetic
for rounding issues, and its precision (e.g. binary32 or binary64) for
subnormal values.

  * •

interval-specific functions: test cases must also target specific functions,
such as the ”union” of two intervals, which returns the convex hull of the
union, or the radius of the interval [9].

  * •

input and output: I/O functions must be thoroughly tested, in particular
because the most unexpected inputs can be given.

####  1.2.3 Testing IEEE 1788-2015 compliance

We do not detail here the content of the IEEE 1788-2015 standard [12] for
interval arithmetic. We only allude to the features relevant to the
development of test cases.

  * •

As I/O are very precisely defined by the standard, the different possible
inputs and outputs must be tested according to these specifications.

  * •

IEEE 1788-2015 is designed to allow for several mathematical theories to be
the theoretical foundations for the implemented interval arithmetic: each
theory is called a flavor. As the only standardized flavor so far is the set-
based flavor, test cases for the corresponding flavor must be given: they
encompass unbounded intervals as well as empty intervals.

  * •

The standard gives two lists of mathematical functions, the required ones
(such as addition or $\sin$) and the recommended ones (such as log2p1). The
tests should include test cases for all of these functions, but should not
fail in the event a recommended function is not implemented.

  * •

Decorations are a rather unique feature of the IEEE 1788-2015 standard: they
provide some additional information on the history of the computation that
yields the current result. Test cases must include all possible decorations
for the inputs, in order to check that the decoration of the output is
correctly determined.

  * •

The IEEE 1788-2015 standard specifies different accuracies, from tight, where
each result is the tightest representable interval that encloses the exact
result, to valid, where the computed result must simply include the exact one.
For the set-based flavor, an additional mode is the accurate mode, which is
very precisely defined by the standard, but corresponding values require care
to be determined. Test cases must include all of these possibilities, even if
the tested library provides only some of them; its behavior according to the
claimed accuracy must be tested.

####  1.2.4 Tests specific to some libraries

Interval arithmetic libraries may use different representations for the
interval, e.g. by their endpoints as in MPFI [24] or by their midpoint and
radius as in Arb [13]. They may rely on usual floating-point arithmetic as
Octave interval package [10], or arbitrary precision floating-point arithmetic
as MPFI, or exact rational arithmetic as JInterval [21]. Some libraries use
usual floating-point arithmetic but avoid changing the rounding mode, as Filib
[18] or JuliaInterval [25], or as GAOL [8] which can use either rounding-to-
nearest exclusively, or rounding-upwards exclusively. Each specificity must be
tested through several different test cases. Typically, several examples using
largely varied precisions must be tested for a library using arbitrary
precision. Examples which are cornercases for EFT (Error Free Transform) [20,
Chapter 4] must be included to test libraries that are based on rounding-to-
nearest and that use EFT to round outwardly their results.

However, since these tests are very specific to a given library, either their
presence in a database is optional, or running them should not be systematic.

###  1.3 Our contributions

Our contributions can be summarized as follows.

  * •

We introduce a JSON schema to describe interval tests in Section 2. This
allows the tests to be easily integrated in existing libraries.

  * •

We present an extensive database of tests in Section 3. This builds on
previous efforts such as ITF1788 presented in Section 1.4.2, though we expand
it with new test cases.

  * •

In particular, we generate tests whose intervals have hard-to-round floating
point numbers. These cases were generated using the program BaCSeL and are
particularly important to ensure that an interval implementation is robust to
floating-point rounding issues. How these tests have been generated is
detailed in Section 3.1.

  * •

Finally, we run our tests on the IntervalArithmetic.jl. As developed in
Section 4, we show how our newly generated hard-to-round cases can help unveil
hard to detect bugs, making hence our contribution valuable for interval
arithmetic library developers.

###  1.4 Related work

Our work relates to unit testing, that has largely been studied and discussed.
We sum up the main points below, before introducing existing platforms
offering unit testing of interval arithmetic libraries.

####  1.4.1 Unit testing

In the field of software testing, various approaches are possible and are
usually applied one after the other. Unit testing is usually the first step,
it consists in testing each developed function independently of the other
ones.

Of course, unit testing is useful to detect bugs. However, a function that
successfully passes the provided unit tests is not guaranteed to be correct.
An approach to guarantee correctness would be to use formal proof, but this is
not in the scope of this work. We refer to Gappa [5] and [20, Section 13.3],
and Flocq, based on Coq.Interval [3].

Unit testing is usually the first step to test a library. It is then followed
by integration testing that tests whether a whole program behaves properly,
and thus whether the various functions of the program perform properly
together. Relevant tests for this phase of integration testing are benchmarks
such as the one given in [27].

The field of software testing has developed various metrics to quantify the
quality of the tests, such as code coverage or branch coverage [29]. These
metrics are not relevant for our work. Indeed, our tests are oblivious to the
implementation, they only focus on the semantics of the computations: the
image of an interval by the sine function has a mathematical meaning, which is
decorrelated from the actual implementation of the tested sine function.

Finally, we would like to emphasize the conclusions of Ellims, Bridges and
Ince [6] about unit testing: unit tests allow developers to detect, at an
early stage of the development, many errors, for a cost per error, measured in
hours of human work, which is about a tenth of the cost of detecting an error
during the integration phase.

####  1.4.2 Existing platforms

Although interval arithmetic was standardized in 2015, no official test suite
to test for compliance was provided. Traditionally, each library relied on its
own in-house tests. However, Kuliamin [17] lists all kind of input values
needed to test thoroughly a floating-point arithmetic library. For the
exponential function for instance, this amounts to 15,000 test cases. For an
interval-valued function, one gets $15,000\times 15,000$ possible inputs. This
may be extreme, but it illustrates the need to share unit tests, rather than
letting each developer of an interval arithmetic library devise her/his own
tests.

The first effort in this direction was the JInterval P1788 Test Launcher [22]
based on the JInterval library [21] for interval arithmetic in Java. For each
tested library, a wrapper must be written to enable to call the operations and
functions of this library. Wrappers are available for Profil/BIAS,
boost/interval, C-XSC, Filib, libieeep1788, libMoore, MPFI. The launcher reads
tests from plain text files of simple human-readable format and writes the
results computed using tested library and JInterval and their relation into a
plain text report. Test set included in the JInterval P1788 Test Launcher
consists of over 14,000 tests which partly originated from libieeep1788,
Filib, libMoore while the other ones are original.

Another attempt to introduce a systematic testing framework for interval
arithmetic was the Interval Testing Framework 1788 (ITF1788) [14]. Developers
of ITF1788 introduced a domain specific language called ITL (Interval Testing
Language) that could be used to express tests for interval arithmetic. ITF1788
is a Python engine that inputs ITL-files and converts tests into code for the
specified language, test framework and library. However, while being an
elegant language, this ITL solution presents a few drawbacks. First, being an
ad-hoc language, it also requires an ad-hoc parser for parsing the tests into
a structured representation and/or produce code in the target language. To
overcome this, we choose to use the JSON format, which is widely used in
software engineering and comes with available parsers. Secondly, the proposed
language lacked flexibility to test different precision or to test for both
accurate and tight mode; our proposal offers the required flexibility.

After the introduction of ITL, a first attempt to produce an exhaustive test-
suite was presented in [15]. This work both collected tests from existing
interval libraries, namely C-XSC, FILIB, MPFI, libieeep1788, and produced
several new test cases for features introduced by the standard, such as
decorated intervals. In our work, we build on top of this previous attempt, by
porting the existing ITF1788 tests into JSON format and enriching the test
suite with new test cases for the accurate mode, as well as new tests for
different precision and tests for hard-to-round cases.

###  1.5 Parallel or distributed computations

Testing correctness of interval libraries under parallel behavior is
important, especially for libraries that handle floating-point error by
changing the rounding mode. In several common processors, the rounding mode of
the floating processing unit is controlled by a global state in the unit
registers. If not handled with care, manually changing the rounding mode in
parallel codes can lead to race conditions and undefined behavior. To test
library soundness under multithreaded execution, the tests presented here
could be run in parallel. In addition to speeding up the tests execution, this
could also be used to test that correct rounding is preserved under parallel
execution.

Parallelism was also used to generate the hard-to-round cases: the code that
determines hard-to-round cases, called BaCSeL, uses multithreading to explore
in parallel many subintervals. Some timings are given in [26], that justify
the need for large computing power and parallelism.

##  2 A specification language for interval arithmetic testing

We propose to specify the tests for our framework using the JSON format [16].
To quote the introduction to JSON, “JSON (JavaScript Object Notation) is a
lightweight data-interchange format. It is easy for humans to read and write.
It is easy for machines to parse and generate. It is based on a subset of the
JavaScript Programming Language$[\ldots]$ JSON is a text format that is
completely language independent but uses conventions that are familiar to
programmers of the C-family of languages$[\ldots]$” This being a popular
format, it can be easily read seamlessly by several languages, while still
being human-readable. JSON being the most common format for web-development,
this could potentially allow to develop a web-interface to test interval
libraries.

Each test-case is a JSON object with the following fields:

  * •

function: specifies the name of the function-under-test (FUT), must be as
named in the standard;

  * •

precision: the number of bits of precision used;

  * •

input: specifies the input, see below for description of how the input is
specified;

  * •

output: an object with two fields: tight (required), and accurate (optional),
specifying the result for the tight and accurate modes. For functions
returning a number or a boolean (e.g. a function returning the midpoint of the
interval or a comparison), only the tight mode shall be provided. The format
corresponding to each mode is similar to the format used for the input.

Both inputs and outputs must have a field type, which can have one of the
following values: number, interval, boolean or string. For types boolean,
number and string, the only required additional field is val, which must
contain the value. If the type is interval, then the object has two compulsory
fields inf and sup, specifying the lower and upper endpoints of the interval,
respectively. Intervals can have an optional third field, dec, specifying the
decoration.

The following code snippet shows an example of a function.

      {
        "function": "atanh",
        "precision": 23,
        "input": [
          {
            "type": "interval",
            "inf": "-0xf.fe1e00@-1",
            "sup": "0xf.fe1e00@-1"
            "dec": "com"
          }
        ],
        "output": {
          "tight": {
            "type": "interval",
            "inf": "-0x4.305fa0@0",
            "sup": "0x4.305fa0@0",
            "dec": "com"
          },
          "accurate": {
            "type": "interval",
            "inf": "-0x4.306830@0",
            "sup": "0x4.306830@0",
            "dec": "com"
          }
        }
      }

In the above code snippet, the inf and sup entries are hexadecimal numbers
whose exponent is in base 16.

##  3 New test cases

One of our contributions is to devise and to share publicly a suite of unit
test cases, for arithmetic and mathematical operations and functions. Ideally,
the whole set of functions recommended by the IEEE 1788-2015 standard would be
covered. In this work, we share the set of test cases which have already been
gathered within ITF-1788 [15] (which was a collection of test cases from
various existing libraries: C-XSC, FILIB, MPFI, libieeep1788 in particular)
and new ones.

Our goal is to offer and to share test cases in all categories in Section 1.2.
In what follows, we do not develop the first categories, and we focus on hard-
to-round cases because they took most of our time.

So, here are just a few indications about the generation of the other test
cases.

  * •

Easy cases are small, usual values. The choices may be arbitrary, but they
offer the advantage to be easy to read and there are quantities that are easy
to check, such as sign or magnitude; examples include $\sin 0$ or $\exp 0$ or
$\operatorname{cbrt}(-64)$, where $\operatorname{cbrt}$ stands for cubic root.

  * •

Special cases are also included, they cover both exact cases and floating-
point special values. For transcendental functions, there are very few ”exact”
cases, apart from $0$ (as argument or as result): this is true for functions
deriving from the exponential as stated by Lindemann-Weierstrass theorem [19,
28], and for functions deriving from the logarithm as proven by Baker’s
theorem [2] (see also Wikipedia for a more accessible presentation of these
theorems). Special values, for floating-point arithmetic, are $0$ (with its
two signed representations $-0$ and $+0$) and infinities, both as inputs and
as outputs, as in $\exp(-\infty)$ or $\operatorname{atanh}(1)$. Each of these
special values can appear as the left endpoint or right endpoint of an
interval, or both.

  * •

Exceptions. We must then test for instance $\exp[-\infty,-\infty]$: in this
case, the input argument is ill-formed and can be considered either as the
empty set (if the flavor so permits) or as an error. This is an example of an
exceptional behavior. Other examples of exceptional values that must be tested
include NaN as one or both endpoints. We did our best to include all such
possibilities in our test cases.

  * •

Cornercases. The bulk of our contribution in this section consists in the
determination of cornercases, in particular in the determination of values
which are hard to round. Indeed, these cases are difficult to evaluate
accurately and thus a library is most likely to fail. In what follows, we
develop our approach to determine these cornercases.

###  3.1 Generation of the test cases

####  3.1.1 Generation of hard-to-round cases

For a floating-point number $x$, its image by a function $f$: $f(x)$ is
written in the working base $B$ as $f(x)=\pm y_{0}.y_{1}\ldots
y_{p-1}y_{p}\ldots.B^{e}$ for some exponent $e$, with an infinite sequence
$(y_{i})$ of bits (if $B=2$) or digits (if $B=10$). Without loss of
generality, let’s assume that $B=2$. The floating-point number $\hat{y}$ that
represents $f(x)$ in the floating-point format has $p$ bits, it is the
rounding of the value $f(x)$.

If, for instance, the rounding mode is towards $+\infty$ and the approximation
of $f(x)$ using $n\geq p$ bits is $\pm y_{0}.y_{1}\ldots y_{p-1}0\ldots
0.B^{e}$, i.e. $y_{p}=\ldots=y_{n-1}=0$ and the error of the approximation is
bounded by $B^{e-n+1}$ which is the weight of the last bit of this
approximation, then one doesn’t know how to round $f(x)$, indeed there are two
possible candidates: $\pm y_{0}.y_{1}\ldots y_{p-1}.B^{e}$ and
$\pm(y_{0}.y_{1}\ldots y_{p-1}+B^{-p+1}).B^{e}$. When the difference between
$p$ and $n$ is large, one must evaluate $f(x)$ with a large precision $\geq n$
to be able to round $f(x)$ correctly into $\hat{y}$: $x$ is said to be a hard-
to-round case for $f$ and the floating-point format. In what follows, we do
not focus on the hardest-to-round cases, that is, the values of $x$ that
require the largest value of $n$, but only on sufficiently hard-to-round ones.
This problem is well-known to people who produced (by hand, before computers
existed) tables for mathematical functions, hence its name: the Table-Maker
Dilemma, TMD in short.

The hard-to-round cases we propose have been determined using BaCSeL [1], a
program of about 5,000 lines of code in C, written and maintained mostly by P.
Zimmermann. It is primarily intended for a personal use, with a concise
documentation given as a README and an installation that requires GMP, MPFR,
FPLLL and OpenMP. It has been developed since 2020 and it aims at determining
hard-to-round cases for several mathematical functions, over specific
subdomains. It is based on several algorithms that rely on LLL to prune easy
cases and to focus on hard ones. The computing basis is either $2$ or $10$, in
what follows we focus on the basis $2$. The list of available functions
contains $\exp$, $\exp_{2}$, $\log$ and $\log_{2}$, $\sin$, $\cos$, $1/x$,
$\sqrt{x}$, $1/\sqrt{x}$, $\sqrt[3]{x}$, $\operatorname{acos}$ where the
subscript indicates the basis of the exponential or the logarithm:
$\exp_{2}(x)=2^{x}$ and $\log_{2}x=\log x/\log 2$. Each subdomain is split
into smaller intervals, that can be further subdivided to reach the required
accuracy. These smaller intervals are explored in parallel, using OpenMP to
distribute the work. Each function is internally approximated by a Taylor
expansion on each small interval.

More precisely, the BaCSeL determines all hard-to-round cases for a given
function $f$ over an interval $[t_{0},t_{1}]$. The interval is such that
$t_{0}<t_{1}$ and $t_{0}$, $t_{1}$ belong to the same binade, in other words
there exists an integer $n$ such that $2^{n-1}\leq t_{0}<t_{1}\leq 2^{n}$ (and
the interval $[2^{n-1},2^{n})$ is called a binade). It is required that the
image of $[t_{0},t_{1}]$ by $f$ does not overlap several binades. The
”hardness-to-round” is given by the user, it corresponds to the number $n-p$
in our definition of a hard-to-round number. The precision of the input
arguments and of the evaluated results must be specified by the user, as well
as the computing precision needed to evaluate accurately enough the function,
in order to determine whether an argument $x$ is a hard-to-round case; in
particular this computing precision must be larger than $n$. Several other
values must be fed to BaCSeL so that it determines hard-to-round cases. The
short description given here gives a hint of the number of trial-and-errors
required to choose proper intervals and precision, in order to get a positive
and reasonable number of hard-to-round cases. In other words, some expertise
along with computational power are useful to determine hard-to-round cases.

####  3.1.2 Computation of the expected results

Once the input values are listed, either easy ones, or special ones, or
cornercases, then the expected result must be computed. In this work, we rely
on MPFI [24] to evaluate the given function $f$ over the input interval
$\mathbf{x}$. MPFI is designed to return $f_{\mathrm{tightest}}(\mathbf{x})$,
the tightest interval enclosing the mathematical image $f(\mathbf{x})$, with
endpoints that are floating-point numbers at the prescribed precision. It thus
gives the target result for the ”tight” mode of the IEEE 1788-2015 standard.

MPFI is also used to give the largest admissible result for the ”accurate”
mode of the set-based flavor. This mode is defined as follows in the standard
[12]. Let $x$ be a floating-point number and let us denote by
$\mathrm{nextUp}(x)$ the smallest floating-point number strictly larger than
$x$, it can be $+\infty$. Similarly, let us denote by $\mathrm{nextDown}(x)$
the largest floating-point number strictly smaller than $x$, it can be
$-\infty$. Finally, let us denote by nextOut the function that returns, for an
interval $\mathbf{x}=[\underline{x},\bar{x}]$ where $\underline{x}$ and
$\bar{x}$ are floating-point numbers, the interval
$[\mathrm{nextDown}(\underline{x}),\mathrm{nextUp}(\bar{x})]$. The ”accurate”
mode requires that the evaluation of a function $f$ over an interval
$\mathbf{x}$ is enclosed in
$\mathrm{nextOut}(f_{\mathrm{tightest}}(\mathrm{nextOut}(\mathbf{x})))$. MPFI
has been enriched with the $\mathrm{nextOut}$ function, so as to determine
this interval as well.

####  3.1.3 Format for the input and the expected output

In this work, we use the capability of both BaCSel and MPFI to output their
results in hexadecimal format: this enables writing and reading floating-point
numbers without any conversion error from and to radix 10.

###  3.2 Four mathematical functions as illustration: cubic root, exp, sin,
atanh

As the IEEE 1788-2015 standard for interval arithmetic requires over 50
mathematical operations and functions and recommends 16 additional ones for
the set-based flavor, we will not detail the test cases for each function but
we will focus on a few selected ones. The ones we choose here are the cubic
root, the exponential, the sine and the hyperbolic arc-tangent. These test
cases can be found at https://github.com/lucaferranti/IntervalArithmeticTests.

####  3.2.1 Test cases for cbrt

The cubic root is a relatively simple function, as it is algebraic, it is
defined on the whole set of real numbers, and it exhibits a symmetry (indeed
it is odd) which can easily be checked. Furthermore, examples involving the
cubic root are easy to check, at least approximately: first because it is
monotonic, and also as multiplying each endpoint of the result by itself twice
should land close to the original argument. Last, because the cubic root is
algebraic, many inputs yield an exactly representable result. Thus, exact
cases are numerous. Examples where the floating-point endpoints have exponents
which are integer multiples of 3 are easily reduced to the former case.

Figure 1: Graph of the cubic root. (Credit: Krishnavedala - CC BY-SA 4.0.)

In our test cases, we focus on the zone close to $0$, where the function has a
vertical tangent, see Figure 1: since it is the place where the variations are
the fastest, this is the zone where we expect numerous hard-to-round cases to
lie.

####  3.2.2 Test cases for exp

Our other three functions are transcendental ones. The simplest one is the
exponential: it is a very well studied function, it usually comes with a high
quality implementation. It is also a ”simple” function, in the sense that it
is defined everywhere, it is positive, monotonic, it has only one ”exact”
value in the sense that only one floating-point value, namely $0$, has an
image which is exactly representable as a floating-point number, namely $1$.
The domain on which it takes finitely representable values is bounded from
above: for the Binary64 format, $expx$ overflows for $x\geq 710$.

Our test cases include mostly hard-to-round values around $0$, that have been
determined using BaCSeL. This is a first step to assess the quality of the
implementation of the exponential function. Indeed, to evaluate $\exp x$, the
implementation usually involves a range-reduction step, where the argument $x$
is repeatedly divided by $2$ until it becomes $x^{\prime}=x/2^{n}$ close to
$0$. Then $\exp x^{\prime}$ is evaluated accurately. Eventually $\exp x$ is
obtained as $\exp x^{\prime}$ squared $n$ times: this last phase is called
reconstruction.

Future work will consist in adding values:

  * •

close, but smaller than the upper limit of the domain where $\exp$ overflows;

  * •

large in absolute value, and negative, to check whether subnormal results are
properly handled;

  * •

more varied in magnitude, to check the accuracy of the reconstruction phase as
well.

####  3.2.3 Test cases for sin

Our second transcendental function is the sine function. Again, there is only
one ”exact” value: $\sin 0=0$. This function is more difficult to implement,
as it is non monotonic. As it is periodic, the first step to evaluate $\sin x$
is again a range-reduction phase: it consists in translating $x$ by the
appropriate integer $k$ multiple of the period $2\pi$, to obtain
$x^{\prime}=x-2k\pi$ close to $0$. Contrary to what happened for the
exponential, there is no reconstruction phase for the sine function.

The range-reduction can be delicate: as $\pi$ is transcendental, the
implementation can require a very large computing precision to determine $k$
and then to compute accurately $x^{\prime}=x-2k\pi$, especially when $x$ is
large. We test such extremely large arguments $x$ in our test cases. The other
test cases we propose are hard-to-round cases close to $0$: these values are
useful to determine whether the evaluation of the sine is accurate, once the
range-reduction phase is performed.

It is easy to check whether the results are correct, in the ”tight” mode: for
such tiny values, $\sin x\sim x$ applies. It is also easy to check whether the
results are correct in the ”accurate” mode: the input argument has been
enlarged by one ulp at each endpoint, and the tight evaluation with this
enlarged argument is again enlarged by one ulp at each endpoint, the
difference between the tight and the accurate evaluations are 2 ulps at each
endpoint. As the values are given in hexadecimal format, this difference of 2
ulps at each endpoint is easy to check.

Future work will consist in adding more varied input arguments, both in sign
and in magnitude, and also in adding intervals closer to an extremum
$\frac{\pi}{2}+k\pi$, either containing the extremum or not. Some expertise in
approximation theory is needed to identify floating-point arguments that are
very close to an optimum; continued fractions are useful as illustrated in
Chapters 4 and 10 of [20]. Again, all our test cases have been determined
using BaCSeL.

####  3.2.4 Test cases for atanh

Our last transcendental function is the hyperbolic arc-tangent
$\operatorname{atanh}$. We choose this function because, in our experience,
this function is usually implemented with much less accuracy than the two
previous, more usual, ones. Again there is only one ”exact” value”
$\operatorname{atanh}0=0$. This function is rather simple to study: it is odd,
monotonic, defined on a bounded domain $(-1,1)$, with infinite limits at the
boundaries, see Figure 2:

$\lim_{x\rightarrow-1^{+}}\operatorname{atanh}x=-\infty\mbox{ and
}\lim_{x\rightarrow 1^{-}}\operatorname{atanh}x=+\infty.$

Figure 2: Graph of the hyperbolic arc-tangent function.

Our test cases contain hard-to-round cases determined using BaCSeL. We first
had to integrate the hyperbolic arc-tangent function in the list of functions
known by BaCSeL. Our goal was to cover several intervals bounded by successive
powers of $2$ and included in the domain of $\operatorname{atanh}$, including
values close to the boundaries of the domain, however $\operatorname{atanh}$
does not grow very rapidly and we did not find (yet) any value where the
evaluation overflowed.

##  4 Results with the Julia Interval library

The IntervalArithmetic.jl package is written in the Julia programming
language, motivated by Julia’s mathematical-style syntax, near C-performance,
high-level interactive usability, very robust type system with easy to exploit
composability properties. These capabilities of the language permit to have
parametric interval that depend on the underlying floating-point
representation, which can be extended to incorporate the decoration system,
the possibility for different accuracy modes, and flavors. The package is on
its way to be fully conforming with the IEEE 1788-2015 Standard, with most
tests of the ITL test suite successfully passing.

IntervalArithmetic.jl implements interval-types in the form
Interval$\\{$T$\\}$, where T refers to the underlying floating-point
representation. For the tight accuracy mode, the package relies on the CRlibm
library [4] for correctly rounded Binary64, and on the MPFR library [7] for
Binary32 and arbitrary precision. The CRlibm library does not have all
functions that the IEEE 1788-2015 standard requires, for instance,
$\operatorname{atanh}$ is not included. In those cases, we use the MPFR
library, and convert the final result to the floating-point representation of
the input.

The results discussed below use the 1.0-dev branch, and are produced with the
tight mode of accuracy only. For this accuracy mode, we check the identity of
the result produced by the library and the hard-to-round results. Below we do
not report the tests that involve the empty set, since they are already
addressed by the ITL test suite.

All test cases for the $\exp$ and $\sin$ functions pass correctly. For these
functions, the hard-to-round tests we produced include examples for Binary32
and Binary64 numbers only. While these results are positive, it remains to
produce and test hard-to-round cases that involve other precisions than the
usual formats.

For the $\operatorname{atanh}$ and $\operatorname{cbrt}$ functions, the
results show tests that pass and others that fail. In both cases, the tests
that systematically pass correspond to the Binary64 format. This statement is
not trivial since none of these functions is included in the CRlibm library.
The failing tests then involve the Binary32 format, or cases involving
different precision values for extended precision format. The positive aspect
of the failing tests that we have uncovered is that some functions that have
naively been extended, say to Binary32, should be handled more carefully.

##  5 Conclusions and future work

As interval arithmetic aims at producing verified computations that can be
used for rigorous mathematical proofs, standardization and testing is a
crucial topic. Despite the standardization of interval arithmetic in 2015,
little attention has been given to how testing should be implemented to ensure
the library is sound and standard compliant. In this paper, we reviewed the
principal approaches and challenges of testing interval libraries. Building on
previous work, we produced an open source collection of tests for interval
arithmetic libraries. This test suite is released in JSON format, which allows
it to be easily integrated by different libraries in different languages. In
addition to collecting already existing tests, our test suite tackles
challenges that were previously unaddressed. First, our test suite contains
test cases for the accurate mode described in the standard. This helps
libraries developers to determine whether their results are tight, accurate or
simply valid. Second, using the BaCSeL program, we generated hard-to-round
cases for some elementary functions. These are important cornercases to test
the library robustness against floating-point rounding issues.

As future work, we plan to keep expanding our test suite with more hard-to-
round cases both for Binary32 and Binary64, for the whole set of functions
that are either required or recommended by the IEEE 1788-2015 standard for
interval arithmetic. We also plan to expand our test suite to check more
thoroughly specificities of the standard, such as the computation of
decorations, or the handling of unbounded or empty intervals.

We also plan to test other libraries. For instance, we could test the accuracy
of the ARB library, which uses a midpoint-radius representation for the
intervals. Indeed, while interval arithmetic is usually concerned with inf-sup
format, the midpoint-radius format has proven itself an appealing alternative,
as one can choose the appropriate trade-off between accuracy and efficiency.
Algorithms that use the midpoint-radius representation can be fast but at the
expense of the accuracy, and this seems to be the case of ARB, at least for
some operations, or they can be slower and use more complicated formulas to
remain very accurate. It would be interesting to determine whether the results
computed by ARB, once converted to a representation by endpoints, fall into
the tight, accurate or simply valid category.

Finally, midpoint-radius arithmetic has received less attention and there has
not been much effort in standardizing it. Future work will focus on generating
test cases for midpoint-radius arithmetic, which can be another approach to
assess the accuracy of a library as popular as ARB.

### Acknowledgement

Luis Benet acknowledges funding from IG-101122 DGAPA (UNAM) project.

A CC-BY public copyright license has been applied by the authors to the
present document and will be applied to all subsequent versions up to the
Author Accepted Manuscript arising from this submission.

## References

  * [1] BaCSeL. https://gitlab.inria.fr/zimmerma/bacsel
  * [2] Murty MR, Rath P. (2014). Baker’s Theorem. Transcendental Numbers 2014: 95–100.  10.1007/978-1-4939-0832-5_19. ISBN 978-1-4939-0831-8.
  * [3] Boldo S, Melquiond G. Some formal tools for computer arithmetic: Flocq and Gappa. In: 28th IEEE Symposium on Computer Arithmetic (ARITH). ; 2021, pp. 111–114.
  * [4] Daramy-Loirat C, Defour D, de Dinechin F, et al. CR-LIBM A library of correctly rounded elementary functions in double-precision. research report, LIP,; https://hal-ens-lyon.archives-ouvertes.fr/ensl-01529804 2006\.
  * [5] Daumas M, Melquiond G. Certification of bounds on expressions involving rounded operators. ACM Transactions on Mathematical Software 2010, no 37(1), pp. 1–20.
  * [6] Ellims M, Bridges J, Ince D. Unit testing in practice. In: IEEE 15th International Symposium on Software Reliability Engineering. ; 2004, pp. 3–13.
  * [7] Fousse L, Hanrot G, Lefèvre V, Pélissier P, Zimmermann P. MPFR. ACM Transactions on Mathematical Software 2007, no 33(2), article 13.  10.1145/1236463.1236468.
  * [8] Goualard F. Gaol: Not just another interval library. https://sourceforge.net/projects/gaol/ 2005\.
  * [9] Goualard F. How do you compute the midpoint of an interval?. ACM Transactions on Mathematical Software 2014, no 40(2), pp. 1–25.
  * [10] Heimlich O. Interval arithmetic in GNU Octave. In: SWIM 2016, 9th Summer Workshop on Interval Methods. ; 2016; ENS de Lyon, France.
  * [11] IEEE: Institute of Electrical and Electronic Engineers. IEEE Standard for Floating-Point Arithmetic IEEE Std 754-2008 2008, pp. 1–97.
  * [12] IEEE: Institute of Electrical and Electronic Engineers. 1788-2015 - IEEE Standard for Interval Arithmetic. IEEE Std 1788-2015 2015, pp. 1–97.
  * [13] Johansson F. Arb: a C library for ball arithmetic. ACM Communications in Computer Algebra 2013, no 47(4), pp. 166–169.
  * [14] Kiesner M, Nehmeier M, Gudenberg W. vJ. ITF1788: An Interval Testframework for IEEE 1788. Report 495, Department of Computer Science; University of Würzburg: https://www.researchgate.net/publication/278620157_ITF1788_An_Interval_Testframework_for_IEEE_1788 2015\.
  * [15] ITF1788 - Interval Test Framework for IEEE Std 1788-2015. https://github.com/oheim/ITF1788
  * [16] JSON: https://www.json.org/
  * [17] Kuliamin V. Standardization and Testing of Implementations of Mathematical Functions in Floating Point Numbers. Programming and Computer Software 2007, no 33(3), pp. 154–173.
  * [18] Lerch M, Tischler G, Gudenberg W. vJ, Hofschuster W, Krämer W. filib++, a Fast Interval Library Supporting Containment Computations. ACM Transactions on Mathematical Software 2006, no 32(2), pp. 299–324.
  * [19] Lindemann, F. Über die Ludolph’sche Zahl.. Sitzungsberichte der Königlich Preussischen Akademie der Wissenschaften zu Berlin 1882, no 2, pp. 679–682.
  * [20] Muller JM, Brunie N, de Dinechin F, et al. Handbook of Floating-Point Arithmetic (2nd ed.).  Birkhäuser.  2018\.
  * [21] Nadezhin D, Zhilin S. JInterval Library: Principles, Development, and Perspectives. Reliable Computing 2013, no 19(3), pp. 229–247.
  * [22] P1788 Test Launcher (based on JInterval Library). https://github.com/jinterval/jinterval/tree/master/p1788-launcher-java
  * [23] Revol N, Benet L, Ferranti L, Zhilin S. Testing interval arithmetic libraries, including their IEEE-1788 compliance. In: Wryzykowski R. , ed. PPAM 2022: Parallel Processing and Applied Mathematics. LNCS, no 13827, pp. 428-440, 2023. Springer.
  * [24] Revol N, Rouillier F. Motivations for an Arbitrary Precision Interval Arithmetic and the MPFI Library. Reliable Computing 2005, no 11(4), pp. 275–290.  10.1007/s11155-005-6891-y.
  * [25] Sanders D, Benet L, Ferranti L, Richard B, and collaborators. JuliaIntervals/IntervalArithmetic.jl. https://github.com/JuliaIntervals/IntervalArithmetic.jl DOI: 10.5281/zenodo.3336308 (all versions).  10.5281/zenodo.3336308.
  * [26] Sibidanov A, Zimmermann P, Glondu S. The CORE-MATH Project. In: 29th IEEE Symposium on Computer Arithmetic (ARITH). ; 2022.  10.1109/ARITH54963.2022.00014.
  * [27] Tang Z, Ferguson Z, Schneider T, Zorin D, Kamil S, Panozzo D. A Cross-Platform Benchmark for Interval Computation Libraries. In: Wryzykowski R. , ed. PPAM 2022: Parallel Processing and Applied Mathematics. LNCS, no 13827, pp. 415-427, 2023. Springer.
  * [28] Weierstrass, K. Zu Lindemann’s Abhandlung. ”Über die Ludolph’sche Zahl”. Sitzungsberichte der Königlich Preussischen Akademie der Wissenschaften zu Berlin 1885, no 5, pp. 1067–1085.
  * [29] Zhu H, Hall P, May J. Software unit test coverage and adequacy. ACM Computing Surveys 1997, no 29(4), pp. 366–427.


Lagrangian and Hamiltonian Formulation of

Classical Electrodynamics without Potentials

Dan N. Vollick

Irving K. Barber School of Arts and Sciences

University of British Columbia Okanagan

3333 University Way

Kelowna, B.C.

Canada

V1V 1V7

Abstract

In the standard Lagrangian and Hamiltonian approach to Maxwell’s theory the
potentials $A^{\mu}$ are taken as the dynamical variables. In this paper I
take the electric field $\vec{E}$ and the magnetic field $\vec{B}$ as the the
dynamical variables. I find a Lagrangian that gives the dynamical Maxwell
equations and include the constraint equations by using Lagrange multipliers.
In passing to the Hamiltonian one finds that the canonical momenta
$\vec{\Pi}_{E}$ and $\vec{\Pi}_{B}$ are constrained giving 6 second class
constraints at each point in space. Gauss’s law and
$\vec{\nabla}\cdot\vec{B}=0$ can than be added in as additional constraints.
There are now 8 second class constraints, leaving 4 phase space degrees of
freedom. The Dirac bracket is then introduced and is calculated for the field
variables and their conjugate momenta.

##  1 Introduction

In the standard Lagrangian and Hamiltonian approach to Maxwell’s theory the
potentials $A^{\mu}$ are taken as the dynamical variables (see for example [1,
2]). The Lagrangian $L=\frac{1}{4}F^{\mu\nu}F_{\mu\nu}$, where
$F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$, gives Gauss’s law
and Maxwell’s version of Ampere’s law without sources while,
$F_{[\mu\nu,\alpha]}=0$ gives Faraday’s law and $\vec{\nabla}\cdot\vec{B}=0$.
The canonical momenta are given by $\Pi^{\mu}=F^{\mu 0}$, so that one runs
into a problem in the Hamiltonian approach since $\Pi^{0}=0$. This problem can
be dealt with by using Dirac’s approach [3, 4] for systems with constraints.
The consistency of the primary constraint $\Pi^{0}=0$ implies a secondary
constraint $\partial_{k}\Pi^{k}=0$. These two constraints are first class and
together generate a set of transformations that includes the standard gauge
transformations.

In this paper I find a Lagrangian with $\vec{E}$ and $\vec{B}$ as the
dynamical variables that gives the dynamical Maxwell equations. The constraint
equations can be included by adding them to this Lagrangian using Lagrange
multipliers. The canonical momenta $\vec{\Pi}_{E}$ and $\vec{\Pi}_{B}$ are
constrained giving 6 second class constraints at each point in space. Gauss’s
law and $\vec{\nabla}\cdot\vec{B}=0$ can than be added in as additional
constraints. There are now 8 second class constraints at each point in space
leaving 4 phase space degrees of freedom. The Dirac bracket is then introduced
and is calculated for the field variables and their conjugate momenta.

##  2 The Lagrangian Formulation

Maxwell’s equations can be divided into the dynamical equations (with $c=1$)

$\vec{\nabla}\times\vec{E}=-\frac{\partial\vec{B}}{\partial t}\hskip
72.26999pt\vec{\nabla}\times\vec{B}=\frac{\partial\vec{E}}{\partial
t}+\frac{4\pi}{c}\vec{j}$ (1)

and the constraint equations

$\vec{\nabla}\cdot\vec{E}=4\pi\rho\hskip
72.26999pt\vec{\nabla}\cdot\vec{B}=0.$ (2)

In this section I will find a Lagrangian that involves $\vec{E}$ and $\vec{B}$
and their first derivatives that gives the dynamical equations of motion. The
constraint equations can then be included using Lagrange multipliers.

Consider the Lagrangian

$L_{D}=a\vec{E}\cdot\left(\frac{\partial\vec{B}}{\partial
t}\right)+b\vec{B}\cdot\left(\frac{\partial\vec{E}}{\partial
t}\right)+c\vec{E}\cdot\left(\vec{\nabla}\times\vec{B}\right)+d\vec{E}\cdot\left(\vec{\nabla}\times\vec{E}\right)+e\vec{B}\cdot\left(\vec{\nabla}\times\vec{B}\right),$
(3)

where $a,b,c,d$ and $e$ are constants. I have not included a term of the form
$\vec{B}\cdot\left(\vec{\nabla}\times\vec{E}\right)$ since
$\vec{B}\cdot\left(\vec{\nabla}\times\vec{E}\right)=\vec{E}\cdot\left(\vec{\nabla}\times\vec{B}\right)+\vec{\nabla}\cdot(\vec{E}\times\vec{B})$.
The Euler-Lagrange equations ($\mu=0,1,2,3$ and $k=1,2,3$)

$\frac{\partial}{\partial x^{\mu}}\left(\frac{\partial
L_{D}}{\partial(\partial_{\mu}E_{k})}\right)-\frac{\partial L_{D}}{\partial
E_{k}}=0$ (4)

give

$(b-a)\frac{\partial\vec{B}}{\partial
t}-c(\vec{\nabla}\times\vec{B})-2d(\vec{\nabla}\times\vec{E})=0\;.$ (5)

To obtain Faraday’s law we require that

$c=0\hskip 36.135ptand\hskip 36.135pta-b=2d\;,$ (6)

with $a-b\neq 0$. The Euler-Lagrange equations

$\frac{\partial}{\partial x^{\mu}}\left(\frac{\partial
L_{D}}{\partial(\partial_{\mu}B_{k})}\right)-\frac{\partial L_{D}}{\partial
B_{k}}=0$ (7)

give

$(a-b)\frac{\partial\vec{E}}{\partial t}-2e(\vec{\nabla}\times\vec{B})=0\;.$
(8)

To obtain Ampere’s law without sources we require that

$a-b=2e\;.$ (9)

From this we can also conclude that $d=e=\frac{1}{2}(a-b)$. Thus, the
Lagrangian

$L_{D}=f\left[\vec{E}\cdot\left(\vec{\nabla}\times\vec{E}+\frac{\partial\vec{B}}{\partial
t}\right)+\vec{B}\cdot\left(\vec{\nabla}\times\vec{B}-\frac{\partial\vec{E}}{\partial
t}\right)\right]$ (10)

gives the dynamical Maxwell equations for arbitrary $f=\frac{1}{2}(a-b)\neq
0$. I also used the relationship $\vec{E}\cdot\frac{\partial\vec{B}}{\partial
t}=-\vec{B}\cdot\frac{\partial\vec{E}}{\partial
t}+\frac{\partial(\vec{E}\cdot\vec{B})}{\partial t}$ to derive (10).
Throughout the rest of the paper I will take $f=1$.

Sources can be included by adding

$L_{S}=-8\pi\vec{j}\cdot\vec{B}$ (11)

to $L_{D}$.

Now consider the constraint equations. One approach is to add them to to the
dynamical equations outside the action principle. Another approach is to
include them in the action with the use of Lagrange multipliers. However, as
we shall see below, additional constraints have to be imposed on the Lagrange
multipliers to obtain Maxwell’s theory. In this approach

$L_{C}=\lambda_{1}\left(\vec{\nabla}\cdot\vec{E}-4\pi\rho\right)+\lambda_{2}\left(\vec{\nabla}\cdot\vec{B}\right)$
(12)

is added to $L_{D}+L_{S}$, where $\lambda_{1}(\vec{x},t)$ and
$\lambda_{2}(\vec{x},t)$ are Lagrange multipliers (or auxiliary fields). The
addition of this term modifies the dynamical equations of motion. They become

$\vec{\nabla}\times\vec{E}=-\frac{\partial\vec{B}}{\partial
t}+\vec{\nabla}\lambda_{1}$ (13)

and

$\vec{\nabla}\times\vec{B}=\frac{\partial\vec{E}}{\partial
t}+4\pi\vec{j}+\vec{\nabla}\lambda_{2}\;.$ (14)

The additional terms involving $\lambda_{1}$ and $\lambda_{2}$ act as electric
and magnetic current sources for the fields. This approach yields Maxwell’s
equations if and only if these sources vanish. Taking the divergence of both
sides of these equations, and assuming charge conservation, gives

$\nabla^{2}\lambda_{1}=\nabla^{2}\lambda_{2}=0\;.$ (15)

These equations do not fix $\lambda_{1}$ or $\lambda_{2}$, so some additional
conditions are required. We are free to choose conditions that set
$\vec{\nabla}\lambda_{1}=\vec{\nabla}\lambda_{2}=0$, so that the equations of
motion are Maxwell’s equations. One possibility is to require that
$\lambda_{1}$ and $\lambda_{2}$ vanish uniformly as $r\rightarrow\infty$. This
implies that $\lambda_{1}=\lambda_{2}=0$. Another possibility is to require
that $\lambda_{1}$ and $\lambda_{2}$ are bounded everywhere. According to
Liouville’s theorem any harmonic function $\chi(\vec{x})$ that is bounded on
all of $R^{3}$ must be a constant. Since $\lambda_{1}$ and $\lambda_{2}$ are
functions of $\vec{x}$ and $t$ we will have $\lambda_{1}=C_{1}(t)$ and
$\lambda_{2}=C_{2}(t)$ implying that
$\vec{\nabla}\lambda_{1}=\vec{\nabla}\lambda_{2}=0$.

The properties of the Lagrangian and constraints under Lorentz transformations
are discussed in Appendix I. It is shown that the Lagrangian and constraints
can be written as

$L_{\alpha}=\alpha^{\mu}\left[F_{\mu\nu}\partial_{\beta}G^{\nu\beta}-G_{\mu\nu}\partial_{\beta}F^{\nu\beta}+8\pi
G_{\mu\beta}j^{\beta}\right]\;,$ (16)
$\alpha^{\mu}[\partial_{\alpha}F^{\alpha}_{\;\;\mu}+j_{\mu}]=0\;\;\;\;\;\;\;\;\;\;\;and\;\;\;\;\;\;\;\;\;\;\;\alpha^{\mu}[\partial_{\alpha}G^{\alpha}_{\;\;\mu}]=0\;,$
(17)

where $F_{\mu\nu}$ is the electromagnetic field strength tensor,
$G^{\mu\nu}=\frac{1}{2}\epsilon^{\mu\nu\alpha\beta}F_{\alpha\beta}$ is its
dual, $j^{\mu}$ is the current density and $\alpha^{\mu}=(1,\vec{0})$. It is
also shown that $L_{\alpha}$ and the constraints give Maxwell’s equations in
all inertial frames for any constant time-like four vector $\alpha^{\mu}$.
This gives a manifestly covariant formulation of the theory.

There have been previous papers that have discussed electrodynamics without
potentials. The approaches taken in these papers differ from the approach
taken here. In [5] Infeld and Plebanski formulate a Lagrangian approach to
classical electrodynamics without the use of potentials. Their action is

$S=\int[\Omega_{0}(F)+A_{\alpha}F^{\alpha\beta}_{\;\;\;\;\;\;;\beta}]\sqrt{g}d^{4}x$
(18)

where $F=-\frac{1}{4}F^{\alpha\beta}F_{\alpha\beta}$ and the $A_{\alpha}$ are
Lagrange multipliers. The field equations that follow from this action are

$F^{\alpha\beta}_{\;\;\;\;\;\;;\beta}=0$ (19)

and

$\Omega_{0F}F_{\alpha\beta}=A_{\beta,\alpha}-A_{\alpha,\beta}\;,$ (20)

where $\Omega_{0F}=\frac{d\Omega_{0}}{dF}$. These are the field equations
equations for non-linear electrodynamics introduced by Born [6] and by Born
and Infeld [7]. Maxwell’s equations are obtained when $\Omega_{0}=F$. This
approach is quite different from the one taken here. In the approach taken by
Infeld and Plebanski the field equations $\partial_{\alpha}F^{\alpha\mu}=0$,
which contain one constraint equation and and one (vector) dynamical equation,
are enforced by Lagrange multipliers. Also, from (20) it can be seen that the
Lagrange multiplier $A_{\alpha}$ is really the vector potential. In the
approach taken here only the constraint equations are enforced by Lagrange
multipliers and there is no quantity which plays the role of the vector
potential. The Hamiltonian formalism is also not investigated in [5].

In [8] Fiutak and Zukowski take the action to be (in their notation)

$W=-\frac{1}{4}\int f_{\mu\nu}f^{\mu\nu}d^{4}x$ (21)

where $f_{\mu\nu}$ is required to satisfy the “equations of constraint”

$\partial_{\alpha}f_{\beta\gamma}+\partial_{\gamma}f_{\alpha\beta}+\partial_{\beta}f_{\gamma\alpha}=0\;.$
(22)

They are, therefore, assuming half of Maxwell’s equations at this point (one
constraint equation and and one vector dynamical equation). The action is
varied with respect to $f_{\mu\nu}$ giving

$\delta_{0}W=-\frac{1}{2}\int f^{\mu\nu}\delta_{0}f_{\mu\nu}d^{4}x=0$ (23)

The $\delta_{0}f_{\mu\nu}$ are not independent but must satisfy

$\partial_{\alpha}\delta_{0}f_{\beta\gamma}+\partial_{\gamma}\delta_{0}f_{\alpha\beta}+\partial_{\beta}\delta_{0}f_{\gamma\alpha}=0\;.$
(24)

These constraints are taken into account by taking $\delta_{0}f_{\mu\nu}$ to
be of the form

$\delta_{0}f_{\mu\nu}=\partial_{\mu}\psi_{\nu}-\partial_{\nu}\psi_{\nu}$ (25)

Fiutak and Zukowski then state “Note that all we need is the form of the
variation (25), and this implies that the electromagnetic potentials are not
introduced.” (I have changed the equation number to match the number here).
This approach is also mentioned in [9].

In [10] Strazhev and Shkol’nikov develop classical and quantum theories of the
electromagnetic field without potentials. They start with the action (in their
notation)

$W=-\frac{1}{4}\int\Phi_{\mu\nu}^{2}d^{4}x\;\;\;\;\;\;\;\;\;\;\;where\;\;\;\;\;\;\;\;\;\;\;\Phi_{\mu\nu}=\partial_{[\mu}\int
F_{|\rho|\nu]}f_{\rho}(x-x^{\prime})dx^{\prime}\;,$ (26)
$f_{\mu}(x)=\frac{1}{2}\int_{0}^{\infty}[\delta(x-\xi)-\delta(x+\xi)]d\xi_{\mu}$
(27)

and $\xi_{\mu}=\xi_{\mu}(\eta)$ is a spacelike path with $\xi_{\mu}(0)=0$ and
$\xi_{\mu}(\eta)\rightarrow+\infty$ as $\eta\rightarrow+\infty$. The dynamical
variables are taken to be $F_{\mu\nu}(x)$ and they are taken to satisfy
$\partial_{\nu}\tilde{F}_{\mu\nu}=0$ where $\tilde{F}_{\mu\nu}$ is the dual of
$F_{\mu\nu}$. They have, therefore, assumed half of Maxwell’s equations by
imposing this condition on $\tilde{F}_{\mu\nu}$. They then add
$\int\lambda_{\mu}\partial_{\nu}\tilde{F}_{\mu\nu}d^{4}x$ to the action and
state “In order to avoid misunderstandings we stress that the Lagrange
multipliers $\lambda_{\mu}$ are not varied, since the relations (5) are, by
definition, known independently of the variational procedure. Hence the
$\lambda_{\mu}$ cannot be treated as auxiliary dynamical variables.” Here (5)
is the equation $\partial_{\nu}\tilde{F}_{\mu\nu}=0$. The variation of the
action with respect to $F_{\mu\nu}$ using $\partial_{\nu}\tilde{F}_{\mu\nu}=0$
does not give the remaining Maxwell equations but instead gives

$\int
dx^{\prime}[\partial^{{}^{\prime}}_{\nu}F_{\alpha\nu}(x^{\prime})]f_{\beta}(x-x^{\prime})-\partial^{{}^{\prime}}_{\nu}F_{\beta\nu}(x^{\prime})f_{\alpha}(x-x^{\prime})]+\frac{1}{2}\epsilon_{\mu\nu\alpha\beta}\partial_{\mu}\lambda_{\nu}=0\;.$
(28)

This appears to be very different than the remaining set of Maxwell equations:
$\partial_{\nu}F_{\mu\nu}=0$. Maxwell’s equations can be obtained by
differentiating (28) with respect to $x_{\beta}$. However (28) cannot be
expected to be equivalent to Maxwell’s theory, since we need to differentiate
(28) to obtain Maxwell’s equations. For example, it is well known that one can
obtain wave equations for the electric and magnetic field by differentiating
Maxwell’s equations. However, not all solutions to the wave equations will
satisfy Maxwell’s equations. The solution space of these wave equations is
larger than that of Maxwell’s equations. The Hamiltonian formalism is not
investigated in this paper.

In [11] Mandelstam develops a quantum theory of the electromagnetic field
coupled to a complex scalar field without potentials. He considers the
electromagnetic field coupled to a complex scalar field and takes the gauge
invariant variables $\Phi(x,P)$, $\Phi^{*}(x,P)$ and $F_{\mu\nu}(x)$ are taken
as the field variables where

$\Phi(x,P)=\phi(x)\exp\left\\{-ie\int_{-\infty}^{x}d\xi_{\mu}A_{\mu}(\xi)\right\\}$
(29)

and the integral is taken over the spacelike path $P$. The Lagrangian is given
by

$L=-\partial_{\mu}\Phi^{*}\partial_{\mu}\Phi-m^{2}\Phi^{*}\Phi-\frac{1}{4}\left\\{F_{\mu\nu}\right\\}^{2}$
(30)

where $\partial_{\mu}\Phi$ is the gauge invariant derivative of $\Phi$. This
approach differs significantly from the approach taken in this paper.

##  3 The Hamiltonian Formulation

In this section I will consider the free electromagnetic field (i.e. $\rho=0$
and $\vec{j}=0$). The Hamiltonian for the dynamical equations will be derived
first and the constraint equations will then be included. The canonical
momenta are given by

$\vec{\Pi}_{E}=\frac{\partial L_{D}}{\partial\dot{\vec{E}}}=-\vec{B}$ (31)

and

$\vec{\Pi}_{B}=\frac{\partial L_{D}}{\partial\dot{\vec{B}}}=\vec{E}.$ (32)

We therefore have the primary constraints

$\vec{\phi}_{1}=\vec{\Pi}_{E}+\vec{B}\approx 0\hskip 28.90755ptand\hskip
28.90755pt\vec{\phi}_{2}=\vec{\Pi}_{B}-\vec{E}\approx 0\;,$ (33)

where $\approx$ denotes a weak equality which can only be imposed after the
Poisson brackets have been evaluated. These constraints satisfy

$\\{\phi_{1n}(\vec{x}),\phi_{1m}(\vec{y})\\}=\\{\phi_{2n}(\vec{x}),\phi_{2m}(\vec{y})\\}=0$
(34)

and

$\\{\phi_{1n}(\vec{x}),\phi_{2m}(\vec{y})\\}=2\delta_{mn}\delta^{3}(\vec{x}-\vec{y})\;.$
(35)

Thus, these constraints are second class.

The canonical Hamiltonian $H_{C}=\int
d^{3}x[\vec{\Pi}_{E}\cdot\dot{\vec{E}}+\vec{\Pi}_{B}\cdot\dot{\vec{B}}-L_{D}]$
is given by

$H_{C}=\int
d^{3}x\left[\vec{\phi}_{1}\cdot\dot{\vec{E}}+\vec{\phi}_{2}\cdot\dot{\vec{B}}-\vec{E}\cdot\left(\vec{\nabla}\times\vec{E}\right)-\vec{B}\cdot\left(\vec{\nabla}\times\vec{B}\right)\right]$
(36)

and the total Hamiltonian is given by

$H_{T}=-\int
d^{3}x\left[\vec{E}\cdot\left(\vec{\nabla}\times\vec{E}\right)+\vec{B}\cdot\left(\vec{\nabla}\times\vec{B}\right)-\vec{u}_{1}\cdot\vec{\phi}_{1}-\vec{u}_{2}\cdot\vec{\phi}_{2}\right]\;,$
(37)

where $\vec{u}_{1}$ and $\vec{u}_{2}$ are undetermined parameters.

For consistency we require that

$\dot{\vec{\phi}_{1}}=\\{\vec{\phi}_{1},H_{T}\\}\approx 0\hskip
28.90755ptand\hskip
28.90755pt\dot{\vec{\phi}}_{2}=\\{\vec{\phi}_{2},H_{T}\\}\approx 0\;.$ (38)

These two equations give

$\vec{u}_{1}=\vec{\nabla}\times\vec{B}$ (39)

and

$\vec{u}_{2}=-\vec{\nabla}\times\vec{E}\;.$ (40)

Substituting these expressions for $\vec{u}_{1}$ and $\vec{u}_{2}$ into (37)
gives

$H_{T}=\int
d^{3}x\left[\vec{\Pi}_{E}\cdot\left(\vec{\nabla}\times\vec{B}\right)-\vec{\Pi}_{B}\cdot\left(\vec{\nabla}\times\vec{E}\right)\right]\;.$
(41)

It is easy to see that
$\dot{\vec{E}}=\\{\vec{E},H_{T}\\}=\vec{\nabla}\times\vec{B}$ and
$\dot{\vec{B}}=\\{\vec{B},H_{T}\\}=-\vec{\nabla}\times\vec{E}$.

The constraints

$\phi_{3}=\vec{\nabla}\cdot\vec{E}\approx 0\hskip 28.90755ptand\hskip
28.90755pt\phi_{4}=\vec{\nabla}\cdot\vec{B}\approx 0$ (42)

must also be included. It is well known that these constraints are preserved
under time evolution so there are no additional constraints. It is easy to
show that

$\\{\phi_{3}(x),\vec{\phi}_{1}(y)\\}=\vec{\nabla}_{x}\delta^{3}(\vec{x}-\vec{y})\;,\hskip
36.135pt\\{\phi_{3}(x),\vec{\phi}_{2}(y)\\}=0\;,$ (43)

and

$\\{\phi_{4}(x),\vec{\phi}_{1}(y)\\}=0\;,\hskip
36.135pt\\{\phi_{4}(x),\vec{\phi}_{2}(y)\\}=\vec{\nabla}_{x}\delta^{3}(\vec{x}-\vec{y})\;.$
(44)

the constraints $\\{\vec{\chi}_{1},\;\vec{\chi_{2}}\;,\chi_{3}\;,\chi_{4}\\}$
are, therefore, a set of second class constraints.

Systems involving second class constraints can be treated using the Dirac
bracket. Before introducing the Dirac bracket it will be convenient to relabel
the constraints as follows

$\eta_{1}=\phi_{1x},\;\;\;\;\;\;\;\eta_{2}=\phi_{1y},\;\;\;\;\;\;\;\eta_{3}=\phi_{1z},\;\;\;\;\;\;\;\eta_{4}=\phi_{2x}$
(45)

and

$\eta_{5}=\phi_{2y},\;\;\;\;\;\;\;\eta_{6}=\phi_{2z},\;\;\;\;\;\;\;\eta_{7}=\phi_{3},\;\;\;\;\;\;\;\eta_{8}=\phi_{4}$
(46)

and to define the matrix

$\Delta_{m\vec{x},n\vec{y}}=\\{\eta_{m}(\vec{x}),\eta_{n}(\vec{y})\\}.$ (47)

The Dirac bracket is defined as

$\\{A,B\\}_{D}=\\{A,B\\}-\Sigma_{m,n}\int
d^{3}xd^{3}y\\{A,\eta_{m}(\vec{x})\\}C_{m\vec{x},n\vec{y}}\\{\eta_{n}(\vec{y}),B\\}\;,$
(48)

where $C_{m\vec{x};n\vec{y}}$ is the inverse of $\Delta_{m\vec{x},n\vec{y}}$ :

$\Sigma_{k}\int
d^{3}z\;C_{m\vec{x},k\vec{z}}\Delta_{k\vec{z},n\vec{y}}=\delta_{mn}\delta^{3}(\vec{x}-\vec{y})\;.$
(49)

The non-zero elements of $C_{m\vec{x},n\vec{y}}$ are given by (see Appendix
II)

$C_{1\vec{x},4\vec{y}}=-\frac{1}{2}\delta^{3}(\vec{x}-\vec{y})+\frac{1}{8\pi|\vec{x}-\vec{y}|^{3}}-\frac{3(x^{1}-y^{1})^{2}}{8\pi|\vec{x}-\vec{y}|^{5}}\;,$
(50)
$C_{2\vec{x},5\vec{y}}=-\frac{1}{2}\delta^{3}(\vec{x}-\vec{y})+\frac{1}{8\pi|\vec{x}-\vec{y}|^{3}}-\frac{3(x^{2}-y^{2})^{2}}{8\pi|\vec{x}-\vec{y}|^{5}}\;,$
(51)
$C_{3\vec{x},6\vec{y}}=-\frac{1}{2}\delta^{3}(\vec{x}-\vec{y})+\frac{1}{8\pi|\vec{x}-\vec{y}|^{3}}-\frac{3(x^{3}-y^{3})^{2}}{8\pi|\vec{x}-\vec{y}|^{5}}\;,$
(52)
$C_{1\vec{x},5\vec{y}}=C_{2\vec{x},4\vec{y}}=-\frac{3(x^{1}-y^{1})(x^{2}-y^{2})}{8\pi|\vec{x}-\vec{y}|^{5}}\;,$
(53)
$C_{1\vec{x},6\vec{y}}=C_{3\vec{x},4\vec{y}}=-\frac{3(x^{1}-y^{1})(x^{3}-y^{3})}{8\pi|\vec{x}-\vec{y}|^{5}}\;,$
(54)
$C_{2\vec{x},6\vec{y}}=C_{3\vec{x},5\vec{y}}=-\frac{3(x^{2}-y^{2})(x^{3}-y^{3})}{8\pi|\vec{x}-\vec{y}|^{5}}\;,$
(55)
$C_{1\vec{x},7\vec{y}}=C_{4\vec{x},8\vec{y}}=\frac{(x^{1}-y^{1})}{4\pi|\vec{x}-\vec{y}|^{3}}\;,$
(56)
$C_{2\vec{x},7\vec{y}}=C_{5\vec{x},8\vec{y}}=\frac{(x^{2}-y^{2})}{4\pi|\vec{x}-\vec{y}|^{3}}\;,$
(57)
$C_{3\vec{x},7\vec{y}}=C_{6\vec{x},8\vec{y}}=\frac{(x^{3}-y^{3})}{4\pi|\vec{x}-\vec{y}|^{3}}\;,$
(58) $C_{7\vec{x},8\vec{y}}=\frac{1}{2\pi|\vec{x}-\vec{y}|}\;.$ (59)

plus terms related by antisymmetry:
$C_{m\vec{x};n\vec{y}}=-C_{n\vec{y};m\vec{x}}$. The advantage of using Dirac
brackets instead of Poisson brackets is that the constraints can be set to
zero strongly. The Dirac brackets involving $\vec{E}$ and $\vec{B}$ are

$\\{E_{i}(\vec{x}),E_{j}(\vec{y})\\}_{D}=\\{B_{i}(\vec{x}),B_{j}(\vec{y})\\}_{D}=0$
(60)

and

$\\{E_{i}(\vec{x}),\Pi_{Ej}(\vec{y})\\}_{D}=\frac{1}{2}\left[\delta^{3}(\vec{x}-\vec{y})\delta_{ij}+\frac{\partial^{2}}{\partial
x^{i}\partial x^{j}}\left(\frac{1}{4\pi|\vec{x}-\vec{y}|}\right)\right]\;.$
(61)

It is interesting to note that in the standard approach in the Coulomb gauge
one finds [12] that $A_{i}(\vec{x})$ and $\Pi^{j}(\vec{x})$ satisfy the same
Dirac bracket up to a factor of 1/2. The other Dirac brackets involving the
momenta can be found from these using the constraints (i.e.
$\\{E_{i}(\vec{x}),B_{j}(\vec{y})\\}_{D}=-\\{E_{i}(\vec{x}),\Pi_{Ej}(\vec{y})\\}_{D}$).

The constraints $\vec{\nabla}\cdot\vec{E}=0$ and $\vec{\nabla}\cdot\vec{B}=0$
imply that the longitudinal components of $\vec{E}$ and $\vec{B}$ vanish. The
constraints $\vec{\phi}_{1}=0$ and $\vec{\phi}_{2}=0$ then imply that the
longitudinal components of $\vec{\Pi}_{E}$ and $\vec{\Pi}_{B}$ vanish. The
four independent degrees of freedom can therefore be taken to be the
transverse components of $\vec{E}$ and $\vec{\Pi}_{E}$.

One could also have started with $L_{D}+L_{C}$ and obtained similar results
with two additional constraints, $\Pi_{\lambda_{1}}\approx 0$ and
$\Pi_{\lambda_{2}}\approx 0$.

##  4 Conclusion

In this paper I found a Lagrangian that gives the dynamical Maxwell equations
and included the constraint equations using Lagrange multipliers. In the
Hamiltonian formalism one finds that the canonical momenta $\vec{\Pi}_{E}$ and
$\vec{\Pi}_{B}$ are constrained giving 6 second class constraints at each
point in space. Gauss’s law and $\vec{\nabla}\cdot\vec{B}=0$ can than be added
in as additional constraints giving a total of 8 constraints. This leaves 4
independent degrees of freedom, as expected.

## Acknowledgements

This research was supported by the Natural Sciences and Engineering Research
Council of Canada.

##  5 Appendix I

In this appendix I examine the transformation properties of $L=L_{D}+L_{S}$
and the constraints under Lorentz transformations and find covariant
expressions for these quantities. For notational simplicity define

$\vec{a}=\frac{\partial\vec{E}}{\partial
t}-\vec{\nabla}\times\vec{B}+4\pi\vec{j}$ (62)

and

$\vec{b}=\frac{\partial\vec{B}}{\partial t}+\vec{\nabla}\times\vec{E}$ (63)

Now consider the constraints. They can be written as

$\vec{\nabla}\cdot\vec{E}-4\pi\rho=\partial_{\alpha}F^{\alpha}_{\;\;0}+4\pi
j_{0}=0$ (64)

and

$\vec{\nabla}\cdot\vec{B}=\partial_{\alpha}G^{\alpha}_{\;\;0}=0\;,$ (65)

where $F_{\mu\nu}$ is the electromagnetic field strength tensor,
$G^{\mu\nu}=\frac{1}{2}\epsilon^{\mu\nu\alpha\beta}F_{\alpha\beta}$ is its
dual and $j^{\mu}$ is the current density. Consider the constraints in the
inertial frame $\bar{S}$. Under the Lorentz transformation
$\Lambda_{\mu}^{\;\;\nu}$ from $\bar{S}$ to the inertial frame $S$

$\bar{\partial}_{\alpha}\bar{F}^{\alpha}_{\;\;0}+4\pi\bar{j}_{0}=\Lambda_{0}^{\;\;\mu}(\partial_{\alpha}F^{\alpha}_{\;\;\mu}+4\pi
j_{\mu})=0$ (66)

and

$\bar{\partial}_{\alpha}\bar{G}^{\alpha}_{\;\;0}=\Lambda_{0}^{\;\;\mu}\partial_{\alpha}G^{\alpha}_{\;\;\mu}=0\;.$
(67)

The constraints in $S$ are then

$\Lambda_{0}^{\;\;\mu}(\partial_{\alpha}F^{\alpha}_{\;\;\mu}+4\pi
j_{\mu})=0\;\;\;\;\;\;\;\;\;\;\;and\;\;\;\;\;\;\;\;\;\;\;\Lambda_{0}^{\;\;\mu}\partial_{\alpha}G^{\alpha}_{\;\;\mu}=0\;.$
(68)

They can also be written as

$\vec{\nabla}\cdot\vec{E}-4\pi\rho+\vec{\Lambda}\cdot\vec{a}=0$ (69)

and

$\vec{\nabla}\cdot\vec{B}+\vec{\Lambda}\cdot\vec{b}=0\;,$ (70)

where $\vec{\Lambda}_{k}=\frac{\Lambda_{0}^{\;\;k}}{\Lambda_{0}^{\;\;0}}$.

The Lagrangian $L=L_{D}+L_{S}$ can be written as

$L=F_{0\mu}\partial_{\alpha}G^{\mu\alpha}-G_{0\mu}\partial_{\alpha}F^{\mu\alpha}+8\pi
G_{0\mu}j^{\mu}\;.$ (71)

Consider the Lagrangian $\bar{L}=\bar{L}_{d}+\bar{L}_{S}$ in the inertial
frame $\bar{S}$. Under the Lorentz transformation $\Lambda_{\mu}^{\;\;\nu}$
from $\bar{S}$ to the inertial frame $S$

$\bar{F}_{0\mu}\bar{\partial}_{\alpha}\bar{G}^{\mu\alpha}-\bar{G}_{0\mu}\partial_{\alpha}\bar{F}^{\mu\alpha}+8\pi\bar{G}_{0\mu}\bar{j}^{\mu}=\Lambda_{0}^{\;\;\nu}(F_{\nu\mu}\partial_{\alpha}G^{\mu\alpha}-G_{\nu\mu}\partial_{\alpha}F^{\mu\alpha}+8\pi
G_{\nu\mu}j^{\mu})\;.$ (72)

The Lagrangian in $S$ is then given by.

$L=\Lambda_{0}^{\;\;\mu}(F_{\mu\nu}\partial_{\alpha}G^{\nu\alpha}-G_{\mu\nu}\partial_{\alpha}F^{\nu\alpha}+8\pi
G_{\mu\nu}j^{\nu})\;.$ (73)

Expressing this Lagrangian in terms of $\vec{E}$ and $\vec{B}$ gives

$L=\Lambda_{0}^{\;\;0}\left\\{\vec{b}\cdot\vec{E}-\left(\vec{a}+4\pi\vec{j}\right)\cdot\vec{B}+\vec{\Lambda}\cdot\left[\vec{E}(\vec{\nabla}\cdot\vec{B})-\vec{B}(\vec{\nabla}\cdot\vec{E})+\vec{B}\times\vec{b}+\vec{E}\times\left(\vec{a}+4\pi\vec{j}\right)+8\pi\rho\vec{B}\right]\right\\}.$
(74)

The field equations, using the constraints (69) and (70), generated by this
Lagrangian are

$\vec{a}+\vec{\Lambda}\times\vec{b}-(\vec{a}\cdot\vec{\Lambda})\vec{\Lambda}=0$
(75)

and

$\vec{b}-\vec{\Lambda}\times\vec{a}-(\vec{b}\cdot\vec{\Lambda})\vec{\Lambda}=0\;.$
(76)

Taking the inner product of both sides of these equations with $\vec{\Lambda}$
gives

$(1-\Lambda^{2})(\vec{\Lambda}\cdot\vec{a})=0\;\;\;\;\;\;\;\;\;\;\;and\;\;\;\;\;\;\;\;\;\;\;(1-\Lambda^{2})(\vec{\Lambda}\cdot\vec{b})=0\;,$
(77)

where $\Lambda^{2}=\vec{\Lambda}\cdot\vec{\Lambda}$. Thus, for
$\Lambda^{2}\neq 1$, which will be the case for Lorentz transformations, we
have $\vec{\Lambda}\cdot\vec{a}=\vec{\Lambda}\cdot\vec{b}=0$ Substituting (76)
into (75) and using $\vec{\Lambda}\cdot\vec{a}=\vec{\Lambda}\cdot\vec{b}=0$
gives

$\vec{a}=\frac{\partial\vec{E}}{\partial
t}-\vec{\nabla}\times\vec{B}+4\pi\vec{j}=0\;.$ (78)

A similar argument gives

$\vec{b}=\frac{\partial\vec{B}}{\partial t}+\vec{\nabla}\times\vec{E}=0\;.$
(79)

In the above derivation of Maxwell’s equations I used $|\vec{\Lambda}|\neq 1$,
which is equivalent to
$-(\Lambda_{0}^{\;\;0})^{2}+\Sigma_{k}\Lambda_{0}^{\;\;k}\Lambda_{0}^{\;\;k}\neq
0$, and $\Lambda_{0}^{\;\;0}\neq 0$, both of which are true for Lorentz
transformations. It is easy to see that the Lagrangian and constraints will
still produce Maxwell’s equations if $\Lambda_{0}^{\;\;\mu}$ in (68) and (73)
is replaced by any constant four-vector $\alpha^{\mu}$ satisfying
$\alpha^{\mu}\alpha_{\mu}\neq 0$ and $\alpha^{0}\neq 0$. The condition
$\alpha^{\mu}\alpha_{\mu}\neq 0$ constrains $\alpha^{\mu}$ to be non-null. If
$\alpha^{\mu}$ is a space-like four vector there exists an inertial frame in
which $\alpha^{0}=0$. Therefore, if the Lagrangian and constraints are to
produce Maxwell’s theory in all inertial frames of reference $\alpha^{\mu}$
must be a time-like four-vector. Thus, the Lagrangian

$L_{\alpha}=\alpha^{\mu}\left[F_{\mu\nu}\partial_{\beta}G^{\nu\beta}-G_{\mu\nu}\partial_{\beta}F^{\nu\beta}+8\pi
G_{\mu\beta}j^{\beta}\right]$ (80)

and the constraints

$\alpha^{\mu}[\partial_{\beta}F^{\beta}_{\;\;\mu}+4\pi
j_{\mu}]=0\;\;\;\;\;\;\;\;\;\;\;and\;\;\;\;\;\;\;\;\;\;\;\alpha^{\mu}[\partial_{\beta}G^{\beta}_{\;\;\mu}]=0$
(81)

are Lorentz scalars and give Maxwell equations for any constant time-like
four-vector $\alpha^{\mu}$. For simplicity I have imposed the constraints
outside the variational principle so that the Lagrange multipliers do not
appear in the equations for the fields.

##  6 Appendix II

In this appendix I derive a few of the $C_{i\vec{x},j\vec{y}}$ to illustrate
how the calculations are done. I will assume that $C_{n\vec{x},n\vec{y}}=0$
(no sum on $n$) and that $C_{m\vec{x},n\vec{y}}$ is a function of
$(x^{1}-y^{1},x^{2}-y^{2},x^{3}-y^{3})$. These assumptions will be justified
by the success of the calculations. Some of the equations for
$C_{i\vec{x},j\vec{y}}$ are differential equations, so boundary conditions
need to be imposed. I will take these conditions to be

$C_{i\vec{x},j\vec{y}}\rightarrow 0\hskip 21.68121ptas\hskip
21.68121pt|\vec{x}-\vec{y}|\rightarrow\infty\;.$ (82)

Define $P_{i\vec{x},j\vec{y}}$ by

$P_{m\vec{x},n\vec{y}}=\Sigma_{k}\int
d^{3}z\;C_{m\vec{x},k\vec{z}}\Delta_{k\vec{z},n\vec{y}}\;.$ (83)

Consider $m=7$ and $n=1$

$P_{7\vec{x},1\vec{y}}=C_{7\vec{x},4\vec{y}}=0$ (84)

and $m=1$ and $n=4$

$P_{1\vec{x},4\vec{y}}=-\frac{\partial}{\partial
y^{1}}C_{1\vec{x},8\vec{y}}=0\;.$ (85)

This implies that $C_{1\vec{x},8\vec{y}}$ is independent of $x^{1}-y^{1}$. The
boundary condition (82) then gives $C_{1\vec{x},8\vec{y}}=0$.

Now consider

$P_{1\vec{x},1\vec{y}}=-2C_{1\vec{x},4\vec{y}}-\frac{\partial}{\partial
y^{1}}C_{1\vec{x},7\vec{y}}=\delta^{3}(\vec{x}-\vec{y})\;,$ (86)
$P_{1\vec{x},2\vec{y}}=-2C_{1\vec{x},5\vec{y}}-\frac{\partial}{\partial
y^{2}}C_{1\vec{x},7\vec{y}}=0\;,$ (87)
$P_{1\vec{x},3\vec{y}}=-2C_{1\vec{x},6\vec{y}}-\frac{\partial}{\partial
y^{3}}C_{1\vec{x},7\vec{y}}=0\;,$ (88)

and

$P_{1\vec{x},8\vec{y}}=-\frac{\partial}{\partial
y^{1}}C_{1\vec{x},4\vec{y}}-\frac{\partial}{\partial
y^{2}}C_{1\vec{x},5\vec{y}}-\frac{\partial}{\partial
y^{3}}C_{1\vec{x},6\vec{y}}=0\;.$ (89)

Substituting (86), (87) and (88) into (89) gives

$\nabla^{2}_{y}C_{1\vec{x},7\vec{y}}=-\frac{\partial}{\partial
y^{1}}\delta^{3}(\vec{x}-\vec{y})\;.$ (90)

This equation has the solution

$C_{1\vec{x},7\vec{y}}=\frac{(x^{1}-y^{1})}{4\pi|\vec{x}-\vec{y}|^{3}}\;.$
(91)

Equations (87) and (88) now give

$C_{1\vec{x},5\vec{y}}=-\frac{3(x^{1}-y^{1})(x^{2}-y^{2})}{8\pi|\vec{x}-\vec{y}|^{5}}$
(92)

and

$C_{1\vec{x},6\vec{y}}=-\frac{3(x^{1}-y^{1})(x^{3}-y^{3})}{8\pi|\vec{x}-\vec{y}|^{5}}\;.$
(93)

## References

  * [1] J.D. Jackson, Classical electrodynamics (3rd edition, Wiley, 1999).
  * [2] L.D. Landau and E.M. Lifshitz, The classical theory of fields (Pergamon, Oxford, 1962).
  * [3] P.A.M. Dirac, Lectures on Quantum Mechanics (Dover Publications, Mineola, New York, 1964), pp 1-24.
  * [4] P.A.M. Dirac, Can. J. Math. 2, 129 (1950).
  * [5] L. Infeld and J. Plebanski, Proc. Roy. Soc. A222, 224 (1954)
  * [6] M. Born, Proc. Roy. Soc. A143, 410 (1934)
  * [7] M. Born and L. Infeld, Proc. Roy. Soc. A144, 425 (1934)
  * [8] J. Fiutak and M. Zukowski, J. Phys. A14, 3229 (1981)
  * [9] V. I. Ogievetskii and I. V. Polubarinov, Zh. Eksp. Teor. Fiz. 43, 1365 (1962), JETP 16, 969 (1963)
  * [10] V. I. Strazhev and P. L. Shkol’nikov, Sov. Phys. J. 32, 378 (1989)
  * [11] S. Mandelstam, Annals Phys. 19, 1 (1962)
  * [12] Steven Weinberg, The Quantum Theory of Fields: Volume I Foundations (Cambridge University Press 1995), p 348.