formal_lp/old/arith/informal/informal_arith.hl

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