MathgeomGLS

1
{---------------------------------------------------------------------------}
2
{                                                                           }
3
{ File:       Velthuis.BigRationals.pas                                     }
4
{ Function:   A multiple precision rational nubmer implementation, based    }
5
{             on the BigInteger implementation in                           }
6
{             Velthuis.BigIntegers.pas.                                     }
7
{ Language:   Delphi version XE2 or later                                   }
8
{ Author:     Rudy Velthuis                                                 }
9
{ Copyright:  (c) 2016,2017 Rudy Velthuis                                   }
10
{                                                                           }
11
{ License:    Redistribution and use in source and binary forms, with or    }
12
{             without modification, are permitted provided that the         }
13
{             following conditions are met:                                 }
14
{                                                                           }
15
{             * Redistributions of source code must retain the above        }
16
{               copyright notice, this list of conditions and the following }
17
{               disclaimer.                                                 }
18
{             * Redistributions in binary form must reproduce the above     }
19
{               copyright notice, this list of conditions and the following }
20
{               disclaimer in the documentation and/or other materials      }
21
{               provided with the distribution.                             }
22
{                                                                           }
23
{ Disclaimer: THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDER "AS IS"     }
24
{             AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT     }
25
{             LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND     }
26
{             FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO        }
27
{             EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE     }
28
{             FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY,     }
29
{             OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,      }
30
{             PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,     }
31
{             DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED    }
32
{             AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT   }
33
{             LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)        }
34
{             ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF   }
35
{             ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.                    }
36
{                                                                           }
37
{---------------------------------------------------------------------------}
38

39
unit Velthuis.BigRationals experimental;
40

41
//$$RV comments are internal comments by RV, denoting a problem. Should be all removed as soon as problem solved.
42
// TODO: handle extended properly, e.g. in Create(Extended)
43

44
// To implement:
45
// Create(Value: Double; MaxDenominator: Integer); -- see below
46
// AsSingle
47
// AsDouble
48
// AsExtended
49
// AsBigDecimal()
50
// AsBigDecimal(RoundingMode)
51
// AsBigDecimal(Scale, RoundingMode)
52
// AsInteger
53
// AsInt64
54
// Add(BigRational, BigRational)
55
// Add(BigRational, Integer)
56
// Add(BigRational, Int64)
57
// Add(BigRational, BigInteger)
58
// same for Subtract, Multiply and Divide
59
// Negate
60
// PercentageValue ?
61
// Pow(Exponent: Integer)
62
// Pow(Exponent: Int64);
63
// Pow(Exponent: BigInteger);
64
// Pow(Exponent: Double): Double;
65
// My own:
66
// ToMixedString (returns "97/15" as "6 7/15")
67
// IsProper (proper fraction is when num < denom)
68

69
interface
70

71
uses
72
  Velthuis.BigIntegers, Velthuis.BigDecimals, System.Math, System.SysUtils, CompilerAndRTLVersions;
73

74
{$IF CompilerVersion < CompilerVersionDelphiXE8}
75
  {$IF (DEFINED(WIN32) OR DEFINED(CPUX86)) AND NOT DEFINED(CPU32BITS)}
76
    {$DEFINE CPU32BITS}
77
  {$IFEND}
78
  {$IF (DEFINED(WIN64) OR DEFINED(CPUX64)) AND NOT DEFINED(CPU64BITS)}
79
    {$DEFINE CPU64BITS}
80
  {$IFEND}
81
{$IFEND}
82

83
{$IF CompilerVersion >= CompilerVersionDelphiXE3}
84
  {$LEGACYIFEND ON}
85
{$IFEND}
86

87
{$IF CompilerVersion >= CompilerVersionDelphiXE}
88
  {$CODEALIGN 16}
89
  {$ALIGN 16}
90
{$IFEND}
91

92
{$IF CompilerVersion >= CompilerVersionDelphi2010}
93
  {$DEFINE HasClassConstructors}
94
{$IFEND}
95

96
{$IF SizeOf(Extended) > SizeOf(Double)}
97
  {$DEFINE HasExtended}
98
{$IFEND}
99

100
type
101
  PBigRational = ^BigRational;
102

103
  /// <summary>BigRational is a multiple precision rational data type, where each value is expressed as the quotient
104
  ///   of two BigIntegers, the numerator and the denominator.</summary>
105
  /// <remarks><para>BigRationals are simplified, i.e. common factors in numerator and denominator are eliminated.
106
  ///   </para>
107
  ///   <para>The resulting sign is always moved to the numerator. The denominator must always be unsigned.</para>
108
  /// </remarks>
109
  BigRational = record
110
  private
111
    type
112
      // Error code for the Error procedure.
113
      TErrorCode = (ecParse, ecDivByZero, ecConversion, ecInvalidArg, ecZeroDenominator);
114

115
    var
116
      // The numerator (the "top" part of the fraction).
117
      FNumerator: BigInteger;
118
      // The denominator (the "bottom" part of the fraction).
119
      FDenominator: BigInteger;
120

121
    class var
122
      // Field for AlwaysReduce property.
123
      FAlwaysReduce: Boolean;
124

125
    // Returns -1 for negative values, 1 for positive values and 0 for zero.
126
    function GetSign: Integer;
127

128
    // Checks for invalid values, simplifies the fraction by eliminating common factors in numerator and denominator,
129
    // and moves sign to numerator. In other words, turns a BigRational into its canonical form.
130
    procedure Normalize(Forced: Boolean = False);
131

132
    // Raises exception using error code and additional data to decide which and how.
133
    procedure Error(Code: TErrorCode; ErrorInfo: array of const);
134
  public
135
    const
