Iterated quaternion multiplication
This example originates from a real industrial code using quaternions to represent rotations in in 3D space. It is presented in details in the paper Verification of the Functional Behavior of a Floating-Point Program: an Industrial Case Study.
Authors: Claude Marché
Topics: Floating-Point Computations
Tools: Frama-C / Jessie / Gappa
see also the index (by topic, by tool, by reference, by year)
#pragma JessieTerminationPolicy(user)
/*@ axiomatic PowerInt {
logic real power(real x, integer n);
}
*/
#pragma JessieBuiltin(power, "\\real_int_pow")
/*********** Quaternions: types, equality *************/
typedef double quat[4];
typedef double *quat_ptr;
// equality of two quaternions in two different memory states
/*@ predicate quat_eq{L1,L2}(quat_ptr q1, quat_ptr q2) =
@ \at(q1[0],L1) == \at(q2[0],L2) && \at(q1[1],L1) == \at(q2[1],L2)
@ && \at(q1[2],L1) == \at(q2[2],L2) && \at(q1[3],L1) == \at(q2[3],L2);
@*/
/*@ requires \valid(src+(0..3));
@ requires \valid(dst+(0..3));
@ assigns dst[0..3];
@ ensures quat_eq{Here,Here}(dst,src);
@*/
void Quat_copy(const quat src, quat dst) {
dst[0] = src[0];
dst[1] = src[1];
dst[2] = src[2];
dst[3] = src[3];
}
/*** norm ***/
/*@ logic real norm2(real p1, real p2, real p3, real p4) =
@ p1*p1 + p2*p2 + p3*p3 + p4*p4;
@*/
/*@ lemma norm2_pos: \forall real p1,p2,p3,p4;
@ 0.0 <= norm2(p1,p2,p3,p4);
@*/
/* lemma norm2_neg: \forall real p1,p2,p3,p4;
@ norm2(-p1,-p2,-p3,-p4) == norm2(p1,p2,p3,p4);
@*/
/*@ logic real norm_vect(real p1, real p2, real p3, real p4) =
@ \sqrt(norm2(p1,p2,p3,p4));
@*/
/*@ lemma norm_vect_pos: \forall real p1,p2,p3,p4;
@ 0.0 <= norm_vect(p1,p2,p3,p4);
@*/
/*@ logic real quat_norm{L}(quat_ptr q) =
@ norm_vect(q[0],q[1],q[2],q[3]);
@*/
/*@ lemma quat_norm_pos{L}: // needed for Coq proof of the assertion
@ // in the loop
@ \forall quat_ptr q; quat_norm{L}(q) >= 0.0 ;
@*/
/*** distance ***/
/*@ logic real distance2(real p0, real p1, real p2, real p3,
@ real q0, real q1, real q2, real q3) =
@ norm2(q0 - p0, q1 - p1, q2 - p2, q3 - p3) ;
@*/
/*@ lemma distance2_sym :
@ \forall real p0,p1,p2, p3,q0,q1,q2,q3;
@ norm2(q0 - p0, q1 - p1, q2 - p2, q3 - p3) ==
@ norm2(p0 - q0, p1 - q1, p2 - q2, p3 - q3) ;
@*/
/*@ logic real distance_quat_vect{L}(quat_ptr q, real p0, real p1,
@ real p2, real p3) =
@ \sqrt(distance2(q[0], q[1], q[2], q[3], p0, p1, p2, p3)) ;
@*/
/*************** Quaternions: product ******************/
/*@ logic real product1{L}(quat_ptr q1, quat_ptr q2) =
@ q1[0]*q2[0] - q1[1]*q2[1] - q1[2]*q2[2] - q1[3]*q2[3] ;
@
@ logic real product2{L}(quat_ptr q1, quat_ptr q2) =
@ q1[0]*q2[1] + q1[1]*q2[0] + q1[2]*q2[3] - q1[3]*q2[2] ;
@
@ logic real product3{L}(quat_ptr q1, quat_ptr q2) =
@ q1[0]*q2[2] - q1[1]*q2[3] + q1[2]*q2[0] + q1[3]*q2[1] ;
@
@ logic real product4{L}(quat_ptr q1, quat_ptr q2) =
@ q1[0]*q2[3] + q1[1]*q2[2] - q1[2]*q2[1] + q1[3]*q2[0] ;
@*/
/*@ predicate bounded(real x,real k) = -k <= x <= k;
@*/
/*@ lemma bounded_sqr:
@ \forall real x,k; 0.0 <= k && x*x <= k*k ==> bounded(x,k);
@*/
/*@ lemma bounded_norm2:
@ \forall real p1,p2,p3,p4,k;
@ 0.0 <= k && norm2(p1,p2,p3,p4) <= k * k ==>
@ bounded(p1,k) &&
@ bounded(p2,k) &&
@ bounded(p3,k) &&
@ bounded(p4,k);
@*/
/*@ lemma sqrt_le_le_sqr :
@ \forall real x,y; 0.0 <= x && \sqrt(x) <= y ==> x <= y * y;
@*/
/*@ lemma bounded_norm_vect:
@ \forall real p1,p2,p3,p4,k;
@ 0.0 <= k && norm_vect(p1,p2,p3,p4) <= k ==>
@ bounded(p1,k) &&
@ bounded(p2,k) &&
@ bounded(p3,k) &&
@ bounded(p4,k);
@*/
/*@ lemma bounded_norm{L}: // needed for proving Quat_prod
@ \forall quat_ptr q, real k; quat_norm(q) <= k ==>
@ bounded(q[0],k) &&
@ bounded(q[1],k) &&
@ bounded(q[2],k) &&
@ bounded(q[3],k);
@*/
//@ logic real EPS0 = 0x3p-50;
//@ logic real BOUND = 1.125;
/*@ requires \valid(q1+(0..3));
@ requires \valid(q2+(0..3));
@ requires \valid(q+(0..3));
@ requires quat_norm(q1) <= BOUND;
@ requires quat_norm(q2) <= BOUND;
@ assigns q[0..3];
@ ensures
@ distance_quat_vect(q,product1(q1,q2),product2(q1,q2),
@ product3(q1,q2),product4(q1,q2)) <= EPS0;
@*/
void Quat_prod(const quat q1, const quat q2, quat q) {
// asserts needed for Gappa, that can be deduced from lemma
// bounded_norm
//@ assert bounded(q1[0],BOUND);
//@ assert bounded(q1[1],BOUND);
//@ assert bounded(q1[2],BOUND);
//@ assert bounded(q1[3],BOUND);
//@ assert bounded(q2[0],BOUND);
//@ assert bounded(q2[1],BOUND);
//@ assert bounded(q2[2],BOUND);
//@ assert bounded(q2[3],BOUND);
q[0] = q1[0]*q2[0] - q1[1]*q2[1] - q1[2]*q2[2] - q1[3]*q2[3];
q[1] = q1[0]*q2[1] + q1[1]*q2[0] + q1[2]*q2[3] - q1[3]*q2[2];
q[2] = q1[0]*q2[2] - q1[1]*q2[3] + q1[2]*q2[0] + q1[3]*q2[1];
q[3] = q1[0]*q2[3] + q1[1]*q2[2] - q1[2]*q2[1] + q1[3]*q2[0];
}
/*************** Quaternions: random ******************/
//@ logic real EPS1 = 0x3p-53;
/*@ requires \valid(q+(0..3));
@ assigns q[0..3];
@ ensures \abs(quat_norm(q) - 1.0) <= EPS1;
@*/
void random_unit_quat(quat q);
/*************** Main loop and property ******************/
/*@ lemma norm_product{L}:
@ \forall quat_ptr q1,q2;
@ \let p1 = product1(q1,q2);
@ \let p2 = product2(q1,q2);
@ \let p3 = product3(q1,q2);
@ \let p4 = product4(q1,q2);
@ norm_vect(p1,p2,p3,p4) == quat_norm(q1) * quat_norm(q2);
@*/
//@ logic real EPS = EPS0 + EPS1;
//@ logic integer MAX = 35000000000000; // 35 x 10^{12}
// Two lemmas needed for proving pre-condition of call to Quat_prod
/*@ lemma pow_eps_max_int:
@ power(1.0 + EPS, MAX) <= BOUND;
@
@ lemma pow_eps2_max_int:
@ power(1.0 + EPS, 10000000000) <= BOUND;
@
@ lemma power_monotonic:
@ \forall integer n,m, real x;
@ 0 <= n <= m && 1.0 <= x ==> power(x,n) <= power(x,m);
@*/
/*@ lemma triangle_inequality:
@ \forall real p1,p2,p3,p4,q1,q2,q3,q4;
@ norm_vect(p1+q1,p2+q2,p3+q3,p4+q4) <=
@ norm_vect(p1,p2,p3,p4) + norm_vect(q1,q2,q3,q4);
@*/
/*@ lemma norm_distance_inequality:
@ \forall real p1,p2,p3,p4,q1,q2,q3,q4;
@ \sqrt(norm2(p1,p2,p3,p4)) <=
@ \sqrt(distance2(p1,p2,p3,p4,q1,q2,q3,q4)) + \sqrt(norm2(q1,q2,q3,q4));
@*/
int test1(void) {
quat current, next, incr;
random_unit_quat(current);
//@ ghost long long n = 1; // integer would be better !
/*@ loop invariant
@ 0 <= n <= MAX;
@ loop invariant
@ quat_norm(current+0) <= power(1.0 + EPS,n);
@*/
while (1) {
random_unit_quat(incr);
Quat_prod(current,incr,next);
/*@ assert // we restate the post-condition of random_unit_quat
@ // just to simplify the Coq proof below
@ quat_norm(incr+0) <= 1.0 + EPS1 ;
@*/
/*@ assert // we restate the loop invariant after the product
@ // for the same reason
@ quat_norm(current+0) <= power(1.0 + EPS,n) ;
@*/
/*@ assert // Coq is needed for this one, because ATPs
@ // are too weak with multiplication
@ quat_norm(current+0) * quat_norm(incr+0) <=
@ power(1.0 + EPS,n) * (1.0 + EPS1) ;
@*/
Quat_copy(next, current);
//@ ghost n++ ;
//@ assert n <= MAX;
}
return 0;
}
/*i
Local Variables:
compile-command: "frama-c -jessie quat2.c"
End:
*/
download ZIP archive