Flint C API#

Pure C flint objects and functions#

Rounded floating point interval arithmetic and simple math functions in C

struct flint#

Rounded floating point interval with tracked value

There are three tracked values a lower and upper bound of the interval as well as a tracked valuse that acts exactly like a 64 bit float (c double) for easy cast back to a float.

double a#: The lower bound

double b#: The upper bound

double v#: The tracked value

Pre-defined constants in flint types#

These are the real number mathicatical constants as flint types. They the lower and upper boundaries are two consecutive floating point numbers that bound the actual constant, and the tracked value is the boundary value that is closer to the constant’s value.

Todo

Add in the other constants defined in the GNU math.h header file.

FLINT_2PI#

FLINT_PI#

FLINT_PI_2#

Casting to and from flints#

The following rules are followed when casting to flints:

For integers, if the integer can be exactly represented as a 64-bit float, then the lower, upper, and tracked values are all set to the tracked values
For intergers, if they are larger than can be represented as a 64-bit float, then the tracked value is set the the closest floating point value, and the lower and upper values are raised and lowered by 1 ULP.
For floating point values, it will be assumed that the intended value is NOT exactly represetned by the input. The tracked value shall be set to the input and the upper and lower bounds are moved outware by 1 ULP of the input type.

When casting from flints, the tracked value is returns as a 64 bit floating point (c double) which can then be cast using standard C casting rules.

static inline flint int_to_flint(long long l)#

static inline flint float_to_flint(float f)#

static inline flint double_to_flint(double f)#

static inline double flint_to_double(flint f)#

Floating point special value queries#

There are several special values in the IEEE-754 floating point representation including +/- infinity and not-a-number (NaN). Since those values might be fed as inputs into various functions that have been chained together, numpy requires four functions to check for those values: nonzero, isnan, isinf, and isfinite.

For the flint types:

nonzero : the interval does NOT contain zero. isnan : neither boundary or tracked value are a NaN. isinf : one of the two boundaries is infinite. isfinite : neither boundary is infinite.

static inline int flint_nonzero(flint f)#

static inline int flint_isnan(flint f)#

static inline int flint_isinf(flint f)#

static inline int flint_isfinite(flint f)#

Comparisons#

For comparisons we have to worry about equality as well as a sense or order (greater than or less than). Since flints are intervals, with an lower and upper boundary we can use standard interval order where we define equality as the intervals overlapping, and less that or greater than as an intervals not overlapping at all and lying completely to the left or right along the number line.

Note

Although the flints have a partial order, they do not have a strict total order. Specifically the equality of values is NOT transitive.

\[a = b, \quad\text{and}\quad b = c, \quad\text{but}\quad a \ne c.\]

The following functions implement all comparisons ==, !=, >, >=, <, and <=. The return value is 1 if it evalutes true, 0 if false.

static inline int flint_eq(flint f1, flint f2)#

static inline int flint_ne(flint f1, flint f2)#

static inline int flint_le(flint f1, flint f2)#

static inline int flint_lt(flint f1, flint f2)#

static inline int flint_ge(flint f1, flint f2)#

static inline int flint_gt(flint f1, flint f2)#

Arithmetic#

The IEEE-754 standard for floating point operations is as follows:

Assume that all input operand are exact
The result should be calculated exactly. If the result can not be exactly represented by a floating point value it should be rounded to a neighboring floating point value.

In this work, we assume that the flint contains the exact value somewhere within the interval, whose boundaries are given by floating point numbers. We can then carry out exact interval arithmetic and then round lower and upper the boundarys of the result down and up to gaurantee that the result lies somewhere in the new rounded interval.

Unary operations#

There are two unary operators defined: 1. plus + 2. negation -

static inline flint flint_positive(flint f)#

static inline flint flint_negative(flint f)#

Binary operations#

There are 4 binary operations defined and their in-place versions:

Addition `+` and `+=`#

static inline flint flint_add(flint f1, flint f2)#

static inline void flint_inplace_add(flint *f1, flint f2)#

static inline flint flint_scalar_add(double s, flint f)#

static inline flint flint_add_scalar(flint f, double s)#

static inline void flint_inplace_add_scalar(flint *f, double s)#

Subtraction `-` and `-=`#

static inline flint flint_subtract(flint f1, flint f2)#

static inline void flint_inplace_subtract(flint *f1, flint f2)#

static inline flint flint_scalar_subtract(double s, flint f)#

static inline flint flint_subtract_scalar(flint f, double s)#

static inline void flint_inplace_subtract_scalar(flint *f, double s)#

Multiplication `-` and `-=`#

static inline flint flint_multiply(flint f1, flint f2)#

static inline void flint_inplace_multiply(flint *f1, flint f2)#

static inline flint flint_scalar_multiply(double s, flint f)#

static inline flint flint_multiply_scalar(flint f, double s)#

static inline void flint_inplace_multiply_scalar(flint *f, double s)#

Division `/` and `/=`#

static inline flint flint_divide(flint f1, flint f2)#

static inline void flint_inplace_divide(flint *f1, flint f2)#

static inline flint flint_scalar_divide(double s, flint f)#

static inline flint flint_divide_scalar(flint f, double s)#

