Programming in Isabelle/HOL

.

ATTENTION: this tutorial is NOT updated to recent versions of Isabelle

.

Isabelle/jEdit is appropriate for programming in Standard ML as well as Isabelle/HOL's executable kernel is appropriate for programming plus verifying the programmed functions and generating efficient code.

These pages shall support students at RISC Linz to enter an advanced level of programming, in particular programming algorithms in Computer Algebra (CA). Presently these pages collect experiences from the feasibility study on implementation of the greatest common divisor (GCD) for multivariate polynomials. The experiences are recorded in a fairly ad-hoc manner and thus are likely to be improved in later revisions.

Programming in Standard ML (SML)
Standard ML is the implementation language of Isabelle and thus of ISAC. SML is acknowledged as one of the languages most appropriate for the development of theorem provers. Isabelle works best with PolyML, an implementation of SML featuring parallel execution on multi-core machines and featuring state-of-the-art prover IDEs, for instance Isabelle/jEdit. Actually, the latter is recommended as an Integrated Developmein Environment (IDE) for development in SML (without using Isabelle's prover itself).

Start programming in SML
Mathematicians commonly are not educated in functional programming, but experiences from other programming languages easily carry over to programming in SML. For instance, this is an algorithm for calculating the GCD of univariate polynomials (p.93 in ) described in a very common style:

Fig. TODO: copy of p.93 in Franz Winkler. Polynomial algorithms

The above algorithm can be translated in a straight forward manner to SML as follows, see lines 263..317 in the initial Version of GCD_Poly.thy:

01 fun GCD_MOD (a: unipoly) (b: unipoly) = 02  if a = [0] orelse b = [0] then [1] else 03    if a = b then a else 04      let 05 (*1*)  val d = abs (Integer.gcd (lc a) (lc b)); 06        val M = 2*d*landau_mignotte_bound a b; 07         fun GOTO2 a b d M p = (*==============================*) 08          let 09 (*2*)      val p = find_next_prime_not_divide d p 10             val c_p = normirt_polynom ( mod_gcd_p a b p) p 11             val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p 12             fun GOTO3 a b d M p g_p = (*........................*) 13 (*3*)        if (deg g_p) = 0 then [1] else 14                let 15                  val P = p 16                   val g = g_p 17                  fun WHILE a b d M P g p = (*--*) 18 (*4*)              if P > M then g else 19                      let 20                        val p = find_next_prime_not_divide d p 21                         val c_p = normiert_polynom ( mod_gcd_p a b p) p 22                         val g_p = poly_centr((c_p %* (d mod p)) poly_mod p) p 23                       in 24                         if deg g_p < deg g 25                         then GOTO3 a b d M p g_p 26                        else 27                          if deg g_p = deg g 28                           then 29                            let 30                              val g = CRA_2_poly (P,p)(g,g_p) 31                              val P = P*p 32                              val g = poly_centr(g poly_mod P)P 33                            in WHILE a b d M P g p end 34                          else WHILE a b d M P g p 35                       end (*--*) 36                  val g = WHILE a b d M P g p (* << 1.Mal -*) 37                in g end (*....................................*) 38            val g = GOTO3 a b d M p g_p (* << 1.Mal .........*) 39 (*5*)      val g = pp g 40           in 41             if (g %|% a) andalso (g %|% b) 42             then g 43             else GOTO2 a b d M p 44           end (*==============================================*) 45        val g = GOTO2 a b d M (*p*)1(* << 1. Mal =============*) 46      in g end;

In functional programming GOTOs and WHILE loops are replaced by function calls. Using these as above (marked by (*<<--*)) leads to a above one-to-one translation from the original description in Winkler's book to and SML program, where the labels for the GOTOs are marked by comments (*2*) and (*3*) respectively.

This simple kind of translation should not blind us to the fact that functional programming is an art of its own with particular strengths (for instance, getting parallel execution for free). There is a wealth of textbooks introducing to this art; Isabelle extends and completes the Standard ML Basis Library with a very elegant library.ML of its own.

This library is worth to study as a good starting point for programming in SML in particular and for functional programming in general. Using the functions in this library as much as possible is mandatory for the ISAC project.

Development of Computer Algebra (CA) in ISAC
CA shall be developed by use of Isabelle's function package and code generator as shown in the sequel. Since these Isabelle/HOL technologies are more advanced than SML and not yet as mature as SML, it might be advisable to do an intermediate step and implement in SML first. For instance, given the SML code

type monom = (int * int list); type poly = monom list; fun mult_monom ((c1, es1): monom) ((c2, es2): monom) = (c1 * c2, map2 Integer.add es1 es2): monom;

one can define multiplication of polynomials, for instance p1 = 2x + 3xy^2 and p2 = 2x - 3xy^2 which are encoded as

val p1 = [(2, [1,0]), (3, [1,2])]; val p2 = [(2, [1,0]), (-3, [1,2])];

and build up the code by stepwise transforming p1, p2 from values (as defined above) to arguments of functions, for instance starting by

map (mult_monom (hd p1)) p2; map ... fold ...

This style of programming, i.e. interactively building code bottom-up, is not supported by Isabelle/Isar, but very useful.

The first work in CA in ISAC at RISC has been doneby Diana Meindl implementing an algorithm for the greatest common divisor (GCD) of multivariate polynomials; this algorithm is required for cancellation of fractions and many other applications in CA. Her work is documented as guidelines for further work in her  thesis and on these web pages. The first step was translation from Winkler's mathematics textbook to SML as described above.

This step can be resumed in the ISAC repository at changeset 7aea961f5b50 (just copy&paste this hex-number into the input field top left in your browser); the above SML code is copied from GCD_Poly.thy. Because GCD_Poly.thy is part of the code under development, the original version translated from Winkler has been copied to a separate file gcd_poly_winkler.sml and kept unchanged in ISAC's tests. The respective files are as follows (where '' indicate the root directory of an Isabelle/ISAC distribution):

So, from changeset 7aea961f5b50 onwards changes can be traced in ISAC's repository from gcd_poly_winkler.sml via GCD_Poly.thy. In changeset 7aea961f5b50 the code in both, gcd_poly_winkler.sml and GCD_Poly.thy, were equal, and then GCD_Poly.thy started to change to the current state. The changes involved adaption of identifiers to Isabelle's coding standards and improvement of the whole code. This improvement allowed for simple translation from SML to Isabelle's function package -- a translation which started with changeset 1c906991244b "Winkler --> Isabelle", a process which is discussed in the subsequent sections.

Programming with Isabelle's function package (FP)
Programming in Isabelle/HOL is more advanced than programming in SML, because the former


 * enforces proof of termination for all functions (where most proofs are done automatically)
 * creates theorems from the function definition as a prerequisite to verify correctness
 * creates theorems required by the code generator for translating into efficient code.

This section records the first approach of translating from SML (in GCD_Poly.thy) to the FP (resulting in GCD_POly_FP.thy), so the section is likely to be changed alongside further experiences.

Getting started
Starting with 'hands-on' experience is well supported by tutorials included in the Isabelle distribution:

Defining Recursive Functions in Isabelle/HOL Programming and Proving in Isabelle/HOL: §1 and §2 for the first approach.

For interactive experiments you can copy code snippets into a separate theories, for instance Try_Functions.thy (don't forget to take the same identifier for the file and for the head in the file): theory Try_Functions imports Main begin end

You have to observe some details when copying von PDF to Isabelle/Isar. In the Isabelle distribution there are also theories with examples; relevant for our purpose are

Examples of function definitions Examples and regression tests for automated termination proofs

Interactive experiments are most efficient when copying code from these theories, and from GCD_Poly_FP.thy to separate files, for instance: theory Try_GCD_Poly_FP imports "/src/Tools/isac/Knowledge/GCD_Poly_FP" begin end In this separate file one can use the same identifiers for function definitions again (which are allowed only uniquely within one and the same theory). Since you are supposed to have Isabelle already installed, you find the above theories in the distribution, too ('' addresses the root of your installation): Examples of function definitions: /src/HOL/ex/Fundefs.thy Examples and regression tests for automated termination proofs: /src/HOL/ex/Termination.thy

Depending on the working style, background theory is helpful sooner or later. Sources of information beyond the above mentioned are cited at the point in these web pages, where they are required first.

From SML to Isabelle/HOL
The transition from programming in SML to programming in Isabelle/HOL is straightforward such that starting with trials is possible before looking to below.

From fun to definition .. function
The tutorial Defining Recursive Functions in Isabelle/HOL immediately introduces 'fun'. However, 'fun' is too powerful for simple cases; these are the means to define functions in increasing power (look up GCD_Poly_FP.thy for examples):


 * definition: without recursion. Actually, all definitions werde done by 'fun' in GCD-Poly_FP.thy down to the function LANDAU_MIGNOTTE_bound, a function which caused troubles with 'fun' and 'function'. Only afterwards 'definition' was used from the beginning of GCD-Poly_FP.thy.


 * primrec: for primitive recursion. Just try 'primrec', and if you get errors like ' Primrec definition error: constructor missing ', then go on to the next ...

Unfinished subgoals: (a, 1, <): 1. ⋀r m. 0 < m ⟹ r mod m ≠ 1 ⟹ False 2. ⋀r m rinv. rinv < m ⟹ r mod m ≠ 1 ⟹ ¦rinv + 1¦ < ¦rinv¦ : Could not find lexicographic termination order.
 * fun: for cases of recursion where termination can be proved automatically. Go on to 'function' if you get errors and an output like

function ... :  equations : by pat_completeness auto termination sorry
 * function: produces proof obligations which express completeness and compatibility of patterns. The combination of the methods 'pat_completeness' and 'auto' is used to solve these proof obligations. In the first approach we postpone termination proofs (by termination sorry) until Programming_in_Isabelle/HOL:

A surprising variety of function definitions can be automatically completed like that: for instance, the total of 42 functions defining the gcd of univariate polynomials in GCD_Poly_FP.thy are implemented by:

The table is significant insofar as the algorithms are typical for Computer Algebra and insofar as GCD_Poly_FP.thy strictly uses the weakest kind of definition possible ('weakest' with respect to the order given above).

From SML code to Isabelle/HOL code
Comparison between SML in /src/Tools/isac/Knowledge/GCD_Poly.thy and Isabelle in /src/Tools/isac/Knowledge/GCD_Poly_FP.thy shows how to write function bodies. Actually, it is just 'copy and paste' from SML, in general.

The similarities between SML and Isabelle/HOL are in the focus of Isabelle development. List functions are used most frequently, so for instance, compare the '(*lists*)' section in Pure/library.ML with HOL/List.thy --- there are almost no differences. But in HOL/List.thy there are lots of helpful lemmas about these functions. It is worth to regularly watch the extensions and improvements during development in these files and others.

Take the most appropriate types
Isabelle has a more elaborated type system than ML. For instance the ML function

ML {* fun chinese_remainder (r1, m1) (r2, m2) = (r1 mod m1) + (((r2 - (r1 mod m1)) * (mod_inv m1 m2)) mod m2) * m1; chinese_remainder (17, 9) (3, 4) *} has the signature val chinese_remainder = fn: int * int -> int * int -> int val it = 35: int

and is implemented in Isabelle more concisely with a different signature as

definition chinese_remainder :: "int ⇒ nat ⇒ int ⇒ nat ⇒ nat" where "chinese_remainder r1 m1 r2 m2 =     ((nat (r1 mod (int m1))) +        (nat (((r2 - (r1 mod (int m1))) * (int (mod_inv m1 m2))) mod (int m2))) * m1)"

Later on, in the termination proofs, the type nat is more convenient than the type int. However, the type conversions require some experience to handle. During the process of translation from ML to Isabelle the following trials were helpful:

declare\[\[show_types\]\] value "nat (17 mod (int 9))" --" = 8" value "      (3 - (17 mod (int 9)))" --" = -5::int" value "      (3 - (17 mod (int 9)))" --" = -5::int" value "                               (int (mod_inv 9 4))"          --" = 1::int" value "     ((3 - (17 mod (int 9))) * (int (mod_inv 9 4)))"           --" = -5::int" value "(    ((3 - (17 mod (int 9))) * (int (mod_inv 9 4))) mod 4)"      --" = 3::int" value "(    ((3 - (17 mod (int 9))) * (int (mod_inv 9 4))) mod 4)  * 9"   --" = 27::int" value "(nat (((3 - (17 mod (int 9))) * (int (mod_inv 9 4))) mod 4)) * 9"  --" = 27::nat" (*                        |-|          |---|                    |--|    |-|               ||              ||             ||        |--|*)

And this is the way to see details within an already defined ML function

ML {* fun chinese_remainder (r1, m1) (r2, m2) = (writeln ("chinese_remainder" ^ ", r1 = " ^ PolyML.makestring r1 ^ ", m1 = " ^ PolyML.makestring m1 ^ ", r2 = " ^ PolyML.makestring r2 ^ ", m2 = " ^ PolyML.makestring m2 ^ ",      r1 mod m1 = " ^ PolyML.makestring (r1 mod m1) ^ ",      r2 - (r1 mod m1) mod m1 = " ^ PolyML.makestring (r2 - (r1 mod m1)) ^ ",      (mod_inv m1 m2) = " ^ PolyML.makestring (mod_inv m1 m2) ^ ",      ((r2 - (r1 mod m1)) * (mod_inv m1 m2)) = " ^ PolyML.makestring (((r2 - (r1 mod m1)) * (mod_inv m1 m2))) ^ ",      ((...) mod m2) * m1 = " ^ PolyML.makestring (((r2 - (r1 mod m1)) * (mod_inv m1 m2) mod m2) * m1) ^ ",      res = " ^ PolyML.makestring ((r1 mod m1) +      ((r2 - (r1 mod m1)) * (mod_inv m1 m2) mod m2) * m1) );   (r1 mod m1) +      (((r2 - (r1 mod m1)) * (mod_inv m1 m2)) mod m2) * m1 *}

while this is a convenient way to build up a function in ML:

val (r1, m1, r2, m2) = (17, 9, 3, 4); r1 mod (int m1));       (r2 - (r1 mod (int m1)));        (r2 - (r1 mod (int m1)));                             (mod_inv m1 m2);       ((r2 - (r1 mod m1)) * (mod_inv m1 m2)); (     ((r2 - (r1 mod m1)) * (mod_inv m1 m2)) mod m2); (     ((r2 - (r1 mod m1)) * (mod_inv m1 m2)) mod m2)  * m1;

The function package (FP) at work
Going a 'hands-on' approach we investigate how the FP works on an example, the first function definition in GCD_Poly_FP.thy, which requires extra effort, 'make_primes'.

For your trials load GCD_Poly_FP.thy into Isabelle/jEdit and go to this function definition:

function make_primes :: "nat list ⇒ nat ⇒ nat ⇒ nat list" where "make_primes ps last_p n =    (if (List.nth ps (List.length ps - 1)) < n     then       if is_prime ps (last_p + 2)        then make_primes (ps @ [(last_p + 2)]) (last_p + 2) n       else make_primes  ps (last_p + 2) n     else ps)" by pat_completeness auto termination sorry

This function definition automatically creates a series of theorems required later -- just copy these lines into GCD_Poly_FP.thy, click on them and watch the output window:

thm "make_primes.simps" thm "make_primes.induct" thm "make_primes.cases" thm "make_primes.psimps" thm "make_primes.pinduct"

And this function definition is already executable (by rewriting!), which can be tried by

value "make_primes [2, 3] 3 4" value "make_primes [2, 3] 3 5" value "make_primes [2, 3] 3 6" :

where clicking on the respective line shows the result (the "value" of execution) in the output window. Then click on the function definition and observe the information in the output window:

proof (prove): step 0 goal (2 subgoals): 1. ⋀P x. (⋀ps last_p n. x = (ps, last_p, n) ⟹ P) ⟹ P 2. ⋀ps last_p n psa last_pa na. (ps, last_p, n) = (psa, last_pa, na) ⟹ (if ps ! (length ps - 1) < n        then if is_prime ps (last_p + 2) then make_primes_sumC (ps @ [last_p + 2], last_p + 2, n)              else make_primes_sumC (ps, last_p + 2, n)         else ps) = (if psa ! (length psa - 1) < na        then if is_prime psa (last_pa + 2) then make_primes_sumC (psa @ [last_pa + 2], last_pa + 2, na)              else make_primes_sumC (psa, last_pa + 2, na)         else psa)

In order to observe details the provers can be called differently

function make_primes :: "nat list ⇒ nat ⇒ nat ⇒ nat list" where "make_primes ps last_p n =    (if (List.nth ps (List.length ps - 1)) < n     then       if is_prime ps (last_p + 2)        then make_primes (ps @ [(last_p + 2)]) (last_p + 2) n       else make_primes  ps (last_p + 2) n     else ps)" '''apply pat_completeness '''apply auto '''done termination sorry

Now clicking on apply pat_completeness resolves the first goal (which encodes 'pattern completeness', see the footnote on p.9 of the tutorial).

Clicking on apply auto resolves the second goal, which is a trivial example for 'pattern compatibility' (trivial, because this function is defined only by one pattern on the left-hand-side of the defining equation). Nevertheless, auto runs into problems with this goal. termination sorry postpones the respective proof, which will be addressed in.

In order to inspect auto when running into problems in we need to declare the simplifier to trace its activity:

function make_primes :: "nat list ⇒ nat ⇒ nat ⇒ nat list" where "make_primes ps last_p n =    (if (List.nth ps (List.length ps - 1)) < n     then       if is_prime ps (last_p + 2)        then make_primes (ps @ [(last_p + 2)]) (last_p + 2) n       else make_primes  ps (last_p + 2) n     else ps)" apply pat_completeness '''using \[\[simp_trace_depth_limit = 99\]\] '''using \[\[simp_trace = true\]\] apply auto '''using \[\[simp_trace = false\]\] done termination sorry

Note that \[\[ above is enforced by the wiki, so please, remove the backslashes and only keep the brackets. The trace you get is not easily readable, a known problem which is under construction for improvement already. The tracing for the simplifier can also be declared outside function definitions and proof sequences (in 'theory mode') as follows declare \[\[simp_trace_depth_limit=99\]\] declare \[\[simp_trace=true\]\]

Proving 'pattern compatibility'
'Pattern completeness' was proved automatically for all functions so far. However, proving 'pattern compatibility' requires special handling of auto, which in turn calls the simplifier. Isabelle's simplifier attacks equality by descending into the function definitions. For CA this is not the appropriate strategy: mod, div or gcd cause error "linarith_split_limit exceeded (current value is 9)" and the mathematical algorithms tend to make the simplifier never come back in proofs for 'pattern compatibility'. So the simplifier needs to be adjusted to these kinds of proofs.

Below the 'hands-on approach' is continued. Since the examples are from a working version of GCD_Poly_FP.thy, the situation from the initial approach needs to be restored, in order to show the respective proofs not working, and then to demonstrate how to make them work.

Remove looping rules from the simplifier
The first function definition in GCD_Poly_FP.thy, which requires extra effort, is 'make_primes'. First we re-establish the situation of the first approach, when the troubles appeared -- we re-introduce harmful rewrite rules is_prime.simps (which are removed as shown in below) by "declare is_prime.simps [simp add]":

'''declare is_prime.simps [simp add] function make_primes :: "nat list ⇒ nat ⇒ nat ⇒ nat list" where "make_primes ps last_p n =    (if (List.nth ps (List.length ps - 1)) < n     then       if is_prime ps (last_p + 2)        then make_primes (ps @ [(last_p + 2)]) (last_p + 2) n       else make_primes  ps (last_p + 2) n     else ps)" apply pat_completeness '''using \[\[simp_trace_depth_limit = 99\]\] '''using \[\[simp_trace = true\]\] apply auto '''using \[\[simp_trace = false\]\] done termination sorry

Now Isabelle 'hangs' on apply auto and clicking on this line creates a huge amount of output, ending with the line : : Tracing paused. Stop, or continue with next 100, 1000, 10000 messages?

In order to get the information required, this trace is too short; so select "1000" three times (and not 10000, in Isabelle2013 this still freezes jEdit). We get the following trace (':' indicates (many) further lines inbetween):

[1]SIMPLIFIER INVOKED ON THE FOLLOWING TERM: ⋀ps last_p n psa last_pa na. (ps, last_p, n) = (psa, last_pa, na) ⟹ (if ps ! (length ps - 1) < n then if is_prime ps (last_p + 2) then make_primes_sumC (ps @ [last_p + 2], last_p + 2, n) else make_primes_sumC (ps, last_p + 2, n) else ps) = (if psa ! (length psa - 1) < na then if is_prime psa (last_pa + 2) then make_primes_sumC (psa @ [last_pa + 2], last_pa + 2, na) else make_primes_sumC (psa, last_pa + 2, na) else psa) [1]Applying instance of rewrite rule "Product_Type.prod.inject": (?a1, ?b1) = (?a'1, ?b'1) ≡ ?a1 = ?a'1 ∧ ?b1 = ?b'1 [1]Rewriting: : :  [3]SIMPLIFIER INVOKED ON THE FOLLOWING TERM: is_prime ps (last_p + 2) ≡ ?c1 [3]Applying instance of rewrite rule "??.unknown": ps ≡ psa [3]Rewriting: ps ≡ psa [3]Applying instance of rewrite rule "??.unknown": last_p ≡ last_pa [3]Rewriting: last_p ≡ last_pa ''[3]Applying instance of rewrite rule "??.unknown": ''is_prime ?primes1 ?n1 ≡ if 0 < length ?primes1 then if ?n1 mod hd ?primes1 = 0 then False else is_prime (tl ?primes1) ?n1 else True ''[3]Rewriting: ''is_prime psa (last_pa + 2) ≡ if 0 < length psa then if (last_pa + 2) mod hd psa = 0 then False else is_prime (tl psa) (last_pa + 2) else True ''[3]Applying congruence rule: ''0 < length psa ≡ ?c2 ⟹ ''if 0 < length psa then if (last_pa + 2) mod hd psa = 0 then False else is_prime (tl psa) (last_pa + 2) else True ≡ ''if ?c2 then if (last_pa + 2) mod hd psa = 0 then False else is_prime (tl psa) (last_pa + 2) else True : :

Now this tracing output requires closer investigation. In Isabelle2013 the output buffer still has no search tools; so copy the contents of the output window to an editor of your choice. Closer investigation reveals, that is_prime occurs again and again while expanding arguments; this is a selection:

[3]Rewriting: is_prime psa (last_pa + 2) ≡ if 0 < length psa then if (last_pa + 2) mod hd psa = 0 then False else is_prime (tl psa) (last_pa + 2) else True : [3]Rewriting: is_prime psa (last_pa + 2) ≡ if 0 < length psa then if (last_pa + 2) mod hd psa = 0 then False else is_prime (tl psa) (last_pa + 2) else True : [4]Rewriting: is_prime (tl psa) (last_pa + 2) ≡ if 0 < length (tl psa) then if (last_pa + 2) mod hd (tl psa) = 0 then False else is_prime (tl (tl psa)) (last_pa + 2) else True : [6]Rewriting: is_prime (tl psa) (last_pa + 2) ≡ if 0 < length (tl psa) then if (last_pa + 2) mod hd (tl psa) = 0 then False else is_prime (tl (tl psa)) (last_pa + 2) else True : [6]Rewriting: is_prime (tl (tl psa)) (last_pa + 2) ≡ if 0 < length (tl (tl psa)) then if (last_pa + 2) mod hd (tl (tl psa)) = 0 then False else is_prime (tl (tl (tl psa))) (last_pa + 2) else True : [7]Rewriting: is_prime (tl (tl (tl psa))) (Suc (Suc last_pa)) ≡ ... : :

This kind of descent into recursion is not required for proving pattern compatibility, and the solution is to prevent the descent by deleting the respective rule from the simplifier:

function make_primes :: "nat list ⇒ nat ⇒ nat ⇒ nat list" where "make_primes ps last_p n =    (if (List.nth ps (List.length ps - 1)) < n     then       if is_prime ps (last_p + 2)        then make_primes (ps @ [(last_p + 2)]) (last_p + 2) n       else make_primes  ps (last_p + 2) n     else ps)" apply pat_completeness apply (auto simp del: is_prime.simps) done termination sorry

Exclude certain rules from simplification
The following function seems to cause the same problem as above, but it requires a different solution. For the purpose of investigation we re-establish the situation of the first approach, when the troubles appeared -- we re-introduce the harmful rewrite rules:

declare is_prime.simps [simp add] declare make_primes.simps [simp add] declare primes_upto.simps [simp add]

These declarations put above the following function definition make by pat_completeness auto hang:

function next_prime_not_dvd :: "nat ⇒ nat ⇒ nat" (infixl "next'_prime'_not'_dvd" 70) where "p next_prime_not_dvd n =     (let       ps = primes_upto (p + 1);       nxt = List.nth ps (List.length ps - 1)     in       if n mod nxt ≠ 0       then nxt       else nxt next_prime_not_dvd n)" by pat_completeness auto termination sorry

This simple version creates the following lines in the obligation-prover output in analogy to the previous case:

[1]Applying instance of rewrite rule "??.unknown": primes_upto ?n1 ≡ if ?n1 < 3 then [2] else make_primes [2, 3] 3 ?n1 [1]Rewriting: primes_upto (pa + 1) ≡ if pa + 1 < 3 then [2] else make_primes [2, 3] 3 (pa + 1) : : [4]Applying instance of rewrite rule "??.unknown": make_primes ?ps1 ?last_p1 ?n1 ≡ if ?ps1 ! (length ?ps1 - 1) < ?n1 then if is_prime ?ps1 (?last_p1 + 2) then make_primes (?ps1 @ [?last_p1 + 2]) (?last_p1 + 2) ?n1 else make_primes ?ps1 (?last_p1 + 2) ?n1 else ?ps1 [4]Rewriting: make_primes [2, 3] 3 (pa + 1) ≡ if [2, 3] ! (length [2, 3] - 1) < pa + 1 then if is_prime [2, 3] (3 + 2) then make_primes ([2, 3] @ [3 + 2]) (3 + 2) (pa + 1) else make_primes [2, 3] (3 + 2) (pa + 1) else [2, 3] : [5]Applying instance of rewrite rule "??.unknown": make_primes ?ps1 ?last_p1 ?n1 ≡ if ?ps1 ! (length ?ps1 - 1) < ?n1 then if is_prime ?ps1 (?last_p1 + 2) then make_primes (?ps1 @ [?last_p1 + 2]) (?last_p1 + 2) ?n1 else make_primes ?ps1 (?last_p1 + 2) ?n1 else ?ps1 [5]Rewriting: make_primes [2, 3, 5, 7] 9 (pa + 1) ≡ if [2, 3, 5, 7] ! (length [2, 3, 5, 7] - 1) < pa + 1 then if is_prime [2, 3, 5, 7] (9 + 2) then make_primes ([2, 3, 5, 7] @ [9 + 2]) (9 + 2) (pa + 1) else make_primes [2, 3, 5, 7] (9 + 2) (pa + 1) else [2, 3, 5, 7] : :

Now following the recommendation for the previous case does not succeed:

:

However, this dynamic exclusion does not prevent auto to call these rules and to unfold them inappropriately:

[1]SIMPLIFIER INVOKED ON THE FOLLOWING TERM: ⋀p n pa na. (p, n) = (pa, na) ⟹ (let ps = primes_upto (p + 1); nxt = ps ! (length ps - 1) in if n mod nxt ≠ 0 then nxt else next_prime_not_dvd_sumC (nxt, n)) = (let ps = primes_upto (pa + 1); nxt = ps ! (length ps - 1) in if na mod nxt ≠ 0 then nxt else next_prime_not_dvd_sumC (nxt, na)) [1]Applying instance of rewrite rule "Product_Type.prod.inject": (?a1, ?b1) = (?a'1, ?b'1) ≡ ?a1 = ?a'1 ∧ ?b1 = ?b'1 [1]Rewriting: : :  [1]Applying instance of rewrite rule "??.unknown": primes_upto ?n1 ≡ if ?n1 < 3 then [2] else make_primes [2, 3] 3 ?n1 [1]Rewriting: primes_upto (Suc n1) ≡ if Suc n1 < 3 then [2] else make_primes [2, 3] 3 (Suc n1) : : [4]Applying instance of rewrite rule "??.unknown": make_primes ?ps1 ?last_p1 ?n1 ≡ if ?ps1 ! (length ?ps1 - 1) < ?n1 then if is_prime ?ps1 (?last_p1 + 2) then make_primes (?ps1 @ [?last_p1 + 2]) (?last_p1 + 2) ?n1 else make_primes ?ps1 (?last_p1 + 2) ?n1 else ?ps1 [4]Rewriting: make_primes [2, 3] 3 (Suc n1) ≡ if [2, 3] ! (length [2, 3] - 1) < Suc n1 then if is_prime [2, 3] (3 + 2) then make_primes ([2, 3] @ [3 + 2]) (3 + 2) (Suc n1) else make_primes [2, 3] (3 + 2) (Suc n1) else [2, 3] : : [6]Applying instance of rewrite rule "??.unknown": is_prime ?ps1 ?n1 ≡ if 0 < length ?ps1 then if ?n1 mod hd ?ps1 = 0 then False else is_prime (tl ?ps1) ?n1 else True [6]Rewriting: is_prime [2, 3] 5 ≡ if 0 < length [2, 3] then if 5 mod hd [2, 3] = 0 then False else is_prime (tl [2, 3]) 5 else True : :

This case is solved by never including certain rule in any rule set: declare such exclusion after the respective function definition: fun is_prime :: "nat list ⇒ nat ⇒ bool" where : declare is_prime.simps [simp del] -- "make_primes, next_prime_not_dvd"

function make_primes :: "nat list ⇒ nat ⇒ nat ⇒ nat list" where : declare make_primes.simps [simp del] -- "next_prime_not_dvd"

fun primes_upto :: "nat ⇒ nat list" where : declare primes_upto.simps [simp del] -- "next_prime_not_dvd"

With these three declarations the function package resolves next_prime_not_dvd. The other ''declare primes_XXX.simps [simp del] serve the same purpose for other functions: The respective relations are noted in comments for the first occurrence. For instance,

by auto --"simp del: is_prime.simps, make_primes.simps, primes_upto.simps <-- declare"

means: The proof for this function requires the following declarations (which are located close to the respective function definitions):

declare is_prime.simps [simp add] declare make_primes.simps [simp add] declare primes_upto.simps [simp add]

Observe impossible patterns
In interesting example is the following function, because it allows to understand, what 'pattern compatibility' means (again, the leading declare re-establish the initial state causing troubles):

function dvd_up :: "unipoly ⇒ unipoly ⇒ bool" (infixl "dvd'_up" 70) where "[d] dvd_up [p] = ((¦d¦ ≤ ¦p¦) ∧ (p mod d = 0))" | "ds dvd_up ps =    (let        ds = drop_lc0_up ds; ps = drop_lc0_up ps;       d000 = (List.replicate (List.length ps - List.length ds) 0) @ ds;       quot = (lcoeff_up ps) div2 (lcoeff_up d000);       rest = drop_lc0_up (ps %-% (d000 %* quot))     in        if rest = [] then True else          if quot ≠ 0 & List.length ds ≤ List.length rest then ds %|% rest else False)" apply pat_completeness apply auto done termination sorry

At apply auto there are the following proof obligations for 'pattern compatibility'

proof (prove): step 1 goal (3 subgoals): 1. ⋀d p da pa. ([d], [p]) = ([da], [pa]) ⟹ (¦d¦ ≤ ¦p¦ ∧ p mod d = 0) = (¦da¦ ≤ ¦pa¦ ∧ pa mod da = 0) 2. ⋀d p ds ps.       ([d], [p]) = (ds, ps) ⟹ (¦d¦ ≤ ¦p¦ ∧ p mod d = 0) = (let ds = drop_lc0_up ds; ps = drop_lc0_up ps; d000 = replicate (length ps - length ds) 0 @ ds;        ...  3. ⋀ds ps dsa psa.        (ds, ps) = (dsa, psa) ⟹        (let ds = drop_lc0_up ds; ps = drop_lc0_up ps; d000 = replicate (length ps - length ds) 0 @ ds; ... =       (let ds = drop_lc0_up dsa; ps = drop_lc0_up psa; d000 = replicate (length ps - length ds) 0 @ ds;         ...

The most instructive subgoal is no.2. which we shall call 'mixed case' and which means: calling the two different definitions with the same arguments [d] and [p], then the two definitions must be the same, i.e.

"[d] dvd_up [p] =   (let       ds = drop_lc0_up [d]; ps = drop_lc0_up [p];      d000 = (List.replicate (List.length ps - List.length ds) 0) @ ds;      quot = (lcoeff_up ps) div2 (lcoeff_up d000);      rest = drop_lc0_up (ps %-% (d000 %* quot))    in       if rest = [] then True else         if quot ≠ 0 & List.length ds ≤ List.length rest then ds %|% rest else False)"

with the replacement ([d], [p]) = (ds, ps) must be the same as (¦d¦ ≤ ¦p¦ ∧ p mod d = 0). Actually, it can be shown, that the recursion in the above definition collapses to [d] dvd_up [p - d] -- which agrees with an essential property of mod. However, the proof obligation encodes the recursive call as dvd_up_sumC (ds, rest), i.e. dvd_up_sumC ([d], [p - d]), the external interface to the function package knows nothing about dvd_up_sumC and thus

dvd_up_sumC ([d], [p - d]) = (¦d¦ ≤ ¦p¦ ∧ p mod d = 0)

cannot be proved. The method to tackle cases like this one is (sequential):

function (sequential) dvd_up :: "unipoly ⇒ unipoly ⇒ bool" (infixl "dvd'_up" 70) where "[d] dvd_up [p] = ((¦d¦ ≤ ¦p¦) ∧ (p mod d = 0))" | "ds dvd_up ps =    (let        ds = drop_lc0_up ds; ps = drop_lc0_up ps;       d000 = (List.replicate (List.length ps - List.length ds) 0) @ ds;       quot = (lcoeff_up ps) div2 (lcoeff_up d000);       rest = drop_lc0_up (ps %-% (d000 %* quot))     in        if rest = [] then True else          if quot ≠ 0 & List.length ds ≤ List.length rest then ds %|% rest else False)" apply pat_completeness :

(sequential) automatically creates all different possibilities for patterns (see p.3 of the tutorial), for this function there are 15 possibilities, which are tackled one by one:

apply simp apply simp defer (*a "mixed" obligation*) apply simp defer (*a "mixed" obligation*) apply simp defer (*a "mixed" obligation*) apply simp apply simp apply simp (* > 1 sec *) apply simp apply simp (* > 1 sec *) apply simp defer (*a "mixed" obligation*) apply simp (* > 1 sec *) defer (*a "mixed" obligation*) apply (rule ex_falso_1) apply simp apply (rule ex_falso_2) apply simp apply (rule ex_falso_3) apply simp apply (rule ex_falso_1) apply simp done (* > 1 sec IMPROVED BY declare simp del: centr_up_def normalise.simps mod_up_gcd.simps lcoeff_up.simps*)

The defer shift the 'mixed' subgoals to the tail, which finally is:

proof (prove): step 17 goal (4 subgoals): 1. ⋀d p ds v vb vc. ([d], [p]) = (ds, v # vb # vc) ⟹ ... 2. ⋀ps v vb vc psa. ([], ps) = (v # vb # vc, psa) ⟹ ... 3. ⋀ds dsa v vb vc. (ds, []) = (dsa, v # vb # vc) ⟹ ... 4. ⋀d p v vb vc ps. ([d], [p]) = (v # vb # vc, ps) ⟹ ...

A closer look shows, that all the assumptions of the 'mixed' subgoals are false; for instance, in the 1.subgoal the assumption is false, because [p] = v # vb # vc is impossible (# is the list constructor). So there are lemmas like

lemma ex_falso_1: "([d], [p]) = (v # vb # vc, ps) ⟹ QUODLIBET" by simp

which are used at the bottom of the calls to the provers:

:

Using Isabelle's code generator (CG)
"Programming in Isabelle/HOL" as described in the previous section has an appealing advantage over programming in SML: The function package translates the (functional and "functional-logic" ) "programs" into higher-order rewrite systems, and from this intermediate language automated code generation into SML, Scala and Haskell is possible.

The ISAC project focuses generation of SML code. As soon as and  are accomplished, this code can be integrated into the Isabelle distribution.

But already before these proof tasks are accomplished, code generation is possible.

TODO
TODO

TODO
TODO

TODO
TODO

Prove termination
As mentioned above, proving termination can be separated from implementing functions (by 'definition', 'primrec', 'fun' or 'function'), because executable code (executable by rewriting in 'value "..."') is generated by the function package without proof of termination. Only 'function' requires an explicit proof.

It has also been mentioned, that the function package provides a high degree of automation: The 42 functions defining the univariate case in GCD_Poly_FP.thy require only 11 times 'function', all others state termination automatically. Among the 11 'function' definitions "termination by lexicographic_order" causes

So, the default lexicographic order is not very successful in these 11 visible cases (there are 11 invisible successes by 'fun', however); probably this is due to specifics of Computer Algebra. Interestingly, the only automatically successful "termination by lexicographic_order" is established for the final function gcd_up which calls all the other functions and thus could be considered the most complicated case.

TODO
TODO

function modiTEST :: "int ⇒ int ⇒ nat ⇒ nat ⇒ nat" where "modiTEST r rold m rinv =     (if m ≤ rinv       then 0 else         if r mod (int m) = 1          then rinv          else modiTEST (rold * ((int rinv) + 1)) rold m (rinv + 1))" by pat_completeness auto termination (*proof 1. wf ?R  2. ⋀r rold m rinv.        ¬ m ≤ rinv ⟹        r mod int m ≠ 1 ⟹ ((rold * (int rinv + 1), rold, m, rinv + 1), r, rold, m, rinv) ∈ ?R *) apply (relation "measure (λ(r, rold, m, rinv). m - rinv)") (* 1. wf (measure (λ(r, rold, m, rinv). m - rinv))  2. ⋀r rold m rinv.        ¬ m ≤ rinv ⟹        r mod int m ≠ 1 ⟹        ((rold * (int rinv + 1), rold, m, rinv + 1), r, rold, m, rinv)        ∈ measure (λ(r, rold, m, rinv). m - rinv) *) by auto

TODO
TODO

see how elements are addressed in tuples, for instance in a triple: value "fst (a,b,c)"      --"a" value "snd (a,b,c)" value "fst (snd (a,b,c))" --"b" value "snd (snd (a,b,c))" --"c"

Now one can read the Measures in the error message: Measures: 1) λp. size (snd (snd p))  2) λp. size (fst (snd p)) 3) λp. size (snd p)  4) λp. list_size size (fst p)  5) λp. length (fst p)   6) size

p is the triple of arguments (ps, last_p, n) and Measure 2) λp. size (fst (snd p)) concerns last_p, the measure with the largest number of ? in the error message: Result matrix:    1  2  3  4  5  6  a:  <= ?  <= ?  ?  <= b:  <= ?  <= <= <= <=

TODO

function make_primes :: "nat list ⇒ nat ⇒ nat ⇒ nat list" where "make_primes ps last_p n =    (if (List.nth ps (List.length ps - 1)) < n     then       if is_prime ps (last_p + 2)       then make_primes (ps @ [(last_p + 2)]) (last_p + 2) n       else make_primes  ps (last_p + 2) n     else ps)" by pat_completeness auto --"simp del: is_prime.simps <-- declare" termination (*Calls:  a) (ps, last_p, n) ~> (ps @ [last_p + 2], last_p + 2, n)   b) (ps, last_p, n) ~> (ps, last_p + 2, n) *) by (relation "measure (λ(ps, last_p, n). n - last_p)") auto

