kevinclancy / levitate

Utilities for verified floating-point computation in Haskell: Controlled rounding, interval arithmetic, and affine arithmetic.

levitate

Utilities for verified floating-point computation in Haskell: Controlled rounding, interval arithmetic, and affine arithmetic.

Prerequisites

• ghc version >= 7.10.2
• cabal version >= 1.22.6.0

Installation

Enter the following console commands:

git clone https://github.com/kevinclancy/levitate/
cd levitate
cabal sandbox init
cabal install
cabal run tests

Overview

This library is intended for the implementation of experimental numerical verification tools. As such, it is optimized for readability and simplicity rather than performance. Each controlled rounding operator records the current rounding mode, changes to the specified one, performs the operation, and then restores the original rounding mode. It isn't very efficient, but it allows a simple programming interface.

Controlled rounding

The Levitate module provides double-precision floating-point operators with controlled rounding. (+↑) is used for upward-rounded addtion, (+↓) for downward-rounded addition, etc. For a smooth experience using these operators, you will want an editor with unicode support. I use xah math input mode.

In addition to basic arithmetic operators, controlled rounding for the exponent, logarithm, and square root functions is provided.

Interval and affine arithmetic.

The interval and affine arithmetic schemes described in Self-validated Numerical Methods and Applications have been implemented in the Interval and Affine modules, respectively.

The Interval.Interval datatype has a constructor E for empty intervals, and another constructor I for non-empty intervals. The low endpoint of an interval is accessed using the lp function and the high endpoint is accessed using the hp function. Basic arithmetic operators are suffixed with !: (+!) is used for interval addition, (/!) is used for interval division, etc.

levitate's Issues

Question about multiplication

levitate/src/Affine.hs

Lines 226 to 244 in 8dfcf39

	(*@) :: Affine -> Affine -> AffineMonad Affine
	(*@) Empty _ = return Empty
	(*@) _ Empty = return Empty
	(*@) RealLine _ = return RealLine
	(*@) _ RealLine = return RealLine
	(*@) x@(Affine cx _) y@(Affine cy _) = do
	let rx = radius x
	let ry = radius y
	let δ = rx *↑ ry
	let p = (Interval.make cx cx) *! (Interval.make cy cy)
	let δ' = δ +↑ (Interval.radius p)
	aζ <- makeFreshAffine (-(Interval.midpoint p)) []
	aδ <- makeFreshAffine δ' []
	a1 <- cy ·@ x
	a2 <- cx ·@ y
	a3 <- a1 +@ a2
	a4 <- a3 +@ aζ
	a5 <- a4 +@ aδ
	return a5

Hi,

as I understand from this piece of code, for two affines

x = x0 + x1 . e1
y = y0 + y1 . e2

their multiplication are calculated as:

z = x0y0 + y0x0 - x0y0 + y0x1e1 + x0y1e2 + |x1||y1|
  = (x0y0 + |x1||y1|) + y0x1e1 + x0y1e2

But by reading Stolffi reference, it should be:

z = x0y0  + y0x1e1 + x0y1e2 + |x1||y1|e3

In my tests your implementation seems to lead to more accurate results. Do you have any reference for that?