(* :Copyright: Copyright 1990-1999, Wolfram Research, Inc. *)

(*:Mathematica Version: 2.0 *)

(*:Package Version: 1.0 *)

(*:Name: NumberTheory`Rationalize` *)

(*:Author:
	Daniel R. Grayson
*)

(*:Keywords:
	rationalize, real number
*)

(*:Requirements: None. *)

(*:Warnings:
*)

(*:Sources:
*)

(*:Summary:
The built-in function Rationalize gives a rational approximation to a single
real number.  The functions provided by this package give a list of 
rational numbers that approximate a list of real numbers.  The approximation
is simultaneous and preserves the relationships among the given numbers.
*)

BeginPackage["NumberTheory`Rationalize`"]

ProjectiveRationalize::usage = "ProjectiveRationalize[{x0,x1,...,xn}] returns
a list of integers {s0,s1,...,sn} so that the ratios si/sj approximate
the ratios xi/xj well for each pair 0 <= i,j <= n.
ProjectiveRationalize[{x0,x1,...,xn}, prec] allows a tolerance of 10^-prec.
The numbers xi must be real."

AffineRationalize::usage = "AffineRationalize[{x1,...,xn}] returns rational
numbers {u1,...,un} so that ui approximates xi well for each
1 <= i <= n, and so the least common denominator of the ui's is small.
AffineRationalize[{x0,x1,...,xn}, prec] allows a tolerance of 10^-prec.
The numbers xi must be real."

Unprotect[ ProjectiveRationalize, AffineRationalize ]

Begin["NumberTheory`Rationalize`Private`"]

ProjectiveRationalize[x:{__Real}, prec_Integer:0] :=
    Module[{scale, res, len = Length[x]},
	scale = If[prec===0,
			1,  (* tolerance depends on Accuracy[x] *)
		    (* else *)
			scale = Max[Accuracy[x] - prec len/(len-1), 0];
			scale = 10^Round[scale]
		];
	res = Append[IdentityMatrix[len] scale, Round[x 10^Accuracy[x]]];
	res = Drop[LatticeReduce[Transpose[res]], -1];
	res = Minors[Drop[Transpose[res], -1], len-1];
	res = Reverse[Flatten[res]] (-1)^Range[len] / scale^(len-1);
	If[res[[1]] < 0, -res, res]
	]

AffineRationalize[{x__Real}, prec_Integer:0] := 
    Module[{res},
	res = ProjectiveRationalize[{N[1, Accuracy[{x}]], x}, prec];
	res = Drop[res, 1]/res[[1]];
	If[Sign[N[{x}[[1]]]] != Sign[res[[1]]], -res, res]
	]

End[]   (* NumberTheory`Rationalize`Private` *)

Protect[ ProjectiveRationalize, AffineRationalize ]

EndPackage[]   (* NumberTheory`Rationalize` *)

(*:Limitations: None known. *)


(*:Examples:

ProjectiveRationalize[ N[{647367, 475748, 978694, 536475}/Pi] ]

AffineRationalize[ N[{647367, 475748, 978694} / 536475] ]

*)

