formal_ineqs/informal/informal_arith.hl

   1 (* =========================================================== *)
   2 (* Informal arithmetic procedures                              *)
   3 (* Author: Alexey Solovyev                                     *)
   4 (* Date: 2012-10-27                                            *)
   5 (* =========================================================== *)
   6
   7 (* Dependencies *)
   8 needs "misc/misc.hl";;
   9 needs "arith_options.hl";;
  10
  11
  12 (* Natural numbers *)
  13
  14 module type Informal_nat_sig =
  15   sig
  16     type nat
  17         val arith_base : int
  18     val mk_nat : num -> nat
  19     val mk_small_nat : int -> nat
  20     val dest_nat : nat -> num
  21     val eq_nat : nat -> nat -> bool
  22     val suc_nat : nat -> nat
  23     val pre_nat : nat -> nat
  24     val eq0_nat : nat -> bool
  25     val gt0_nat : nat -> bool
  26     val lt_nat : nat -> nat -> bool
  27     val le_nat : nat -> nat -> bool
  28     val add_nat : nat -> nat -> nat
  29     val sub_nat : nat -> nat -> nat
  30     (* If sub_and_le_nat m n = (m - n, true) if n <= m; (n - m, false) if m < n *)
  31     val sub_and_le_nat : nat -> nat -> nat * bool
  32     val mul_nat : nat -> nat -> nat
  33     val div_nat : nat -> nat -> nat
  34     val even_nat : nat -> bool
  35     val odd_nat : nat -> bool
  36
  37     (* normalize_nat m = (n, e) s.t. m = n * base^e, e >= 0 *)
  38     val normalize_nat : nat -> nat * int
  39     val denormalize_nat : nat * int -> nat
  40     (* hi_nat p m = (n, e) s.t. m <= n * base^e and n contains at most p "digits" *)
  41     val hi_nat : int -> nat -> nat * int
  42     val hi_lt_nat : int -> nat -> nat * int
  43     (* lo_nat p m = (n, e) s.t. n * base^e <= m and n contains at most p "digits" *)
  44     val lo_nat : int -> nat -> nat * int
  45   end;;
  46
  47
  48
  49 module Informal_nat : Informal_nat_sig = struct
  50
  51 open Arith_misc;;
  52 open Big_int;;
  53
  54 type nat = big_int;;
  55
  56 let arith_base = !Arith_options.base;;
  57
  58 let mk_nat n =
  59   let result = big_int_of_num n in
  60     if sign_big_int result < 0 then zero_big_int else result;;
  61
  62 let mk_small_nat n =
  63   if n < 0 then zero_big_int else big_int_of_int n;;
  64
  65 let dest_nat = num_of_big_int;;
  66
  67 let eq_nat = eq_big_int;;
  68
  69 let suc_nat = succ_big_int;;
  70
  71 let pre_nat n =
  72   let result = pred_big_int n in
  73     if sign_big_int result < 0 then zero_big_int else result;;
  74
  75 let eq0_nat n = sign_big_int n = 0;;
  76
  77 let gt0_nat n = sign_big_int n > 0;;
  78
  79 let lt_nat = lt_big_int;;
  80
  81 let le_nat = le_big_int;;
  82
  83 let add_nat = add_big_int;;
  84
  85 let sub_nat m n =
  86   let result = sub_big_int m n in
  87     if sign_big_int result < 0 then zero_big_int else result;;
  88
  89 let sub_and_le_nat m n =
  90   let result = sub_big_int m n in
  91     if sign_big_int result >= 0 then (result, true) else (abs_big_int result, false);;
  92
  93 let mul_nat = mult_big_int;;
  94
  95 let div_nat = div_big_int;;
  96
  97 let two_big_int = big_int_of_int 2;;
  98
  99 let even_nat n = sign_big_int (mod_big_int n two_big_int) = 0;;
 100
 101 let odd_nat n = sign_big_int (mod_big_int n two_big_int) > 0;;
 102
 103 (*******************************)
 104 (* num_exp *)
 105
 106 let base_nat = mk_small_nat arith_base;;
 107
 108 (* normalize_nat m = (n, e) s.t. m = n * base^e, e >= 0 *)
 109 let normalize_nat =
 110   let rec normalize n e =
 111     let q, r = quomod_big_int n base_nat in
 112       if sign_big_int r > 0 then
 113         (n, e)
 114       else
 115         normalize q (succ e) in
 116     fun n ->
 117       if sign_big_int n = 0 then (n, 0) else normalize n 0;;
 118
 119
 120 let denormalize_nat (n, e) =
 121   mult_big_int n (power_int_positive_int arith_base e);;
 122
 123
 124 let lo_nat pp =
 125   let max = power_int_positive_int arith_base pp in
 126   let rec lo m e =
 127     if lt_big_int m max then
 128       (m, e)
 129     else
 130       let q = div_big_int m base_nat in
 131         lo q (succ e) in
 132     fun m ->
 133       if sign_big_int m = 0 then
 134         (m, 0)
 135       else
 136         let n1, e1 = lo m 0 in
 137         let n, e2 = normalize_nat n1 in
 138           n, e1 + e2;;
 139
 140
 141 let hi_nat pp =
 142   if pp <= 0 then failwith "hi_nat: pp <= 0" else
 143     let max = power_int_positive_int arith_base pp in
 144     let rec hi m e =
 145       if lt_big_int m max then
 146         (m, e)
 147       else
 148         let q, r = quomod_big_int m base_nat in
 149           if sign_big_int r = 0 then
 150             hi q (succ e)
 151           else
 152             hi (succ_big_int q) (succ e) in
 153       fun m ->
 154         if sign_big_int m = 0 then
 155           (m, 0)
 156         else
 157           let n1, e1 = hi m 0 in
 158           let n, e2 = normalize_nat n1 in
 159             n, e1 + e2;;
 160
 161
 162 let hi_lt_nat pp m =
 163   hi_nat pp (succ_big_int m);;
 164
 165
 166 end;;
 167
 168
 169
 170 (* Floating point numbers *)
 171
 172 module type Informal_float_sig =
 173   sig
 174     type ifloat
 175         val min_exp : int
 176     val mk_float : num -> int -> ifloat
 177     val mk_num_float : num -> ifloat
 178     val mk_small_num_float : int -> ifloat
 179     val dest_float : ifloat -> bool * num * int
 180     val sign_float : ifloat -> bool
 181    (* Compares representations, not numbers themselves *)
 182     val eq_float : ifloat -> ifloat -> bool
 183     val lo_float : int -> ifloat -> ifloat
 184     val hi_float : int -> ifloat -> ifloat
 185     val neg_float : ifloat -> ifloat
 186     val abs_float : ifloat -> ifloat
 187     val lt0_float : ifloat -> bool
 188     val gt0_float : ifloat -> bool
 189     val le0_float : ifloat -> bool
 190     val ge0_float : ifloat -> bool
 191     val lt_float : ifloat -> ifloat -> bool
 192     val le_float : ifloat -> ifloat -> bool
 193     val min_float : ifloat -> ifloat -> ifloat
 194     val max_float : ifloat -> ifloat -> ifloat
 195     val mul_float_eq : ifloat -> ifloat -> ifloat
 196     val mul_float_lo : int -> ifloat -> ifloat -> ifloat
 197     val mul_float_hi : int -> ifloat -> ifloat -> ifloat
 198     val div_float_lo : int -> ifloat -> ifloat -> ifloat
 199     val div_float_hi : int -> ifloat -> ifloat -> ifloat
 200     val add_float_lo : int -> ifloat -> ifloat -> ifloat
 201     val add_float_hi : int -> ifloat -> ifloat -> ifloat
 202     val sub_float_lo : int -> ifloat -> ifloat -> ifloat
 203     val sub_float_hi : int -> ifloat -> ifloat -> ifloat
 204     val sqrt_float_lo : int -> ifloat -> ifloat
 205     val sqrt_float_hi : int -> ifloat -> ifloat
 206   end;;
 207
 208 module Informal_float : Informal_float_sig = struct
 209
 210 open Informal_nat;;
 211
 212 type ifloat = bool * nat * int;;
 213
 214 let min_exp = !Arith_options.min_exp;;
 215
 216 (* Creates a non-negative float *)
 217 let mk_float n e : ifloat = false, mk_nat n, e + min_exp;;
 218
 219 let mk_num_float n = false, mk_nat n, min_exp;;
 220
 221 let mk_small_num_float n = false, mk_small_nat n, min_exp;;
 222
 223 let zero_float = mk_small_num_float 0;;
 224
 225 let dest_float ((s, n, e) : ifloat) = s, dest_nat n, e;;
 226
 227 let sign_float ((s,_,_) : ifloat) = s;;
 228
 229 let eq_float (s1,n1,e1) (s2,n2,e2) = s1 = s2 && eq_nat n1 n2 && e1 = e2;;
 230
 231 let lo_float pp (s,n,e) =
 232   let n1, e1 = if s then hi_nat pp n else lo_nat pp n in
 233     (s, n1, e + e1);;
 234
 235 let hi_float pp (s,n,e) =
 236   let n1, e1 = if s then lo_nat pp n else hi_nat pp n in
 237     (s, n1, e + e1);;
 238
 239 (* Auxiliary num_exp functions *)
 240
 241 let num_exp_add =
 242   let (+) = add_nat in
 243     fun (n1,e1) (n2,e2) ->
 244       if e1 <= e2 then
 245         n1 + denormalize_nat (n2, e2 - e1), e1
 246       else
 247         n2 + denormalize_nat (n1, e1 - e2), e2;;
 248
 249
 250 (* Returns (n,e),true if (n1,e1) >= (n2,e2) and (n,e) = (n1,e1) - (n2,e2)
 251    Returns (n,e),false if (n1,e1) <= (n2,e2) and (n,e) = (n2,e2) - (n1,e1) *)
 252 let num_exp_sub =
 253   let (--) = sub_and_le_nat in
 254     fun (n1,e1) (n2,e2) ->
 255       if e2 <= e1 then
 256         let a = denormalize_nat (n1, e1 - e2) and
 257             b = n2 in
 258         let sub, flag = a -- b in
 259           (sub, e2), flag
 260       else
 261         let a = n1 and
 262             b = denormalize_nat (n2, e2 - e1) in
 263         let sub, flag = a -- b in
 264           (sub, e1), flag;;
 265
 266
 267 let num_exp_le =
 268   let (<=/) = le_nat in
 269     fun (n1,e1) (n2,e2) ->
 270       if e1 <= e2 then
 271         n1 <=/ denormalize_nat (n2, e2 - e1)
 272       else
 273         denormalize_nat (n1, e1 - e2) <=/ n2;;
 274
 275
 276 let num_exp_lt =
 277   let (</) = lt_nat in
 278     fun (n1,e1) (n2,e2) ->
 279       if e1 <= e2 then
 280         n1 </ denormalize_nat (n2, e2 - e1)
 281       else
 282         denormalize_nat (n1, e1 - e2) </ n2;;
 283
 284
 285
 286 (* neg *)
 287
 288 let neg_float (s,n,e) = (not s, n, e);;
 289
 290
 291 (* abs *)
 292
 293 let abs_float (_,n,e) = (false, n, e);;
 294
 295
 296
 297 (* lt0, gt0 *)
 298
 299 let lt0_float (s,n,e) =
 300   if not s then false else gt0_nat n;;
 301
 302 let gt0_float (s,n,e) =
 303   if s then false else gt0_nat n;;
 304
 305
 306 (* le0, ge0 *)
 307
 308 let le0_float (s,n,e) =
 309   if s then true else eq0_nat n;;
 310
 311 let ge0_float (s,n,e) =
 312   if s then eq0_nat n else true;;
 313
 314
 315 (* lt *)
 316
 317 let lt_float (s1,n1,e1) (s2,n2,e2) =
 318   if not s1 then
 319     if s2 then false else num_exp_lt (n1,e1) (n2,e2)
 320   else
 321     if s2 then num_exp_lt (n2,e2) (n1,e1)
 322     else
 323       (* TF *)
 324       if eq0_nat n1 then gt0_nat n2 else true;;
 325
 326
 327 let le_float (s1,n1,e1) (s2,n2,e2) =
 328   if s1 then
 329     if s2 then num_exp_le (n2,e2) (n1,e1) else true
 330   else
 331     if not s2 then num_exp_le (n1,e1) (n2,e2)
 332     else
 333       (* FT *)
 334       if eq0_nat n2 then eq0_nat n1 else false;;
 335
 336
 337
 338 (* min, max *)
 339
 340 let min_float f1 f2 =
 341   if le_float f1 f2 then f1 else f2;;
 342
 343 let max_float f1 f2 =
 344   if le_float f1 f2 then f2 else f1;;
 345
 346
 347 (* mul *)
 348
 349 let badd b1 b2 =
 350   if b1 then not b2 else b2;;
 351
 352
 353 let mul_float_eq (s1,n1,e1) (s2,n2,e2) =
 354   let s = badd s1 s2 in
 355   let n = mul_nat n1 n2 in
 356   let e = e1 + e2 - min_exp in
 357     if e < 0 then
 358       failwith "mul_float_eq: underflow"
 359     else
 360       (s, n, e);;
 361
 362
 363 let mul_float_lo pp (s1,n1,e1) (s2,n2,e2) =
 364   if eq0_nat n1 or eq0_nat n2 then
 365     zero_float
 366   else
 367     let s = badd s1 s2 in
 368     let n' = mul_nat n1 n2 in
 369     let n, e' = if s1 = s2 then lo_nat pp n' else hi_nat pp n' in
 370     let e = e1 + e2 + e' - min_exp in
 371       if e < 0 then
 372         failwith "mul_float_lo: underflow"
 373       else
 374         (s, n, e);;
 375
 376
 377 let mul_float_hi pp (s1,n1,e1) (s2,n2,e2) =
 378   if eq0_nat n1 or eq0_nat n2 then
 379     zero_float
 380   else
 381     let s = badd s1 s2 in
 382     let n' = mul_nat n1 n2 in
 383     let n, e' = if s1 = s2 then hi_nat pp n' else lo_nat pp n' in
 384     let e = e1 + e2 + e' - min_exp in
 385       if e < 0 then
 386         failwith "mul_float_hi: underflow"
 387       else
 388         (s, n, e);;
 389
 390
 391 (* div *)
 392
 393 let div_float_lo pp (s1,n1,e1) (s2,n2,e2) =
 394   if eq0_nat n1 then
 395     zero_float
 396   else
 397     let s = badd s1 s2 in
 398     let k = 2 * pp in
 399     let nn1 = denormalize_nat (n1, k) in
 400     let n' = div_nat nn1 n2 in
 401     let n, e' = if s1 = s2 then lo_nat pp n' else hi_lt_nat pp n' in
 402     let e = min_exp + e' + e1 - e2 - k in
 403       if e < 0 then
 404         failwith "div_float_lo: underflow"
 405       else
 406         (s, n, e);;
 407
 408
 409 let div_float_hi pp (s1,n1,e1) (s2,n2,e2) =
 410   if eq0_nat n1 then
 411     zero_float
 412   else
 413     let s = badd s1 s2 in
 414     let k = 2 * pp in
 415     let nn1 = denormalize_nat (n1, k) in
 416     let n' = div_nat nn1 n2 in
 417     let n, e' = if s1 = s2 then hi_lt_nat pp n' else lo_nat pp n' in
 418     let e = min_exp + e' + e1 - e2 - k in
 419       if e < 0 then
 420         failwith "div_float_hi: underflow"
 421       else
 422         (s, n, e);;
 423
 424
 425 (* add *)
 426
 427 let add_float_lo pp (s1,n1,e1) (s2,n2,e2) =
 428   if s1 = s2 then
 429     let n', e' = num_exp_add (n1,e1) (n2,e2) in
 430     let n, e'' = if s1 then hi_nat pp n' else lo_nat pp n' in
 431       (s1, n, e' + e'')
 432   else
 433     if s1 then
 434       let (n', e'), flag = num_exp_sub (n2,e2) (n1,e1) in
 435         if flag then
 436           let n, e'' = lo_nat pp n' in
 437             (false, n, e' + e'')
 438         else
 439           let n, e'' = hi_nat pp n' in
 440             (true, n, e' + e'')
 441     else
 442       let (n', e'), flag = num_exp_sub (n1,e1) (n2,e2) in
 443         if flag then
 444           let n, e'' = lo_nat pp n' in
 445             (false, n, e' + e'')
 446         else
 447           let n, e'' = hi_nat pp n' in
 448             (true, n, e' + e'');;
 449
 450
 451 let add_float_hi pp (s1,n1,e1) (s2,n2,e2) =
 452   if s1 = s2 then
 453     let n', e' = num_exp_add (n1,e1) (n2,e2) in
 454     let n, e'' = if s1 then lo_nat pp n' else hi_nat pp n' in
 455       (s1, n, e' + e'')
 456   else
 457     if s1 then
 458       let (n', e'), flag = num_exp_sub (n2,e2) (n1,e1) in
 459         if flag then
 460           let n, e'' = hi_nat pp n' in
 461             (false, n, e' + e'')
 462         else
 463           let n, e'' = lo_nat pp n' in
 464             (true, n, e' + e'')
 465     else
 466       let (n', e'), flag = num_exp_sub (n1,e1) (n2,e2) in
 467         if flag then
 468           let n, e'' = hi_nat pp n' in
 469             (false, n, e' + e'')
 470         else
 471           let n, e'' = lo_nat pp n' in
 472             (true, n, e' + e'');;
 473
 474
 475 (* sub *)
 476
 477 let sub_float_lo pp f1 f2 = add_float_lo pp f1 (neg_float f2);;
 478 let sub_float_hi pp f1 f2 = add_float_hi pp f1 (neg_float f2);;
 479
 480
 481 (* sqrt *)
 482
 483 let rec sqrt_float_lo pp (s,n1,e1) =
 484   if s then
 485     failwith "sqrt_float_lo: negative argument"
 486   else
 487     if e1 land 1 = 1 then
 488       sqrt_float_lo pp (s, denormalize_nat (n1, 1), e1 - 1)
 489     else
 490       let p2 = pp * 2 in
 491       let f1' = denormalize_nat (n1, p2) in
 492       let f1 = Big_int.sqrt_big_int (big_int_of_num (dest_nat f1')) in
 493       let n, e' = lo_nat pp (mk_nat (num_of_big_int f1)) in
 494       let e = ((e1 + min_exp) lsr 1) + e' - pp in
 495         if e < 0 then
 496           failwith "sqrt_float_lo: underflow"
 497         else
 498           (s, n, e);;
 499
 500
 501 let rec sqrt_float_hi pp (s,n1,e1) =
 502   if s then
 503     failwith "sqrt_float_hi: negative argument"
 504   else
 505     if e1 land 1 = 1 then
 506       sqrt_float_hi pp (s, denormalize_nat (n1, 1), e1 - 1)
 507     else
 508       let p2 = pp * 2 in
 509       let x = (big_int_of_num o dest_nat o denormalize_nat) (n1, p2) in
 510       let f1' = Big_int.sqrt_big_int x in
 511       let f1 = (mk_nat o num_of_big_int) f1' in
 512       let n, e' =
 513         let ( * ) = Big_int.mult_big_int and
 514             (==) = Big_int.eq_big_int in
 515           hi_nat pp (if f1' * f1' == x then f1 else suc_nat f1) in
 516       let e = ((e1 + min_exp) lsr 1) + e' - pp in
 517         if e < 0 then
 518           failwith "sqrt_float_hi: underflow"
 519         else
 520           (s, n, e);;
 521
 522 end;;
 523
 524
 525 (* Interval arithmetic *)
 526
 527 module type Informal_interval_sig =
 528   sig
 529     type interval
 530     val zero_interval : interval
 531     val one_interval : interval
 532     val two_interval : interval
 533     val mk_interval : Informal_float.ifloat * Informal_float.ifloat -> interval
 534     val mk_num_interval : num -> interval
 535     val mk_small_num_interval : int -> interval
 536     val dest_interval : interval -> Informal_float.ifloat * Informal_float.ifloat
 537     val round_interval : int -> interval -> interval
 538     val neg_interval : interval -> interval
 539     val mul_interval : int -> interval -> interval -> interval
 540     val div_interval : int -> interval -> interval -> interval
 541     val add_interval : int -> interval -> interval -> interval
 542     val sub_interval : int -> interval -> interval -> interval
 543     val sqrt_interval : int -> interval -> interval
 544     val inv_interval : int -> interval -> interval
 545     val pow_interval : int -> int -> interval -> interval
 546     (* Computes max(-lo, hi) *)
 547     val abs_interval : interval -> Informal_float.ifloat
 548   end;;
 549
 550 module Informal_interval : Informal_interval_sig = struct
 551
 552 open Informal_float;;
 553
 554
 555 type interval = ifloat * ifloat;;
 556
 557 let mk_interval (lo,hi) =
 558   if lt_float hi lo then failwith "mk_interval: hi < lo" else (lo,hi);;
 559
 560 let mk_num_interval n =
 561   let f = mk_num_float n in (f, f);;
 562
 563 let mk_small_num_interval n =
 564   let f = mk_small_num_float n in (f, f);;
 565
 566 let zero_interval = mk_small_num_interval 0;;
 567 let one_interval = mk_small_num_interval 1;;
 568 let two_interval = mk_small_num_interval 2;;
 569
 570 let dest_interval ((lo,hi) : interval) = (lo,hi);;
 571
 572 let round_interval pp (lo,hi) = (lo_float pp lo, hi_float pp hi);;
 573
 574 let neg_interval (lo,hi) = (neg_float hi, neg_float lo);;
 575
 576 let abs_interval (lo,hi) = max_float hi (neg_float lo);;
 577
 578 let add_interval pp (lo1,hi1) (lo2,hi2) =
 579   (add_float_lo pp lo1 lo2, add_float_hi pp hi1 hi2);;
 580
 581 let sub_interval pp (lo1,hi1) (lo2,hi2) =
 582   (sub_float_lo pp lo1 hi2, sub_float_hi pp hi1 lo2);;
 583
 584 let sqrt_interval pp (lo,hi) =
 585   if sign_float lo then
 586     failwith "sqrt_interval: negative lower bound"
 587   else
 588     (sqrt_float_lo pp lo, sqrt_float_hi pp hi);;
 589
 590 (* mul *)
 591 let mul_interval pp (l_lo,l_hi) (r_lo,r_hi) =
 592   let s1 = sign_float l_lo and
 593       s2 = sign_float l_hi and
 594       s3 = sign_float r_lo and
 595       s4 = sign_float r_hi in
 596     if s1 <> s2 && not s1 then
 597       zero_interval
 598     else if s3 <> s4 && not s3 then
 599       zero_interval
 600     else
 601       if s1 <> s2 && s3 <> s4 then
 602         let lo1, lo2 = mul_float_lo pp l_lo r_hi, mul_float_lo pp l_hi r_lo and
 603             hi1, hi2 = mul_float_hi pp l_lo r_lo, mul_float_hi pp l_hi r_hi in
 604           (min_float lo1 lo2, max_float hi1 hi2)
 605       else
 606         let lo1, lo2, hi1, hi2 =
 607           if s1 <> s2 then
 608             if not s3 then
 609               l_lo, r_hi, l_hi, r_hi
 610             else
 611               l_hi, r_lo, l_lo, r_lo
 612           else
 613             if s3 <> s4 then
 614               if not s1 then
 615                 l_hi, r_lo, l_hi, r_hi
 616               else
 617                 l_lo, r_hi, l_lo, r_lo
 618             else
 619               if not s1 then
 620                 if not s3 then
 621                   l_lo, r_lo, l_hi, r_hi
 622                 else
 623                   l_hi, r_lo, l_lo, r_hi
 624               else
 625                 if not s3 then
 626                   l_lo, r_hi, l_hi, r_lo
 627                 else
 628                   l_hi, r_hi, l_lo, r_lo in
 629           (mul_float_lo pp lo1 lo2, mul_float_hi pp hi1 hi2);;
 630
 631 (* div *)
 632 let div_interval pp (l_lo,l_hi) (r_lo,r_hi) =
 633   let s1 = sign_float l_lo and
 634       s2 = sign_float l_hi and
 635       s3 = sign_float r_lo and
 636       s4 = sign_float r_hi in
 637     if s1 <> s2 && not s1 then
 638       zero_interval
 639     else if s3 <> s4 then
 640       failwith "div_interval: division by an interval containing 0"
 641     else
 642       let lo1, lo2, hi1, hi2 =
 643         if s1 = s2 then
 644           if not s1 then
 645             if not s3 then
 646               l_lo, r_hi, l_hi, r_lo
 647             else
 648               l_hi, r_hi, l_lo, r_lo
 649           else
 650             if not s3 then
 651               l_lo, r_lo, l_hi, r_hi
 652             else
 653               l_hi, r_lo, l_lo, r_hi
 654         else
 655           if not s3 then
 656             l_lo, r_lo, l_hi, r_lo
 657           else
 658             l_hi, r_hi, l_lo, r_hi in
 659         (div_float_lo pp lo1 lo2, div_float_hi pp hi1 hi2);;
 660
 661 (* inv *)
 662 let inv_interval pp int =
 663   div_interval pp one_interval int;;
 664
 665 (* pow *)
 666 let pow_interval pp n int =
 667   let rec pow n =
 668     if n <= 0 then
 669       one_interval
 670     else if n = 1 then
 671       int
 672     else
 673       let i2 = pow (n - 1) in
 674         mul_interval pp int i2 in
 675     pow n;;
 676
 677
 678 (* Arith_misc.gen_pow (mul_interval pp) one_interval n;; *)
 679
 680 end;;
 681
 682
 683 (* atn *)
 684 module type Informal_atn_sig =
 685   sig
 686     val atn_interval : int -> Informal_interval.interval -> Informal_interval.interval
 687     val acs_interval : int -> Informal_interval.interval -> Informal_interval.interval
 688     val pi_approx_array : Informal_interval.interval array
 689     val pi2_approx_array : Informal_interval.interval array
 690   end;;
 691
 692
 693 module Informal_atn : Informal_atn_sig = struct
 694
 695 open Informal_float;;
 696 open Informal_interval;;
 697
 698
 699 let rec poly_f_interval pp l x =
 700   if length l = 0 then
 701     failwith "poly_f_interval: an empty coefficient list"
 702   else
 703     let first = hd l in
 704       if length l = 1 then
 705         first
 706       else
 707         let r = poly_f_interval pp (tl l) x in
 708           add_interval pp first (mul_interval pp x r);;
 709
 710
 711 let poly_f_even_interval pp l x =
 712   let xx = mul_interval pp x x in
 713     poly_f_interval pp l xx;;
 714
 715
 716 let poly_f_odd_interval pp l x =
 717   let even = poly_f_even_interval pp l x in
 718     mul_interval pp x even;;
 719
 720
 721 let halfatn_interval pp x =
 722   let xx = mul_interval pp x x in
 723   let one_xx = add_interval pp one_interval xx in
 724   let sqrt = sqrt_interval pp one_xx in
 725   let one_sqrt = add_interval pp sqrt one_interval in
 726     div_interval pp x one_sqrt;;
 727
 728
 729 let halfatn4_interval pp x =
 730   (halfatn_interval pp o halfatn_interval pp o halfatn_interval pp o halfatn_interval pp) x;;
 731
 732
 733 (* Computes an interval for 16 * sum(0..n) (halfatn4_co x) *)
 734 let atn_sum_interval =
 735   let interval_16 = mk_small_num_interval 16 in
 736     fun pp l x ->
 737       let halfatn4 = halfatn4_interval pp x in
 738       let poly = poly_f_odd_interval pp l halfatn4 in
 739         mul_interval pp interval_16 poly;;
 740
 741
 742 let atn0_interval pp l eps x =
 743   let sum = atn_sum_interval pp l x in
 744   let a, b = dest_interval sum in
 745   let _, d = dest_interval eps in
 746   let hi = add_float_hi pp b d in
 747   let lo = sub_float_lo pp a d in
 748     mk_interval (lo, hi);;
 749
 750
 751 (* Computes an interval for 2 ^ -(6n + 5) *)
 752 let compute_eps1 pp n =
 753   let pow = pow_interval pp (6 * n + 5) two_interval in
 754     inv_interval pp pow;;
 755
 756
 757 let mk_atn_co_table pp n =
 758   let get_val k =
 759     let l = if (k land 1) = 0 then one_interval else neg_interval (one_interval) in
 760     let r = mk_small_num_interval (2 * k + 1) in
 761       div_interval pp l r in
 762     map get_val (0--n);;
 763
 764 (* Lookup tables *)
 765 let n_of_p pp =
 766   let x = (float_of_int (pp + 1) *. log (float_of_int Informal_nat.arith_base) /. log (2.0) -. 5.0) /. 6.0 in
 767   let n = (int_of_float o ceil) x in
 768     if n < 1 then 1 else n;;
 769
 770 let atn_co_array = Array.init 21 (fun i -> mk_atn_co_table (i + 1) (n_of_p i));;
 771 let eps1_array = Array.init 21 (fun i -> compute_eps1 (i + 1) (n_of_p i));;
 772
 773
 774 let atn_interval pp x =
 775   atn0_interval pp atn_co_array.(pp) eps1_array.(pp) x;;
 776
 777
 778 (* pi approximation *)
 779 let pi_approx_array, pi2_approx_array =
 780   let pp = 20 in
 781   let x = one_interval in
 782   let r1 = atn_interval pp x in
 783   let r2 = mul_interval pp (mk_small_num_interval 4) r1 in
 784   let float_pi = r2 in
 785   let float_pi2 = div_interval pp float_pi two_interval in
 786   let pi_int0 = mk_small_num_interval 0 in
 787   let pi2_int0 = pi_int0 in
 788     Array.init 19 (fun i -> if i = 0 then pi_int0 else round_interval i float_pi),
 789   Array.init 19 (fun i -> if i = 0 then pi2_int0 else round_interval i float_pi2);;
 790
 791
 792 (* acs *)
 793 let acs0_interval pp l eps1 x =
 794   let int1 = sub_interval pp one_interval (mul_interval pp x x) in
 795   let int2 = div_interval pp x (sqrt_interval pp int1) in
 796   let atn_int = atn0_interval pp l eps1 int2 in
 797     sub_interval pp pi2_approx_array.(pp + 1) atn_int;;
 798
 799
 800 let acs_interval pp x =
 801   acs0_interval pp atn_co_array.(pp) eps1_array.(pp) x;;
 802
 803
 804 end;;
 805