diff options
Diffstat (limited to 'otherlibs/num/ratio.ml')
| -rw-r--r-- | otherlibs/num/ratio.ml | 619 |
1 files changed, 0 insertions, 619 deletions
diff --git a/otherlibs/num/ratio.ml b/otherlibs/num/ratio.ml deleted file mode 100644 index 04f9c5e8ff..0000000000 --- a/otherlibs/num/ratio.ml +++ /dev/null @@ -1,619 +0,0 @@ -(**************************************************************************) -(* *) -(* OCaml *) -(* *) -(* Valerie Menissier-Morain, projet Cristal, INRIA Rocquencourt *) -(* *) -(* Copyright 1996 Institut National de Recherche en Informatique et *) -(* en Automatique. *) -(* *) -(* All rights reserved. This file is distributed under the terms of *) -(* the GNU Lesser General Public License version 2.1, with the *) -(* special exception on linking described in the file LICENSE. *) -(* *) -(**************************************************************************) - -open Int_misc -open Nat -open Big_int -open Arith_flags - -(* Definition of the type ratio : - Conventions : - - the denominator is always a positive number - - the sign of n/0 is the sign of n -These convention is automatically respected when a ratio is created with -the create_ratio primitive -*) - -type ratio = { mutable numerator : big_int; - mutable denominator : big_int; - mutable normalized : bool} - -let failwith_zero name = - let s = "infinite or undefined rational number" in - failwith (if String.length name = 0 then s else name ^ " " ^ s) - -let numerator_ratio r = r.numerator -and denominator_ratio r = r.denominator - -let null_denominator r = sign_big_int r.denominator = 0 - -let verify_null_denominator r = - if sign_big_int r.denominator = 0 - then (if !error_when_null_denominator_flag - then (failwith_zero "") - else true) - else false - -let sign_ratio r = sign_big_int r.numerator - -(* Physical normalization of rational numbers *) -(* 1/0, 0/0 and -1/0 are the normalized forms for n/0 numbers *) -let normalize_ratio r = - if r.normalized then r - else if verify_null_denominator r then begin - r.numerator <- big_int_of_int (sign_big_int r.numerator); - r.normalized <- true; - r - end else begin - let p = gcd_big_int r.numerator r.denominator in - if eq_big_int p unit_big_int - then begin - r.normalized <- true; r - end else begin - r.numerator <- div_big_int (r.numerator) p; - r.denominator <- div_big_int (r.denominator) p; - r.normalized <- true; r - end - end - -let cautious_normalize_ratio r = - if (!normalize_ratio_flag) then (normalize_ratio r) else r - -let cautious_normalize_ratio_when_printing r = - if (!normalize_ratio_when_printing_flag) then (normalize_ratio r) else r - -let create_ratio bi1 bi2 = - match sign_big_int bi2 with - -1 -> cautious_normalize_ratio - { numerator = minus_big_int bi1; - denominator = minus_big_int bi2; - normalized = false } - | 0 -> if !error_when_null_denominator_flag - then (failwith_zero "create_ratio") - else cautious_normalize_ratio - { numerator = bi1; denominator = bi2; normalized = false } - | _ -> cautious_normalize_ratio - { numerator = bi1; denominator = bi2; normalized = false } - -let create_normalized_ratio bi1 bi2 = - match sign_big_int bi2 with - -1 -> { numerator = minus_big_int bi1; - denominator = minus_big_int bi2; - normalized = true } -| 0 -> if !error_when_null_denominator_flag - then failwith_zero "create_normalized_ratio" - else { numerator = bi1; denominator = bi2; normalized = true } -| _ -> { numerator = bi1; denominator = bi2; normalized = true } - -let is_normalized_ratio r = r.normalized - -let report_sign_ratio r bi = - if sign_ratio r = -1 - then minus_big_int bi - else bi - -let abs_ratio r = - { numerator = abs_big_int r.numerator; - denominator = r.denominator; - normalized = r.normalized } - -let is_integer_ratio r = - eq_big_int ((normalize_ratio r).denominator) unit_big_int - -(* Operations on rational numbers *) - -let add_ratio r1 r2 = - if !normalize_ratio_flag then begin - let p = gcd_big_int ((normalize_ratio r1).denominator) - ((normalize_ratio r2).denominator) in - if eq_big_int p unit_big_int then - {numerator = add_big_int (mult_big_int (r1.numerator) r2.denominator) - (mult_big_int (r2.numerator) r1.denominator); - denominator = mult_big_int (r1.denominator) r2.denominator; - normalized = true} - else begin - let d1 = div_big_int (r1.denominator) p - and d2 = div_big_int (r2.denominator) p in - let n = add_big_int (mult_big_int (r1.numerator) d2) - (mult_big_int d1 r2.numerator) in - let p' = gcd_big_int n p in - { numerator = div_big_int n p'; - denominator = mult_big_int d1 (div_big_int (r2.denominator) p'); - normalized = true } - end - end else - { numerator = add_big_int (mult_big_int (r1.numerator) r2.denominator) - (mult_big_int (r1.denominator) r2.numerator); - denominator = mult_big_int (r1.denominator) r2.denominator; - normalized = false } - -let minus_ratio r = - { numerator = minus_big_int (r.numerator); - denominator = r.denominator; - normalized = r.normalized } - -let add_int_ratio i r = - ignore (cautious_normalize_ratio r); - { numerator = add_big_int (mult_int_big_int i r.denominator) r.numerator; - denominator = r.denominator; - normalized = r.normalized } - -let add_big_int_ratio bi r = - ignore (cautious_normalize_ratio r); - { numerator = add_big_int (mult_big_int bi r.denominator) r.numerator ; - denominator = r.denominator; - normalized = r.normalized } - -let sub_ratio r1 r2 = add_ratio r1 (minus_ratio r2) - -let mult_ratio r1 r2 = - if !normalize_ratio_flag then begin - let p1 = gcd_big_int ((normalize_ratio r1).numerator) - ((normalize_ratio r2).denominator) - and p2 = gcd_big_int (r2.numerator) r1.denominator in - let (n1, d2) = - if eq_big_int p1 unit_big_int - then (r1.numerator, r2.denominator) - else (div_big_int (r1.numerator) p1, div_big_int (r2.denominator) p1) - and (n2, d1) = - if eq_big_int p2 unit_big_int - then (r2.numerator, r1.denominator) - else (div_big_int r2.numerator p2, div_big_int r1.denominator p2) in - { numerator = mult_big_int n1 n2; - denominator = mult_big_int d1 d2; - normalized = true } - end else - { numerator = mult_big_int (r1.numerator) r2.numerator; - denominator = mult_big_int (r1.denominator) r2.denominator; - normalized = false } - -let mult_int_ratio i r = - if !normalize_ratio_flag then - begin - let p = gcd_big_int ((normalize_ratio r).denominator) (big_int_of_int i) in - if eq_big_int p unit_big_int - then { numerator = mult_big_int (big_int_of_int i) r.numerator; - denominator = r.denominator; - normalized = true } - else { numerator = mult_big_int (div_big_int (big_int_of_int i) p) - r.numerator; - denominator = div_big_int (r.denominator) p; - normalized = true } - end - else - { numerator = mult_int_big_int i r.numerator; - denominator = r.denominator; - normalized = false } - -let mult_big_int_ratio bi r = - if !normalize_ratio_flag then - begin - let p = gcd_big_int ((normalize_ratio r).denominator) bi in - if eq_big_int p unit_big_int - then { numerator = mult_big_int bi r.numerator; - denominator = r.denominator; - normalized = true } - else { numerator = mult_big_int (div_big_int bi p) r.numerator; - denominator = div_big_int (r.denominator) p; - normalized = true } - end - else - { numerator = mult_big_int bi r.numerator; - denominator = r.denominator; - normalized = false } - -let square_ratio r = - ignore (cautious_normalize_ratio r); - { numerator = square_big_int r.numerator; - denominator = square_big_int r.denominator; - normalized = r.normalized } - -let inverse_ratio r = - if !error_when_null_denominator_flag && (sign_big_int r.numerator) = 0 - then failwith_zero "inverse_ratio" - else {numerator = report_sign_ratio r r.denominator; - denominator = abs_big_int r.numerator; - normalized = r.normalized} - -let div_ratio r1 r2 = - mult_ratio r1 (inverse_ratio r2) - -(* Integer part of a rational number *) -(* Odd function *) -let integer_ratio r = - if null_denominator r then failwith_zero "integer_ratio" - else if sign_ratio r = 0 then zero_big_int - else report_sign_ratio r (div_big_int (abs_big_int r.numerator) - (abs_big_int r.denominator)) - -(* Floor of a rational number *) -(* Always less or equal to r *) -let floor_ratio r = - ignore (verify_null_denominator r); - div_big_int (r.numerator) r.denominator - -(* Round of a rational number *) -(* Odd function, 1/2 -> 1 *) -let round_ratio r = - ignore (verify_null_denominator r); - let abs_num = abs_big_int r.numerator in - let bi = div_big_int abs_num r.denominator in - report_sign_ratio r - (if sign_big_int - (sub_big_int - (mult_int_big_int - 2 - (sub_big_int abs_num (mult_big_int (r.denominator) bi))) - r.denominator) = -1 - then bi - else succ_big_int bi) - -let ceiling_ratio r = - if (is_integer_ratio r) - then r.numerator - else succ_big_int (floor_ratio r) - - -(* Comparison operators on rational numbers *) -let eq_ratio r1 r2 = - ignore (normalize_ratio r1); - ignore (normalize_ratio r2); - eq_big_int (r1.numerator) r2.numerator && - eq_big_int (r1.denominator) r2.denominator - -let compare_ratio r1 r2 = - if verify_null_denominator r1 then - let sign_num_r1 = sign_big_int r1.numerator in - if (verify_null_denominator r2) - then - let sign_num_r2 = sign_big_int r2.numerator in - if sign_num_r1 = 1 && sign_num_r2 = -1 then 1 - else if sign_num_r1 = -1 && sign_num_r2 = 1 then -1 - else 0 - else sign_num_r1 - else if verify_null_denominator r2 then - -(sign_big_int r2.numerator) - else match compare_int (sign_big_int r1.numerator) - (sign_big_int r2.numerator) with - 1 -> 1 - | -1 -> -1 - | _ -> if eq_big_int (r1.denominator) r2.denominator - then compare_big_int (r1.numerator) r2.numerator - else compare_big_int - (mult_big_int (r1.numerator) r2.denominator) - (mult_big_int (r1.denominator) r2.numerator) - - -let lt_ratio r1 r2 = compare_ratio r1 r2 < 0 -and le_ratio r1 r2 = compare_ratio r1 r2 <= 0 -and gt_ratio r1 r2 = compare_ratio r1 r2 > 0 -and ge_ratio r1 r2 = compare_ratio r1 r2 >= 0 - -let max_ratio r1 r2 = if lt_ratio r1 r2 then r2 else r1 -and min_ratio r1 r2 = if gt_ratio r1 r2 then r2 else r1 - -let eq_big_int_ratio bi r = - (is_integer_ratio r) && eq_big_int bi r.numerator - -let compare_big_int_ratio bi r = - ignore (normalize_ratio r); - if (verify_null_denominator r) - then -(sign_big_int r.numerator) - else compare_big_int (mult_big_int bi r.denominator) r.numerator - -let lt_big_int_ratio bi r = compare_big_int_ratio bi r < 0 -and le_big_int_ratio bi r = compare_big_int_ratio bi r <= 0 -and gt_big_int_ratio bi r = compare_big_int_ratio bi r > 0 -and ge_big_int_ratio bi r = compare_big_int_ratio bi r >= 0 - -(* Coercions *) - -(* Coercions with type int *) -let int_of_ratio r = - if ((is_integer_ratio r) && (is_int_big_int r.numerator)) - then (int_of_big_int r.numerator) - else failwith "integer argument required" - -and ratio_of_int i = - { numerator = big_int_of_int i; - denominator = unit_big_int; - normalized = true } - -(* Coercions with type nat *) -let ratio_of_nat nat = - { numerator = big_int_of_nat nat; - denominator = unit_big_int; - normalized = true } - -and nat_of_ratio r = - ignore (normalize_ratio r); - if not (is_integer_ratio r) then - failwith "nat_of_ratio" - else if sign_big_int r.numerator > -1 then - nat_of_big_int (r.numerator) - else failwith "nat_of_ratio" - -(* Coercions with type big_int *) -let ratio_of_big_int bi = - { numerator = bi; denominator = unit_big_int; normalized = true } - -and big_int_of_ratio r = - ignore (normalize_ratio r); - if is_integer_ratio r - then r.numerator - else failwith "big_int_of_ratio" - -let div_int_ratio i r = - ignore (verify_null_denominator r); - mult_int_ratio i (inverse_ratio r) - -let div_ratio_int r i = - div_ratio r (ratio_of_int i) - -let div_big_int_ratio bi r = - ignore (verify_null_denominator r); - mult_big_int_ratio bi (inverse_ratio r) - -let div_ratio_big_int r bi = - div_ratio r (ratio_of_big_int bi) - -(* Functions on type string *) -(* giving floating point approximations of rational numbers *) - -(* Compares strings that contains only digits, have the same length, - from index i to index i + l *) -let rec compare_num_string s1 s2 i len = - if i >= len then 0 else - let c1 = int_of_char s1.[i] - and c2 = int_of_char s2.[i] in - match compare_int c1 c2 with - | 0 -> compare_num_string s1 s2 (succ i) len - | c -> c;; - -(* Position of the leading digit of the decimal expansion *) -(* of a strictly positive rational number *) -(* if the decimal expansion of a non null rational r is equal to *) -(* sigma for k=-P to N of r_k*10^k then msd_ratio r = N *) -(* Nota : for a big_int we have msd_ratio = nums_digits_big_int -1 *) - -(* Tests if s has only zeros characters from index i to index lim *) -let rec only_zeros s i lim = - i >= lim || s.[i] == '0' && only_zeros s (succ i) lim;; - -(* Nota : for a big_int we have msd_ratio = nums_digits_big_int -1 *) -let msd_ratio r = - ignore (cautious_normalize_ratio r); - if null_denominator r then failwith_zero "msd_ratio" - else if sign_big_int r.numerator == 0 then 0 - else begin - let str_num = string_of_big_int r.numerator - and str_den = string_of_big_int r.denominator in - let size_num = String.length str_num - and size_den = String.length str_den in - let size_min = min size_num size_den in - let m = size_num - size_den in - let cmp = compare_num_string str_num str_den 0 size_min in - match cmp with - | 1 -> m - | -1 -> pred m - | _ -> - if m >= 0 then m else - if only_zeros str_den size_min size_den then m - else pred m - end -;; - -(* Decimal approximations of rational numbers *) - -(* Approximation with fix decimal point *) -(* This is an odd function and the last digit is round off *) -(* Format integer_part . decimal_part_with_n_digits *) -let approx_ratio_fix n r = - (* Don't need to normalize *) - if (null_denominator r) then failwith_zero "approx_ratio_fix" - else - let sign_r = sign_ratio r in - if sign_r = 0 - then "+0" (* r = 0 *) - else - (* r.numerator and r.denominator are not null numbers - s1 contains one more digit than desired for the round off operation *) - if n >= 0 then begin - let s1 = - string_of_nat - (nat_of_big_int - (div_big_int - (base_power_big_int - 10 (succ n) (abs_big_int r.numerator)) - r.denominator)) in - (* Round up and add 1 in front if needed *) - let s2 = - if round_futur_last_digit (Bytes.unsafe_of_string s1) 0 - (String.length s1) - then "1" ^ s1 - else s1 in - let l2 = String.length s2 - 1 in - (* if s2 without last digit is xxxxyyy with n 'yyy' digits: - <sign> xxxx . yyy - if s2 without last digit is yy with <= n digits: - <sign> 0 . 0yy *) - if l2 > n then begin - let s = Bytes.make (l2 + 2) '0' in - Bytes.set s 0 (if sign_r = -1 then '-' else '+'); - String.blit s2 0 s 1 (l2 - n); - Bytes.set s (l2 - n + 1) '.'; - String.blit s2 (l2 - n) s (l2 - n + 2) n; - Bytes.unsafe_to_string s - end else begin - let s = Bytes.make (n + 3) '0' in - Bytes.set s 0 (if sign_r = -1 then '-' else '+'); - Bytes.set s 2 '.'; - String.blit s2 0 s (n + 3 - l2) l2; - Bytes.unsafe_to_string s - end - end else begin - (* Dubious; what is this code supposed to do? *) - let s = string_of_big_int - (div_big_int - (abs_big_int r.numerator) - (base_power_big_int - 10 (-n) r.denominator)) in - let len = succ (String.length s) in - let s' = Bytes.make len '0' in - Bytes.set s' 0 (if sign_r = -1 then '-' else '+'); - String.blit s 0 s' 1 (pred len); - Bytes.unsafe_to_string s' - end - -(* Number of digits of the decimal representation of an int *) -let num_decimal_digits_int n = - String.length (string_of_int n) - -(* Approximation with floating decimal point *) -(* This is an odd function and the last digit is round off *) -(* Format (+/-)(0. n_first_digits e msd)/(1. n_zeros e (msd+1) *) -let approx_ratio_exp n r = - (* Don't need to normalize *) - if (null_denominator r) then failwith_zero "approx_ratio_exp" - else if n <= 0 then invalid_arg "approx_ratio_exp" - else - let sign_r = sign_ratio r - and i = ref (n + 3) in - if sign_r = 0 then - String.concat "" ["+0."; String.make n '0'; "e0"] - else - let msd = msd_ratio (abs_ratio r) in - let k = n - msd in - let s = - (let nat = nat_of_big_int - (if k < 0 - then - div_big_int (abs_big_int r.numerator) - (base_power_big_int 10 (- k) - r.denominator) - else - div_big_int (base_power_big_int - 10 k (abs_big_int r.numerator)) - r.denominator) in - string_of_nat nat) in - if round_futur_last_digit (Bytes.unsafe_of_string s) 0 - (String.length s) - then - let m = num_decimal_digits_int (succ msd) in - let str = Bytes.make (n + m + 4) '0' in - (String.blit (if sign_r = -1 then "-1." else "+1.") 0 str 0 3); - Bytes.set str !i ('e'); - incr i; - (if m = 0 - then Bytes.set str !i '0' - else String.blit (string_of_int (succ msd)) 0 str !i m); - Bytes.unsafe_to_string str - else - let m = num_decimal_digits_int (succ msd) - and p = n + 3 in - let str = Bytes.make (succ (m + p)) '0' in - (String.blit (if sign_r = -1 then "-0." else "+0.") 0 str 0 3); - (String.blit s 0 str 3 n); - Bytes.set str p 'e'; - (if m = 0 - then Bytes.set str (succ p) '0' - else (String.blit (string_of_int (succ msd)) 0 str (succ p) m)); - Bytes.unsafe_to_string str - -(* String approximation of a rational with a fixed number of significant *) -(* digits printed *) -let float_of_rational_string r = - let s = approx_ratio_exp !floating_precision r in - if String.get s 0 = '+' - then (String.sub s 1 (pred (String.length s))) - else s - -(* Coercions with type string *) -let string_of_ratio r = - ignore (cautious_normalize_ratio_when_printing r); - if !approx_printing_flag - then float_of_rational_string r - else string_of_big_int r.numerator ^ "/" ^ string_of_big_int r.denominator - -(* XL: j'ai puissamment simplifie "ratio_of_string" en virant la notation - scientifique. - | I have strongly simplified "ratio_of_string" by deleting scientific notation -*) - -let ratio_of_string s = - try - let n = String.index s '/' in - create_ratio (sys_big_int_of_string s 0 n) - (sys_big_int_of_string s (n+1) (String.length s - n - 1)) - with Not_found -> - { numerator = big_int_of_string s; - denominator = unit_big_int; - normalized = true } - -(* Coercion with type float *) - -let float_of_ratio r = - let p = r.numerator and q = r.denominator in - (* Special cases 0/0, 0/q and p/0 *) - if sign_big_int q = 0 then begin - match sign_big_int p with - | 0 -> nan - | 1 -> infinity - | -1 -> neg_infinity - | _ -> assert false - end - else if sign_big_int p = 0 then 0.0 - else begin - let np = num_bits_big_int p and nq = num_bits_big_int q in - if np <= 53 && nq <= 53 then - (* p and q convert to floats exactly; use FP division to get the - correctly-rounded result. *) - Int64.to_float (int64_of_big_int p) - /. Int64.to_float (int64_of_big_int q) - else begin - let ap = abs_big_int p in - (* |p| is in [2^(np-1), 2^np) - q is in [2^(nq-1), 2^nq) - hence |p|/q is in (2^(np-nq-1), 2^(np-nq+1)). - We define n such that |p|/q*2^n is in [2^54, 2^56). - >= 2^54 so that the round to odd technique applies. - < 2^56 so that the integral part is representable as an int64. *) - let n = 55 - (np - nq) in - (* Scaling |p|/q by 2^n *) - let (p', q') = - if n >= 0 - then (shift_left_big_int ap n, q) - else (ap, shift_left_big_int q (-n)) in - (* Euclidean division of p' by q' *) - let (quo, rem) = quomod_big_int p' q' in - (* quo is the integral part of |p|/q*2^n - rem/q' is the fractional part. *) - (* Round quo to float *) - let f = round_big_int_to_float quo (sign_big_int rem = 0) in - (* Apply exponent *) - let f = ldexp f (-n) in - (* Apply sign *) - if sign_big_int p < 0 then -. f else f - end - end - - -let power_ratio_positive_int r n = - create_ratio (power_big_int_positive_int (r.numerator) n) - (power_big_int_positive_int (r.denominator) n) - -let power_ratio_positive_big_int r bi = - create_ratio (power_big_int_positive_big_int (r.numerator) bi) - (power_big_int_positive_big_int (r.denominator) bi) |