leaves the subgoals

All make_primes_dom 1. ⋀ps last_p n. ps ! (length ps - Suc 0) < n ⟹ is_prime ps (Suc (Suc last_p)) ⟹ n - Suc (Suc last_p) < n - last_p 2. ⋀ps last_p n. ps ! (length ps - Suc 0) < n ⟹ ¬ is_prime ps (Suc (Suc last_p)) ⟹ n - Suc (Suc last_p) < n - last_p

next trial

by (relation "measure (λ(ps, last_p, n). n + 2 - last_p)") auto

leaves these subgoals

All make_primes_dom 1. ⋀ps last_p n. ps ! (length ps - Suc 0) < n ⟹ is_prime ps (Suc (Suc last_p)) ⟹ n - last_p < Suc (Suc n) - last_p 2. ⋀ps last_p n. ps ! (length ps - Suc 0) < n ⟹ ¬ is_prime ps (Suc (Suc last_p)) ⟹ n - last_p < Suc (Suc n) - last_p

TODO
TODO

TODO
TODO

Verify the implemented functions
So far no deep knowledge about proving in Isabelle was required: Isabelle's function packages supports implementation of functions with a high degree of automation: about 3/4 is done automatically, the remaining 1/4 require explicit termination proofs, which in turn are supported by handy proof tools.

However, ...

TODO

locale gcd_poly = fixes gcd_poly :: int poly => int poly => int poly assumes gcd_poly a b = d    and not (a = [:0:]) and not (b = [:0:]) ==> d dvd_poly a ∧  d dvd_poly b  ∧ (∀ d'. (d' dvd_poly a ∧ d’ dvd_poly b) ==> d’ ≤ d)

TODO