Fractions

[Bart]

From http://mathforum.org/library/drmath/view/51886.html

Code: [Select]

//Source: http://mathforum.org/library/drmath/view/51886.html
Function FloatToFraction(val, Precision: Double): TFraction;

Sorry this code from that page seems to be copyrighted, so I removed it here.


Example:

Code: [Select]

  D := 1.0 + 307/1700;
  Prec := 0.00001;
  F1 := FloatTofraction(D, Prec);
  writeln('F1      = ',F1.ToString);
  writeln('D       = ',D:10:10);
  writeln('ToFloat = ',(F1.Numerator/F1.Denominator):10:10);
  writeln('Diff    = ',Abs(D-(F1.Numerator/F1.Denominator)):10:10);
  writeln('Prec    = ',Prec:10:10);

Ouputs

Code: [Select]

F1      = 438/371
D       = 1.1805882353
ToFloat = 1.1805929919
Diff    = 0.0000047566
Prec    = 0.0000100000

Nice.

Bart

[circular]

You could also do using Bart fractions like that:

Code: [Select]

const fracDenom = 8*9*5*7*11;
var
  value: double;
  frac: TFraction;
begin
  value := 3.141592653589793238;
  frac.Init(round(value*fracDenom),fracDenom);
  frac.Normalize;
  writeln(frac.ToString);
  writeln(frac.Numerator/frac.Denominator);
  readln;
end;   That gives you Pi = 17417/5544.

However you have not much control of the number of digits here. You could use the continued fraction and Bart fractions like that:

Code: [Select]

function ContinuedFraction(AValue: Double; AMaxNumerator, AMaxDenominator: int64): TFraction;
const maxDigits = 99999999999999;
var
  backup: TFraction;
  depth: integer;
  epsilon: double;
begin
  if AMaxNumerator > maxDigits then AMaxNumerator := maxDigits;
  if AMaxDenominator > maxDigits then AMaxDenominator := maxDigits;
  depth := 1;
  epsilon:= 1/AMaxDenominator;
  backup := ContinuedFractionWithDepth(AValue, depth, epsilon);
  while depth < 50 do
  begin
    inc(depth);
    result := ContinuedFractionWithDepth(AValue, depth, epsilon);
    if (abs(result.Numerator) > AMaxNumerator) or
      (result.Denominator > AMaxDenominator) then
    begin
      result := backup;
      exit;
    end;
    backup := result;
  end;
end;

function ContinuedFractionWithDepth(AValue: Double; AMaxDepth: integer; AEpsilon: Double = 1e-15): TFraction;
var fracPart: Double;
  unitFraction: TFraction;
  intPart: int64;
begin
  if AValue < 0 then
  begin
    result := ContinuedFractionWithDepth(-AValue, AMaxDepth, AEpsilon);
    result.Numerator := -result.Numerator;
  end
  else
  begin
    intPart := Trunc(AValue);
    result.Init(intPart,1);
    if AMaxDepth > 1 then
    begin
      unitFraction.Init(1,1);
      fracPart := AValue-intPart;
      if fracPart > AEpsilon then
        result += unitFraction / ContinuedFractionWithDepth(1/fracPart, AMaxDepth-1, AEpsilon);
    end;
  end;
end;

const maxNb = 999;
var
  value: double;
  frac: TFraction;
begin
  value := 3.141592653589793238;
  frac:= ContinuedFraction(value, maxNb,maxNb);
  writeln(frac.ToString);
  writeln(frac.Numerator/frac.Denominator);
  readln;
end;

With maxNb = 99, this gives you Pi = 22/7
with maxNb = 999, this gives you Pi = 355/113

Note that is not optimum because 355/113 is a better approximation.

[wp]

This is what I found (looks like the algorithm Bart uses), yields pi as 355/113. Still some work required for controlling the number of digits - for this circular's method is best (and very elegant!).

Code: [Select]

procedure FracApprox(AValue: Double; out ANumerator, ADenominator: Integer);
// "Stern-Brocot-Tree"
// http://stackoverflow.com/questions/5124743/algorithm-for-simplifying-decimal-to-fractions
const
  EPS = 0.000001;
var
  n: Integer;
  lower_n, lower_d, upper_n, upper_d, middle_n, middle_d: Integer;
  isNeg: Boolean;
