(* Copyright (C) <2012> This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program. If not, see . *) open Big_int let precision = 100 let base = power_int_positive_int 10 precision type t = { mantisse:big_int; exponent:int; } let big_float_of_int i = {mantisse=big_int_of_int i;exponent=0} let is_neg ta = sign_big_int ta.mantisse < 0 let rec compare ta tb = if sign_big_int ta.mantisse = 0 && compare_big_int ta.mantisse tb.mantisse = 0 then 0 else if ta.exponent < tb.exponent then - (compare tb ta) else if ta.exponent > tb.exponent + precision then 1 else compare_big_int tb.mantisse (mult_big_int ta.mantisse (power_int_positive_int 10 (ta.exponent - tb.exponent))) let string_of_big_float t = let is_neg = sign_big_int t.mantisse < 0 in let str = string_of_big_int t.mantisse in let index = if is_neg then 2 else 1 in let fpart = (String.sub str 0 index) in fpart ^ "."^(Str.global_replace (Str.regexp ("0*$")) "" (String.sub str index (String.length str - index)))^(let exponent = t.exponent + (String.length str - index) in if exponent <> 0 then "e"^(string_of_int exponent) else "") let zero = big_float_of_int 0 let one = big_float_of_int 1 let add ta tb = let ta,tb = if ta.exponent > tb.exponent then ta,tb else tb,ta in if compare_big_int ta.mantisse zero_big_int = 0 then tb else if compare_big_int tb.mantisse zero_big_int = 0 then ta else if ta.exponent > tb.exponent + precision then ta else let mantisse = add_big_int (mult_big_int ta.mantisse (power_int_positive_int 10 (ta.exponent-tb.exponent))) tb.mantisse in let str = string_of_big_int mantisse in let numstr = String.length str in if numstr > precision then let valeur = tb.exponent + numstr - precision in {mantisse=big_int_of_string (String.sub str 0 precision);exponent=valeur} else {mantisse=mantisse;exponent=tb.exponent} let minus ta = {mantisse=minus_big_int ta.mantisse;exponent=ta.exponent} let sub ta tb = add ta (minus tb) let mult ta tb = let ta,tb = if ta.exponent > tb.exponent then ta,tb else tb,ta in if compare_big_int ta.mantisse zero_big_int = 0 then zero else if compare_big_int tb.mantisse zero_big_int = 0 then zero else let mantisse = mult_big_int (mult_big_int ta.mantisse (power_int_positive_int 10 (ta.exponent-tb.exponent))) (tb.mantisse) in let str = string_of_big_int mantisse in let numstr = String.length str in let pr = min precision numstr in {mantisse=big_int_of_string (String.sub str 0 pr);exponent=tb.exponent + tb.exponent + numstr - pr} let normalize ta = let size = String.length (string_of_big_int ta.mantisse) - (if is_neg ta then 1 else 0) in {mantisse= mult_big_int ta.mantisse (power_int_positive_int 10 (precision-size));exponent = ta.exponent - precision + size} let inv ta = let ta = normalize ta in let res = string_of_big_int (div_big_int (power_int_positive_int 10 (precision + precision - 1)) ta.mantisse) in let reste = if String.length res > precision + (if is_neg ta then 1 else 0) then 1 else 0 in let reste = reste + (if is_neg ta then 1 else 0) in {mantisse=big_int_of_string (String.sub res 0 precision);exponent= 1 - ta.exponent - precision - precision + reste} let div ta tb = mult ta (inv tb) let two = add one one let rec power t i = match i with |0 -> one |1 -> t |_ -> let sur2 = power t (i/2) in if i mod 2 = 0 then mult sur2 sur2 else mult t (mult sur2 sur2) let rec fact ta = if is_neg ta then raise (Invalid_argument "fact") else if compare ta zero = 0 then one else mult ta (fact (sub ta one)) let exp t = let res = ref zero in let new_val = ref one in let power = ref one in let fact = ref one in let i = ref 0 in while compare !res !new_val <> 0 do incr i; res := !new_val; power := mult !power t; fact := mult !fact (big_float_of_int !i); let div = div !power !fact in new_val := add !res div; done; !new_val let e = exp one let log t = let res = ref zero in let z = div (sub t one) (add t one) in let new_val = ref z in let power = ref z in let z2 = mult z z in let i = ref 1 in while compare !res !new_val <> 0 do i := !i + 2; res := !new_val; let bigi = big_float_of_int !i in power := mult !power z2; let div = div !power bigi in new_val := add !res div; done; mult (add one one) !new_val let log10 = log (big_float_of_int 10) let log t = let t = normalize t in let a = {mantisse=t.mantisse;exponent= - precision + 1} in let lt = log a in add lt (mult log10 (big_float_of_int (t.exponent+precision-1))) let pi = let res = ref one in let new_val = ref zero in let i = ref (-1) in let n396 = big_float_of_int 396 in let n1103 = big_float_of_int 1103 in let n26390 = big_float_of_int 26390 in while compare !res !new_val <> 0 do res := !new_val; incr i; let k = big_float_of_int !i in let k4 = fact (big_float_of_int (4 * !i)) in let kfact = fact k in let lin = add n1103 (mult k n26390) in let den = power n396 (4 * !i) in new_val := add !res (div (mult k4 lin) (mult kfact den)); done; let n9801 = big_float_of_int 9801 in let n2sqrt2 = mult two (exp (div (log two) two)) in inv (div (mult n2sqrt2 !new_val) n9801) let sin t = let res = ref zero in let new_val = ref t in let power = ref t in let t2 = mult t t in let min = ref false in let bigi = ref one in let i = ref 1 in while compare !res !new_val <> 0 do min := not !min; i := !i + 2; res := !new_val; bigi := mult !bigi (mult (big_float_of_int !i) (big_float_of_int (!i-1))); power := mult !power t2; let div = (if !min then minus else (fun x -> x)) (div !power !bigi) in new_val := add !res div; done; !new_val let ( =. ) ta tb = compare ta tb = 0 let ( *. ) ta tb = mult ta tb let ( /. ) ta tb = div ta tb let ( +. ) ta tb = add ta tb let ( -. ) ta tb = sub ta tb let ( **. ) ta i = power ta i let display ta = print_endline (string_of_big_float ta)