The doubles of the C code

The C code of ratpoints holds the ends of its search intervals as doubles: the intervals of the caller (args->domain), the default interval [-h, h] of the library and the ends of the intervals that the Sturm step keeps. Apart from the estimate of the work in the choice of the moduli (numerators_for, Ratpoints.numeratorsFor, where the ends enter the arithmetic of Lean's Float by their bits, Ratpoints.Double.toFloat), it does no floating-point arithmetic on them: it converts longs to doubles, compares doubles, takes them apart into mantissa and exponent (frexp), scales them by powers of 2 (ldexp), steps to the next double (nextafter), converts doubles to longs and unsigned longs, and prints them, as printf of glibc does under the conversions %f and %.pf (Ratpoints.Double.printFixed); the report of the library prints doubles of its cost model too (Ratpoints.Double.ofFloat). The model holds a double as its IEEE 754 binary64 bits (Ratpoints.Double), with these operations as functions on the bits, each doing what the operation of the C language or the C library does to a double on the platforms of ratpoints (IEEE 754 binary64 with rounding to nearest, ties to even, as gcc and glibc implement them): the comparisons, in which a NaN is unordered and -0.0 equals 0.0, the negation and fabs, which act on the sign bit, the conversions, frexp, ldexp and nextafter towards ±HUGE_VAL, but for a signalling NaN, which glibc's frexp, ldexp and nextafter return quieted and the model as it is (ratpoints passes them no NaN). The exact value of a finite double is m 2 ^ e for the integers of Ratpoints.Double.toDyadic, and a dyadic number is rounded to a double, to nearest with ties to even, by Ratpoints.Double.ofDyadic, a fraction by Ratpoints.Double.ofRatio (for strtod, which the program uses).

Main definitions

• Ratpoints.Double: a C double, by its bits.
• Ratpoints.Double.isNaN, Ratpoints.Double.isInf, Ratpoints.Double.isFinite: the NaNs, the infinities and the finite doubles.
• Ratpoints.Double.lt, Ratpoints.Double.le, Ratpoints.Double.beq: the comparisons <, <= and == of the C language on doubles.
• Ratpoints.Double.neg, Ratpoints.Double.fabs: the negation and the absolute value.
• Ratpoints.Double.toDyadic, Ratpoints.Double.ofDyadic: the exact value of a finite double, and the double nearest to a dyadic number.
• Ratpoints.Double.toRat, Ratpoints.Double.ofRatio: the exact value of a finite double as a rational number, and the double nearest to a fraction.
• Ratpoints.Double.assemble, Ratpoints.Double.ofSignificand: the double of a significand and the exponent of its last place, the last step of both roundings.
• Ratpoints.Double.ofInt64, Ratpoints.Double.toInt64, Ratpoints.Double.toUInt64: the conversions between longs and doubles.
• Ratpoints.Double.ldexp, Ratpoints.Double.frexp, Ratpoints.Double.nextUp, Ratpoints.Double.nextDown: the functions of the C library.
• Ratpoints.Double.toFloat, Ratpoints.Double.ofFloat: a double as a value of Lean's Float, for the cost model of the choice of the moduli, and back, for the report of the library.
• Ratpoints.Double.printFixed: a double as printf prints it under %.pf.

structure Ratpoints.Double : Type

A C double, by its IEEE 754 binary64 bits: the sign bit (bit 63), the biased exponent (bits 52 to 62) and the fraction (bits 0 to 51). Equality of the bits is not the == of the C language, under which -0.0 equals 0.0 and a NaN equals nothing; that is Ratpoints.Double.beq.

• bits : UInt64

  The 64 bits of the double.

@[instance_reducible]

instance Ratpoints.instInhabitedDouble : Inhabited Double

Equations

• Ratpoints.instInhabitedDouble = { default := Ratpoints.instInhabitedDouble.default }

def Ratpoints.instReprDouble.repr : Double → Nat → Std.Format

Equations

• Ratpoints.instReprDouble.repr x✝ prec✝ = Std.Format.bracket "{ " (Std.Format.nil ++ Std.Format.text "bits" ++ Std.Format.text " := " ++ (Std.Format.nest 8 (repr x✝.bits)).group) " }"

@[instance_reducible]

instance Ratpoints.instReprDouble : Repr Double

Equations

• Ratpoints.instReprDouble = { reprPrec := Ratpoints.instReprDouble.repr }

def Ratpoints.Double.mantDig : Nat

DBL_MANT_DIG of float.h, the number of bits of the significand of a double, 53: a fact of the format, not a parameter of ratpoints.

Equations

• Ratpoints.Double.mantDig = 53

def Ratpoints.Double.zero : Double

0.0.

Equations

• Ratpoints.Double.zero = { bits := 0 }

def Ratpoints.Double.huge : Double

HUGE_VAL of math.h, which is +∞ for the doubles of IEEE 754.

Equations

• Ratpoints.Double.huge = { bits := 9218868437227405312 }

def Ratpoints.Double.twoPow63 : Double

The double 2 ^ 63, the literal 9223372036854775808.0 of the C code.

Equations

• Ratpoints.Double.twoPow63 = { bits := 4890909195324358656 }

def Ratpoints.Double.minusOne : Double

The double -1.0, the literal -1.0 with which main.c sets the tuning fields of the arguments that are doubles, meaning "the default" (Ratpoints.SieveOptions).

Equations

• Ratpoints.Double.minusOne = { bits := 13830554455654793216 }

def Ratpoints.Double.toFloat (x : Double) : Float

The double x as a value of Lean's Float, the IEEE 754 binary64 type with the arithmetic of the platform, by its bits: the doubles that the caller gives the cost model of the choice of the moduli enter its arithmetic so, the tuning values of the options of the sieve (Ratpoints.SieveOptions) and the ends of the search intervals (Ratpoints.numeratorsFor).

Equations

• x.toFloat = Float.ofBits x.bits

def Ratpoints.Double.ofFloat (x : Float) : Double

The value x of Lean's Float as a double, by its bits: the doubles of the cost model of the choice of the moduli that the report of the library prints (Ratpoints.Report.choice, Ratpoints.Report.moduli).

Equations

• Ratpoints.Double.ofFloat x = { bits := x.toBits }

def Ratpoints.Double.signBit (x : Double) : Bool

The sign bit.

Equations

• x.signBit = (x.bits >>> 63 == 1)

def Ratpoints.Double.expField (x : Double) : Nat

The biased exponent, from 0 (zeros and subnormals) to 2047 (infinities and NaNs).

Equations

• x.expField = (x.bits >>> 52 &&& 2047).toNat

def Ratpoints.Double.fracField (x : Double) : Nat

The fraction, the 52 low bits.

Equations

• x.fracField = (x.bits &&& 4503599627370495).toNat

def Ratpoints.Double.isNaN (x : Double) : Bool

Whether x is a NaN, as isnan of math.h decides it: the largest exponent and a fraction that is not zero.

Equations

• x.isNaN = (x.expField == 2047 && x.fracField != 0)

def Ratpoints.Double.isInf (x : Double) : Bool

Whether x is +∞ or -∞, as isinf of math.h decides it: the largest exponent and the fraction zero.

Equations

• x.isInf = (x.expField == 2047 && x.fracField == 0)

def Ratpoints.Double.isZero (x : Double) : Bool

Whether x is 0.0 or -0.0.

Equations

• x.isZero = (x.bits &&& 9223372036854775807 == 0)

def Ratpoints.Double.isFinite (x : Double) : Bool

Whether x is finite, as isfinite of math.h decides it: an exponent below the largest, so neither an infinity nor a NaN.

Equations

• x.isFinite = (x.expField != 2047)

def Ratpoints.Double.orderValue (x : Double) : Int

The key of a double that is not a NaN in the order of the values: the bits without the sign bit, negated when the sign bit is set, which increase with the value for doubles of either sign, the infinities included, and are 0 for both zeros.

Equations

• x.orderValue = if x.signBit = true then -↑(x.bits &&& 9223372036854775807).toNat else ↑(x.bits &&& 9223372036854775807).toNat

def Ratpoints.Double.lt (x y : Double) : Bool

x < y of the C language on doubles: false when either is a NaN, and otherwise the order of the values, in which -0.0 and 0.0 are equal.

Equations

• x.lt y = (!x.isNaN && !y.isNaN && decide (x.orderValue < y.orderValue))

def Ratpoints.Double.le (x y : Double) : Bool

x <= y of the C language on doubles: false when either is a NaN, and otherwise the order of the values, in which -0.0 and 0.0 are equal.

Equations

• x.le y = (!x.isNaN && !y.isNaN && decide (x.orderValue ≤ y.orderValue))

def Ratpoints.Double.beq (x y : Double) : Bool

x == y of the C language on doubles: false when either is a NaN, and otherwise the equality of the values, under which -0.0 and 0.0 are equal.

Equations

• x.beq y = (!x.isNaN && !y.isNaN && x.orderValue == y.orderValue)

def Ratpoints.Double.neg (x : Double) : Double

The negation -x of the C language on a double: the sign bit flipped, also for a zero, an infinity and a NaN.

Equations

• x.neg = { bits := x.bits ^^^ 9223372036854775808 }

def Ratpoints.Double.fabs (x : Double) : Double

fabs of math.h: the sign bit cleared.

Equations

• x.fabs = { bits := x.bits &&& 9223372036854775807 }

def Ratpoints.Double.toDyadic (x : Double) : Int × Int

The exact value m 2 ^ e of a finite double x, as the pair (m, e): for a normal x the significand 2 ^ 52 + f of the fraction f with the exponent less 1075, for a subnormal x or a zero the fraction with e = -1074, and m negated when the sign bit is set (m = 0 for both zeros). Not meant for infinities and NaNs.

Equations

• x.toDyadic = match if x.expField = 0 then (x.fracField, -1074) else (x.fracField + 2 ^ 52, ↑x.expField - 1075) with | (m, e) => (if x.signBit = true then -↑m else ↑m, e)

def Ratpoints.Double.toRat (x : Double) : Rat

The exact value of a finite double x as a rational number, m 2 ^ e for the pair (m, e) of Ratpoints.Double.toDyadic, as mpq_set_d of GMP sets it: 0 for both zeros. Not meant for infinities and NaNs.

Equations

• x.toRat = match x.toDyadic with | (m, e) => if 0 ≤ e then ↑(m * 2 ^ e.toNat) else mkRat m (2 ^ (-e).toNat)

def Ratpoints.Double.roundShift (m : Nat) (k : Int) : Nat

m / 2 ^ k rounded to the nearest integer, ties to even, and m 2 ^ (-k) for k ≤ 0.

Equations

• One or more equations did not get rendered due to their size.

def Ratpoints.Double.assemble (neg : Bool) (M : Nat) (q : Int) : Double

The double with the sign - for neg that IEEE 754 encodes by a significand M < 2 ^ 53 and the exponent q of its last place: a subnormal or a zero with the fraction M for M < 2 ^ 52 (whose last place is 2 ^ (-1074)), an infinity when the exponent field q + 1075 would reach 2047, and otherwise a normal double with the exponent field q + 1075 and the fraction M - 2 ^ 52.

Equations

• One or more equations did not get rendered due to their size.

def Ratpoints.Double.ofSignificand (neg : Bool) (M : Nat) (q : Int) : Double

The double with the sign - for neg of a rounded significand M ≤ 2 ^ 53 and the exponent q of its last place, the last step of Ratpoints.Double.ofDyadic and Ratpoints.Double.ofRatio: the significand 2 ^ 53, which rounding up can reach, carried into the next exponent as 2 ^ 52, then the bits of Ratpoints.Double.assemble.

Equations

• Ratpoints.Double.ofSignificand neg M q = match if M = 2 ^ 53 then (2 ^ 52, q + 1) else (M, q) with | (M, q) => Ratpoints.Double.assemble neg M q

def Ratpoints.Double.ofDyadic (neg : Bool) (m : Nat) (e : Int) : Double

The double nearest to ±m 2 ^ e, the sign - for neg, with ties to even, as the operations of IEEE 754 round their exact results in the rounding mode to nearest: a signed zero for m = 0; otherwise the last place 2 ^ q of the result, q = max (E - 52) (-1074) for the exponent E with 2 ^ E ≤ m 2 ^ e < 2 ^ (E + 1); M, the value divided by 2 ^ q and rounded, carried into the next exponent at 2 ^ 53; a subnormal or a zero below 2 ^ 52; and an infinity beyond the largest finite double (Ratpoints.Double.ofSignificand). The proofs show only that a finite result for m ≠ 0 is a faithful rounding: the value lies strictly between the neighbours of the result (Ratpoints.Double.lt_nextUp_and_nextDown_lt_ofDyadic).

Equations

• One or more equations did not get rendered due to their size.

def Ratpoints.Double.ofRatio (neg : Bool) (n d : Nat) : Double

The double nearest to ±n / d for d ≠ 0, the sign - for neg, with ties to even, as IEEE 754 rounds an exact result in the rounding mode to nearest: a signed zero for n = 0; otherwise the exponent E with 2 ^ E ≤ n / d < 2 ^ (E + 1), which the lengths of n and d in bits fix up to one, the last place 2 ^ q of the result, q = max (E - 52) (-1074), the quotient n / (d 2 ^ q) rounded to the nearest integer M (r the remainder), carried into the next exponent at 2 ^ 53, a subnormal or a zero below 2 ^ 52, and an infinity beyond the largest finite double (Ratpoints.Double.ofSignificand). The proofs show only that a finite result for n ≠ 0 and d ≠ 0 is a faithful rounding: the value lies strictly between the neighbours of the result (Ratpoints.Double.lt_nextUp_and_nextDown_lt_ofRatio).

Equations

• One or more equations did not get rendered due to their size.

def Ratpoints.Double.ofInt64 (n : Int64) : Double

The conversion (double)n of a long n to a double: the double nearest to n, ties to even, as gcc converts it on the platforms of ratpoints, exact for |n| ≤ 2 ^ 53, and ±2 ^ 63 for n = ±(2 ^ 63 - 1).

Equations

• Ratpoints.Double.ofInt64 n = Ratpoints.Double.ofDyadic (decide (n < 0)) n.toInt.natAbs 0

def Ratpoints.Double.toInt64 (x : Double) : Int64

The conversion (long)x of a finite double x with |x| < 2 ^ 63 to a long: the value truncated towards zero, as the C language defines it for the doubles whose truncation is a long (for the others the conversion is undefined, and the value here is of no meaning).

Equations

• x.toInt64 = match x.toDyadic with | (m, e) => have a := if 0 ≤ e then m.natAbs <<< e.toNat else m.natAbs >>> (-e).toNat; Int64.ofInt (if m < 0 then -↑a else ↑a)

def Ratpoints.Double.toUInt64 (x : Double) : UInt64

The conversion (unsigned long)x of a finite double 0 ≤ x < 2 ^ 64 to an unsigned long: the value truncated towards zero, as the C language defines it for those doubles (for the others the conversion is undefined, and the value here is of no meaning).

Equations

• x.toUInt64 = match x.toDyadic with | (m, e) => (if 0 ≤ e then m.natAbs <<< e.toNat else m.natAbs >>> (-e).toNat).toUInt64

def Ratpoints.Double.ldexp (x : Double) (n : Int) : Double

ldexp(x, n) of math.h, x 2 ^ n: a zero, an infinity and a NaN as they are, as glibc returns them, but for a signalling NaN, which glibc returns quieted (by x + x; ratpoints passes it no NaN); otherwise the exact value of x times 2 ^ n rounded to a double (Ratpoints.Double.ofDyadic), which is exact unless it leaves the normal range, an infinity on overflow and rounded to nearest in the subnormal range, as glibc computes it.

Equations

• x.ldexp n = if (x.isNaN || x.isInf || x.isZero) = true then x else match x.toDyadic with | (m, e) => Ratpoints.Double.ofDyadic x.signBit m.natAbs (e + n)

def Ratpoints.Double.frexp (x : Double) : Double × Int

frexp(x, &ex) of math.h: the pair (f, ex) with x = f 2 ^ ex and 1/2 ≤ |f| < 1, f of the sign of x, for a finite x ≠ 0, also a subnormal one; (x, 0) for a zero, an infinity and a NaN, as glibc returns them, but for a signalling NaN, which glibc returns quieted (by x + x; ratpoints passes it no NaN).

Equations

• One or more equations did not get rendered due to their size.

def Ratpoints.Double.nextUp (x : Double) : Double

nextafter(x, HUGE_VAL) of math.h: the least double above x; the least subnormal for both zeros, +∞ for +∞ and the largest finite double, a NaN as it is (glibc returns a signalling NaN quieted, by x + y; ratpoints passes it no NaN), and for a negative x the double towards zero, -0.0 above -2 ^ (-1074).

Equations

• One or more equations did not get rendered due to their size.

def Ratpoints.Double.nextDown (x : Double) : Double

nextafter(x, -HUGE_VAL) of math.h: the greatest double below x; the least subnormal negated for both zeros, -∞ for -∞ and the least finite double, a NaN as it is (glibc returns a signalling NaN quieted, by x + y; ratpoints passes it no NaN), and for a positive x the double towards zero, 0.0 below 2 ^ (-1074).

Equations

• One or more equations did not get rendered due to their size.

def Ratpoints.Double.printFixed (p : Nat) (x : Double) : String

What printf of glibc prints for the double x under the conversion %.pf, of which %f and %lf are %.6f: the exact value of x rounded to p places after the point, ties to even, in full however large it is, without the point for p = 0, after a - whenever the sign bit is set, also for -0.0 and for a negative value that rounds to 0; inf for an infinity and nan for a NaN, after a - when the sign bit is set.

Equations

• One or more equations did not get rendered due to their size.

Examples

The operations on the edge cases that matter (both zeros, the subnormals, the largest finite doubles, the infinities, a NaN, 2 ^ 53 and 2 ^ 63 with their neighbours, ties of the rounding), against what a C program compiled with gcc and linked with glibc on the platform of ratpoints prints for them, as bits.