begin
  isNeg := AValue < 0;
  if isNeg then
    AValue := -AValue;

  n := Floor(AValue);
  AValue := AValue - n;

  if AValue < EPS then
  begin
    ANumerator := n;
    ADenominator := 1;
  end;

  if (1 - EPS < AValue) then
  begin
    ANumerator := n+1;
    ADenominator := 1;
  end;

  // Lower fraction is 0/1
  lower_n := 0;
  lower_d := 1;

  // Upper fraction is 1/1
  upper_n := 1;
  upper_d := 1;

  while true do
  begin
    // middle fraction is (lower_n + upper_n) / (lower_d + upper_d)
    middle_n := lower_n + upper_n;
    middle_d := lower_d + upper_d;

    // AValue + EPS < middle
    if middle_d * (AValue + EPS) < middle_n then
    begin
      // middle is our new upper
      upper_n := middle_n;
      upper_d := middle_d;
    end else
    // middle < AValue - EPS
    if middle_n < (AValue - EPS) * middle_d then
    begin
      // moddle is our new lower
      lower_n := middle_n;
      lower_d := middle_d;
    end else
    // middle is our best fraction
    begin
      ANumerator := n * middle_d + middle_n;
      if isNeg then ANumerator := - ANumerator;
      ADenominator := middle_d;
      exit;
    end;
  end;
end;

[Bart]

Yet another implementation:

Code: [Select]

//uses method of Continued fractions
function FloatToFrac(F: Double; Precision: Double): TFraction;
var
  H1, H2, K1, K2, A, tmp: Int64;
  B, diff, test: Double;
begin
  H1 := 1;
  H2 := 0;
  K1 := 0;
  K2 := 1;
  b := F;
  repeat
    A := Round(Floor(b));
    tmp := H1;
    H1 := (a * H1) + H2;
    H2 := tmp;
    tmp := K1;
    K1 := (a * K1) + K2;
    K2 := tmp;
    B := 1 / (B - A);
    test := H1 / K1;
    diff := Abs(test - F);
  until (diff < Precision);
  Result.Numerator := H1;
  Result.Denominator := K1;
end;

Tested with Pi and 0.0000001 precsion:
Gives a slightly different result than the first oe I posted:

Code: [Select]

C:\Users\Bart\LazarusProjecten\ConsoleProjecten\fractions>fractest
FloatToFraction:
F1      = 104348/33215
D       = 3.1415926536
ToFloat = 3.1415926539
Diff    = 0.0000000003
Prec    = 0.0000001000
FloatToFrac:
F1      = 103993/33102
D       = 3.1415926536
ToFloat = 3.1415926530
Diff    = 0.0000000006
Prec    = 0.0000001000

Bart

[Mike.Cornflake]

Definitely one of those times that I wish this forum had "like" functionality!

You guys rock :-)

[typo]

Why not a LazUtils' FloatToFraction function?

[wp]

