100/bernoulli.ml

   1 (* ========================================================================= *)
   2 (* Bernoulli numbers and polynomials; sum of kth powers.                     *)
   3 (* ========================================================================= *)
   4
   5 needs "Library/binomial.ml";;
   6 needs "Library/analysis.ml";;
   7 needs "Library/transc.ml";;
   8
   9 prioritize_real();;
  10
  11 (* ------------------------------------------------------------------------- *)
  12 (* A couple of basic lemmas about new-style sums.                            *)
  13 (* ------------------------------------------------------------------------- *)
  14
  15 let SUM_DIFFS = prove
  16  (`!a m n. m <= n + 1 ==> sum(m..n) (\i. a(i + 1) - a(i)) = a(n + 1) - a(m)`,
  17   GEN_TAC THEN GEN_TAC THEN INDUCT_TAC THEN
  18   REWRITE_TAC[SUM_CLAUSES_NUMSEG] THENL
  19    [REWRITE_TAC[ARITH_RULE `m <= 0 + 1 <=> m = 0 \/ m = 1`] THEN
  20     STRIP_TAC THEN ASM_REWRITE_TAC[ARITH; ADD_CLAUSES; REAL_SUB_REFL];
  21     SIMP_TAC[ARITH_RULE `m <= SUC n + 1 <=> m <= n + 1 \/ m = SUC n + 1`] THEN
  22     STRIP_TAC THEN ASM_SIMP_TAC[ADD1] THENL [REAL_ARITH_TAC; ALL_TAC] THEN
  23     REWRITE_TAC[REAL_SUB_REFL; ARITH_RULE `~((n + 1) + 1 <= n + 1)`] THEN
  24     MATCH_MP_TAC SUM_TRIV_NUMSEG THEN ARITH_TAC]);;
  25
  26 let DIFF_SUM = prove
  27  (`!f f' a b.
  28         (!k. a <= k /\ k <= b ==> ((\x. f x k) diffl f'(k)) x)
  29         ==> ((\x. sum(a..b) (f x)) diffl (sum(a..b) f')) x`,
  30   REPLICATE_TAC 3 GEN_TAC THEN INDUCT_TAC THEN
  31   REWRITE_TAC[SUM_CLAUSES_NUMSEG] THEN COND_CASES_TAC THEN
  32   ASM_SIMP_TAC[ARITH; DIFF_CONST; SUM_TRIV_NUMSEG;
  33                ARITH_RULE `~(a <= SUC b) ==> b < a`] THEN
  34   DISCH_TAC THEN MATCH_MP_TAC DIFF_ADD THEN
  35   ASM_SIMP_TAC[LE_REFL; ARITH_RULE `k <= b ==> k <= SUC b`]);;
  36
  37 (* ------------------------------------------------------------------------- *)
  38 (* Bernoulli numbers.                                                        *)
  39 (* ------------------------------------------------------------------------- *)
  40
  41 let bernoulli = define
  42  `(bernoulli 0 = &1) /\
  43   (!n. bernoulli(SUC n) =
  44        --sum(0..n) (\j. &(binom(n + 2,j)) * bernoulli j) / (&n + &2))`;;
  45
  46 (* ------------------------------------------------------------------------- *)
  47 (* A slightly tidier-looking form of the recurrence.                         *)
  48 (* ------------------------------------------------------------------------- *)
  49
  50 let BERNOULLI = prove
  51  (`!n. sum(0..n) (\j. &(binom(n + 1,j)) * bernoulli j) =
  52        if n = 0 then &1 else &0`,
  53   INDUCT_TAC THEN
  54   REWRITE_TAC[bernoulli; SUM_CLAUSES_NUMSEG; GSYM ADD1; ADD_CLAUSES; binom;
  55               REAL_MUL_LID; LE_0; NOT_SUC] THEN
  56   SIMP_TAC[BINOM_LT; ARITH_RULE `n < SUC n`; BINOM_REFL; REAL_ADD_LID] THEN
  57   REWRITE_TAC[ADD_CLAUSES] THEN REWRITE_TAC[GSYM REAL_OF_NUM_ADD] THEN
  58   REWRITE_TAC[ARITH_RULE `SUC(SUC n) = n + 2`] THEN
  59   MATCH_MP_TAC(REAL_FIELD `x = &n + &2 ==> s + x * --s / (&n + &2) = &0`) THEN
  60   REWRITE_TAC[ADD1; BINOM_TOP_STEP_REAL; ARITH_RULE `~(n = n + 1)`] THEN
  61   REWRITE_TAC[BINOM_REFL] THEN REAL_ARITH_TAC);;
  62
  63 (* ------------------------------------------------------------------------- *)
  64 (* Bernoulli polynomials.                                                    *)
  65 (* ------------------------------------------------------------------------- *)
  66
  67 let bernpoly = new_definition
  68  `bernpoly n x = sum(0..n) (\k. &(binom(n,k)) * bernoulli k * x pow (n - k))`;;
  69
  70 (* ------------------------------------------------------------------------- *)
  71 (* The key derivative recurrence.                                            *)
  72 (* ------------------------------------------------------------------------- *)
  73
  74 let DIFF_BERNPOLY = prove
  75  (`!n x. ((bernpoly (SUC n)) diffl (&(SUC n) * bernpoly n x)) x`,
  76   REPEAT GEN_TAC THEN
  77   GEN_REWRITE_TAC (RATOR_CONV o LAND_CONV) [GSYM ETA_AX] THEN
  78   REWRITE_TAC[bernpoly; SUM_CLAUSES_NUMSEG; LE_0] THEN
  79   GEN_REWRITE_TAC LAND_CONV [GSYM REAL_ADD_RID] THEN
  80   MATCH_MP_TAC DIFF_ADD THEN REWRITE_TAC[SUB_REFL; real_pow; DIFF_CONST] THEN
  81   REWRITE_TAC[GSYM SUM_LMUL] THEN MATCH_MP_TAC DIFF_SUM THEN
  82   REPEAT STRIP_TAC THEN REWRITE_TAC[ADD1; BINOM_TOP_STEP_REAL] THEN
  83   DIFF_TAC THEN ASM_SIMP_TAC[ARITH_RULE `k <= n ==> ~(k = n + 1)`] THEN
  84   REWRITE_TAC[REAL_MUL_LZERO; REAL_ADD_LID] THEN
  85   ASM_SIMP_TAC[ARITH_RULE `k <= n ==> (n + 1) - k - 1 = n - k`] THEN
  86   ASM_SIMP_TAC[GSYM REAL_OF_NUM_SUB; ARITH_RULE `k <= n ==> k <= n + 1`] THEN
  87   UNDISCH_TAC `k <= n:num` THEN
  88   REWRITE_TAC[GSYM REAL_OF_NUM_ADD; GSYM REAL_OF_NUM_LE] THEN
  89   ABBREV_TAC `z = x pow (n - k)` THEN CONV_TAC REAL_FIELD);;
  90
  91 (* ------------------------------------------------------------------------- *)
  92 (* Hence the key stepping recurrence.                                        *)
  93 (* ------------------------------------------------------------------------- *)
  94
  95 let INTEGRALS_EQ = prove
  96  (`!f g. (!x. ((\x. f(x) - g(x)) diffl &0) x) /\ f(&0) = g(&0)
  97          ==> !x. f(x) = g(x)`,
  98   REPEAT STRIP_TAC THEN
  99   MP_TAC(SPECL [`\x:real. f(x) - g(x)`; `x:real`; `&0`] DIFF_ISCONST_ALL) THEN
 100   ASM_REWRITE_TAC[] THEN REAL_ARITH_TAC);;
 101
 102 let RECURRENCE_BERNPOLY = prove
 103  (`!n x. bernpoly n (x + &1) - bernpoly n x = &n * x pow (n - 1)`,
 104   INDUCT_TAC THENL
 105    [REWRITE_TAC[bernpoly; SUM_SING_NUMSEG; REAL_SUB_REFL; SUB_REFL;
 106                 real_pow; REAL_MUL_LZERO];
 107     ALL_TAC] THEN
 108   MATCH_MP_TAC INTEGRALS_EQ THEN CONJ_TAC THENL
 109    [X_GEN_TAC `x:real` THEN FIRST_X_ASSUM(MP_TAC o SPEC `x:real`) THEN
 110     ONCE_REWRITE_TAC[GSYM REAL_SUB_0] THEN
 111     DISCH_THEN(MP_TAC o AP_TERM `(*) (&(SUC n))`) THEN
 112     REWRITE_TAC[REAL_MUL_RZERO] THEN DISCH_THEN(SUBST1_TAC o SYM) THEN
 113     REWRITE_TAC[REAL_SUB_LDISTRIB] THEN
 114     REPEAT(MATCH_MP_TAC DIFF_SUB THEN CONJ_TAC) THEN
 115     SIMP_TAC[SUC_SUB1; DIFF_CMUL; DIFF_POW; DIFF_BERNPOLY; ETA_AX] THEN
 116     GEN_REWRITE_TAC LAND_CONV [GSYM REAL_MUL_RID] THEN
 117     MATCH_MP_TAC DIFF_CHAIN THEN REWRITE_TAC[DIFF_BERNPOLY] THEN
 118     DIFF_TAC THEN REAL_ARITH_TAC;
 119     ALL_TAC] THEN
 120   REWRITE_TAC[bernpoly; GSYM SUM_SUB_NUMSEG] THEN
 121   REWRITE_TAC[REAL_ADD_LID; REAL_POW_ONE; GSYM REAL_SUB_LDISTRIB] THEN
 122   REWRITE_TAC[SUM_CLAUSES_NUMSEG; LE_0; SUB_REFL; real_pow] THEN
 123   REWRITE_TAC[REAL_SUB_REFL; REAL_MUL_RZERO; REAL_ADD_RID] THEN
 124   SIMP_TAC[ARITH_RULE `i <= n ==> SUC n - i = SUC(n - i)`] THEN
 125   REWRITE_TAC[real_pow; REAL_MUL_LZERO; REAL_SUB_RZERO; REAL_MUL_RID] THEN
 126   REWRITE_TAC[BERNOULLI; ADD1] THEN
 127   COND_CASES_TAC THEN ASM_REWRITE_TAC[ARITH; real_pow; REAL_MUL_LID] THEN
 128   CONV_TAC SYM_CONV THEN REWRITE_TAC[REAL_ENTIRE; REAL_POW_EQ_0] THEN
 129   ASM_REWRITE_TAC[ADD_SUB]);;
 130
 131 (* ------------------------------------------------------------------------- *)
 132 (* Hence we get the main result.                                             *)
 133 (* ------------------------------------------------------------------------- *)
 134
 135 let SUM_OF_POWERS = prove
 136  (`!n. sum(0..n) (\k. &k pow m) =
 137         (bernpoly(SUC m) (&n + &1) - bernpoly(SUC m) (&0)) / (&m + &1)`,
 138   GEN_TAC THEN ASM_SIMP_TAC[REAL_EQ_RDIV_EQ; REAL_ARITH `&0 < &n + &1`] THEN
 139   ONCE_REWRITE_TAC[GSYM REAL_MUL_SYM] THEN
 140   REWRITE_TAC[GSYM SUM_LMUL] THEN MATCH_MP_TAC EQ_TRANS THEN EXISTS_TAC
 141    `sum(0..n) (\i. bernpoly (SUC m) (&(i + 1)) - bernpoly (SUC m) (&i))` THEN
 142   CONJ_TAC THENL
 143    [REWRITE_TAC[RECURRENCE_BERNPOLY; GSYM REAL_OF_NUM_ADD] THEN
 144     REWRITE_TAC[GSYM REAL_OF_NUM_SUC; SUC_SUB1];
 145     SIMP_TAC[SUM_DIFFS; LE_0] THEN REWRITE_TAC[REAL_OF_NUM_ADD]]);;
 146
 147 (* ------------------------------------------------------------------------- *)
 148 (* Now explicit computations of the various terms on specific instances.     *)
 149 (* ------------------------------------------------------------------------- *)
 150
 151 let SUM_CONV =
 152   let pth = prove
 153    (`sum(0..0) f = f 0 /\ sum(0..SUC n) f = sum(0..n) f + f(SUC n)`,
 154     SIMP_TAC[SUM_CLAUSES_NUMSEG; LE_0]) in
 155   let econv_0 = GEN_REWRITE_CONV I [CONJUNCT1 pth]
 156   and econv_1 = GEN_REWRITE_CONV I [CONJUNCT2 pth] in
 157   let rec sconv tm =
 158     (econv_0 ORELSEC
 159      (LAND_CONV(RAND_CONV num_CONV) THENC econv_1 THENC
 160       COMB2_CONV (RAND_CONV sconv) (RAND_CONV NUM_SUC_CONV))) tm in
 161   sconv;;
 162
 163 let BINOM_CONV =
 164   let pth = prove
 165    (`a * b * x = FACT c ==> x = (FACT c) DIV (a * b)`,
 166     REPEAT STRIP_TAC THEN CONV_TAC SYM_CONV THEN
 167     MATCH_MP_TAC DIV_UNIQ THEN EXISTS_TAC `0` THEN CONJ_TAC THENL
 168      [POP_ASSUM MP_TAC THEN ARITH_TAC;
 169       POP_ASSUM MP_TAC THEN ONCE_REWRITE_TAC[GSYM CONTRAPOS_THM] THEN
 170       SIMP_TAC[LT_NZ; MULT_ASSOC; MULT_CLAUSES] THEN
 171       MESON_TAC[LT_NZ; FACT_LT]]) in
 172   let match_pth = MATCH_MP pth
 173   and binom_tm = `binom` in
 174   fun tm ->
 175     let bop,lr = dest_comb tm in
 176     if bop <> binom_tm then failwith "BINOM_CONV" else
 177     let l,r = dest_pair lr in
 178     let n = dest_numeral l and k = dest_numeral r in
 179     if n </ k then
 180       let th = SPECL [l;r] BINOM_LT in
 181       MP th (EQT_ELIM(NUM_LT_CONV(lhand(concl th))))
 182     else
 183       let d = n -/ k in
 184       let th1 = match_pth(SPECL [mk_numeral d; r] BINOM_FACT) in
 185       CONV_RULE NUM_REDUCE_CONV th1;;
 186
 187 let BERNOULLIS =
 188   let th_0,th_1 = CONJ_PAIR bernoulli
 189   and b_tm = `bernoulli` in
 190   let conv_1 = GEN_REWRITE_CONV I [th_1] in
 191   let rec bconv n =
 192     if n <= 0 then [th_0] else
 193     let bths = bconv (n - 1)
 194     and tm = mk_comb(b_tm,mk_small_numeral n) in
 195     (RAND_CONV num_CONV THENC conv_1 THENC
 196      LAND_CONV(RAND_CONV SUM_CONV) THENC
 197      ONCE_DEPTH_CONV BETA_CONV THENC
 198      DEPTH_CONV(NUM_RED_CONV ORELSEC BINOM_CONV) THENC
 199      GEN_REWRITE_CONV ONCE_DEPTH_CONV bths THENC
 200      REAL_RAT_REDUCE_CONV) tm :: bths in
 201   bconv;;
 202
 203 let BERNOULLI_CONV =
 204   let b_tm = `bernoulli` in
 205   fun tm -> let op,n = dest_comb tm in
 206             if op <> b_tm or not(is_numeral n) then failwith "BERNOULLI_CONV"
 207             else hd(BERNOULLIS(dest_small_numeral n));;
 208
 209 let BERNPOLY_CONV =
 210   let conv_1 =
 211     REWR_CONV bernpoly THENC SUM_CONV THENC
 212     TOP_DEPTH_CONV BETA_CONV THENC NUM_REDUCE_CONV
 213   and conv_3 =
 214     ONCE_DEPTH_CONV BINOM_CONV THENC REAL_POLY_CONV in
 215   fun tm ->
 216     let n = dest_small_numeral(lhand tm) in
 217     let conv_2 = GEN_REWRITE_CONV ONCE_DEPTH_CONV (BERNOULLIS n) in
 218     (conv_1 THENC conv_2 THENC conv_3) tm;;
 219
 220 let SOP_CONV =
 221   let pth = prove
 222    (`sum(0..n) (\k. &k pow m) =
 223         (\p. (p(&n + &1) - p(&0)) / (&m + &1))
 224         (\x. bernpoly (SUC m) x)`,
 225     REWRITE_TAC[SUM_OF_POWERS]) in
 226   let conv_0 = REWR_CONV pth in
 227   REWR_CONV pth THENC
 228   RAND_CONV(ABS_CONV(LAND_CONV NUM_SUC_CONV THENC BERNPOLY_CONV)) THENC
 229   TOP_DEPTH_CONV BETA_CONV THENC
 230   REAL_POLY_CONV;;
 231
 232 let SOP_NUM_CONV =
 233   let pth = prove
 234    (`sum(0..n) (\k. &k pow p) = &m ==> nsum(0..n) (\k. k EXP p) = m`,
 235     REWRITE_TAC[REAL_OF_NUM_POW; GSYM REAL_OF_NUM_SUM_NUMSEG;
 236                 REAL_OF_NUM_EQ]) in
 237   let rule_1 = PART_MATCH (lhs o rand) pth in
 238   fun tm ->
 239     let th1 = rule_1 tm in
 240     let th2 = SOP_CONV(lhs(lhand(concl th1))) in
 241     MATCH_MP th1 th2;;
 242
 243 (* ------------------------------------------------------------------------- *)
 244 (* The example Bernoulli bragged about.                                      *)
 245 (* ------------------------------------------------------------------------- *)
 246
 247 time SOP_NUM_CONV `nsum(0..1000) (\k. k EXP 10)`;;
 248
 249 (* ------------------------------------------------------------------------- *)
 250 (* The general formulas for moderate powers.                                 *)
 251 (* ------------------------------------------------------------------------- *)
 252
 253 time SOP_CONV `sum(0..n) (\k. &k pow 0)`;;
 254 time SOP_CONV `sum(0..n) (\k. &k pow 1)`;;
 255 time SOP_CONV `sum(0..n) (\k. &k pow 2)`;;
 256 time SOP_CONV `sum(0..n) (\k. &k pow 3)`;;
 257 time SOP_CONV `sum(0..n) (\k. &k pow 4)`;;
 258 time SOP_CONV `sum(0..n) (\k. &k pow 5)`;;
 259 time SOP_CONV `sum(0..n) (\k. &k pow 6)`;;
 260 time SOP_CONV `sum(0..n) (\k. &k pow 7)`;;
 261 time SOP_CONV `sum(0..n) (\k. &k pow 8)`;;
 262 time SOP_CONV `sum(0..n) (\k. &k pow 9)`;;
 263 time SOP_CONV `sum(0..n) (\k. &k pow 10)`;;
 264 time SOP_CONV `sum(0..n) (\k. &k pow 11)`;;
 265 time SOP_CONV `sum(0..n) (\k. &k pow 12)`;;
 266 time SOP_CONV `sum(0..n) (\k. &k pow 13)`;;
 267 time SOP_CONV `sum(0..n) (\k. &k pow 14)`;;
 268 time SOP_CONV `sum(0..n) (\k. &k pow 15)`;;
 269 time SOP_CONV `sum(0..n) (\k. &k pow 16)`;;
 270 time SOP_CONV `sum(0..n) (\k. &k pow 17)`;;
 271 time SOP_CONV `sum(0..n) (\k. &k pow 18)`;;
 272 time SOP_CONV `sum(0..n) (\k. &k pow 19)`;;
 273 time SOP_CONV `sum(0..n) (\k. &k pow 20)`;;
 274 time SOP_CONV `sum(0..n) (\k. &k pow 21)`;;