static inline void flint_inplace_divide_scalar(flint *f, double s)#

Math functions#

In addition to arithmetic, we need to make sure that we can apply most of the elementary functions to the flints and preserve our guarantee that the result interval included the exact value of the function applied to any value in the input interval. There are few parts to consider for this: how accurate is the result for the elementary function on floating point values? And how do we convert a function for a floating point number into a function for a flint object?

To address the first issue: how accurate are the elementary function on floating point numbers, we have to turn to the documentation for the math libraries. The only one I was able to find was the Errors in Math Functions section of the gnu libc manual. The answer turned ot to be more complicated that I anticipated, depending on the hardware architecture as well as the individual functions. But for a quick look through seemed to indicate that for most functions for the arm 64bit and x86 32 and 64 bit, the double precision versions of most math functions had either 1 or 2 ULP of accuracy.

Now lets discuss a naive implementation of the math functions for a flint object. If the function is monotonic,

\[a > b \to f(a) > f(b),\]

then the endpoints of the boundaries of the input interval will be become the boundaries of the output interval. So for a function double func(double x) we can write the flint version as

flint flint_func(flint x) {
    flint res;
    res.a = nextafter(nextafter(func(x.a), -INFINITY), -INFINITY);
    res.b = nextafter(nextafter(func(x.b), INFINITY) INFINITY);
    res.v = func(x.v);
    return res;
}

This works for all monotonically increasing functions whose domain is the full real number line, and will work with trivial changes for monotonically decreasing functions. However, what happens when the input interval contains an extremum of the function? What happens if the domain is NOT the full number line. We will need some special consideration for each case.

The first case is what happens when the function is not monotoic. In most cases, there will be a single extremum, such as the minium in the hypot function as either value goes through zero. If the interval contains that extremum, then the upper or lower boundary of the output flint need to be reset to use the extremum as the boundary AND that particular boundary value does not need to be expanded. Finally some care need to be taken in the case of sine and cosine, that oscillat with multiple extrememum.

The second case is what happens when the function’s domain does not cover the entire number line. The standard response is for a C function to return NaN for an input ouside of it’s accepted domain. However for flints, we can be a little more robust. If the interval spans the edge of the domain, then the function would return NaN for one boundary value and a non NaN number for the other. In that case we will assume that the value’s we’re interested in are only the ‘real’ (non NaN) values and replace the NaN’s with the corresponding values at the edge of the domain. To give an explicity example, consider the C sqrt function. lets say we have a flint whose lower value is slightly below zero (1.5e-16) and upper value is above zero (1.0e-8). The endpoints would evaluate to NaN and 1.0e-4, but we can replace the lower endpoint with 0.0 instead of NaN. The tracked value is either kept as the direct result or replaced with the value at the edge of the domain if it would be NaN.

Note

It is because of the special behavior depending on the interval spaning extremums or domain boundaries that make the evaluation of these functions for flints significantly slower compared to using floating point value directly

The following functions defined in the c99 math.h header file have flint implementation with the following signature:

static inline flint flint_FUNCNAME(flint fa, ...)#

Functions#

power absolute sqrt cbrt hypot exp exp2 expm1 log log10 log2 log1p erf erfc sin cos tan asin acos atan atan2 sinh cosh tanh asinh acosh atanh

FLINT_MONOTONIC(fname)#

Python C Extension objects and functions#

Python C-API for flints

struct PyFlint#

The python flint object

flint obval#: The flint c struct

static PyTypeObject *PyFlint_Type_Ptr#

A pointer to the PyTypeObject for the PyFlints

Because of the programming standard or declaring everything static outside of enclosing functions, it is neccessary to use this pointer when you need to access the PyTypeObject* for the PyFlint instead of the direct address &PyFlint_Type.

static int NPY_FLINT#: The integer enumeration for the custom ‘flint’ type in numpy.

static int import_flint(void)#

Import the c api for numpy-flint python module

Returns:: 0 on success -1 on failure

static inline int PyFlint_Check(PyObject *ob)#

Check if an object is a flint

Parameters:

ob – The PyObject to check

Returns:

1 if the object is a flint, 0 otherwise

static inline PyObject *PyFlint_FromFlint(flint f)#

Create a PyFlint object from a c flint struct.

Parameters:

f – The c flint struct

Returns:

A new PyFlint object that contains a copy of f

static inline flint PyFlint_AsFlint(PyObject *ob)#

Get a C flint from a PyFlint

Parameters:

ob – The PyFlint object

Returns:

The coresponding flint value

static inline double PyFlint_AsDouble(PyObject *ob)#

Get a C double from a PyFlint

Parameters:

ob – The PyFlint object

Returns:

The coresponding flint center value

Flint C API#

Pure C flint objects and functions#

Pre-defined constants in flint types#

Casting to and from flints#

Floating point special value queries#

Comparisons#

Arithmetic#

Unary operations#

Binary operations#

Addition + and +=#

Subtraction - and -=#

Multiplication - and -=#

Division / and /=#

Math functions#

Functions#

Python C Extension objects and functions#

Addition `+` and `+=`#

Subtraction `-` and `-=`#

Multiplication `-` and `-=`#

Division `/` and `/=`#