136
      MinEpsilon = 5e-10;
137

138
    class var
139
      /// <summary>Predefined BigRational value 0.</summary>
140
      Zero: BigRational;
141

142
      /// <summary>Predefined BigRational value 1/10.</summary>
143
      OneTenth: BigRational;
144

145
      /// <summary>Predefined BigRational value 1/4.</summary>
146
      OneFourth: BigRational;
147

148
      /// <summary>Predefined BigRational value 1/3.</summary>
149
      OneThird: BigRational;
150

151
      /// <summary>Predefined BigRational value 1/2.</summary>
152
      OneHalf: BigRational;
153

154
      /// <summary>Predefined BigRational value 2/3.</summary>
155
      TwoThirds: BigRational;
156

157
      /// <summary>Predefined BigRational value 3/4.</summary>
158
      ThreeFourths: BigRational;
159

160
      /// <summary>Predefined BigRational value 1.</summary>
161
      One: BigRational;
162

163
      /// <summary>Predefined BigRational value 10.</summary>
164
      Ten: BigRational;
165

166
{$IFDEF HasClassConstructors}
167
    class constructor Initialize;
168
{$ENDIF}
169

170
    // -- Constructors --
171

172
    /// <summary>Creates a new BigRational with the given values as numerator and denominator, respectively. The
173
    ///  sign is adjusted thus, that the denominator is positive.</summary>
174
    constructor Create(const Numerator, Denominator: BigInteger); overload;
175

176
    /// <summary>Creates a new BigRational from the given value, with a denominator of 1.</summary>
177
    constructor Create(const Value: BigInteger); overload;
178

179
    /// <summary>Creates a new BigRational from the given value, with a denominator of 1.</summary>
180
    constructor Create(const Value: Integer); overload;
181

182
    /// <summary>Creates a new BigRational from the given values.</summary>
183
    constructor Create(const Numerator: Integer; const Denominator: Cardinal); overload;
184

185
    /// <summary>Creates a new BigRational with the exact same value as the given Double value.</summary>
186
    /// <exception cref="EInvalidArgument">EInvalidArgument is raised if the Double represents a positive or
187
    ///   negative infinity or a NaN.</exception>
188
    /// <remarks>Note that this is an exact conversion, taking into account the exact bit representation of the
189
    ///   double. This means that the numerator and denominator values can be rather big.</remarks>
190
    constructor Create(const Value: Double); overload;
191

192
    /// <summary>Creates a new BigRational with the same value as the given Double value. For this, it uses
193
    ///  at most MaxIterations and the error is at most Epsilon.</summary>
194
    /// <exception cref="EInvalidArgument">EInvalidArgument is raised if the Double represents a positive or
195
    ///  negative infinity or a NaN.</exception>
196
    /// <exception cref="EOverflow">EOverflow is raised when the Double represents a value that is not representable
197
    ///  as a fraction of Integers.</exception>
198
    /// <remarks><para>Note that this conversion creates enumerator and denominator that are at most Integer
199
    ///  values.</para>
200
    /// <para>This uses a conversion to a finite continued fraction combined with the unfolding to a simple fraction
201
    ///  in one loop, steadily refining the fraction. MaxIterations and Epsilon determine when the loop ends.
202
    //   MaxIterations governs the maximum number of iterations in the loop, and Epsilon governs the maximum difference
203
    //   between the double and the fraction.</para>
204
    /// <para>Typical values are MaxIterations = 15 and Epsilon = 4e-10 (which is close to 1/MaxInt).</para></remarks>
205
    constructor Create(Value, Epsilon: Double; MaxIterations: Integer); overload;
206

207
    /// <summary>Creates a new BigRational with the best matching value, with a denominator value of at most
208
    ///  MaxDenominator</summary>
209
    constructor Create(F: Double; MaxDenominator: Cardinal); overload;
210

211
    //{$IFDEF HasExtended}
212
//    /// <summary>Creates a new BigRational with the exact same value as the given Extended value.</summary>
213
//    /// <exception cref="EInvalidArgument">EInvalidArgument is raised if the Extended represents a positive or
214
//    ///   negative infinity or a NaN.</exception>
215
//    /// <remarks>Note that this is an exact conversion, taking into account the exact bit representation of the
216
//    ///   Extended. This means that the numerator and denominator values can be rather big.</remarks>
217
//    constructor Create(const Value: Extended); overload;
218
//{$ENDIF}
219

220
    /// <summary>Creates a new BigRational from the given value, with a denominator of 1.</summary>
221
    constructor Create(const Value: Int64); overload;
222

223
    /// <summary>Creates a new BigRational from the given values.</summary>
224
    constructor Create(const Numerator: Int64; const Denominator: UInt64); overload;
225

226
    /// <summary>Creates a new BigRational from the given string, which should be in the format
227
    ///   [optional sign] + numerator + '/' + denominator.</summary>
228
    /// <remarks>Example input: '-1234/5678'.</remarks>
229
    constructor Create(const Value: string); overload;
230

231
    /// <summary>Creates a new BigRational from the given BigDecimal.</summary>
232
    constructor Create(const Value: BigDecimal); overload;
233

234
    // TODO: Create(const Value: Double; MaxDenominator: Integer); overload; // gets best approximation
235
    // TODO: Create(const Value: Double; MaxEpsilon: Double); overload; // gets best approximation
236

237

238
    // -- Mathematical operators and functions --
239

240
    /// <summary>Adds two BigRationals and returns the sum. Simplifies the result.</summary>
241
    class function Add(const Left, Right: BigRational): BigRational; static;
242

243
    /// <summary>Adds two BigRationals and returns the sum. Simplifies the result.</summary>
