Difference between revisions of "Ball Arithmetic"

Introduction

Ball[1] and Interval Arithmetic[2] both refer to the same idea of calculating with arbitrary precision inexact numbers — that is, with a range of floating point values instead of a single representative. This allows such numbers to reflect their inexactness explicitly.

"The most common use is to keep track of and handle rounding errors directly during the calculation and of uncertainties in the knowledge of the exact values of physical and technical parameters. The latter often arise from measurement errors and tolerances for components or due to limits on computational accuracy. Interval arithmetic also helps find reliable and guaranteed solutions to equations and optimization problems."[2]

The difference between the two is in how each approach represents its values: Interval Arithmetic uses the two ends of the interval [a, b] whereas Ball Arithmetic uses a Ball with a Midpoint and Radius. Each has its advantages. Interval Arithmetic is implemented by several libraries such as MPFI[3] and Ball Arithmetic is implemented by ARB[1].

This feature uses the ARB[1] library to implement Ball Arithmetic.

Environment

In the examples below, unless otherwise specified, the precision of multi-precision floating-point numbers is set to 53 bits via ⎕FPC←53 which is the same precision as standard double-precision numbers. The Printing Precision is set to 17 digits via ⎕PP←17 so as to display the maximum number of significant digits at 53 bits of precision. The Comparison Tolerance is set via ⎕CT←1E¯10.

Calculating With Ball Arithmetic

If you have ever entered an expression such as

```      ÷3
0.3333333333333333
```

undoubtedly, you are aware that the result is inexact. Moreover, as you incorporate into your code the many other functions that produce inexact numbers, you should be wondering "How does the accumulation of inexactness affect my results?".

The answer to that question lies in the province of Numerical Analysis where it is described as "the study of algorithms that use numerical approximation". In general, applying the techniques of Numerical Analysis often involves time-consuming analysis. Fortunately, Ball Arithmetic is an easy to use alternative to that. That is, instead of analyzing the algorithm, run the algorithm with all calculations done with Ball Arithmetic.

Here is the same expression in Ball Arithmetic (use Alt-‘q’ or Ctrl-‘q’, depending on your keyboard layout, to enter the ± symbol):

```      ÷3±
0.3333333333333333±5.6E¯17
```

Essentially, each Ball Arithmetic number carries with it an explicit Radius that describes the interval which contains the true value. Because Ball Arithmetic data propagates throughout a calculation and is never demoted to a lower datatype, its advantage is that after each operation the resulting Ball reliably includes the true answer. Across an entire algorithm, the final result is a Ball that indicates explicitly how precise it is, thus eliminating the guesswork as to the effect on the algorithm of the inexactness of the input and any other inexact constants used.

For example,

```      sqrt←{⍺←= ⋄ {0.5×⍵+⍺÷⍵}⍣⍺⍨⍵} ⍝ Calculate square root using Newton's method
√2±
1.414213562373095±2.2E¯16
sqrt 2±
1.414213562373095±1.1E¯15
10 sqrt 2±
1.414213562373095±2.8E¯15
100 sqrt 2±
1.414213562373095±3.3E¯14
1000 sqrt 2±
1.414213562373095±3.3E¯13
10000 sqrt 2±
1.414213562373095±3.5E¯12
```

As you can see, after a few iterations, the accuracy of the result doesn't improve. However, because it is an iterative algorithm, the precision worsens as the number of iterations increases. Only by using Ball Arithmetic can you see just how the precision of the result is affected.

To demonstrate with an extreme example, consider the following function:

```      f←{⍪⍵-3√(⍵*3)-⍵*2}
f 10*3×⍳5
0.33344450621439137
0.3333334452472627
0.33333444595336914
0.3348388671875
2.25
```

Without knowing much about this (ill-conditioned) algorithm, you have no reason to suspect these results. As it turns out, there's a mathematical proof that as the argument increases to infinity, the algorithm converges to ÷3. Only when you calculate with Ball Arithmetic do you see that there is a problem:

```      f 10*3×⍳5±
0.3334445062153009±2.2E¯12
0.3333334461785853±4.3E¯9
0.333337664604187±0.0000062
0.3385009765625±0.0085
1.25±5
```

