Commit 1869fcc2 authored by Heiko Becker's avatar Heiko Becker
Browse files

Fix bug in HOL4 certificates and add more testcases for formal development

parent c6af90f5
INCLUDES = $(HOLDIR)/examples/machine-code/hoare-triple ../
INCLUDES = $(HOLDIR)/examples/machine-code/hoare-triple ../ ../Infra
OPTIONS = QUIT_ON_FAILURE
ifdef POLY
......
......@@ -101,22 +101,26 @@ object CertificatePhase extends DaisyPhase {
if (prover == "coq")
"Require Import Daisy.CertificateChecker."
else if (prover == "hol4")
"open preamble abbrevsTheory CertificateCheckerTheory;\nopen simpLib realTheory realLib RealArith;\n\n" +
"open preamble AbbrevsTheory CertificateCheckerTheory;\nopen simpLib realTheory realLib RealArith;\n\n" +
"val _ = new_theory \"" + fname +"\";\n\n"
else
"needs \"CertificateChecker.hl\";;\nneeds \"Infra/convs.hl\";;"
private def coqVariable(vname:Identifier, id:Int, reporter:Reporter) :(String, String) =
{
val theExpr = s"Definition ExpVar$vname :exp Q := Var Q $id.\n"
(theExpr,s"ExpVar$vname")
}
{
//FIXME: Ugly Hack to get disjoint names for multiple function encodings with same variable names:
val freshId = nextFreshVariable()
val theExpr = s"Definition ExpVar$vname$freshId :exp Q := Var Q $id.\n"
(theExpr,s"ExpVar$vname$freshId")
}
private def hol4Variable(vname:Identifier, id:Int, reporter:Reporter) :(String, String) =
{
val theExpr = s"val ExpVar$vname = Define `ExpVar$vname:(real exp) = Var $id`;\n"
(theExpr, s"ExpVar$vname")
}
{
//FIXME: Ugly Hack to get disjoint names for multiple function encodings with same variable names:
val freshId = nextFreshVariable()
val theExpr = s"val ExpVar$vname${freshId}_def = Define `ExpVar$vname$freshId:(real exp) = Var $id`;\n"
(theExpr, s"ExpVar$vname$freshId")
}
private def holLightVariable(vname:Identifier, id:Int, reporter:Reporter) :(String, String) =
{
......@@ -127,18 +131,22 @@ object CertificatePhase extends DaisyPhase {
private def coqConstant(r:RealLiteral, id:Int, reporter:Reporter) :(String, String) =
r match {
case RealLiteral(v) =>
//FIXME: Ugly Hack to get disjoint names for multiple function encodings with same variable names:
val freshId = nextConstantId()
val rationalStr = v.toFractionString
val coqRational = rationalStr.replace('/','#')
val theExpr = s"Definition ExpCst$id :exp Q := Const ($coqRational).\n"
(theExpr, s"ExpCst$id")
val theExpr = s"Definition ExpCst$id$freshId :exp Q := Const ($coqRational).\n"
(theExpr, s"ExpCst$id$freshId")
}
private def hol4Constant(r:RealLiteral, id:Int, reporter:Reporter) :(String, String) =
r match {
case RealLiteral(v) =>
//FIXME: Ugly Hack to get disjoint names for multiple function encodings with same variable names:
val freshId = nextConstantId()
val rationalStr = v.toFractionString
val theExpr = s"val ExpCst$id = Define `ExpCst$id:(real exp) = Const ($rationalStr)`;\n"
(theExpr, s"ExpCst$id")
val theExpr = s"val ExpCst$id${freshId}_def = Define `ExpCst$id$freshId:(real exp) = Const ($rationalStr)`;\n"
(theExpr, s"ExpCst$id$freshId")
}
private def holLightConstant(r:RealLiteral, id:Int, reporter:Reporter) :(String, String) =
......@@ -501,7 +509,7 @@ object CertificatePhase extends DaisyPhase {
val errY = errorMap(Variable (y)).toFractionString.replace("/","#")
val nameY = expressionNames(Variable(y))
conditional (
s"( expEqBool e Var ${identifierNums(y)} )",
s"( expEqBool e (Var Q ${identifierNums(y)}) )",
"(" + intvY + ", " + errY + ")",
gFun)
......@@ -572,10 +580,10 @@ object CertificatePhase extends DaisyPhase {
rFun)
case x @ Let (y,exp, g) =>
val expFun = coqAbsEnv (exp, errorMap, rangeMap, cont)
val gFun = coqAbsEnv (g, errorMap, rangeMap, expFun)
val intvY = coqInterval((rangeMap(Variable(y)).xlo, rangeMap(Variable(y)).xhi))
val errY = errorMap(Variable (y)).toFractionString.replace("/","#")
val expFun = hol4AbsEnv (exp, errorMap, rangeMap, cont)
val gFun = hol4AbsEnv (g, errorMap, rangeMap, expFun)
val intvY = hol4Interval((rangeMap(Variable(y)).xlo, rangeMap(Variable(y)).xhi))
val errY = errorMap(Variable (y)).toFractionString
val nameY = expressionNames(Variable(y))
conditional (
s"( e = Var ${identifierNums(y)} )",
......@@ -667,7 +675,7 @@ object CertificatePhase extends DaisyPhase {
analysisResultName + " " + precondName + " = true.\n" +
"Proof.\n cbv; auto.\nQed."
else if (prover == "hol4")
"val _ = store_thm (\"ErrorBound_${fName}_Sound\",\n ``CertificateCheckerCmd " +
"val _ = store_thm (\"" + s"ErrorBound_${fName}_Sound" + "\",\n ``CertificateCheckerCmd " +
lastGenName + " " + analysisResultName + " " + precondName + " = T``,\n EVAL_TAC);"
else
"DAISY_CONV CertificateChecker thePrecondition theRewrites `CertificateChecker " + lastGenName + " " + analysisResultName + " " + precondName + "`;;"
......
import daisy.lang._
import Real._
/**
Equation and initial ranges from:
L. Zhang, Y. Zhang, and W. Zhou. Tradeoff between Approximation Accuracy and
Complexity for Range Analysis using Affine Arithmetic.
*/
object Bsplines {
def bspline0(u: Real): Real = {
require(0 <= u && u <= 1)
(1 - u) * (1 - u) * (1 - u) / 6.0
} ensuring (res => 0 <= res && res <= 0.17 && res +/- 1e-15)
// proven in paper: [-0.05, 0.17]
def bspline1(u: Real): Real = {
require(0 <= u && u <= 1)
(3 * u*u*u - 6 * u*u + 4) / 6.0
} ensuring (res => 0.16 <= res && res <= 0.7 && res +/- 1e-15)
// in paper [-0.05, 0.98]
def bspline2(u: Real): Real = {
require(0 <= u && u <= 1)
(-3 * u*u*u + 3*u*u + 3*u + 1) / 6.0
} ensuring (res => 0.16 <= res && res <= 0.7 && res +/- 1e-15)
// in paper [-0.02, 0.89]
def bspline3(u: Real): Real = {
require(0 <= u && u <= 1)
-u*u*u / 6.0
} ensuring (res => -0.17 <= res && res <= 0.0 && res +/- 1e-15)
// in paper [-0.17, 0.05]
}
import daisy.lang._
import Real._
object Doppler {
def doppler(u: Real, v: Real, T: Real): Real = {
require(-100.0 <= u && u <= 100 && 20 <= v && v <= 20000 && -30 <= T && T <= 50)
val t1 = 331.4 + 0.6 * T
(- (t1) *v) / ((t1 + u)*(t1 + u))
}
}
\ No newline at end of file
import daisy.lang._
import Real._
object JetEngine {
def jetEngine(x1: Real, x2: Real): Real = {
require(-5 <= x1 && x1 <= 5 && -20 <= x2 && x2 <= 5)
val t = (3*x1*x1 + 2*x2 - x1)
x1 + ((2*x1*(t/(x1*x1 + 1))*
(t/(x1*x1 + 1) - 3) + x1*x1*(4*(t/(x1*x1 + 1))-6))*
(x1*x1 + 1) + 3*x1*x1*(t/(x1*x1 + 1)) + x1*x1*x1 + x1 +
3*((3*x1*x1 + 2*x2 -x1)/(x1*x1 + 1)))
}
}
\ No newline at end of file
import daisy.lang._
import Real._
//Unrolled loops
object Pendulum {
def pendulumT1(t_0: Real, w_0: Real): Real = {
require(-2 <= t_0 && t_0 <= 2 && -5 <= w_0 && w_0 <= 5)
val h: Real = 0.01
val L: Real = 2.0
val m: Real = 1.5
val g: Real = 9.80665
val k1t = w_0
val k1w = -g/L * (t_0 - t_0 * t_0 * t_0/6 + t_0 * t_0 * t_0 * t_0 * t_0/120)
val k2t = w_0 + h/2*k1w
val arg = t_0 + h/2*k1t
val k2w = -g/L * (arg - arg * arg * arg/6 + arg * arg * arg * arg * arg/120)
t_0 + h*k2t
}
def pendulumW1(t_0: Real, w_0: Real): Real = {
require(-2 <= t_0 && t_0 <= 2 && -5 <= w_0 && w_0 <= 5)
val h: Real = 0.01
val L: Real = 2.0
val m: Real = 1.5
val g: Real = 9.80665
val k1t = w_0
val k1w = -g/L * (t_0 - t_0 * t_0 * t_0/6 + t_0 * t_0 * t_0 * t_0 * t_0/120)
val k2t = w_0 + h/2*k1w
val arg = t_0 + h/2*k1t
val k2w = -g/L * (arg - arg * arg * arg/6 + arg * arg * arg * arg * arg/120)
w_0 + h*k2w
}
// TODO: takes forever, maybe because of Rationals getting too big?
def pendulumT3(t_0: Real, w_0: Real): Real = {
require(-2 <= t_0 && t_0 <= 2 && -5 <= w_0 && w_0 <= 5)
val h: Real = 0.01
val L: Real = 2.0
val m: Real = 1.5
val g: Real = 9.80665
// loop 1
val k1t_1 = w_0
val k1w_1 = -g/L * (t_0 - t_0 * t_0 * t_0/6 + t_0 * t_0 * t_0 * t_0 * t_0/120)
val k2t_1 = w_0 + h/2*k1w_1
val arg_1 = t_0 + h/2*k1t_1
val k2w_1 = -g/L * (arg_1 - arg_1 * arg_1 * arg_1/6 + arg_1 * arg_1 * arg_1 * arg_1 * arg_1/120)
val t_1 = t_0 + h*k2t_1
val w_1 = w_0 + h*k2w_1
// loop 2
val k1t_2 = w_1
val k1w_2 = -g/L * (t_1 - t_1 * t_1 * t_1/6 + t_1 * t_1 * t_1 * t_1 * t_1/120)
val k2t_2 = w_1 + h/2*k1w_2
val arg_2 = t_1 + h/2*k1t_2
val k2w_2 = -g/L * (arg_2 - arg_2 * arg_2 * arg_2/6 + arg_2 * arg_2 * arg_2 * arg_2 * arg_2/120)
val t_2 = t_1 + h*k2t_2
val w_2 = w_1 + h*k2w_2
// loop 3
val k1t_3 = w_2
val k1w_3 = -g/L * (t_2 - t_2 * t_2 * t_2/6 + t_2 * t_2 * t_2 * t_2 * t_2/120)
val k2t_3 = w_2 + h/2*k1w_3
val arg_3 = t_2 + h/2*k1t_3
val k2w_3 = -g/L * (arg_3 - arg_3 * arg_3 * arg_3/6 + arg_3 * arg_3 * arg_3 * arg_3 * arg_3/120)
t_2 + h*k2t_3
}
def pendulumW3(t_0: Real, w_0: Real): Real = {
require(-2 <= t_0 && t_0 <= 2 && -5 <= w_0 && w_0 <= 5)
val h: Real = 0.01
val L: Real = 2.0
val m: Real = 1.5
val g: Real = 9.80665
// loop 1
val k1t_1 = w_0
val k1w_1 = -g/L * (t_0 - t_0 * t_0 * t_0/6 + t_0 * t_0 * t_0 * t_0 * t_0/120)
val k2t_1 = w_0 + h/2*k1w_1
val arg_1 = t_0 + h/2*k1t_1
val k2w_1 = -g/L * (arg_1 - arg_1 * arg_1 * arg_1/6 + arg_1 * arg_1 * arg_1 * arg_1 * arg_1/120)
val t_1 = t_0 + h*k2t_1
val w_1 = w_0 + h*k2w_1
// loop 2
val k1t_2 = w_1
val k1w_2 = -g/L * (t_1 - t_1 * t_1 * t_1/6 + t_1 * t_1 * t_1 * t_1 * t_1/120)
val k2t_2 = w_1 + h/2*k1w_2
val arg_2 = t_1 + h/2*k1t_2
val k2w_2 = -g/L * (arg_2 - arg_2 * arg_2 * arg_2/6 + arg_2 * arg_2 * arg_2 * arg_2 * arg_2/120)
val t_2 = t_1 + h*k2t_2
val w_2 = w_1 + h*k2w_2
// loop 3
val k1t_3 = w_2
val k1w_3 = -g/L * (t_2 - t_2 * t_2 * t_2/6 + t_2 * t_2 * t_2 * t_2 * t_2/120)
val k2t_3 = w_2 + h/2*k1w_3
val arg_3 = t_2 + h/2*k1t_3
val k2w_3 = -g/L * (arg_3 - arg_3 * arg_3 * arg_3/6 + arg_3 * arg_3 * arg_3 * arg_3 * arg_3/120)
w_2 + h*k2w_3
}
def pendulumT5(t_0: Real, w_0: Real): Real = {
require(-2 <= t_0 && t_0 <= 2 && -5 <= w_0 && w_0 <= 5)
val h: Real = 0.01
val L: Real = 2.0
val m: Real = 1.5
val g: Real = 9.80665
// loop 1
val k1t_1 = w_0
val k1w_1 = -g/L * (t_0 - t_0 * t_0 * t_0/6 + t_0 * t_0 * t_0 * t_0 * t_0/120)
val k2t_1 = w_0 + h/2*k1w_1
val arg_1 = t_0 + h/2*k1t_1
val k2w_1 = -g/L * (arg_1 - arg_1 * arg_1 * arg_1/6 + arg_1 * arg_1 * arg_1 * arg_1 * arg_1/120)
val t_1 = t_0 + h*k2t_1
val w_1 = w_0 + h*k2w_1
// loop 2
val k1t_2 = w_1
val k1w_2 = -g/L * (t_1 - t_1 * t_1 * t_1/6 + t_1 * t_1 * t_1 * t_1 * t_1/120)
val k2t_2 = w_1 + h/2*k1w_2
val arg_2 = t_1 + h/2*k1t_2
val k2w_2 = -g/L * (arg_2 - arg_2 * arg_2 * arg_2/6 + arg_2 * arg_2 * arg_2 * arg_2 * arg_2/120)
val t_2 = t_1 + h*k2t_2
val w_2 = w_1 + h*k2w_2
// loop 3
val k1t_3 = w_2
val k1w_3 = -g/L * (t_2 - t_2 * t_2 * t_2/6 + t_2 * t_2 * t_2 * t_2 * t_2/120)
val k2t_3 = w_2 + h/2*k1w_3
val arg_3 = t_2 + h/2*k1t_3
val k2w_3 = -g/L * (arg_3 - arg_3 * arg_3 * arg_3/6 + arg_3 * arg_3 * arg_3 * arg_3 * arg_3/120)
val t_3 = t_2 + h*k2t_3
val w_3 = w_2 + h*k2w_3
// loop 4
val k1t_4 = w_3
val k1w_4 = -g/L * (t_3 - t_3 * t_3 * t_3/6 + t_3 * t_3 * t_3 * t_3 * t_3/120)
val k2t_4 = w_3 + h/2*k1w_4
val arg_4 = t_3 + h/2*k1t_4
val k2w_4 = -g/L * (arg_4 - arg_4 * arg_4 * arg_4/6 + arg_4 * arg_4 * arg_4 * arg_4 * arg_4/120)
val t_4 = t_3 + h*k2t_4
val w_4 = w_3 + h*k2w_4
// loop 5
val k1t_5 = w_4
val k1w_5 = -g/L * (t_4 - t_4 * t_4 * t_4/6 + t_4 * t_4 * t_4 * t_4 * t_4/120)
val k2t_5 = w_4 + h/2*k1w_5
val arg_5 = t_4 + h/2*k1t_5
val k2w_5 = -g/L * (arg_5 - arg_5 * arg_5 * arg_5/6 + arg_5 * arg_5 * arg_5 * arg_5 * arg_5/120)
t_4 + h*k2t_5
}
def pendulumW5(t_0: Real, w_0: Real): Real = {
require(-2 <= t_0 && t_0 <= 2 && -5 <= w_0 && w_0 <= 5)
val h: Real = 0.01
val L: Real = 2.0
val m: Real = 1.5
val g: Real = 9.80665
// loop 1
val k1t_1 = w_0
val k1w_1 = -g/L * (t_0 - t_0 * t_0 * t_0/6 + t_0 * t_0 * t_0 * t_0 * t_0/120)
val k2t_1 = w_0 + h/2*k1w_1
val arg_1 = t_0 + h/2*k1t_1
val k2w_1 = -g/L * (arg_1 - arg_1 * arg_1 * arg_1/6 + arg_1 * arg_1 * arg_1 * arg_1 * arg_1/120)
val t_1 = t_0 + h*k2t_1
val w_1 = w_0 + h*k2w_1
// loop 2
val k1t_2 = w_1
val k1w_2 = -g/L * (t_1 - t_1 * t_1 * t_1/6 + t_1 * t_1 * t_1 * t_1 * t_1/120)
val k2t_2 = w_1 + h/2*k1w_2
val arg_2 = t_1 + h/2*k1t_2
val k2w_2 = -g/L * (arg_2 - arg_2 * arg_2 * arg_2/6 + arg_2 * arg_2 * arg_2 * arg_2 * arg_2/120)
val t_2 = t_1 + h*k2t_2
val w_2 = w_1 + h*k2w_2
// loop 3
val k1t_3 = w_2
val k1w_3 = -g/L * (t_2 - t_2 * t_2 * t_2/6 + t_2 * t_2 * t_2 * t_2 * t_2/120)
val k2t_3 = w_2 + h/2*k1w_3
val arg_3 = t_2 + h/2*k1t_3
val k2w_3 = -g/L * (arg_3 - arg_3 * arg_3 * arg_3/6 + arg_3 * arg_3 * arg_3 * arg_3 * arg_3/120)
val t_3 = t_2 + h*k2t_3
val w_3 = w_2 + h*k2w_3
// loop 4
val k1t_4 = w_3
val k1w_4 = -g/L * (t_3 - t_3 * t_3 * t_3/6 + t_3 * t_3 * t_3 * t_3 * t_3/120)
val k2t_4 = w_3 + h/2*k1w_4
val arg_4 = t_3 + h/2*k1t_4
val k2w_4 = -g/L * (arg_4 - arg_4 * arg_4 * arg_4/6 + arg_4 * arg_4 * arg_4 * arg_4 * arg_4/120)
val t_4 = t_3 + h*k2t_4
val w_4 = w_3 + h*k2w_4
// loop 5
val k1t_5 = w_4
val k1w_5 = -g/L * (t_4 - t_4 * t_4 * t_4/6 + t_4 * t_4 * t_4 * t_4 * t_4/120)
val k2t_5 = w_4 + h/2*k1w_5
val arg_5 = t_4 + h/2*k1t_5
val k2w_5 = -g/L * (arg_5 - arg_5 * arg_5 * arg_5/6 + arg_5 * arg_5 * arg_5 * arg_5 * arg_5/120)
w_4 + h*k2w_5
}
}
import daisy.lang._
import Real._
object RigidBody {
//x1, x2, x3: < 1, 16 11>, [-15, 15]
def out1(x1: Real, x2: Real, x3: Real): Real = {
require(-15 <= x1 && x1 <= 15 && -15 <= x2 && x2 <= 15 && -15 <= x3 && x3 <= 15)
( -1 * x1)*x2 - 2*x2*x3 - x1 - x3
def rigidBody1(x1: Real, x2: Real, x3: Real): Real = {
require(-15.0 <= x1 && x1 <= 15 && -15.0 <= x2 && x2 <= 15.0 && -15.0 <= x3 && x3 <= 15)
-x1*x2 - 2*x2*x3 - x1 - x3
}
// apparently this needs 17 bits, or we've not been accurate enough
def out2(x1: Real, x2: Real, x3: Real): Real = {
require(-15 <= x1 && x1 <= 15 && -15 <= x2 && x2 <= 15 && -15 <= x3 && x3 <= 15)
2*x1*x2*x3 + 3*x3*x3 - x2*x1*x2*x3 + 3*x3*x3 - x2
def rigidBody2(x1: Real, x2: Real, x3: Real): Real = {
require(-15.0 <= x1 && x1 <= 15 && -15.0 <= x2 && x2 <= 15.0 &&
-15.0 <= x3 && x3 <= 15)
2*(x1*x2*x3) + (3*x3*x3) - x2*(x1*x2*x3) + (3*x3*x3) - x2
}
}
}
\ No newline at end of file
import daisy.lang._
import Real._
object Science {
def verhulst(r: Real, K: Real, x: Real): Real = {
require(r >= 4.0 && r <= 4.0 && K >= 1.11 && K <= 1.11 && 0.1 <= x && x <= 0.3)
(r*x) / (1 + (x/K))
}
def predatorPrey(r: Real, K: Real, x: Real): Real = {
require(r >= 4.0 && r <= 4.0 && K >= 1.11 && K <= 1.11 && 0.1 <= x && x <= 0.3)
(r*x*x) / (1 + (x/K)*(x/K))
}
def carbonGas(T: Real, a: Real, b: Real, N: Real, p: Real, V: Real): Real = {
require(T >= 300 && T <= 300 && a >= 0.401 && a <= 0.401 && b >= 42.7e-6 && b <= 42.7e-6 && N >= 1000 && N <= 1000 &&
p >= 3.5e7 && p <= 3.5e7 && 0.1 < V && V < 0.5)
val k: Real = 1.3806503e-23
(p + a * (N / V) * (N / V)) * (V - N * b) - k * N * T
}
}
\ No newline at end of file
import daisy.lang._
import Real._
/*
The famous formula for computing the area of a triangle.
*/
object Triangle {
def triangleUnstable(a: Real, b: Real, c: Real): Real = {
// TODO: this precondition can only be handled with SMT, and only
// if we add the usage of the additional constraints
//require(1 < a && a < 9 && 1 < b && b < 9 && 1 < c && c < 9 &&
// a + b > c + 0.000001 && a + c > b + 0.000001 && b + c > a + 0.000001)
require(4 <= a && a <= 5 && 4 <= b && b <= 4 && 4 <= c && c <= 5)
val s = (a + b + c)/2.0
sqrt(s * (s - a) * (s - b) * (s - c))
} ensuring (res => res >= 0.0 && res +/- 2e-9)
// def triangleKahan(aa: Real, bb: Real, cc: Real): Real = {
// var a = aa
// var b = bb
// var c = cc
// if(b < a) {
// val t = a
// if(c < b) {
// a = c; c = t
// }
// else {
// if(c < a) {
// a = b; b = c; c = t
// }
// else {
// a = b; b = t
// }
// }
// }
// else if(c < b) {
// val t = c; c = b;
// if(c < a) {
// b = a; a = t
// }
// else {
// b = t
// }
// }
// sqrt((a+(b+c)) * (c-(a-b)) * (c+(a-b)) * (a+(b-c))) / 4.0
// }
}
\ No newline at end of file
import daisy.lang._
import Real._
object Turbine {
def turbine1(v: Real, w: Real, r: Real): Real = {
require(-4.5 <= v && v <= -0.3 && 0.4 <= w && w <= 0.9 && 3.8 <= r && r <= 7.8)
3 + 2/(r*r) - 0.125*(3-2*v)*(w*w*r*r)/(1-v) - 4.5
}
def turbine2(v: Real, w: Real, r: Real): Real = {
require(-4.5 <= v && v <= -0.3 && 0.4 <= w && w <= 0.9 && 3.8 <= r && r <= 7.8)
6*v - 0.5 * v * (w*w*r*r) / (1-v) - 2.5
}
def turbine3(v: Real, w: Real, r: Real): Real = {
require(-4.5 <= v && v <= -0.3 && 0.4 <= w && w <= 0.9 && 3.8 <= r && r <= 7.8)
3 - 2/(r*r) - 0.125 * (1+2*v) * (w*w*r*r) / (1-v) - 0.5
}
}
\ No newline at end of file