@circular: Your function hangs for the value 0.0317 (Probably the check "if fracPart > 0 then"  should be modified to something like "if fracPart > eps then" (eps = 1e-9).

@bart: I modified your last function such that it could be added to fpspreadsheet (add max denominator digits, remove dependency on your fractions unit because fpspreadsheet will perform calculations based on the raw values, not the approximated fractions). What I've seen so far, it returns the fractions that Excel produces:

Code: Text  [Select][+][-]

1. procedure FloatToFrac(F: Double; Precision: Double; AMaxDenomDigits: Integer;
2.   out Numerator, Denominator: Integer);
3. var
4.   H1, H2, K1, K2, A, tmp: Integer;
5.   B, diff, test: Double;
6.   prev_H1, prev_K1: Integer;
7. begin
8.   H1 := 1;
9.   H2 := 0;
10.   K1 := 0;
11.   K2 := 1;
12.   b := F;
13.   prev_H1 := H1;
14.   prev_K1 := K1;
15.   repeat
16.     A := Round(Floor(b));
17.     tmp := H1;
18.     H1 := (a * H1) + H2;
19.     H2 := tmp;
20.     tmp := K1;
21.     K1 := (a * K1) + K2;
22.     K2 := tmp;
23.     B := 1 / (B - A);
24.     test := H1 / K1;
25.     diff := Abs(test - F);
26.     if length(IntToStr(K1)) > AMaxDenomDigits then
27.     begin
28.       H1 := prev_H1;
29.       K1 := prev_K1;
30.       break;
31.     end;
32.     prev_H1 := H1;
33.     prev_K1 := K1;
34.   until (diff < Precision);
35.   Numerator := H1;
36.   Denominator := K1;
37. end;
38.  

[circular]

@circular: Your function hangs for the value 0.0317 (Probably the check "if fracPart > 0 then"  should be modified to something like "if fracPart > eps then" (eps = 1e-9).

Hmm, it does not hang on my computer. But yes, I imagine some epsilon would help.

How about this so that you can choose the exact maximum for numerator and denominator, and the precision is computed for you:

function FracApprox(AValue: Double; AMaxNumerator, AMaxDenominator: Int64): TFraction;
// "Stern-Brocot-Tree"
// http://stackoverflow.com/questions/5124743/algorithm-for-simplifying-decimal-to-fractions
var
  n: Integer;
  lower_n, lower_d, upper_n, upper_d, middle_n, middle_d: Int64;
  isNeg: Boolean;
  backup, newResult: TFraction;
  EPS: Double;
begin
  EPS := 0.01/AMaxDenominator;

  isNeg := AValue < 0;
  if isNeg then
    AValue := -AValue;

  n := Trunc(AValue);
  newResult.Init(round(AValue),1);
  if isNeg then newResult.Numerator := -newResult.Numerator;
  backup := newResult;

  AValue := AValue - n;

  // Lower fraction is 0/1
  lower_n := 0;
  lower_d := 1;

  // Upper fraction is 1/1
  upper_n := 1;
  upper_d := 1;

  while true do
  begin
    if abs(newResult.Numerator/newResult.Denominator - n - AValue) <
       abs(backup.Numerator/backup.Denominator - n - AValue) then
      backup := newResult;

    // middle fraction is (lower_n + upper_n) / (lower_d + upper_d)
    middle_n := lower_n + upper_n;
    middle_d := lower_d + upper_d;
    newResult.Init(n * middle_d + middle_n, middle_d);
    newResult.Normalize;
    if (newResult.Numerator > AMaxNumerator) or
     (newResult.Denominator > AMaxDenominator) then
    begin
      result := backup;
      exit;
    end;

    if isNeg then newResult.Numerator := -newResult.Numerator;

    // AValue + EPS < middle
    if middle_d * (AValue + EPS) < middle_n then
    begin
      // middle is our new upper
      upper_n := middle_n;
      upper_d := middle_d;
    end else
    // middle < AValue - EPS
    if middle_n < (AValue - EPS) * middle_d then
    begin
      // moddle is our new lower
      lower_n := middle_n;
      lower_d := middle_d;
    end else
    // middle is our best fraction
    begin
      result := newResult;
      exit;
    end;
  end;
end;

22/7
 3.1428571428571429E+0000
355/113
 3.1415929203539823E+0000
355/113
 3.1415929203539823E+0000
75948/24175
 3.1415925542916236E+0000
100798/32085
 3.1415926445379461E+0000

[wp]

Nice

[circular]

I guess I understand the bug in my ContinuedFraction function. The thing is that using a recursive function is not really a good idea because it makes it complicated to compare successive results. So I have changed the code, so that each depth is tried. That's not the most optimised, but it works. Now Pi gives 355/113 for the format 999/999.

[typo]

My compiler stops on the line:

Code: [Select]

newResult.Init(round(AValue),1);

It does not find Init member of TFraction class or record.

[howardpc]

You added

Code: [Select]

{$modeswitch advancedrecords}
to the unit?

[typo]

I cannot find any declaration of TFraction which would need this.

[typo]

A small improvement to Canonicalize:

Code: [Select]

Function Canonicalize(q1:Fraction):Fraction;
var
    GCD :integer;
    q2:Fraction;
begin
    GCD := GreatestCommonDivisor(q1.num, q1.den);
    q2.num := q1.num div GCD;
    q2.den := q1.den div GCD;
    Result:=q2;

    if ((Result.num < 0) and (Result.den < 0))
    or ((Result.den < 0) and (Result.num >= 0)) then
    begin
      Result.num := - Result.num;
      Result.den := - Result.den;
    end;
end;       


[HappyLarry]

@typo
Your improvement to Canonicalize isn't needed - all negatives (either in numerator or denominatorar or both) are taken care of in CreateFraction (which has to be called before Canonicalize)

Code: [Select]

function CreateFraction(n, d: Int64): Fraction;
var
     q:Fraction;
begin
      if (d = 0) then raise EZeroDivide.Create('denominator cannot be 0');
      if (d < 0) then
      begin
           q.num := n*-1;
           q.den := d*-1;
      end
      else
      begin
           q.num := n;
           q.den := d;
      end;
      Result := q;
end;
Besides which, in your if statement

    if ((Result.num < 0) and (Result.den < 0))
    or ((Result.den < 0) and (Result.num >= 0)) then

is exactly the same as
If Result.den < 0.