(** Projet Euler problem #290 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? It can be easily computed using memoization (or, equivalently, dynamic programming) once the problem is generalized as follows: How many integers 0 <= n < 10^m are such that sd(n) = sd(137n + a) + b? 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. Memoization must be added if one wants to solve the initial problem in a reasonable amount of time. (See for instance fib_memo.mlw for an example of memoized function.) *) module Euler290 use import int.Int use import ref.Ref use import int.EuclideanDivision as E use import int.Power use int.NumOf function sum_digits int : int axiom Sum_digits_def: forall n : int. sum_digits n = if n <= 0 then 0 else sum_digits (div n 10) + mod n 10 (* the number of n st 0 <= n mod 10 < c and sd(n) = sd(137n+a)+b *) predicate p (a b c n: int) = 0 <= mod n 10 < c /\ sum_digits n = sum_digits (137 * n + a) + b function numof (a b c: int) (l u: int) : int = NumOf.numof (\ n. p a b c n) l u function solution(a b m: int) : int = numof a b 10 0 (power 10 m) (* shortcut for the number of n st n mod 10 = c and ... *) 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 (* helper lemmas *) lemma Base: forall a b: int. 0 <= a -> sum_digits a + b = 0 -> p a b 10 0 lemma Empty: forall a b x y: int. numof a b 0 x y = 0 lemma Induc: forall a b c: int,m:int. 0 <= a -> 0 <= c < 10 -> m > 0 -> let x = 137 * c + a in let a' = div x 10 in let b' = mod x 10 + b - c in solution a' b' (m-1) = num_of_modc a b c 0 (power 10 m) (* code *) use import int.ComputerDivision let rec sd (n: int) : int requires { n >= 0 } ensures { result = sum_digits n } variant { n } = if n = 0 then 0 else sd (div n 10) + mod n 10 (* f(m,a,b) = the number of 0 <= n < 10^m such that digitsum(137n + a) + b = digitsum(n). *) let rec f (m a b: int) : int requires { 0 <= m /\ 0 <= a } ensures { result = solution a b m } variant { m } = if m = 0 then begin (* only n = 0 is possible *) assert { power 10 m = 1 }; assert { numof a b 10 0 (power 10 m) = if sum_digits a + b = 0 then 1 else 0 }; if sd a + b = 0 then 1 else 0 end else begin let sum = ref 0 in let ghost p = power 10 m in for c = 0 to 9 do (* count the n st n mod 10 = c *) invariant { !sum = numof a b c 0 p } 'L: let x = 137 * c + a in let q = div x 10 in let r = mod x 10 in let b' = (r+b-c) in sum := !sum + f (m-1) q b'; assert { q = E.div x 10 && r = E.mod x 10 && !sum - !(at sum 'L) = num_of_modc a b c 0 p } done; !sum end end