Dafny
lijiahai
2022-03-25
目录

定义、证明、算法正确性的案例研究GCD

# 定义、证明、算法正确性的案例研究:GCD

K. Rustan M. Leino Manuscript KRML 279, 22 June 2021

摘要 本文的目的是展示一个程序开发的示例，介绍支持程序规范的定义，陈述和证明那些定义的引理，并使用引理来证明程序的正确性。以欧几里得的计算最大公约数的减法算法为例。

# 问题描述

让我们指定并验证一个算法来计算两个数的最大公约数(GCD)。在规范中，我们将引入一个函数，它的定义“显然是正确的”。我们不会使用这个函数来计算GCD，因为如果直接编译的话，“明显正确”的定义会给出非常低效的代码。相反，我们将使用Euclid算法来计算“明显正确”函数定义的值。我们将证明这个算法确实计算出了这个值。

本质上，我们有

function Gcd(x: pos, y: pos): pos

method EuclidGcd(x: pos, y: pos) returns (gcd: pos)
  ensures gcd == Gcd(x, y)
1
2
3
4

其中pos表示正整数的类型。

# 正整数

我们所做的一切都与正整数有关。Dafny为自然数(即非负整数)而非正整数构建了一种类型。我们可以在Dafny中使用子集类型来定义它们:

type pos = x | 1 <= x  // error: cannot find witness to show type is inhabited
1

Dafny想知道这种类型是否有居民，但它自己不知道。对于我们的示例来说，这无关紧要，但是我们确实需要处理我们得到的错误。为此，我们提供了witness条款:

type pos = x | 1 <= x witness 1
1

如果我们真的不愿意显示显示该类型为非空的witness，我们可以写witness *，这导致Dafny将pos类型视为可能为空。对于我们的例子，你可以选择其中一种，但因为提供一个实际的证人很容易，所以我们就这么做。

在续集中，当我指的是正整数时，我会说number。

# 因素

一个数的除数是它的因数。我们定义了一个谓词，它说明了数字p是数字x的因数是什么意思:

predicate IsFactor(p: pos, x: pos) {
  exists q :: p * q == x
}
1
2
3

换句话说，p是x的因子，如果存在一个被乘数q，那么x就是p * q的乘积。

为了讨论一个数字的所有因子，我们引入了一个函数factors，我们使用集合理解来定义它。一个简单的定义是:

function Factors(x: pos): set<pos> {
  set p: pos | IsFactor(p, x)  // error: set constructed must be finite
}
1
2
3

在Dafny中，set表示一个有限的集合(对于可能的无限集合，使用isset)。在这种情况下，Dafny并没有立即发现这个理解会生成一个有限集。幸运的是，我们可以简单地在理解中添加另一个连词，让Dafny看到集合是有限的:

function Factors(x: pos): set<pos> {
  set p: pos | p <= x && IsFactor(p, x)
}
1
2
3

在添加这个连词时，我们可能会犯一个错误，因为新的集合可能没有包含我们想要的所有因素。我们的结合点p <= x当然看起来很简单，但为什么不证明加上它不会意外遗漏任何因子呢?我们可以证明这个集合的元素与可能无限集相同:

lemma FactorsHasAllFactors(x: pos)
  ensures forall n :: n in Factors(x) <==> n in iset p: pos | IsFactor(p, x)
{
}
1
2
3
4

引理的证明是在引理主体中给出的(也就是说，在引理规范后面的一对花括号之间)。在这种情况下，证明是空的，因为Dafny自动地证明了引理，而不需要我们提供任何进一步的帮助。

在离开因子的定义之前，让我们说明和证明两个简单的引理。这些引理可以作为对我们定义的检查，它们在以后的开发中也会很有帮助。

lemma FactorsContains1(x: pos)
  ensures 1 in Factors(x)
{
  assert 1 * x == x;
}

lemma FactorsContainsSelf(x: pos)
  ensures x in Factors(x)
{
  assert x * 1 == x;
}
1
2
3
4
5
6
7
8
9
10
11

为了证明一个数字n(这里是1或x)在集合Factors(x)中，我们需要确定n满足集合理解的条件(在Factors(x)主体中)。连词n <= x被自动证明，但IsFactor(n, x)不是。根据IsFactor的定义，我们需要证明n * q == x的被乘数q的存在性。这种证明通常包括证明证人，这就是上面两个引理中的断言陈述所做的。由这些断言，验证者完成了引理的证明。

