## 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:
*/