Computation of the trajectory of an object submitted to gravity
This example contains code from an embedded safety-critical software, which program computes the speed of a device submitted to gravitational acceleration and drag from the atmosphere around it.
This program was used as challenge example in the article Automating the Verification of Floating-Point Programs by ClĂ©ment Fumex, Claude MarchĂ©, and Yannick Moy, published at VSTTE 2017 conference. This case study was initially provided by Florian Schanda from Altran, UK.
The Ada files contain multiple variants of the example in increasing order of difficulty:
- simple_trajectory_overflow.ads and simple_trajectory_overflow.adb Simple version of the computation, which only updates the speed by not the distance. Only absence of run-time errors (including overflows) is proved.
- simple_trajectory.ads and simple_trajectory.adb Simple version of the computation, which only updates the speed by not the distance. Both absence of run-time errors (including overflows) and safe bounds on the computed speed are proved.
- complex_trajectory_overflow.ads and complex_trajectory_overflow.adb More complex version of the computation, which updates both the speed and the distance. Only absence of run-time errors (including overflows) is proved, but this depends on the proof of safe bounds on the computed speed.
- complex_trajectory.ads and complex_trajectory.adb More complex version of the computation, which updates both the speed and the distance. Both absence of run-time errors (including overflows) and safe bounds on the computed speed and distance are expressed as properties, but the safe bounds on distance are not proved yet.
Authors: Clément Fumex / Yannick Moy
Topics: Floating-Point Computations
Tools: SPARK 2014
References: ProofInUse joint laboratory
see also the index (by topic, by tool, by reference, by year)
download ZIP archive
The full code is in the ZIP archive above. Here is below the specification and the code for the simple trajectory.
File : simple_trajectory.ads
package Simple_Trajectory with SPARK_Mode is type Float64 is digits 15; type Frame is range 0 .. 25_000; Frame_Length : constant := 1.0 / 60.0; G : constant := 9.807; -- gravitational acceleration subtype Ratio_T is Float64 range -1.0 .. 1.0; subtype Drag_T is Float64 range -64.0 .. 64.0; Bound : constant Float64 := Drag_T'Last + Float64'Ceiling (G * Frame_Length) with Ghost; function Low_Bound (N : Frame) return Float64 is (Float64 (N) * (- Bound)) with Ghost; function High_Bound (N : Frame) return Float64 is (Float64 (N) * Bound) with Ghost; function In_Bounds (V : Float64) return Boolean is (V in - Bound * Float64 (Frame'Last) .. Bound * Float64 (Frame'Last)) with Ghost; function Invariant (N : Frame; Speed : Float64) return Boolean is (Speed in Low_Bound (N) .. High_Bound (N)) with Ghost; procedure Compute_Speed (N : Frame; Factor : Ratio_T; Drag : Drag_T; Old_Speed : Float64; New_Speed : out Float64) with Global => null, Pre => N < Frame'Last and then Invariant (N, Old_Speed), Post => Invariant (N + 1, New_Speed); end Simple_Trajectory;
File : simple_trajectory.adb
with Interfaces; use Interfaces; package body Simple_Trajectory with SPARK_Mode is procedure Compute_Speed (N : Frame; Factor : Ratio_T; Drag : Drag_T; Old_Speed : Float64; New_Speed : out Float64) is Delta_Speed : constant Float64 := Drag + Factor * G * Frame_Length; function T (N : Integer) return Float64 is (Float64 (N)) with Ghost; N_Bv : constant Unsigned_16 := Unsigned_16(N) with Ghost; begin New_Speed := Old_Speed + Delta_Speed; -- Bound all quantities involved with constants pragma Assert (Delta_Speed in -Bound .. Bound); pragma Assert (In_Bounds (High_Bound(N))); pragma Assert (In_Bounds (Low_Bound(N))); -- Intermediate assertions to bound New_Speed pragma Assert (Float64(N_Bv) * Bound + Bound = (Float64(N_Bv) + 1.0) * Bound); pragma Assert (Float64(N) * Bound + Bound = (Float64(N) + 1.0) * Bound); pragma Assert (Float64(N) * (-Bound) - Bound = (Float64(N) + 1.0) * (-Bound)); -- Relate the property on Float64(N) + 1.0 to the same property on Float64(N+1) pragma Assert (T(1) = 1.0); pragma Assert (Float64(N) + 1.0 = Float64(N + 1)); pragma Assert (New_Speed >= Float64 (N) * (-Bound) - Bound); pragma Assert (New_Speed >= Float64 (N + 1) * (-Bound)); pragma Assert (New_Speed <= Float64 (N) * Bound + Bound); pragma Assert (New_Speed <= Float64 (N + 1) * Bound); end Compute_Speed; end Simple_Trajectory;