Wiki Agenda Contact Version française

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 to appear in Science of Computer Programming


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