As you can see, the precision worsens quickly to the point where the result is meaningless. The last line even says that the answer lies somewhere between ¯3.75 and 6.25. In other words, that answer is so imprecise as to be inaccurate! At higher powers of ten, the answers only worsen.

The good news with Ball Arithmetic is that if you need greater precision, just increase the value of ⎕FPC. For example,

```      ⎕FPC←128 ⋄ ⎕PP←40
sqrt 2±
1.41421356237309504880168872420969807857±5.9E¯39
```

In the case of an ill-conditioned algorithm, we just push the problem off a bit, but not very far — the problem is that the algorithm requires infinite precision:

```      f 10*3×⍳10±
0.333444506214021971364848444464298690637±3.4E¯35
0.333333444444506172880658466392407138179±1.1E¯31
0.333333333444444444506172839632668381652±1.6E¯28
0.333333333333444444444444547458328123363±2.2E¯25
0.333333333333333444444619474517508476812±3.2E¯22
0.333333333333333333552432522356445687706±1.3E¯19
0.333333333333333491771410805881714622956±4.5E¯16
0.333333333333417414223731611855328083038±4.5E¯13
0.33333333363771089352667331695556640625±4.5E¯10
0.3333334587514400482177734375±2.6E¯7
```

While your floating point code might not be so ill-conditioned, the only way to be sure of its results is one of

• Hire a Numerical Analyst
• Become a Numerical Analyst
• Use Ball Arithmetic

Ball Arithmetic Numbers

A Ball Arithmetic number consist of two parts: a Midpoint and Radius. The Midpoint value represents the middle of the Ball and the Radius represents the distance around the Ball of possible values. There is no way to tell which value in the Ball is the correct answer — you should assume that all values are potential candidates.

The precision at which the Midpoint and Radius are stored depends upon the current value of ⎕FPC. For example, the following display

```      ⎕FPC←53
÷3±
0.3333333333333333±5.6E¯17
```

indicates that the number ÷3± has a Midpoint of 0.3333333333333333 at 53-bit precision and a Radius of 5.6E¯17. The Midpoint is stored as a variable precision number, and the Radius is stored as a low-precision number.

Alternatively,

```      ⎕FPC←128 ⋄ ⎕PP←40
÷3±
0.333333333333333333333333333333333333332±1.5E¯39
```

shows the corresponding value at a higher precision and correspondingly narrower Radius.

A Ball describes a range of variable-precision floating point numbers, and so the expression 23±1 represents the set of Real numbers between 22 and 24, inclusive, not the three integers 22, 23, 24.

In general, suffix notation on Ball Arithmetic constants should be sufficient for almost all purposes. That is, the Radius is already correctly calculated for constants such as 23± whose Radius is zero, and 1.2± whose Radius at 53 bits of precision is on the order of 2*¯53 (because of the inexactness of the Midpoint 1.2) which displays as 1.2±2.2E¯16. Although the input notation supports it, the need to override the default Radius should be quite rare.

New Functions and Operator

Functions To Convert To And From Ball Arithmetic Numbers

A Ball Arithmetic constant number may be created through notation (as seen above). Otherwise, a Ball Arithmetic number also may be created from its constituent parts under program control with a new primitive function Contract, and separated into its constituent parts with a new primitive function Distract.

Contract

Contract is a monadic function (≤R) whose right argument is a Real array. It may be used to create a Ball Arithmetic number from the one- or two-elements of its Midpoint and (optional) Radius. The Axis Operator may specify any length 1 or 2 coordinate along which the conversion takes place. It converts the first element along the specified dimension to the Midpoint of the result and (if present) the second element to its Radius. If the second element is not present (i.e., length 1 coordinate), the Radius used is calculated from the precision of the Midpoint.

```      ≤23
23
≤1.2
1.2±1.1E¯16
≤1.3 1E¯20
1.3±1E¯20
≤1r2 1E¯20
0.5±1E¯20
≤○1
3.141592653589793±1.1E¯16
≤(○1) 1E¯20
3.141592653589793±1E¯20
```

Distract