# 集合的最大值

为了讨论最大公约数，我们需要一个函数来挑选集合中最大的数。一种有点声明性的方法是使用such-that结构。特别地，对于集合s, let-such-that表达式

var x :| x in s && forall y :: y in s ==> y <= x;
x
1
2

表示将x绑定到满足条件x在s &&中的所有y:: y在s ==> y <= x的值，然后返回表达式x的值。这个条件说x在集合s中，并且在集合s中，x是最大的。

使用such-that结构有一个证明义务，即满足给定条件的值存在。如果我们要求s为非空，那么x in s条件很容易满足，但它需要更多的工作来说服验证者x的值满足量词。为此，我们将定义一个引理。我们将引理命名为MaxExists，然后我们可以这样写我们的函数Max:

function Max(s: set<pos>): pos
  requires s != {}
{
  MaxExists(s);
  var x :| x in s && forall y :: y in s ==> y <= x;
  x
}

lemma MaxExists(s: set<pos>)
  requires s != {}
  ensures exists x :: x in s && forall y :: y in s ==> y <= x
1
2
3
4
5
6
7
8
9
10
11

Dafny使用引理调用MaxExists(s)来建立后续表达式的格式良好性。顺便说一下，注意Max(以及引理MaxExists)有一个前置条件s !={}(关键字requires)。这意味着函数(以及引理)只能在非空集合中调用。

那么我们如何证明MaxExists呢?证明这样一个x存在的最直接的方法是计算一个满足所需性质的x。我们将引入另一个计算最大值的函数，称为FindMax，并在MaxExists引理的证明中使用它。函数FindMax将被递归实现。

lemma MaxExists(s: set<pos>)
  requires s != {}
  ensures exists x :: x in s && forall y :: y in s ==> y <= x
{
  var x := FindMax(s);
}
 
function FindMax(s: set<pos>): pos
  requires s != {}
  ensures max in s && forall y :: y in s ==> y <= FindMax(s)
1
2
3
4
5
6
7
8
9
10

我们现在不是在兜圈子吗?是的，在某些方面，我们让生活变得比必要的更困难。如果我们有FindMax，我们不需要Max，然后我们也不需要引理MaxExists。事实上，我们可以只编写和使用FindMax，而不引入Max或MaxExists。但在这个例子中，我希望主要的定义尽可能清晰，而不考虑如何计算。从这个意义上说，Max的主体比我们将要为FindMax编写的主体更具声明性。

以下是FindMax的完整定义:

function FindMax(s: set<pos>): (max: pos)
  requires s != {}
  ensures max in s && forall y :: y in s ==> y <= max
{
  var x :| x in s;
  if s == {x} then
    x
  else
    var s' := s - {x};
    assert s == s' + {x};
    var y := FindMax(s');
    if x < y then y else x
}
1
2
3
4
5
6
7
8
9
10
11
12
13

当函数的后置条件想要提到函数的结果值时，你可以使用函数本身，给出参数:FindMax(s)。我在上面第一次介绍FindMax时就这样做了。在完整的定义中，我展示了另一种方法，即为结果值引入一个名称:max。该名称只能在函数的后置条件中使用。很多时候，为结果引入这样的名称会导致更短的规范。

# GCD

有了我们定义的函数，现在就可以定义GCD了。取x的因子与y的因子，与之相交，得到它们的公因式，取其最大值:

function Gcd(x: pos, y: pos): pos {
  var common := Factors(x) * Factors(y);
  Max(common)  // error: common must be nonempty
}
1
2
3
4

对于这个简单的定义，验证者报告了一个前提条件的违反，因为它无法证明common满足Max的前提条件。我们知道公因数是非空的，因为我们知道1是任意两个数x和y的公因数。为了引起验证者的注意，我们写了一个断言:

function Gcd(x: pos, y: pos): pos {
  var common := Factors(x) * Factors(y);
  assert 1 in common;  // error: assertion violation
  Max(common)
}
1
2
3
4
5

唉，验证者不能证明这个断言。但我们可以看到，这一主张的存在足以消除先决条件的违反。现在我们来证明这个断言。这就是我们使用前面介绍的FactorsContains1引理的地方。对该引理的两次调用将证明断言，在程序文本中最好的捕获方法是将assert更改为assert by，并在by块中给出断言的证明:

function Gcd(x: pos, y: pos): pos {
  var common := Factors(x) * Factors(y);
  assert 1 in common by {
    FactorsContains1(x);
    FactorsContains1(y);
  }
  Max(common)
}
1
2
3
4
5
6
7
8

够了!我们现在已经给出了Gcd的一个格式良好的定义。

# GCD的性质

我们将证明我们的Gcd函数的三个属性——如果你愿意，可以称之为健全检查。(我们还需要第四个属性，稍后再介绍。)

作为第一个完整性检查，我们期望Gcd(x, y)返回一个同时是x和y因子的数字。此外，在所有同时是x和y因数的数字中，Gcd(x, y)的收益应该是最大的。

lemma AboutGcd(x: pos, y: pos)
  ensures IsFactor(Gcd(x, y), x)
  ensures IsFactor(Gcd(x, y), y)
  ensures forall p: pos :: IsFactor(p, x) && IsFactor(p, y) ==> p <= Gcd(x, y)
1
2
3
4

这个引理的前两个后置条件是自动证明的，但第三个不是。我们如何证明一个全称量词(也就是forall表达式)是成立的?我们使用达夫尼的forall语句。当用于证明时，forall命题对应于逻辑中的“普遍引入”规则。这个规则说的是"如果你想要证明对于所有的x:: P(x)，那么你所需要做的就是任意选择一个x，然后对那个x证明P(x)。

我们像这样引入forall语句:

  forall p: pos | IsFactor(p, x) && IsFactor(p, y)
    ensures p <= Gcd(x, y)
1
2

为了证明它，我们只需要提出这样一个事实，即p，既是x和y的因子，又是x和y的因子的交集。验证者就能够完成证明。

lemma AboutGcd(x: pos, y: pos)
  ensures IsFactor(Gcd(x, y), x)
  ensures IsFactor(Gcd(x, y), y)
  ensures forall p: pos :: IsFactor(p, x) && IsFactor(p, y) ==> p <= Gcd(x, y)
{
  forall p: pos | IsFactor(p, x) && IsFactor(p, y)
    ensures p <= Gcd(x, y)
  {
    assert p in Factors(x) * Factors(y);
  }
}
1
2
3
4
5
6
7
8
9
10
11

Dafny验证器经常需要这样的属性帮助。要证明它们，就把它们写成断言。换句话说，验证者知道集合交集的这个性质，但它没有足够的创意将这个性质引入证明中。通过断言该属性，我们要求验证者确认该属性(它能够这样做)，然后在其余的证明中使用该属性(在本例中，这就完成了证明)。

作为第二个完整性检查，我们证明了Gcd是对称的。

lemma GcdSymmetric(x: pos, y: pos)
  ensures Gcd(x, y) == Gcd(y, x)
{
  assert Factors(x) * Factors(y) == Factors(y) * Factors(x);
}
1
2
3
4
5

这个证明归结到集合交点是对称的这个事实，我们把它写成引理来引起验证者的注意。

作为第三个健全的检验，我们证明Gcd是幂等的。也就是说，如果你给它相同的参数两次，它会返回那个参数。

lemma GcdIdempotent(x: pos)
  ensures Gcd(x, x) == x
{
  FactorsContainsSelf(x);
  assert x in Factors(x) * Factors(x);
}
1
2
3
4
5
6

这个性质的证明可以归结为集合交点是幂等的，以及一个数是它自己的因子之一的性质。

# 欧几里德算法

欧几里得的求取两个数的GCD的减法算法是将两个数中的较大数反复减去较小的数，直到它们相等。每一个这样的减法都保留了GCD——一个我们需要证明的不变量——两个相等的数的GCD就是那个数——这个数是我们通过上面的“GCD幂等”引理建立的。

该算法具有循环不变量和幂等引理，其结果如下:

method EuclidGcd(X: pos, Y: pos) returns (gcd: pos)
  ensures gcd == Gcd(X, Y)
{
  var x, y := X, Y;
  while
    invariant Gcd(x, y) == Gcd(X, Y)  // error: invariant not maintained
    decreases x + y
  {
    case x < y =>
      y := y - x;
    case y < x =>
      x := x - y;
  }
  GcdIdempotent(x);
  return x;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

此方法使用while - case循环。(如果您熟悉Dijkstra的保护命令[1 (opens new window)]，这是do-od循环。)这个循环的每次迭代都会选择一个case来执行。所选的case必须是一个其守护条件求值为true的对象(如果几个case的守护条件求值为true，则循环在这些case之间任意选择)。如果没有这样的保护条件，则循环停止迭代。EuclidGcd中的循环当然可以是一个普通的while x != y循环，但while - case循环提供的两种情况的对称性使其美观。

除了一个循环不变式外，循环还声明了一个终止度量(关键字decreases)。证明循环终止归结为证明每次迭代使终止度量的值减小(按照Dafny内置的基于良好基础的整数顺序)。

上面的EuclidGcd方法没有验证，因为验证器无法证明每次迭代都保持循环不变。为此，我们需要我在上面提到的GCD的第四个性质:

lemma GcdSubtract(x: pos, y: pos)
  requires x < y
  ensures Gcd(x, y) == Gcd(x, y - x)
1
2
3

利用这个引理和GCD的对称性，我们可以完成GCD的证明:

method EuclidGcd(X: pos, Y: pos) returns (gcd: pos)
  ensures gcd == Gcd(X, Y)
{
  var x, y := X, Y;
  while
    invariant Gcd(x, y) == Gcd(X, Y)
    decreases x + y
  {
    case x < y =>
      GcdSubtract(x, y);
      y := y - x;
    case y < x =>
      calc {
        Gcd(x, y);
      ==  { GcdSymmetric(x, y); }
        Gcd(y, x);
      ==  { GcdSubtract(y, x); }
        Gcd(y, x - y);
      ==  { GcdSymmetric(y, x - y); }
        Gcd(x - y, y);
      }
      x := x - y;
  }
  GcdIdempotent(x);
  return x;
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

这个版本在循环的第一个分支中添加了对GcdSubtract的调用。在循环的第二个分支中，证明计算使用保持等式的步骤将表达式Gcd(x, y)转换为Gcd(x - y, y)。步骤中给出的提示求助于GcdSubtract和GcdSymmetric引理。

# GCD减性质

在这个案例中，GcdSubtract的证明比其他任何定义和引理都要复杂。

证明首先引入Gcd(x, y)的名称:

  var p := Gcd(x, y);
1

由Gcd的定义可知，p是x和y的因子，我们可以证明p也是y - x的因子:

  assert IsFactor(p, y - x) by {
    var a :| p * a == x;
    var b :| p * b == y;
    calc {
      y - x;
    ==
      p * b - p * a;
    ==
      p * (b - a);
    }
  }
1
2
3
4
5
6
7
8
9
10
11

为了证明IsFactor(p, y - x)，我们对IsFactor的定义告诉我们存在的被乘数引入a和b名称(因为p是x和y的因子)。用基本的算术步骤进行简单的计算，我们就可以得到p可以乘上另一个数(即b - A)得到y - x。

因为p同时是x和y - x的因子，所以它是x和y - x的公因式。我们写了两行代码来确保验证器使用了这个属性，用集合交集来表达:

  var common := Factors(x) * Factors(y - x);
  assert p in common;
1
2

最后，我们需要证明p是这个公因式的最大。我们使用forall语句来声明这个属性:

  forall q | q in common
    ensures q <= p
1
2

为了证明这个性质，我们填充了forall语句体。对于表示“公共性”集合中的任意数字q，我们分别将生成x和y - x乘积的被乘数命名为:

  {
    var a :| q * a == x;
    var b :| q * b == y - x;
1
2
3

使用简单的算术步骤，我们可以用一个证明计算来确定q也是y的因子:

    assert IsFactor(q, y) by {
      calc {
        y;
      ==
        x + (y - x);
      ==
        q * a + q * b;
      ==
        q * (a + b);
      }
    }
1
2
3
4
5
6
7
8
9
10
11

所以，因为q既是x又是y的因子，所以Gcd(x, y)的定义告诉我们q <= Gcd(x, y)通过给出另一个关于集合交集的提示:

    assert q in Factors(x) * Factors(y);
  }
1
2

验证者完成验证。

# 更多的对称

虽然我们现在已经有了GCD算法的完整证明，但您的美感可能会因为我们在两种情况下提供证明的方式的不对称性而受到影响。既然while - case循环为我们提供了这两种情况的对称表述，如果我们也能使这两种情况的证明更加相似就好了。

有几种方法可以改善这种情况。一种是重构第二个case的证明计算到它自己的引理中。然后，每个“案例”都有一行证明。

为了好玩，让我来描述另一个“技巧”，让这两种情况(不完全对称，但至少)更相似。诀窍在于使(已经不对称的)GcdSubtract引理也把参数转换为Gcd。我们将其改写为:

lemma GcdSubtract(x: pos, y: pos)
  requires x < y
  ensures Gcd(y, x) == Gcd(x, y - x)
{
  GcdSymmetric(x, y);
  // ... the proof continues as before
}
1
2
3
4
5
6
7

注意后置条件的左边现在是Gcd(y, x)，而不是Gcd(x, y)就像我们在这个引理的第一个版本中一样。这个证明所需要的唯一改变就是诉诸于Gcd的对称性，我们可以通过引理内部的一个引理来做到这一点。这给了我们一个重新表述的“GcdSubtract”引理的证明。

通过这种重新表述，我们可以简化EuclidGcd的第二个“情形”，以使第一个“情形”更加复杂为代价。本质上，我们把一个引理调用从第二种情况移到第一种情况，所以不是有1和3个引理调用在两种情况下，我们会有2和2。

    case x < y =>
      GcdSubtract(x, y);
      GcdSymmetric(y, x);
      y := y - x;
    case y < x =>
      GcdSymmetric(x - y, y);
      GcdSubtract(y, x);
      x := x - y;
1
2
3
4
5
6
7
8

它不是完全对称的，但也许你还是喜欢它。或者你可能会在另一种情况下记住这个技巧，当鞋子非常合适的时候。如果没有别的，您可以坚持使用我们在上面开发的第一个完整的证明。

# 主方法

如果证明本身不满足你，你仍然想看到算法的运行，你可以写一个Main方法，编译并运行程序。(在命令行中使用dafny工具的/compile:3选项是一种简单的方法。它将验证然后运行程序。)

下面是一个Main的例子:

method Main() {
  Test(15, 9);
  Test(14, 22);
  Test(371, 1);
  Test(1, 2);
  Test(1, 1);
  Test(13, 13);
  Test(60, 60);
}

method Test(x: pos, y: pos) {
  var gcd := EuclidGcd(x, y);
  print x, " gcd ", y, "  =  ", gcd, "\n";
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14

# 结论

这个案例研究展示了如何定义一个感兴趣的领域(这里是数字因子，导致了GCD的定义)，陈述和证明关于这些定义的一些引理，然后在一个小程序的证明中使用这些引理。

该程序，包括所有引理和其他与定义相关的证明义务，只需要少于3秒的Dafny验证器验证。您可以在Dafny测试套件[3 (opens new window)]中找到整个程序。

欧几里德的GCD算法是一个常见的教科书例子。它在不同的验证器中以不同的形式被证明。例如，TLA+教程以这个程序为例[0 (opens new window)]。它假设了我们在这里证明的GCD的数学性质。Why3程序库包含欧几里德GCD算法的一个版本，它在每一步中使用模而不是减法，这样可以减少迭代次数[2 (opens new window)]。

# 致谢

我感谢Reto Kramer提出这个问题作为一个有用的案例研究。

# 参考文献

[0]TLA proof system. https://tla.msr-inria.inria.fr/tlaps/content/Documentation/Tutorial/The_example.html (opens new window). 🔎 (opens new window)

[1]Edsger W. Dijkstra. A Discipline of Programming. Prentice Hall, Englewood Cliffs, NJ, 1976. 🔎 (opens new window)

[2]Jean-Christophe Filliâtre and Claude Marché. Greatest common divisor, using the euclidean algorithm. http://toccata.lri.fr/gallery/gcd.en.html (opens new window). 🔎 (opens new window)

[3]K. Rustan M. Leino. gcd.dfy. https://github.com/dafny-lang/dafny/blob/master/Test/dafny4/gcd.dfy (opens new window), June 2021. 🔎 (opens new window)

编辑 (opens new window)
上次更新: 2022/03/26, 14:38:56
自动调用引理
各种推导式

← 各种推导式→

最近更新
01
寻找最大和最小数
04-06
02
斐波那契数列
04-06
03
线性查询
04-06
更多文章>
Theme by Vdoing | Copyright © 2022-2022 Li Jiahai | Dafny Community | 2022