244
    class operator Add(const Left, Right: BigRational): BigRational;
245

246
    /// <summary>Subtracts two BigRationals and returns the difference. Simplifies the result.</summary>
247
    class function Subtract(const Left, Right: BigRational): BigRational; static;
248

249
    /// <summary>Subtracts two BigRationals and returns the difference. Simplifies the result.</summary>
250
    class operator Subtract(const Left, Right: BigRational): BigRational;
251

252
    /// <summary>Multiplies two BigRationals and returns the product. Simplifies the result.</summary>
253
    class function Multiply(const Left, Right: BigRational): BigRational; static;
254

255
    /// <summary>Multiplies two BigRationals and returns the product. Simplifies the result.</summary>
256
    class operator Multiply(const Left, Right: BigRational): BigRational;
257

258
    /// <summary>Divides two BigRationals by multiplying Left by the reciprocal of Right. Returns the quotient.
259
    ///   Simplifies the result.</summary>
260
    class function Divide(const Left, Right: BigRational): BigRational; static;
261

262
    /// <summary>Divides two BigRationals by multiplying Left by the reciprocal of Right. Returns the quotient.
263
    ///   Simplifies the result.</summary>
264
    class operator Divide(const Left, Right: BigRational): BigRational;
265

266
    /// <summary>Divides two BigRationals returning a BigInteger result.</summary>
267
    class operator IntDivide(const Left, Right: BigRational): BigInteger;
268

269
    /// <summary>Divides two BigRationals returning the remainder. Simplifies the result.</summary>
270
    class operator Modulus(const Left, Right: BigRational): BigRational;
271

272
    /// <summary>Divides two BigRationals returning the remainder. Simplifies the result.</summary>
273
    class function Remainder(const Left, Right: BigRational): BigRational; static;
274

275
    /// <summary>Divides two BigRationals and returns the quotient and remainder.</summary>
276
    class procedure DivMod(const Left, Right: BigRational; var Quotient: BigInteger; var Remainder: BigRational); static;
277

278
    /// <summary>Returns the negation of the given BigRational value.</summary>
279
    class operator Negative(const Value: BigRational): BigRational;
280

281
    // Instance functions
282

283
    function IsNegative: Boolean; inline;
284
    function IsPositive: Boolean; inline;
285
    function IsZero: Boolean; inline;
286

287
    /// <summary>Returns the absolute value of the current BigRational.</summary>
288
    function Abs: BigRational;
289

290
    /// <summary>Returns the negation of the current BigRational.</summary>
291
    function Negate: BigRational;
292

293
    /// <summary>Returns the multiplicative inverse of the current BigRational value by swapping numerator and
294
    ///   denominator.</summary>
295
    function Reciprocal: BigRational;
296

297
    /// <summary>Reduces numerator and denominator to their smallest values representing the same ratio.</summary>
298
    function Reduce: BigRational;
299

300

301
    // -- Comparison and relational operators --
302

303
    /// <summary>Compares two BigRationals. Returns -1 if Left < Right, 1 if Left > Right and
304
    ///   0 if Left = Right.</summary>
305
    class function Compare(const Left, Right: BigRational): Integer; static;
306

307
    /// <summary>Returns True only if Left < Right.</summary>
308
    class operator LessThan(const Left, Right: BigRational): Boolean;
309

310
    /// <summary>Returns True only if Left <= Right.</summary>
311
    class operator LessThanOrEqual(const Left, Right: BigRational): Boolean;
312

313
    /// <summary>Returns True only if Left = Right.</summary>
314
    class operator Equal(const Left, Right: BigRational): Boolean;
315

316
    /// <summary>Returns True only if Left >= Right.</summary>
317
    class operator GreaterThanOrEqual(const Left, Right: BigRational): Boolean;
318

319
    /// <summary>Returns True only if Left > Right.</summary>
320
    class operator GreaterThan(const Left, Right: BigRational): Boolean;
321

322
    /// <summary>Returns True only if Left <> Right.</summary>
323
    class operator NotEqual(const Left, Right: BigRational): Boolean;
324

325

326
    // --- Conversion operators --
327

328
    /// <summary>Converts the string to a BigRational.</summary>
329
    /// <exception cref="EConvertError"></exception>
330
    class operator Implicit(const Value: string): BigRational;
331

332
    /// <summary>Explicitly converts the BigRational to a string (using ToString).</summary>
333
    class operator Explicit(const Value: BigRational): string;
334

335
    /// <summary>Converts the integer to a BigRational.</summary>
336
    class operator Implicit(const Value: Integer): BigRational;
337

338
    /// <summary>Converts the BigRational to an Integer. If necessary, excess top bits are cut off.</summary>
339
    class operator Explicit(const Value: BigRational): Integer;
340

341
    /// <summary>Converts the Int64 to a BigRational.</summary>
342
    class operator Implicit(const Value: Int64): BigRational;
343

344
    /// <summary>Converts the BigRational to an Int64. If necessary, excess top bits are cut off.</summary>
345
    class operator Explicit(const Value: BigRational): Int64;
346

347
    /// <summary>Converts the given Double to a BigRational.</summary>
348
    /// <exception cref="EInvalidArgument">Raises an EInvalidArgument exception if Value represents
349
    ///   (+/-)infinity or NaN.</exception>
350
    class operator Implicit(const Value: Double): BigRational;
351

352
    /// <summary>Converts the given BigRational to a Double.</summary>
353
    class operator Explicit(const Value: BigRational): Double;
354

355
    /// <summary>Converts the given BigDecimal to a BigRational.</summary>
356
    /// <remarks>This conversion is exact.</remarks>