Distract is a monadic function (≥R) whose right argument is a Real array. It may be used to separate out a Ball Arithmetic's Midpoint and Radius into two distinct values. The Axis Operator may specify the place where the new length 2 coordinate is to be inserted. The pair of values along this new coordinate are the Midpoint and Radius of the corresponding number from the right argument. If the number can be represented exactly, the Radius is zero.

```      ⎕PP←10
≥1.2
1.2 1.1102230246251565E¯16
≥1.2x
6r5 0
≥1.2v
1.2 1.110223024625157E¯16
≥1.2±
1.2 2.220446049250313E¯16
≥1.5
1.5 0
≥1.5x
3r2 0
≥1.5v
1.5 1.110223024625157E¯16
≥1.5±
1.5 0
a←(○1±)×⍳4
a,'|',≥a
3.141592654±4.4E¯16 |  3.141592654 4.440892148E¯16
6.283185307±8.9E¯16 |  6.283185307 8.881784296E¯16
9.424777961±1.3E¯15 |  9.424777961 1.332267643E¯15
12.56637061±1.8E¯15  | 12.56637061  1.776356859E¯15
```

These numbers can be represented exactly:

```      ≥⍳4±
1 0
2 0
3 0
4 0
≥[2]⍳4±
1 0
2 0
3 0
4 0
≥[1]⍳4±
1 2 3 4
0 0 0 0
```

Identities

The Contract of the Distract of a number both along the same axis, returns the same number.
The Distract of the Contract of a number both along the same length two axis, returns the same number.

Operator To Provide Arbitrary Precision Specific Derived Functions

The DoubleTilde monadic operator (, Alt-'Q' or Alt-shift-'q') takes a primitive function as its left operand and produces an arbitrary-precision-specific derived function.

This operator should be considered preliminary.

Conversion Between Balls and Intervals

At the moment, the monadic derived function is defined for two primitive functions.

Because Ball Arithmetic has a close relationship with Interval Arithmetic, it may be desirable to convert between the two forms: Endpoint to Midpoint-Radius and its inverse. The resulting derived function is not scalar.

For example,

```      x←0±∞ 1±2
≥≈x
¯∞ ∞
¯1 3
≥[1]≈x
¯∞ ¯1
∞  3
≥[2]≈x
¯∞ ∞
¯1 3
≤≈≥≈x
0±∞ 1±2
≤≈≥[2]≈x
0±∞ 1±2
≤[2]≈≥≈x
0±∞ 1±2
≤[2]≈≥[2]≈x
0±∞ 1±2
≤[1]≈≥[1]≈x
0±∞ 1±2
```

Union and Intersection of Balls

At the moment, the dyadic derived function is defined for two primitive functions.

Because Ball arithmetic numbers represent a range of floating point numbers, at times it is desirable to produce their Union or Intersection. The resulting derived function is scalar.

For example,

```      x←0.50.5      ⍝ The interval [0,1]
x∪≈x+2         ⍝ The interval [0,1] ∪ [2,3] ←→ [0,3]
1.51.5
x∪≈x+1         ⍝ The interval [0,1] ∪ [1,2] ←→ [0,2]
11
x∩≈x+1         ⍝ The interval [0,1] ∩ [1,2] ←→ [1,1]
1
```

If two Balls have no points in common, the result is a NaN ():

```      x∩≈x+1.1       ⍝ The interval [0,1] ∩ [1.1,2.1] ←→ ∅
∅
```

Comparison Functions

This might be the place to define new derived functions such as <≈ which can return a non-Boolean result so as to cover the more complicated relationship of one Ball with another with respect to overlap, etc.

Hypercomplex Ball Arithmetic

Ball Arithmetic also extends to all Hypercomplex numbers (Complex, Quaternion, and Octonion) such as 1.2±i2.2 which represents a Complex Ball Arithmetic number whose Real part is 1.2± and whose Imaginary part is 2.2±. Note that the display of that number (at 53 bits of precision)

```      1.2±i2.2
1.2±2.2E¯16i2.2±4.4E¯16
```

demonstrates the inexactness of both floating point coefficients.

Invalid Results

