# Computations with inexact rings

## Classical interval or ball arithmetic

We provide a rigorous way to work with elements of $\mathbf{R}$ or $\mathbf{C}$ using an interval representation. More precisely, instead of working with real numbers, we work with subsets $x \subseteq \mathbf{R}$. It does not really matter how they look like, but at the moment these are implemented as balls $x = \mathrm{B}_r(x_0) \subseteq \mathbf{R}$ with a (rational) number $x_0$ and a (small) rational number $r > 0$. These rational numbers have some specific form, but in the end they are still just rational numbers. When I speak of "for all $x \subseteq \mathbf{R}$", I always mean the $x$ which can be written down in our model.

### The first principle (inclusion principle)

Here is what is promised (Fredrik calls this the *inclusion principle*)

> If $f \colon \mathbf{R} \to S$ is a function and $A_f$ an algorithm realizing $f$, then $f(x) \subseteq A_f(x)$ for all $x \subseteq \mathbf{R}$.

There is actually no guarantee on how large $A_f(x)$ is (compare to $f(x))$.

This has some funny resp. surprising consequences:
- If $f(x) = \lceil x \rceil$ is the usual ceiling, then $f(\mathrm{B}_{r}(0)) = \{0, 1\}$ for all $r > 0$. In particular, we will always have that $A_f(x)$ is a ball containing $0$ and $1$. In practice

julia> c = ceil(mathbf{R}("0 +/- 0.0001"))
[+/- 1.01]

julia> midpoint(c)
0.50000000000000000000

julia> radius(c)
[0.50000000093132257462 +/- 4.53e-21]

- Operations can fail. For example, the partial function $f \colon \mathbf{R} \to \mathbf{R}, \ x_0 \mapsto x_0^{-1}$ will fail for $x \subseteq \mathbf{R}$ containing $0$.

- The promise applies to the case where $f$ is a predicate, which yields the following:

> If $f \colon \mathbf{R} \to \{0, 1\}$ is a predicate and $A_f$ an algorithm realizing $f$, then $A_f(x)$ holds, if and only if $f(x_0)$ holds for all $x_0 \in x$.

This forces us to work with three-valued logic, although we never mention it anywhere.

The following code illustrates this:

if x < y
  # do this
else
  # we cannot assume that y >= x
end

### The second principle (convergence)

Here is the second promise (although technically, Fredrik does not promise this, but we always assume (and "check") that it holds).

> If $f \colon \mathbf{R} \to S$ is a (continuous) function realized by an algorithm $A_f$, and $(\mathrm{B}_{r_i}(x_i))_{i \geq 1}$ is a sequence of balls with $x_i \to x_\infty$ and $r_i \to 0$, then $A_f(\mathrm{B}_{r_i}(x_i))$ is a sequence of balls $(B_{s_i}(y_i))_{i \geq 1}$ with $y_i \to f(x_\infty)$.

This is why ball arithmetic is useful. In our applications, we almost always start with exact input, do some computation on the real side. If we increase the precision when going from the exact world to the (approximated) real world, we get better approximations for the real information.

### Notes on the implementation

- There is a notion of "exact elements", which are the ones with $x = \{x_0\}$ being a singleton. Those are implemented as the balls with radius $r = 0$. For example, $0$ and $1$ can be exact.
- The set of balls has a true one and zero, which act as multiplicative respectively additive identity elements.
### Remark

The same applies to matrices and polynomials. Because of the inclusion principle and the fact that most ring operations are not continuous, things get weird:
- There is no useful notion of gcd. The gcd of $x \pm 0.1$ and $x$ must be a ball containing $1$ and $x$, for example $(0.5 \pm 0.5)x + (0.5 \pm 0.5)$
- There is no useful notion of rank or solving for underdetermined systems.

Other things are possible, although the answer can be "input does not have enough precision":
- computing interval containing the roots of a polynomial
- isolating the roots of a polynomial
- solving linear systems with unique solutions

## $p$-adic numbers

$p$-numbers will be modelled as cosets $x = x_0 + O(p^n) = x_0 + p^n \mathbf{Z}_p$ which are exactly the balls $\mathrm{B}_{1/p^n}(x_0)$ with respect to the $p$-adic absolute value.

We can and should have a model which supports the inclusion principle and convergence (if possible)
### Current model

Our current model is a mixed of interval and fixed precision arithmetic. The biggest downside is our choice of `==`. In Oscarland we consider $x = x_0 + O(p^n)$ and $y = y_0 + O(p^m)$ "equal", if $x \equiv y \bmod p^{\min(n, m)}$. In other words, if $\mathrm{B}_{1/p^n}(x_0) \cap \mathrm{B}_{1/p^m}(y_0) \neq \emptyset$ . It is immediate, that this is not transitive:

julia> Qp = PadicField(2);

julia> x = Qp(3, precision = 3)
2^0 + 2^1 + O(2^3)

julia> y = Qp(1, precision = 1)
2^0 + O(2^1)

julia> z = Qp(7, precision = 3)
2^0 + 2^1 + 2^2 + O(2^3)

julia> x == y && y == z
true

julia> x == z
false

This is the single most irritating things about our current implementation and makes it *useless* unless one really really knows what one is doing. It does not help, that other algorithms/methods are utterly underspecified.

Let's call this "equality" the "weak, intransitive" equality.

#### Consequences of `==` being the the weak equality

As soon as specifications of methods involve `==`, the inclusion principle is not guaranteed and in fact often violated. Examples include:
- Specifying a solve function by saying that `x = solve(A, b)` satisfies `A*x = b` .
- Specifying `xgcd` by saying that `g, u, v = gcdx(f, h)` satisfies something (what exactly?).

### Magma model (free precision rings and fields)

This Magma model is very much like ball arithmetic, but distinguishes between true mathematical equality and "intransitive, weak" equality.

| | | |
| -------------------------- | -------------------------------------------------------------------------------------- | ------------------------------------------------------------------------ |
| Representation of elements | $x = x_0 + p^n \mathbf{Z}_p = \mathrm{B}_{1/p^n}(x_0)$ | |
| Equality (`eq`) | set theoretic | Consistent. |
| One | ```One(Q2); 1 + O(2^20)``` | This is not a multiplicative identity. |
| IsOne(One()) | `false` | Consistent. |
| Zero | ```> Zero(Q2); 0``` | This is a exact (with radius $0$). (Although the documentation claims that the zero is inexact) |
| IsZero(Zero()) | `true`<br> But note that `IsZero(x - x) eq false` (if `x` is not exact) | Consistent. |
| IsWeaklyZero() | `IsZero(x - x) eq true` | |
| IsWeaklyEqual() | this is the intransitive relation | |
| Valuation() | `Valutation(O(p^n))` gives `n` | I don't think this is consistent. |
| | | |
| (X)GCD | `g, u, v = XGCD(f, h)` does satisfy<br>`IsWeaklyEqual(u * f + v * h, g)`, but not `eq` | It is not specified what exactly is returned. |
| | | |
### Sage model

Same as ours (`==` is the weak, intransitive equality), with two exceptions:
- relative/absolute precision is always capped
- there exists an *exact* zero with valuation `+Infinity`

## Final comments

- I don't think our current model makes any sense (nor does the Sage one)
- I want the Magma model (except for the `valuation(O(p^n))` behaviour) including exact zero and one.
- I think there is a need for the "weak equality", but it should be carefully documented and or made clear that this is now used. I'd rather have an error for `xgcd(f, g)` tell me to use `xgcd_weak(f, g)` or whatever the name is.