357
    class operator Implicit(const Value: BigDecimal): BigRational;
358

359
    /// <summary>Converts the given BigRational to a BigDecimal, using the default BigDecimal precision and
360
    ///   rounding mode.</summary>
361
    class operator Implicit(const Value: BigRational): BigDecimal;
362

363
    /// <summary>Converts the current BigRational to a string. This string can be used as valid input for
364
    ///  conversion.</summary>
365
    function ToString: string;
366

367
    function Parse(const S: string): BigRational;
368
    function TryParse(const S: string; out Value: BigRational): Boolean;
369

370

371
    // -- Properties --
372

373
    /// <summary>The numerator of the fraction represented by the current BigRational.</summary>
374
    property Numerator: BigInteger read FNumerator;
375

376
    /// <summary>The denominator of the fraction represented by the current BigRational.</summary>
377
    property Denominator: BigInteger read FDenominator;
378

379
    /// <summary>The sign of the fraction represented by the current BigRational.</summary>
380
    property Sign: Integer read GetSign;
381

382
    /// <summary>If AlwaysReduce is set (default), fractions like 2/10 are reduced to 1/5.
383
    ///   If it is not set, you can manually call Reduce on any BigRational.</summary>
384
    class property AlwaysReduce: Boolean read FAlwaysReduce write FAlwaysReduce;
385

386
  end;
387

388
implementation
389

390
uses
391
  Velthuis.FloatUtils, Velthuis.StrConsts;
392

393
{ BigRational }
394

395
function BigRational.Abs: BigRational;
396
begin
397
  if Self.FNumerator.IsNegative then
398
    Result := Self.Negate
399
  else
400
    Result := Self;
401
end;
402

403
class operator BigRational.Add(const Left, Right: BigRational): BigRational;
404
begin
405
  if Left.FDenominator = Right.FDenominator then
406
  begin
407
    Result.FNumerator := Left.FNumerator + Right.FNumerator;
408
    Result.FDenominator := Left.FDenominator;
409
  end
410
  else
411
  begin
412
    Result.FDenominator := Left.FDenominator * Right.FDenominator;
413
    Result.FNumerator := Left.FNumerator * Right.FDenominator + Right.FNumerator * Left.FDenominator;
414
  end;
415
  Result.Normalize;
416
end;
417

418
class function BigRational.Add(const Left, Right: BigRational): BigRational;
419
begin
420
  Result := Left + Right;
421
end;
422

423
constructor BigRational.Create(const Value: Integer);
424
begin
425
  FNumerator := BigInteger(Value);
426
  FDenominator := BigInteger.One;
427
end;
428

429
constructor BigRational.Create(const Numerator: Integer; const Denominator: Cardinal);
430
begin
431
  FNumerator := BigInteger(Numerator);
432
  FDenominator := BigInteger(Denominator);
433
  Normalize;
434
end;
435

436
constructor BigRational.Create(const Numerator, Denominator: BigInteger);
437
begin
438
  FNumerator := Numerator;
439
  FDenominator := Denominator;
440
  Normalize;
441
end;
442

443
constructor BigRational.Create(const Value: BigInteger);
444
begin
445
  FNumerator := Value;
446
  FDenominator := BigInteger.One;
447
end;
448

449
(*
450
   Test with e.g.
451

452
   0.750000000   1 / 1   3 / 4   3 / 4   3 / 4   3 / 4
453
   0.518518000   1 / 1   1 / 2   14 / 27   14 / 27   14 / 27
454
   0.905405400   1 / 1   9 / 10   67 / 74   67 / 74   67 / 74
455
   0.142857143   0 / 1   1 / 7   1 / 7   1 / 7   1 / 7
456
   3.141592654   3 / 1   22 / 7   22 / 7   355 / 113   355 / 113
457
   2.718281828   3 / 1   19 / 7   193 / 71   1457 / 536   25946 / 9545
458
  -0.423310825   0 / 1  -3 / 7  -11 / 26  -69 / 163  -1253 / 2960
459
  31.415926536   31 / 1   157 / 5   377 / 12   3550 / 113   208696 / 6643
460
   0.000000000
461
*)
462

463
constructor BigRational.Create(F: Double; MaxDenominator: Cardinal);
464

465
// https://rosettacode.org/wiki/Convert_decimal_number_to_rational#C
466

467
var
468
  A: Int64;
469
  H, K: array[0..2] of Int64;
470
  X, D, N: Int64;
471
  I: Integer;
472
  LNegative: Boolean;
473
  MustBreak: Boolean;
474
begin
475
  H[0] := 0; H[1] := 1; H[2] := 0;
476
  K[0] := 1; K[1] := 0; K[2] := 0;
477
  N := 1;
478

479
  if MaxDenominator <= 1 then
480
  begin
481
    FDenominator := 1;
482
    FNumerator := Trunc(F);
483
    Exit;
484
  end;
485

486
  LNegative := F < 0;
487
  if LNegative then
488
    F := -F;
489

490
  while (F <> System.Trunc(F)) and (N < (High(Int64) shr 1)) do
491
  begin
492
    N := N shl 1;
493
    F := F * 2;
494
  end;
495
  D := System.Trunc(F);
496

497
  // Continued fraction and check denominator each step
498
  for I := 0 to 63 do
499
  begin
500
    MustBreak := False;
501
    if N <> 0 then
502
      A := D div N
503
    else
504
      A := 0;
505

506
    if (I <> 0) and (A = 0) then
507
      Break;
508

509
    if N = 0 then
510
      Break;
511

512
    X := D;
513
    D := N;
514
    N := X mod N;
515

516
    X := A;
517
    if K[1] * A + K[0] >= MaxDenominator then
