sum_of_digits.mlw 2.94 KB
Newer Older
Jean-Christophe Filliatre's avatar
Jean-Christophe Filliatre committed
1

2 3
(**
   Projet Euler problem #290
Jean-Christophe Filliatre's avatar
Jean-Christophe Filliatre committed
4

5 6
     How many integers 0 <= n < 10^18 have the property that the sum
     of the digits of n equals the sum of digits of 137n?
Jean-Christophe Filliatre's avatar
Jean-Christophe Filliatre committed
7

8
   It can be easily computed using memoization (or, equivalently, dynamic
Jean-Christophe Filliatre's avatar
Jean-Christophe Filliatre committed
9 10
   programming) once the problem is generalized as follows:
   How many integers 0 <= n < 10^m are such that sd(n) = sd(137n + a) + b?
11

Jean-Christophe Filliatre's avatar
Jean-Christophe Filliatre committed
12 13
   In the following, we prove the correctness of a recursive function f
   which takes m, a, and b as arguments and returns the number of such n.
14 15
   Memoization must be added if one wants to solve the initial problem
   in a reasonable amount of time.
Jean-Christophe Filliatre's avatar
Jean-Christophe Filliatre committed
16
   (See for instance fib_memo.mlw for an example of memoized function.) *)
17

18 19
module Euler290

20
  use import int.Int
21
  use import ref.Ref
22
  use import int.EuclideanDivision as E
23
  use import int.Power
24
  use int.NumOf
25

Andrei Paskevich's avatar
Andrei Paskevich committed
26
  function sum_digits int : int
27

28
  axiom Sum_digits_def: forall n : int. sum_digits n =
29
    if n <= 0 then 0 else sum_digits (div n 10) + mod n 10
30

31 32
  (* the number of n st 0 <= n mod 10 < c and sd(n) = sd(137n+a)+b *)

33
  predicate p (a b c n: int) =
Andrei Paskevich's avatar
Andrei Paskevich committed
34
    0 <= mod n 10 < c /\ sum_digits n = sum_digits (137 * n + a) + b
35

36 37
  function numof (a b c: int) (l u: int) : int =
    NumOf.numof (\ n. p a b c n) l u
38

39
  function solution(a b m: int) : int = numof a b 10 0 (power 10 m)
40

41
  (* shortcut for the number of n st n mod 10 = c and ... *)
42

43 44
  function num_of_modc (a b c: int) (x y: int) : int =
    numof a b (c+1) x y - numof a b c x y
45 46

  (* helper lemmas *)
47

48
  lemma Base:
49
    forall a b: int. 0 <= a -> sum_digits a + b = 0 -> p a b 10 0
50 51

  lemma Empty:
52
    forall a b x y: int. numof a b 0 x y = 0
53

54
  lemma Induc:
55
    forall a b c: int,m:int. 0 <= a -> 0 <= c < 10 -> m > 0 ->
56 57 58
    let x = 137 * c + a in
    let a' = div x 10 in
    let b' = mod x 10 + b - c in
59
    solution a' b' (m-1) = num_of_modc a b c 0 (power 10 m)
60

61 62
  (* code *)

63
  use import int.ComputerDivision
64

65 66 67 68
  let rec sd (n: int) : int
    requires { n >= 0 }
    ensures  { result = sum_digits n }
    variant  { n }
69
  = if n = 0 then 0 else sd (div n 10) + mod n 10
Jean-Christophe Filliatre's avatar
Jean-Christophe Filliatre committed
70 71

  (* f(m,a,b) = the number of 0 <= n < 10^m such that
72 73
                digitsum(137n + a) + b = digitsum(n). *)
  let rec f (m a b: int) : int
74
    requires { 0 <= m /\ 0 <= a }
75 76
    ensures  { result = solution a b m }
    variant  { m }
77
  = if m = 0 then begin
Jean-Christophe Filliatre's avatar
Jean-Christophe Filliatre committed
78
      (* only n = 0 is possible *)
79 80
      assert { power 10 m = 1 };
      assert { numof a b 10 0 (power 10 m) = if sum_digits a + b = 0 then 1 else 0 };
Jean-Christophe Filliatre's avatar
Jean-Christophe Filliatre committed
81 82 83
      if sd a + b = 0 then 1 else 0
    end else begin
      let sum = ref 0 in
84
      let ghost p = power 10 m in
85
      for c = 0 to 9 do (* count the n st n mod 10 = c *)
86
        invariant { !sum = numof a b c 0 p }
87
        'L:
88
        let x = 137 * c + a in
Jean-Christophe Filliatre's avatar
Jean-Christophe Filliatre committed
89 90
        let q = div x 10 in
        let r = mod x 10 in
91 92 93
        let b' = (r+b-c) in
        sum := !sum + f (m-1) q b';
        assert { q = E.div x 10 && r = E.mod x 10 &&
94
          !sum - !(at sum 'L) = num_of_modc a b c 0 p }
Jean-Christophe Filliatre's avatar
Jean-Christophe Filliatre committed
95 96 97
      done;
      !sum
    end
98

99
end