Not all Balls are valid arguments to all primitive functions. Sometimes, the result cannot be expressed in a Ball, such as

```      √0±1
∅
```

To see why, split apart the argument into two parts: ¯0.5±0.5 and 0.5±0.5 and apply the Square Root function to the separate parts. This yields a result that consists of two quite disjoint Balls: the square root of the Ball ¯0.5±0.5 is 0i0.5±0.5, whereas the square root of the Ball 0.5±0.5 is 0.5±0.5 . As the two parts of the result are orthogonal to each other, they cannot be expressed as a single Midpoint-Radius Ball, and so the result is a NaN (= Not-A-Number = ).

Ball Arithmetic works best when given small Radii, that is, small relative to the Midpoint. Calculations on large Radii numbers such as 0i0.5±0.5 don't always produce useful results. For example,

```      ⎕←x←0i0.5±0.5*2
¯0.25±0.75
```

which follows from the standard way of multiplying Complex numbers[4]. Unfortunately in this case, because the Radius is so large that it includes negative numbers on the low end, the Square Root of this particular Squared interval is not representable, that is, it is a NaN:

```      √x
0i∅
```

These calculations work better using smaller Radii, such as

```      ⎕←y←0i0.5±0.5E¯10*2
¯0.25±5E¯11
√y
0i0.5±5E¯11
```

where Square Root now returns a meaningful result.

Mixed Notation

The notation for Ball Arithmetic may be combined with other Point Notations in many ways, such as

```      1<_p/>1±
3.141592653589793±8.9E¯16
1x1±
2.718281828459045±4.5E¯16
```

For more details about combining different forms of Point Notation, see Point Notation and in particular the section on Mixed Notation.

Implementation

Datataypes

This feature adds four new datatypes to the language: one each for the four Division Algebras. Note that marking any one of the Hypercomplex coefficients as a Ball automatically marks them all as Balls.

 Real 23± or 1.2± Complex 23i1.2± Quaternion 23i1.2j0.5k¯3.14± Octonion 23i1.2j0.5k¯3.14¯3ij19jk¯3kl6±

Primitive Functions

As described above, this feature adds two new monadic primitive functions Contract (≤R) and its inverse Distract (≥R), along with one new monadic primitive operator DoubleTilde f≈ R and L f≈ R.

All primitive and system functions and operators are extended to work with Ball Arithmetic, so you can test out this feature on all algorithms, with the exception of trying to find the Factorial of a multiple-precision/Ball Arithmetic Hypercomplex number, in which case you are temporarily out of luck.

Workspace Version Number

For the moment, if you save a workspace which contains a Ball Arithmetic number in a variable, the workspace version number is 0.03. Otherwise, the version number is 0.02. Workspaces with Ball Arithmetic Numbers may not be loaded by an earlier version of the interpreter.

Pending Design Decisions

• The direct comparison functions (<≤=≥>≠≡≢) and indirect ones (⍳⍸∊⍷⍒⍋) are currently implemented to compare the Midpoints of the arguments ignoring the Radius — this may change.
• The Floor function and its related functions (Ceiling, Residue, GCD, LCM, and Encode) are currently implemented to execute on the Midpoints of the arguments ignoring the Radius — this may change.
• We need new primitives to describe the different ways in which Balls can intersect/overlap, possibly through a new monadic operator (perhaps ≈) on the existing direct and indirect comparison and Floor-related functions.
• Certain primitives are defined on limited ranges, and sometimes multiple ranges where it might use different algorithms for different ranges (e.g. L!R). There are open questions about how to handle a Ball that spans two or more of these ranges as to how to calculate the result.
• Infinity and NaN have not been given any special treatment as yet. The results are whatever ARB has chosen which we may need to override.

This feature is available as an Alpha release — it may be downloaded from the Alpha release directory.

Acknowledgements

• The Arb team has contributed their Ball Arithmetic[1] code which made this feature possible.
• David Rabenhorst has contributed many very useful ideas and many hours debugging the implementation.

References

1. NARS2000 Wiki, Ball Arithmetic
2. Wikipedia, Interval Arithmetic
3. MPFI
4. [1]