518
    begin
519
      X := (MaxDenominator - K[0]) div K[1];
520
      if (X * 2 >= A) or (K[1] >= MaxDenominator) then
521
        MustBreak := True
522
      else
523
        Break;
524
    end;
525
    H[2] := X * H[1] + H[0]; H[0] := H[1]; H[1] := H[2];
526
    K[2] := X * K[1] + K[0]; K[0] := K[1]; K[1] := K[2];
527
    if MustBreak then
528
      Break;
529
  end;
530
  FDenominator := K[1];
531
  if LNegative then
532
    FNumerator := -H[1]
533
  else
534
    FNumerator := H[1];
535
end;
536

537
constructor BigRational.Create(const Value: Double);
538
var
539
  Exponent: Integer;
540
  Mantissa: Int64;
541
begin
542
  if IsInfinite(Value) or IsNaN(Value) then
543
    Error(ecInvalidArg, ['Double']);
544

545
  if Value = 0.0 then
546
  begin
547
    FNumerator := BigInteger.Zero;
548
    FDenominator := BigInteger.One;
549
    Exit;
550
  end;
551

552
  Exponent := GetExponent(Value);
553
  Mantissa := GetMantissa(Value);
554

555
  if IsDenormal(Value) then
556
  begin
557
    FNumerator := Mantissa;
558
    FDenominator := BigInteger.Zero.FlipBit(1023 + 52);
559
  end
560
  else
561
  begin
562
    FDenominator := BigInteger.Zero.FlipBit(52);
563
    FNumerator := Mantissa;
564
    if Exponent < 0 then
565
      FDenominator := FDenominator shl -Exponent
566
    else
567
      FNumerator := FNumerator shl Exponent;
568
  end;
569
  if Value < 0 then
570
    FNumerator := -FNumerator;
571
  Normalize;
572
end;
573

574
constructor BigRational.Create(Value, Epsilon: Double; MaxIterations: Integer);
575

576
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
577
//                                                                                                                  //
578
//  One way to solve this is to convert the Value into a (finite) continous fraction, i.e. a fraction in the form   //
579
//  a0 + 1/(a1 + 1/(a2 + 1/(a3 + ...))) or, often also written as [a0; a1, a2, a3 ...].                             //
580
//  If we reverse this using numerator and denominator, we can easily construct a suitable normal fraction.         //
581
//                                                                                                                  //
582
//  This constructor does something similar, except that the loop to construct a continued fraction and the loop    //
583
//  to decompose the continued fraction into a simple fraction are folded together. This repeatedly refines the     //
584
//  fraction, until the maximum number of iterations is reached, or the difference between fraction and Value is    //
585
//  less than the given epsilon, whichever comes first.                                                             //
586
//                                                                                                                  //
587
//  This uses numerators and denominators as much as possible in the Integer range. A tested version using          //
588
//  BigDecimals and BigIntegers was far too precise and did not stop early enough.                                  //
589
//                                                                                                                  //
590
//  There is a MinEpsilon constant of value 5e-10, which is slightly above 1/MaxInt. Smaller values could cause     //
591
//  invalid  results.                                                                                               //
592
//                                                                                                                  //
593
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
594

595
(*
596
function ConvertDoubleToRational(Value, Epsilon: Double; MaxIterations: Integer; var FNumerator, FDenominator: Int64): Double; overload;
597
*)
598
var
599
  LNegative: Boolean;
600
  LQuot: Double;
601
  LError: Double;
602
  I: Integer;
603
  NewRatio, Ratio: Double;
604
  Rest, Rest0, Inverse, Inverse0: Double;
605
  K, H: array[0..2] of Int64;
606
  A: Int64;
607
  LNum, LDenom: BigInteger;
608
begin
609
  LQuot := 0.0;
610

611
  if IsInfinite(Value) or IsNaN(Value) then
612
    Error(ecInvalidArg, ['Double']);
613

614
  LNegative := Value < 0;
615
  if LNegative then
616
    Value := -Value;
617

618
  if Value > MaxInt then
619
  begin
620
    // reduce Value to range [0..MaxInt)
621
    LQuot := Int(Value / MaxInt) * MaxInt;
622
    Value := Value - LQuot;
623
    // Hmmm... turn it into a suitable range
624
    // FloatMod(Value, MaxInt, Quot, Rem);
625
    // Quot := Int(Value / MaxInt);
626
    // Value := Value - Quot * MaxInt;
627
    // work on rem --> value, remember quot;
628
    // raise EOverflow.CreateFmt('Value %g cannot be converted to an integer ratio', [Value]);
629
  end;
630

631
  if Value = 0.0 then
632
  begin
633
    FNumerator := BigInteger.Zero;
634
    FDenominator := BigInteger.One;
635
    Exit;
636
  end;
637

638
  if Epsilon = 0 then
639
    Epsilon := 5e-10;
640

641
  I := 0;
642
  while I < MaxIterations do
643
  begin
644
    if I = 0 then
645
    begin
646
      A := Trunc(Value);
647
      Rest := Value - A;
648
      if Rest = 0.0 then
649
      begin
650
        FNumerator := A;
651
        FDenominator := 1;
652
        Exit;
653
      end;
654
      Inverse := 1.0 / Rest;
655
      K[2] := A;
656
      H[2] := 1;
657
      NewRatio := K[2] / H[2];
658
      LError := System.Abs(Value - NewRatio);
659
    end
660
    else
661
    begin
662
      A := Trunc(Inverse0);
663
      Rest := Inverse0 - A;
664
      if Rest = 0.0 then
665
        Rest := Epsilon;
666
      Inverse := 1.0 / Rest;
