Commit a297acf9 authored by Antonin Dudermel's avatar Antonin Dudermel
Browse files

msb version 2

parent e929c9c9
......@@ -166,17 +166,9 @@ By case disjunction:
\(F^{(n - m + t)}(\bot) \leq F^{(t)}(\top)\), particularily for \(t = m\)
\end{itemize}\qed
\section*{Copyright}
\section*{Aknowledgement}
\label{sec:copy}
\begin{itemize}
\item The picture of the \figref{fig:zybo} is taken from Digilent Website:\\
https://reference.digilentinc.com/\_media/reference/programmable-logic/zybo/zybo-2.png
\item The \figref{fig:bell-impl} is made by Emma Neiss
\item Figures \ref{fig:difffix},\ref{fig:sfix} are made by Florent de Dinechin
\end{itemize}
Figures \ref{fig:bell-impl}, \ref{fig:difffix},\ref{fig:sfix} are made by
Florent de Dinechin
\end{document}
%%% Local IspellDict: en
......@@ -5,7 +5,7 @@
\label{sec:bounding-signals}
The first problem one has to handle with fixed-point arithmetic is overflow. As
Fixed-point numbers behave exactly like integers, on \(\sfix(0,-3)\) we have
Fixed-point numbers are scaled integers, on \(\sfix(0,-3)\) we have
\(0.875 +_{(0,-3)} 0.125 = -1\). To avoid this, it is necessary and sufficient
to find for each signal a big enough MSB. Formally, we want
......@@ -41,7 +41,7 @@ operation \(\op \in \Ops\), a new operation \(\op\) on intervals:
X \op Y: \forall x\in X, y\in Y, x\op y \in X\op Y
\end{equation*}
Note that once again, there is a minimal interval wich is
Note that once again, there is a minimal interval which is
\begin{equation}
\Cvx(\ens{x \op y, x\in X, y\in Y})\label{eq:inter}
\end{equation}
......@@ -51,7 +51,7 @@ where \(\Cvx\) is the the convex closure.
As this interval can sometimes be very costly to compute, relaxed versions can
be chosen for more efficiency. For instance, defining \(\cos([a, b]) = [0, 1]\)
for any interval \(X\) is perfectly correct. It is however very important that
every operation on intervals remains increasing for \(\subseteq\) {\it ie} for
every operation on intervals remains increasing for \(\subseteq\) {\it pie} for
\(X \subseteq X', Y \subseteq Y'\) we must have
\( X \diamond Y \subseteq X' \diamond Y'\).
......@@ -59,8 +59,8 @@ Many tools already exist for this purpose, such as MPFI\cite{inriaforgempfi}, a
C library which was the reference for the IEEE, or Gappa\cite{gappa}, an
automatic proof generator for arithmetic properties. However, a design choice at
the core of Faust is that the compiler must be as easy as possible to install,
and therfore contain as few dependancies as possible\footnote{Being a language
for musicians and soundmakers, but also a pedagogical tool for discovering
and therefore contain as few dependencies as possible\footnote{Being a language
for musicians and sound-makers, but also a pedagogical tool for discovering
digital music, Faust is targeting mainly people with few computer
knowledge}. Fortunately, there is already a small built-in interval arithmetic
handling in the Faust compiler that can derive intervals for expressions (see
......@@ -74,11 +74,11 @@ library.
{graph}{0.5}{../code/plus-sig.png}
Yet, interval arithmetic is not miraculous. A first
drawback is that interval arithmetic does not handle corellations very
drawback is that interval arithmetic does not handle correlations very
well. Consider for instance the expression \(x - x\), for a variable
\(x \in [-1,1]\). The wanted interval is obviously \([0, 0]\) as \(x - x =
0\). But, following from the equation \eqref{eq:inter}, we have
\(X - X = [-1, 1] - [-1, 1] = [-2, 2]\). Similarily, we want \(x * x\) to be
\(X - X = [-1, 1] - [-1, 1] = [-2, 2]\). Similarly, we want \(x * x\) to be
positive, but \(X * X = [-1, 1]\).
A solution to solve these problems is to rewrite arithmetic expressions to find
......@@ -86,14 +86,15 @@ easier expressions that will be more interval-friendly, so \(x - x\) can be
transformed into \(0\) and \(x * x\) into \(x^2\), both of which are trivial in
interval arithmetic. This is already done by Faust, because each expression is
rewritten into a normal form. But as this normal form is arbitrary, there is no
guarantee that Faust normal form will be the more interval friendly. A simple
correlation example on which the normal form fails is the fractional part.
guarantee that Faust normal form will be the most interval friendly version. A
simple correlation example on which the normal form fails is the fractional
part.
\begin{figure}
\centering
\includegraphics[width=0.5\textwidth]{../code/frac-block.png}
\centering
\smallcode{../code/frac.dsp}\\
\includegraphics[width=0.3\textwidth]{../code/frac-block.png}
\caption{Fractional part, as written in the standard library}
\smallcode{../code/frac.dsp}
\end{figure}
\lstset{language=fortran}
......@@ -137,7 +138,7 @@ Results:
z in [-1, 1]
\end{lstlisting}
The trick here is to use the canonical form of the degree 2 polynom, from which
The trick here is to use the canonical form of the degree 2 polynomial, from which
the extremum is obvious. \(x (1 - x) = - (x - 0.5)^2 + 0.25\). We can give this
hint to Gappa using the following syntax
......@@ -159,17 +160,16 @@ The bell also raises an instructive example with the implementation of
\(s(t) =\cos(\omega t)e^{-t/\tau}\) using a biquad. Looking at the formula, it
is obvious that \(|s(t)| \leq 1\), but except from very naive and costly
implementations, this formula does not appear in code and for the biquad, the
programmer implicitely relies on the second-order differential equation verified
by the formula to implement it efficiently. This is one of the main problems in
program analysis. The programmers don't tell to the computer what to do, but
only how to do it. Knuth developped this idea with his paradigm of literate
programing, arguing that the comments should be more important than the code in
itself\footnote{See \TeX{} example
\url{https://www-cs-faculty.stanford.edu/~knuth/programs/tex.web}}.
programmer implicitly relies on the CCDE verified by the formula to implement
it efficiently. This is one of the main problems in program analysis. The
programmers don't tell to the computer what to do, but only how to do it. Knuth
developed this idea with his paradigm of literate programming, arguing that the
comments should be more important than the code in itself\footnote{See \TeX{}
example \url{https://www-cs-faculty.stanford.edu/~knuth/programs/tex.web}}.
\lstset{language=c}
\subsection{The (possibly Turing-Complete) problem with Feedback Loops}
\subsection{The (possibly undecidable) problem with Feedback Loops}
\label{sec:problem-with-loops}
\subsubsection{Lattices in compilation}
......@@ -179,11 +179,11 @@ Another problem with interval arithmetic appears when we apply it to programs
and not only to expressions. Indeed, program have loops, so the interval
attached to a variable may vary from one iteration to another, and even worse:
the value of a variable may depend on its value at a previous iteration. To
handle this, compilers apply a more general approach based on
lattices\cite{aho2007compilersprinciples}.
handle this, compilers apply a an \emph{abstract analysis} based on
lattices\cite{CousotCousot77-1}.
\begin{defin}[Lattice]
Given partially-ordered set \((L, <)\) is a \emph{lattice} if for any pair
A partially-ordered set \((L, <)\) is a \emph{lattice} if for any pair
\(a, b \in L\), \(\ens{a,b}\) has
\begin{itemize}
\item a least upper bound called \emph{join} and denoted \(a \vee b\)
......@@ -195,20 +195,19 @@ lattices\cite{aho2007compilersprinciples}.
\end{defin}
Real interval with the relation \(\subseteq\) define a bounded lattice with
the two operations \(\vee \wedge\) defined for all \(I, J\) by
the two operations \(\vee \wedge\) defined for all reals \(a,b,c,d\) by
\begin{itemize}
\item \(I \vee J = \Cvx(I \bigcup J)\)
\item \(I \wedge J = I \bigcap J\)
\item \([a, b] \vee [c, d] = [\min(a,c), \max(b,d)]\)
\item \([a, b] \wedge [c, d] = [\max(a,c), \min(b,d)] \)
\item \(\top = \IntR\)
\item \(\bot = \emptyset\)
\end{itemize}
As the cartesian product of two lattices remains a lattice, we can extend the
As the Cartesian product of two lattices remains a lattice, we can extend the
lattice structure to the product of an arbitrary number \(n\) of intervals (a
cuboid) on \(\Real^n\).
In signal processing, the equivalent of loops are feedback loops. They enable to
define a signal \(s(t)\) depending on its value at a previous time \(s(t-1)\),
just like recurrent sequences. To define the sequence characterized by
......@@ -217,11 +216,6 @@ just like recurrent sequences. To define the sequence characterized by
\figref{fig:ramp}. Feedback loops are also the core of digital linear
time-invariant filters, such as the biquads used in the bell.
Furthermore, as every signal has value zero before the time zero, we have to
consider the lattice of intervals containing zero, with \(\bot = [0]\). More
details can be found in annex \ref{sec:lattice-signals}.
To transpose interval arithmetic into signals, we need to attach to a signal
\(s(t)\) an interval \(S'(t)\) for each time \(t\), and we need to ensure that
\(S'(t)\) contains \(s(t)\). To derive a single interval \(\bar{S'}\) from the
......@@ -232,9 +226,9 @@ To transpose interval arithmetic into signals, we need to attach to a signal
To compute the \(S'(t)\), we could simply transform the recurrence relation
\[
\recSig{i}{s(i) = 0}{s(i+1) = f(s'(i)\dots s'(i-k))}
\recSig{i}{s(i) = 0}{s(i+1) = f(s(i)\dots s(i-k))}
\text{into}
\recSig{i}{S(i) = 0}{S(i+1) = f(S'(i)\dots S'(i-k))}
\recSig{i}{S'(i) = 0}{S'(i+1) = f(S'(i)\dots S'(i-k))}
\]
Note that for any function \(f\) on signals, the function \(f\) on interval is
......@@ -248,14 +242,13 @@ all the \(S'(t)\) to compute \(\bar{S}'\). Using instead the relation
\[\recSig{i}{S(i) = [0]}{S(i+1) = S(i) \vee f(S(i)\dots S(i))}\]
can lead to overapproxmations, but ensure that the \(S(t)\) are increasing and
can lead to over-approximations, but ensure that the \(S(t)\) are increasing and
reduce the recurrence relation to an order 1 recurrence and we keep the fact
that \(s(t) \in S(t)\) for any \(t\). Furthermore, writing
\(F: I \to f(I\dots I)\), we can rewrite the recurrence relation as
\[\recSig{i}{S(i) = \bot}{S(i+1) = S(i)\vee F(S(i))}\]
and then by induction
\begin{equation}
......@@ -263,18 +256,17 @@ and then by induction
\end{equation}
with \(F^{(t)}\) denotes \(t\) successive applications of \(F\), and \(F^{(0)}\)
is the identity. This last formula is interesting because it appears in Kleene
fixed-point theorem as the least fixpoint of \(F\). However, \(F\) being not
necessarily continuous, we have to content ourself with a weaker version.
is the identity. This last definition of \(S(t)\) is called Kleene Sequence, as
it appears in Kleene's Theorem (we present here a weaker version)
\begin{theo}\label{th:up}
On a bounded lattice \((L, \leq, \bot, \top)\) for any increasing function
\(F: L\to L\), if \(F\) has a post-fixpoint \(Y\) \ie{} such that
\(F: L\to L\), for any post-fixpoint \(Y\) of \(F\) \ie{} such that
\(F(Y)\leq Y\), then
\[\bigvee_{t\in\Nat} F^{(t)}(\bot) \leq Y\]
\end{theo}
So, by finding a fixpoint on \(F\), we are able to bound \(S\).
So, by finding a post-fixpoint on \(F\), we are able to bound \(S\).
To find such post-fixpoint, a first solution is to compute every \(S(t)\) up to
a certain point and to be patient. If we find \(t_0\) such that
......@@ -308,11 +300,12 @@ On the example of the integer sawtooth
This kind of approach is not very efficient on intervals. On the previous
example, replacing 10 by 1000000 leads to much slower execution time. Worse, on
an unbounded interval, for instance the ramp \figref{fig:ramp}, this algorithm does
not finish. This is due to the fact that the interval lattice unbounded
an unbounded interval, for instance the ramp \figref{fig:ramp}, this algorithm
does not finish. This is due to the fact that the interval lattice unbounded
increasing chains
\([0, 0] \subseteq [0,1] \subseteq [0, 2]\dots \)\footnote{This is however
untrue on the floating-point interval lattice, more details in annex \ref{sec:annex}}
untrue on the floating-point interval lattice, more details in appendix
\ref{sec:annex}}
This is not a surprise. In a classical programming language, it is trivially
Turing-complete to derive an interval from a variable: adding to a program a
......@@ -380,17 +373,17 @@ have \(F: I \mapsto I + [1, 1] \) and the following trace
\]
This is interesting as we don't lose the sign of the signal, compared to the
naïve widening, which would have returned \(\IntR\) as a whole.
naive widening, which would have returned \(\IntR\) as a whole.
An other improvement was to refine the widening interval into something less
categorical than infinity. This came from an intuition from Faust team
concerning the integer sawtooth (see again fig \figref{fig:rec}) with
\(F:I\mapsto (I + [1, 1] \mod [10, 10])\). They noticed that in this case
\(F(\IntR) = [0, 10]\) is smaller than \(\IntR\), which means that the output of
the sawtooth is always in \([0, 10]\), regardless what was its value on the
previous step. My role was to formalize this and to implement it in order to
improve the compiler. The answer comes once again from lattice theory, it is simply the
upperbound part of Kleene Theorem.
categorical than infinity, using the dual of widening: narrowing. This came from
an intuition from Faust team concerning the integer sawtooth (see again fig
\figref{fig:rec}) with \(F:I\mapsto (I + [1, 1] \mod [10, 10])\). They noticed
that in this case \(F(\IntR) = [0, 10]\) is smaller than \(\IntR\), which means
that the output of the sawtooth is always in \([0, 10]\), regardless what was
its value on the previous step. My role was to formalize this and to implement
it in order to improve the compiler. Taking the dual of the above lattice
theorem \ref{th:up}
\begin{theo}\label{th:down}
On a bounded lattice \((L, \leq, \bot, \top)\) for any increasing function
......@@ -399,12 +392,12 @@ upperbound part of Kleene Theorem.
\end{theo}
Particularily, \(\bigvee_{t\in\Nat} F^{(t)}(\bot) \leq F^{(0)}(\top)\), that is
Particularly, \(\bigvee_{t\in\Nat} F^{(t)}(\bot) \leq F^{(0)}(\top)\), that is
the classical widening, and \(\bigvee_{t\in\Nat} F^{(t)}(\bot) \leq F(\top)\),
which is the case discussed by the team. Being able to automatically devise
upperbounds on feedback loops using modulo or max is very convenient, as those
structures are used in Faust to loop over the values of an array, or initialize
the \(N\) first value of a signal using an array.
upperbounds on feedback loops using modulo or max is very convenient, as this
kind of feedback loops are used in Faust to loop over the values of an array, or
initialize the \(N\) first value of a signal using an array.
\subsubsection{Limitations}
\label{sec:limitations}
......@@ -414,16 +407,10 @@ filter described by a simple order-1 recurrence relation in figure
\figref{fig:smooth}. With a bounded input \(|x(t)| \leq 1\), we have the following
inequation
\begin{figure}
\centering
\includegraphics[width=0.5\textwidth]{../code/smooth-block.png}
\caption{A simple yet useful recursive Faust program}
\label{fig:smooth}
\smallcode{../code/smooth.dsp} \(y(t) = x(t) + \frac{1}{2} y(t-1)\)
\end{figure}
\[\recSig{t}{y(t) = 0}{|y(t)| \leq 1 + \frac{1}{2}|y(t-1)|}\]
\faustexfig{LTI filter smoothing the input}{\label{fig:smooth}}
{0.4}{../code/smooth-block.png}
{0.5}{../code/smooth-sig-old.png}
{\(y(t) = x(t) + \frac{1}{2} y(t-1)\)}{../code/plus.dsp}
which we can unroll for \(t\geq 0\) into
......@@ -468,16 +455,16 @@ know that this sum is converging.
There exists algorithms determining the WCPG of arbitrary
filters\cite{volkova2017reliableimplementation}, but they require a lot of
dependancies (mainly Sollya\cite{ChevillardJoldesLauter2010}, which relies on
dependencies (mainly Sollya\cite{ChevillardJoldesLauter2010}, which relies on
MPFR, MPFI, FPLLL libraries).
\subsection{The Temptating Patch: Unsafe Handmade Annotations}
\subsection{The Tempting Patch: Unsafe Handmade Annotations}
\label{sec:user-annotations}
As we have seen in the two previous parts, bounding signals using intervals
requires often to use small tricks such as expression rewriting in an
interval-friendly faishon (\(x - x^2, ~\sin(\omega t)e^{-t/\tau}\)) or computing
infinite sums (\(y(t) = x(t) + 0.5 y(t-1)\)). Since now, the approach used in
often requires to use small tricks such as expression rewriting in an
interval-friendly fashion (\(x - x^2, ~\sin(\omega t)e^{-t/\tau}\)) or computing
infinite sums (\(y(t) = x(t) + 0.5 y(t-1)\)). Until now, the approach used in
Faust was to use the minimum and maximum operators to ensure bounds. On the
example of the smoother, we could have defined an other signal
\(y(t) = \min(\max(x(t) + y(t-1), -2), 2)\). Having added this, we now get
......@@ -487,8 +474,8 @@ mathematically equivalent. An other advantage is that if we went wrong in our
bounds, the minimum and the maximum would have created saturation, which is the
common thing to do in signal processing to deal with overflow.
But using \(\min\) and \(\max\) has one major drawback: it generates code to
compute useless \(\min\) and \(\max\), code which can be time- and
But using \(\min\) and \(\max\) has one major drawback: it generates code
to compute useless \(\min\) and \(\max\), code which can be time- and
space-consuming in hardware. This is the reason why the Faust team decided to
add primitives in order to annotate signals with intervals which are not
computed but given by the user at its own risks. If the interval is too small,
......@@ -496,20 +483,20 @@ the compiled code offers no guarantee concerning its behaviour in case of
overflow. For most users however, writing native Faust programs is very
rare. The Faust compiler is shipped with a very useful and versatile standard
library\footnote{around 24000 lines of code} that implements nearly all the
basical needs of music makers, acousticians, etc\dots It would then be
sufficient to annotate the standard library where it is needed with carfully
basic needs of music makers, acousticians, etc\dots It would then be
sufficient to annotate the standard library where it is needed with carefully
mathematically verified bounds in order to have useful interval computations in
the compiler.
My third contribution to Faust was to propose and implement a version of
interval annotations. The requirement specifications were quite narrow. On the
WCPG example, we want to bound the output signal of the filter by values
depending on the input signal. There is also the problem of the functions, where
output bounds may depend on the parameters of the function. For instance, in the
smoother of the standard lib, there is a function implementing a filter for the
recurrence equation \(y(t) = x(t) + \alpha y(t-1)\) with WCPG
interval annotations. On the WCPG example, we want to bound the output signal of
the filter by values depending on the input signal bounds. There is also the
problem of the functions, where output bounds may depend on the parameters of
the function. For instance, in the smoother of the standard lib, there is a
function implementing a filter for the recurrence equation
\(y(t) = x(t) + \alpha y(t-1)\) we want to bound the output with the WCPG
\((1-\alpha)^{-1}\) for any \(\alpha < 1\). Note that for \(\alpha \geq 1\), the
filter is unstable, so we also need to express bounds that are conditionned on
filter is unstable, so we also need to express bounds that are conditioned on
booleans. These two usecases directly eliminate the very simple idea of having a
signal annotated by numbers.
......@@ -537,7 +524,10 @@ This might seem quite verbose, but recurrent patterns such as%
input and output the second signal with the interval of the first one can be
implemented in Faust standard library.
To sum up everything, devising the MSB of the fixed-point implementation of a
\subsection{Summing up}
\label{sec:summing-up}
Devising the MSB of the fixed-point implementation of a
signal is just upperbounding its value. Upperbounding variable is very common in
compilation, it can partially be achieved with a lattice-based approach using
interval arithmetic along with widening techniques. Both of them can be
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment