examples_in_progress/sum_of_digits.mlw

   1
   2 (**
   3    Projet Euler problem #290
   4
   5      How many integers 0 <= n < 10^18 have the property that the sum
   6      of the digits of n equals the sum of digits of 137n?
   7
   8    It can be easily computed using memoization (or, equivalently, dynamic
   9    programming) once the problem is generalized as follows:
  10    How many integers 0 <= n < 10^m are such that sd(n) = sd(137n + a) + b?
  11
  12    In the following, we prove the correctness of a recursive function f
  13    which takes m, a, and b as arguments and returns the number of such n.
  14    Memoization must be added if one wants to solve the initial problem
  15    in a reasonable amount of time.
  16    (See for instance fib_memo.mlw for an example of memoized function.) *)
  17
  18 module Euler290
  19
  20   use int.Int
  21   use ref.Ref
  22   use import int.EuclideanDivision as E
  23   use int.Power
  24   use int.NumOf
  25
  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   (* the number of n st 0 <= n mod 10 < c and sd(n) = sd(137n+a)+b *)
  32
  33   predicate p (a b c n: int) =
  34     0 <= mod n 10 < c /\ sum_digits n = sum_digits (137 * n + a) + b
  35
  36   function numof (a b c: int) (l u: int) : int =
  37     NumOf.numof (fun 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   function num_of_modc (a b c: int) (x y: int) : int =
  44     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     let x = 137 * c + a in
  57     let a' = div x 10 in
  58     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   (* code *)
  62
  63   use int.ComputerDivision
  64
  65   let rec sd (n: int) : int
  66     requires { n >= 0 }
  67     ensures  { result = sum_digits n }
  68     variant  { n }
  69   = if n = 0 then 0 else sd (div n 10) + mod n 10
  70
  71   (* f(m,a,b) = the number of 0 <= n < 10^m such that
  72                 digitsum(137n + a) + b = digitsum(n). *)
  73   let rec f (m a b: int) : int
  74     requires { 0 <= m /\ 0 <= a }
  75     ensures  { result = solution a b m }
  76     variant  { m }
  77   = if m = 0 then begin
  78       (* only n = 0 is possible *)
  79       assert { power 10 m = 1 };
  80       assert { numof a b 10 0 (power 10 m) = if sum_digits a + b = 0 then 1 else 0 };
  81       if sd a + b = 0 then 1 else 0
  82     end else begin
  83       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         label L in
  88         let x = 137 * c + a in
  89         let q = div x 10 in
  90         let r = mod x 10 in
  91         let b' = (r+b-c) in
  92         sum := !sum + f (m-1) q b';
  93         assert { q = E.div x 10 && r = E.mod x 10 &&
  94           !sum - !(sum at L) = num_of_modc a b c 0 p }
  95       done;
  96       !sum
  97     end
  98
  99 end