667
      if I = 1 then
668
      begin
669
        K[2] := A * K[1] + 1;
670
        H[2] := A * H[1];
671
      end
672
      else
673
      begin
674
        if (A > MaxInt) or (K[1] > MaxInt) or (H[1] > MaxInt) then
675
          Break;
676
        K[2] := A * K[1] + K[0];
677
        H[2] := A * H[1] + H[0];
678
      end;
679
      NewRatio := K[2] / H[2];
680
      LError := System.Abs(NewRatio - Ratio);
681
    end;
682
    if LError < Epsilon then
683
      Break
684
    else
685
      Ratio := NewRatio;
686
    Inc(I);
687
    K[0] := K[1];
688
    K[1] := K[2];
689
    H[0] := H[1];
690
    H[1] := H[2];
691
    Inverse0 := Inverse;
692
    Rest0 := Rest;
693
    if Inverse0 > MaxInt then
694
      Break;
695
  end;
696
  LNum := K[1];
697
  LDenom := H[1];
698
  if LQuot > 0.0 then
699
    LNum := LNum + BigInteger(LQuot) * LDenom;
700
  if LNegative then
701
    LNum := -LNum;
702

703
  FNumerator := LNum;
704
  FDenominator := LDenom;
705
  Normalize;
706
end;
707

708
class function BigRational.Compare(const Left, Right: BigRational): Integer;
709
begin
710
  if Left.FDenominator = Right.FDenominator then
711
    Result := BigInteger.Compare(Left.FNumerator, Right.FNumerator)
712
  else
713
    Result := BigInteger.Compare(Left.FNumerator * Right.FDenominator,
714
                                 Right.FNumerator * Left.FDenominator);
715
end;
716

717
constructor BigRational.Create(const Numerator: Int64; const Denominator: UInt64);
718
begin
719
  FNumerator := BigInteger(Numerator);
720
  FDenominator := BigInteger(Denominator);
721
  Normalize;
722
end;
723

724
constructor BigRational.Create(const Value: Int64);
725
begin
726
  FNumerator := BigInteger(Value);
727
  FDenominator := BigInteger.One;
728
end;
729

730
constructor BigRational.Create(const Value: string);
731
begin
732
  Self := Parse(Value);
733
end;
734

735
constructor BigRational.Create(const Value: BigDecimal);
736
var
737
  Num, Denom: BigInteger;
738
  Scale: Integer;
739
begin
740
  Num := Value.UnscaledValue;
741
  Scale := Value.Scale;
742
  Denom := BigInteger.One;
743

744
  if Scale < 0 then
745
    Num := Num * BigInteger.Pow(BigInteger.Ten, -Scale)
746
  else if Scale > 0 then
747
    Denom := BigInteger.Pow(BigInteger.Ten, Scale);
748

749
  Create(Num, Denom);
750
end;
751

752
//{$IFDEF HasExtended}
753
//constructor BigRational.Create(const Value: Extended);
754
//var
755
//  D: Double;
756
//begin
757
//  D := Value;
758
//  Create(D);
759
//end;
760
//{$ENDIF HasExtended}
761

762
class function BigRational.Divide(const Left, Right: BigRational): BigRational;
763
begin
764
  Result := Left / Right;
765
end;
766

767
class operator BigRational.Divide(const Left, Right: BigRational): BigRational;
768
begin
769
  if Left.FDenominator = Right.FDenominator then
770
  begin
771
    Result.FNumerator := Left.FNumerator;
772
    Result.FDenominator := Right.FNumerator;
773
  end
774
  else
775
  begin
776
    Result.FNumerator := Left.FNumerator * Right.FDenominator;
777
    Result.FDenominator := Left.FDenominator * Right.FNumerator;
778
  end;
779
  Result.Normalize;
780
end;
781

782
class procedure BigRational.DivMod(const Left, Right: BigRational;
783
  var Quotient: BigInteger; var Remainder: BigRational);
784
var
785
  AD, BC: BigInteger;
786
begin
787
  if Left.FDenominator = Right.FDenominator then
788
  begin
789
    AD := Left.FNumerator;
790
    BC := Right.FNumerator;
791
  end
792
  else
793
  begin
794
    AD := Left.FNumerator * Right.FDenominator;
795
    BC := Left.FDenominator * Right.FNumerator;
796
  end;
797

798
  Quotient := AD div BC;
799

800
  Remainder.FNumerator := AD - Quotient * BC;
801
  Remainder.FDenominator := Left.FDenominator * Right.FDenominator;
802
  Remainder.Normalize;
803
end;
804

805
class operator BigRational.Equal(const Left, Right: BigRational): Boolean;
806
begin
807
  Result := Compare(Left, Right) = 0;
808
end;
809

810
procedure BigRational.Error(Code: TErrorCode; ErrorInfo: array of const);
811
begin
812
  case Code of
813
    ecParse:
814
      raise EConvertError.CreateFmt(SErrorParsingFmt, ErrorInfo);
815
    ecDivByZero:
816
      raise EDivByZero.Create(SDivisionByZero);
817
    ecConversion:
818
      raise EConvertError.CreateFmt(SConversionFailedFmt, ErrorInfo);
819
    ecInvalidArg:
820
      raise EInvalidArgument.CreateFmt(SInvalidArgumentFmt, ErrorInfo);
821
    ecZeroDenominator:
822
      raise EDivByZero.Create(SZeroDenominator);
823
  end;
824
end;
825

826
class operator BigRational.Explicit(const Value: BigRational): Integer;
827
begin
828
  if Value.FDenominator = BigInteger.One then
829
    Result := Integer(Value.FNumerator)
