Ball Arithmetic

Introduction

Ball and Interval Arithmetic 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."

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 and Ball Arithmetic is implemented by ARB.

This feature uses the ARB 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 is stored depends upon the current value of ⎕FPC — the precision of the Radius is fixed at 30 bits. 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 an arbitrary precision number; the mantissa of the Radius is stored as a fixed-precision (30 bit) number along with an exponent that may be arbitrarily large or small.

Alternatively,

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

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

The approximate number of digits of precision for a given value of ⎕FPC is about 10⍟2*⎕FPC. The value of ⎕FPC for a given required precision is about ⌈2⍟10*digits.

The value of ⎕FPC should be calculated so as to be sufficient to represent the Midpoint of all ball numbers being used. Otherwise, some results might be misleading.

For example,

2v*⎕FPC←128 3.402823669E38 ⍟ ⎕← * 1±1e¯1000 2.718281828±1.2E¯38 1±7.3E¯39

What happened — the logarithm of an exponential should return the original argument? In this case, it didn't because ⎕FPC was too small!

2v*⎕FPC←4096 1.044388881E1233 ⍟ ⎕← * 1±1e¯1000 2.718281828±2.7E¯1000 1±1E¯1000

A Ball describes a range of arbitrary-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 Real number Midpoint and (optional) Real number 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 

The domain of this function is limited to Real numbers. Hypercomplex numbers signal a DOMAIN ERROR.

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

The domain of this function is limited to Real numbers. Hypercomplex numbers signal a DOMAIN ERROR. To see the Midpoint-Radius of a Hypercomplex number, first Dilate it so as to separate out the coefficients and then apply Distract to extract the Midpoint-Radius as in ≥>R.

For example,

      x←2±0i0±1
      ≥ x
DOMAIN ERROR
      ≥ x
      ∧
      ≥> x
2 0
0 1
      <≤≥> x
2±0i0±1

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

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 

The domain of these functions is limited to Real numbers. Hypercomplex numbers signal a DOMAIN ERROR. To see the Endpoints of a Hypercomplex number, first Dilate it so as to separate out the coefficients and then apply ≥≈ to extract the Endpoints as in ≥≈>R.

For example,

      ≥≈ 2±0i0±1
DOMAIN ERROR
      ≥≈ 2±0i0±1
      ∧
      ≥≈>2±0i0±1
 2 2
¯1 1
      <≤≈≥≈>2±0i0±1
2±0i0±1

Union and Intersection of Balls

The left operand of the DoubleTilde operator is defined for the Union and Intersection functions (∪∩) as follows:

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

For example,

      x←0.5±0.5      ⍝ The interval [0,1]
      x∪≈x+2         ⍝ The interval [0,1] ∪ [2,3] ←→ [0,3]
1.5±1.5
      x∪≈x+1         ⍝ The interval [0,1] ∪ [1,2] ←→ [0,2]
1±1 
      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] ←→ ∅
∅

The identity element for ∩≈ is 0±∞.

The derived functions are scalar.

Comparison Functions

The left operand of the DoubleTilde operator is defined for the six comparison functions (<≤=≥>≠) as follows:

For f as one of the above six functions, Lf≈R is TRUE iff for all t∊L and u∊R, then   t f u   is TRUE. These comparisons are all done with ⎕CT←0 because, in a sense, Ball comparisons carry their own Fuzz.

For example,

      (L R)←4±1 6±1
      L<≈R
0
      L≤≈R
1

The derived functions are scalar.

Containing Functions

The left operand of the DoubleTilde operator is defined for two containing functions (⊂⊃) as follows:

The expression L⊂≈R is TRUE iff for all t∊L it is TRUE that t∊R and that the containment is strict — that is, none of the endpoints match.

The expression L⊃≈R is TRUE iff R⊂≈L is TRUE.

For example,

      (L R)←6±4 6±2
      L⊃≈R
1
      R⊂≈L
1

The derived functions are scalar.

Positional Relationships

For L and R as finite intervals with non-zero Radii (i.e., multi-valued), here are the possible Relative Positions between the two Intervals, where L0, L1 and R0, R1 are the endpoints of L and R, respectively — (L0 L1)←≥≈L and (R0 R1)←≥≈R:

Relative Positions	Endpoint Relationships	Defining Expressions		Sample Data L   R	Distinct Ball {⍺∘≈⍵}
├──L──┤  ├──R──┤	L0 L1≡R0 R1	L=≈R		6±4 6±4	0
├──L──┤  ├──R─┤	(L0=R0)∧L1>R1			6±4 5±3	1
├──L─┤  ├──R──┤	(L0=R0)∧L1<R1			6±4 7±5	¯1
├──L──┤   ├──R─┤	(L0<R0)∧L1=R1			6±4 7±3	2
├──L─┤  ├──R──┤	(L0>R0)∧L1=R1			6±4 5±5	¯2
├──L──┤   ├─R─┤	(L0<R0)∧L1>R1	L⊃≈R		6±4 6±2	3
├─L─┤  ├──R──┤	(L0>R0)∧L1<R1	L⊂≈R		6±4 6±6	¯3
├──L──┤   ├──R──┤	(L0<R0)∧(L1<R1)∧L1>R0			6±4 8±4	4
├──L──┤  ├──R──┤	(L0>R0)∧(L1>R1)∧L0<R1			6±4 4±4	¯4
├──L──┼──R──┤	L1=R0	(L≤≈R)>L<≈R	L=∀R	6±4 14±4	5
├──R──┼──L──┤	L0=R1	(L≥≈R)>L>≈R	L=∃R	6±4 0±2	¯5
├──L──┤ ├──R──┤	L1<R0	L<≈R	L<∀R	6±4 17±5	6
├──R──┤ ├──L──┤	L0>R1	L>≈R	L>∃R	6±4 0±1	¯6

The above table can be interpreted as follows:

• The first column lists all possible Relative Positions two finite non-zero Radii (i.e., multi-valued) Balls may assume.
• The second column lists the Endpoint Relationship that defines the corresponding Relative Position.
• The third column lists a Defining Expression using the DoubleTilde operator that is TRUE for the corresponding Relative Position and FALSE for all of the other Relative Positions in different rows.
• The fourth column lists a Defining Expression using the Always and Possibly operators proposed by David Rabenhorst.
• The fifth column contains Sample Data that defines the corresponding Relative Position.
• The sixth column gives the value of the Distinct Ball function {⍺∘≈⍵} when executed on the Sample Data.

Note the table is organized such that apart from the first row (whose Endpoint Relationship is Symmetric in L and R) the rows occur in six pairs where the Endpoint Relationship for each row in the pair can be obtained from the other by swapping L and R. The pairs of rows are ordered roughly by the amount of overlap between the two Balls.

As you can see, we are missing Defining Expressions for several of the Relative Positions. I was hoping that this table would be filled in with more detail so it would serve as a measure of whether or not we had a complete set of functions/operators/hyperators which could be used to tell apart all of the Relative Positions without having to resort to the Distinct Ball function {⍺∘≈⍵} constructed explicitly for this purpose. Please add to this table if you find a missing expression, or wish to propose whole new functions, operators, and/or hyperators.

If you would like to try a Defining Expression on the Sample Data, collect the data into a 13 by 2 matrix, put your Defining Expression into (say) an Anonymous Function and reduce the same data matrix along the last coordinate — if you have a successful function, you should see a 13-element Boolean vector with a single one in it. For example,

      sd←6±4,⍪6±4 5±3 7±5 7±3 5±5 6±2 6±6 8±4 4±4 14±4 0±2 17±5 0±1
      {(⍺≤≈⍵)>⍺<≈⍵}/sd
0 0 0 0 0 0 0 0 0 1 0 0 0 
      {(⍺≥≈⍵)>⍺>≈⍵}/sd
0 0 0 0 0 0 0 0 0 0 1 0 0 
      {⍺<≈⍵}/sd
0 0 0 0 0 0 0 0 0 0 0 1 0 
      {⍺>≈⍵}/sd
0 0 0 0 0 0 0 0 0 0 0 0 1

Finally, the following example shows that the Sample Data are all distinct:

      {⍺∘≈⍵}/sd
0 1 ¯1 2 ¯2 3 ¯3 4 ¯4 5 ¯5 6 ¯6

Hypercomplex Ball Arithmetic

Real 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 Real and Imaginary coefficients.

Small v. Large Radii

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,

      ⎕PP←8             ⍝ To reduce fractional clutter
      x←0i0.5±0.5
      √x
0.4267767±0.57i0±1.4

Unfortunately in this case, because the Radius is so large it overwhelms the Midpoint, as can be seen when we Square the Square Root:

      (√x)*2
0.18213835±2.8i0±2.8
      x=(√x)*2
0

These calculations work better using smaller Radii, such as

      y←0i0.5±0.5E¯30
      √y
0.5±2.5E¯31i0.5±7.5E¯31
      (√y)*2
0±1E¯30i0.5±1E¯30 
      y=(√y)*2
1

where Square Root now returns a meaningful result, and its Square returns the original value.

Mixed Notation

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

      1p1±
3.141592653589793±4.4E¯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 cannot 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.

Download

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 code which made this feature possible.
• David Rabenhorst has contributed many very useful ideas and many hours debugging the implementation.