830
  else
831
    Result := Integer(Value.FNumerator div Value.FDenominator);
832
end;
833

834
class operator BigRational.Explicit(const Value: BigRational): string;
835
begin
836
  Result := Value.ToString;
837
end;
838

839
class operator BigRational.Explicit(const Value: BigRational): Double;
840
begin
841
  if Value.FDenominator = BigInteger.One then
842
    Result := Value.FNumerator.AsDouble
843
  else
844
    Result := Value.FNumerator.AsDouble / Value.FDenominator.AsDouble;
845
end;
846

847
class operator BigRational.Explicit(const Value: BigRational): Int64;
848
begin
849
  if Value.FDenominator = BigInteger.One then
850
    Result := Int64(Value.FNumerator)
851
  else
852
    Result := Int64(Value.FNumerator div Value.FDenominator);
853
end;
854

855
function BigRational.GetSign: Integer;
856
begin
857
  if FNumerator.IsZero then
858
    Result := 0
859
  else if FNumerator.IsNegative then
860
    Result := -1
861
  else
862
    Result := 1;
863
end;
864

865
class operator BigRational.GreaterThan(const Left, Right: BigRational): Boolean;
866
begin
867
  Result := Compare(Left, Right) > 0;
868
end;
869

870
class operator BigRational.GreaterThanOrEqual(const Left, Right: BigRational): Boolean;
871
begin
872
  Result := Compare(Left, Right) >= 0;
873
end;
874

875
class operator BigRational.Implicit(const Value: Integer): BigRational;
876
begin
877
  Result.FNumerator := BigInteger(Value);
878
  Result.FDenominator := BigInteger.One;
879
end;
880

881
class operator BigRational.Implicit(const Value: string): BigRational;
882
begin
883
  Result.Create(Value);
884
end;
885

886
class operator BigRational.Implicit(const Value: Double): BigRational;
887
begin
888
  Result.Create(Value);
889
end;
890

891
class operator BigRational.Implicit(const Value: Int64): BigRational;
892
begin
893
  Result.FNumerator := BigInteger(Value);
894
  Result.FDenominator := BigInteger.One;
895
end;
896

897
class operator BigRational.Implicit(const Value: BigDecimal): BigRational;
898
begin
899
  if Value.Scale = 0 then
900
  begin
901
    Result.FNumerator := Value.UnscaledValue;
902
    Result.FDenominator := BigInteger.One
903
  end
904
  else if Value.Scale > 0 then
905
  begin
906
    Result.FNumerator := Value.UnscaledValue;
907
    Result.FDenominator := BigInteger.Pow(BigInteger.Ten, Value.Scale);
908
  end
909
  else
910
  begin
911
    Result.FNumerator := Value.UnscaledValue * BigInteger.Pow(BigInteger.Ten, -Value.Scale);
912
    Result.FDenominator := BigInteger.One;
913
  end;
914
  Result.Normalize;
915
end;
916

917
class operator BigRational.Implicit(const Value: BigRational): BigDecimal;
918
begin
919
  Result := BigDecimal(Value.FNumerator) / Value.FDenominator;
920
end;
921

922
{$IFDEF HasClassConstructors}
923
class constructor BigRational.Initialize;
924
{$ELSE}
925
procedure Init;
926
{$ENDIF}
927
begin
928
  BigRational.FAlwaysReduce := True;
929
  BigRational.Zero := BigRational.Create(BigInteger.Zero, BigInteger.One);
930
  BigRational.OneTenth := BigRational.Create(BigInteger.One, BigInteger.Ten);
931
  BigRational.OneFourth := BigRational.Create(BigInteger.One, BigInteger(4));
932
  BigRational.OneThird := BigRational.Create(BigInteger.One, BigInteger(3));
933
  BigRational.OneHalf := BigRational.Create(BigInteger.One, BigInteger(2));
934
  BigRational.TwoThirds := BigRational.OneThird + BigRational.OneThird;
935
  BigRational.ThreeFourths := BigRational.OneHalf + BigRational.OneFourth;
936
  BigRational.One := BigRational.Create(BigInteger.One, BigInteger.One);
937
  BigRational.Ten := BigRational.Create(BigInteger.Ten, BigInteger.One);
938
end;
939

940
class operator BigRational.IntDivide(const Left, Right: BigRational): BigInteger;
941
begin
942
  Result := (Left.FNumerator * Right.FDenominator) div (Left.FDenominator * Right.FNumerator);
943
end;
944

945
function BigRational.IsNegative: Boolean;
946
begin
947
  Result := Self.FNumerator.IsNegative;
948
end;
949

950
function BigRational.IsPositive: Boolean;
951
begin
952
  Result := Self.FNumerator.IsPositive;
953
end;
954

955
function BigRational.IsZero: Boolean;
956
begin
957
  Result := Self.FNumerator.IsZero;
958
end;
959

960
function BigRational.Reciprocal: BigRational;
961
begin
962
  if FNumerator.IsZero then
963
    Error(ecDivByZero, []);
964
  if Self.FNumerator.IsNegative then
965
  begin
966
    Result.FNumerator := -Self.FDenominator;
967
    Result.FDenominator := -Self.FNumerator;
968
  end
969
  else
970
  begin
971
    Result.FNumerator := Self.FDenominator;
972
    Result.FDenominator := Self.FNumerator;
973
  end;
974
end;
975

976
function BigRational.Reduce: BigRational;
977
begin
978
  Result := Self;
979
  Result.Normalize(True);
980
end;
981

982
class function BigRational.Remainder(const Left, Right: BigRational): BigRational;
983
begin
984
  Result := Left mod Right;
985
end;
986

987
class operator BigRational.LessThan(const Left, Right: BigRational): Boolean;
988
begin
989
  Result := Compare(Left, Right) < 0;
990
end;
991

992
class operator BigRational.LessThanOrEqual(const Left, Right: BigRational): Boolean;
993
begin
994
  Result := Compare(Left, Right) <= 0;
995
end;
996

997
class operator BigRational.Modulus(const Left, Right: BigRational): BigRational;
998
var
999
  AD, BC: BigInteger;
1000
begin
1001
  // Result := Left - (Left div Right) * Right;
1002
  // This can be elaborated:
1003
  // X := A/B mod C/D -->
1004
  // X := A/B - (AD div BC) * C/D; -->
1005
  // X := AD/BD - (AD div BC) * BC/BD; -->
1006
  // X := [AD - (AD div BC) * BC] / BD;
1007

1008
  //$$RV: Can this be simplified?
1009
  AD := Left.FNumerator * Right.FDenominator;
1010
  BC := Left.FDenominator * Right.FNumerator;
1011
  Result.FNumerator := AD - (AD div BC) * BC;
1012
  Result.FDenominator := Left.FDenominator * Right.FDenominator;
1013
  Result.Normalize;
1014
end;
1015

1016
class operator BigRational.Multiply(const Left, Right: BigRational): BigRational;
1017
begin
1018
  Result.FNumerator := Left.FNumerator * Right.FNumerator;
1019
  Result.FDenominator := Left.FDenominator * Right.FDenominator;
1020
  Result.Normalize;
1021
end;
1022

1023
function BigRational.Negate: BigRational;
1024
begin
1025
  Result.FDenominator := Self.FDenominator;
1026
  Result.FNumerator := -Self.FNumerator;
1027
end;
1028

1029
class operator BigRational.Negative(const Value: BigRational): BigRational;
1030
begin
1031
  Result.FDenominator := Value.FDenominator;
1032
  Result.FNumerator := -Value.FNumerator;
1033
end;
1034

1035
procedure BigRational.Normalize(Forced: Boolean = False);
1036
var
1037
  GCD: BigInteger;
1038
begin
1039
  if FDenominator.IsZero then
1040
    Error(ecZeroDenominator, []);
1041

1042
  if FDenominator = BigInteger.One then
1043
    Exit;
1044

1045
  if FDenominator.IsNegative then
1046
  begin
1047
    FNumerator := -FNumerator;
1048
    FDenominator := -FDenominator;
1049
  end;
1050

1051
  if (FAlwaysReduce or Forced) then
1052
  begin
1053
    GCD := BigInteger.Abs(BigInteger.GreatestCommonDivisor(FNumerator, FDenominator));
1054

1055
    // TODO: See if this can be simplified by shifting common low zero bits away first
1056
    if GCD > BigInteger.One then
1057
    begin
1058
      FNumerator := FNumerator div GCD;
1059
      FDenominator := FDenominator div GCD;
1060
    end;
1061
  end;
1062
end;
1063

1064
class operator BigRational.NotEqual(const Left, Right: BigRational): Boolean;
1065
begin
1066
  Result := Compare(Left, Right) <> 0;
1067
end;
1068

1069
function BigRational.Parse(const S: string): BigRational;
1070
begin
1071
  if not TryParse(S, Result) then
1072
    Error(ecParse, [S, 'BigRational']);
1073
end;
1074

1075
class function BigRational.Multiply(const Left, Right: BigRational): BigRational;
1076
begin
1077
  Result := Left * Right;
1078
end;
1079

1080
class operator BigRational.Subtract(const Left, Right: BigRational): BigRational;
1081
begin
1082
  if Left.FDenominator = Right.FDenominator then
1083
  begin
1084
    Result.FNumerator := Left.FNumerator - Right.FNumerator;
1085
    Result.FDenominator := Left.FDenominator;
1086
  end
1087
  else
1088
  begin
1089
    Result.FDenominator := Left.FDenominator * Right.FDenominator;
1090
    Result.FNumerator := Left.FNumerator * Right.FDenominator - Right.FNumerator * Left.FDenominator;
1091
  end;
1092
  Result.Normalize;
1093
end;
1094

1095
class function BigRational.Subtract(const Left, Right: BigRational): BigRational;
1096
begin
1097
  Result := Left - Right;
1098
end;
1099

1100
function BigRational.ToString: string;
1101
begin
1102
  if FDenominator = BigInteger.One then
1103
    Result := FNumerator.ToString
1104
  else
1105
    Result := FNumerator.ToString + '/' + FDenominator.ToString;
1106
end;
1107

1108

1109
function BigRational.TryParse(const S: string; out Value: BigRational): Boolean;
1110
var
1111
  LSlash: Integer;
1112
  Num, Denom: BigInteger;
1113
begin
1114
  if S = '' then
1115
    Exit(False);
1116

1117
  Value := Zero;
1118
  LSlash := Pos('/', S);
1119
  if LSlash < 1 then
1120
  begin
1121
    Result := BigInteger.TryParse(Trim(S), Num);
1122
    Denom := BigInteger.One;
1123
  end
1124
  else
1125
  begin
1126
    Result := BigInteger.TryParse(Trim(Copy(S, 1, LSlash - 1)), Num);
1127
    Result := BigInteger.TryParse(Trim(Copy(S, LSlash + 1, MaxInt)), Denom) or Result;
1128
  end;
1129
  if Result then
1130
    Value := BigRational.Create(Num, Denom);
1131
end;
1132

1133
{$IFNDEF HasClassConstructors}
1134
initialization
1135
  Init;
1136
{$ENDIF}
1137

1138
end.
1139