Libraries · Optimization · Uncategorized

High-Precision Floating-Point Types for Delphi

Erik van Bilsen

If Single and Double precision floating-point types don’t cut it for you, then the Neslib.MultiPrecision library may be just the thing. It provides up to 4 times the precision, while still being pretty fast compared to arbitrary precision libraries.

Neslib.MultiPrecision

Neslib.MultiPrecision is a personal project that adds two additional floating-point types for use with Delphi. These offer a precision that is up to 4 times larger than that of the Double type.

It is build on top of the QD 2.3.22 and works on:

  • Windows (32-bit and 64-bit)
  • MacOS (64-bit)
  • iOS (64-bit, no simulator)
  • Android (32-bit and 64-bit)

The algorithms used for these types were developed by David H. Bailey, Yozo Hida and Xiaoye S. Li. Take a look the the QD.pdf file in the repository if you are interested in the details.

This library has no dependencies on other libraries. There are also no run-time dependencies. The underlying C/C++ QD library is linked into the executable using object files or static libraries.

You can find the library on GitHub in my Neslib.MultiPrecision repository.

Two new Floating-Point Types

This library defines two additional floating-point types:

  • DoubleDouble: this is a 128-bit type that has twice the precision of the Double type.
  • QuadDouble: this is a 256-bit type that has four times the precision of the Double type.

Both types have the same range as the Double type (about ± 10308) but much higher precision:

TypeExponent BitsMantissa BitsPrecision (decimal digits)
Single9247
Double125316
DoubleDouble1210632
QuadDouble1221264

The names DoubleDouble and QuadDouble come from the underlying QD library. I considered naming these Float128 and Float256 instead, but there are already official IEEE specifications for such (hardware) types. And since these are not compatible with DoubleDouble and QuadDouble, I didn’t want to add confusion.

The underlying QD library does not use emulation for calculations. Instead, it uses some interesting algorithms and the existing floating-point capabilities of the CPU to use either 2 or 4 Double values to increase the precision. As a result, these types are much faster than other arbitrary/high precision math libraries using the same precision. Also, as opposed to many other libraries, these types don’t require dynamic memory allocations, which further helps performance and also reduces memory fragmentation. In addition, many other libraries have (L)GPL (like) licenses, which may not be suitable for your purposes. Neslib.MultiPrecision uses a simplified BSD license, and the underlying QD library uses a variation of the BSD license as well, making them suitable for practically all situations.

Usage

The DoubleDouble and QuadDouble types can be used much the same way as the Double type. They support the usual operators (+, -, *, /, =, <>, <, <=, > and >=) as well as most methods available for the record helpers for the Double type (IsNan, IsInfinity, IsNegativeInfinity, IsPositiveInfinity, ToString, Parse and TryParse). There are methods and operators to convert to and from Double, DoubleDouble, QuadDouble and String.

Since the Delphi language doesn’t allow you to enter high-precision floating-point literals in source code, the value of a DoubleDouble or QuadDouble variable must be initialized in some other way. There are various options (assuming DD is of type DoubleDouble here):

  • Use one of the Init overloads. For example, DD.Init(1.23); to initialize a DoubleDouble from a Double.
  • Use an explicit cast, as in DD := DoubleDouble(1.23);.
  • Use one of the predefined constants, as in DD := DoubleDouble.Pi;.
  • Or perhaps the easiest way is to just assign a string value, as in DD := '3.1415926535897932384626433832795';.

I didn’t add implicit operators, such as DD := 1.23;, since those can lead to unintentional conversions that impact performance.

It’s important that you call MultiPrecisionInit before performing any DoubleDouble/QuadDouble calculations. This prepares the FPU/CPU for high-precision math. You can use MultiPrecisionReset to restore the FPU/CPU to its previous state.

Customization

The default configuration of the library is suitable for most applications. This configuration sacrifices a bit of accuracy for increased speed. If accuracy is more important than speed for your purposes, then you can compile the library with the MP_ACCURATE define. This will make many calculations a bit slower but more accurate.

Mathematical Functions

The underlying QD library (and thus this library) supports a variety of common mathematical functions, such as trigonometric, exponential and logarithmic functions.

In addition, the Neslib.MultiPrecision library adds numerous equivalents of functions found in the System.SysUtils and System.Math units, such as FloatToStr, Floor, Ceil, Min, Max, EnsureRange etc.

Samples

A fun way to demonstrate high-precision math is by calculating the Mandelbrot fractal. As you zoom into the fractal, you need more and more precision. The Samples subdirectory contains a FireMonkey application that generates the Mandelbrot fractal at 4 levels of precision (Single, Double, DoubleDouble and QuadDouble).

The following image shows a tiny section of the fractal using a magnification of one quadrillion (1015, aka one billiard in Europe) and Double precision:

You can clearly see that the Double type doesn’t provide enough precision at this magnification level. The DoubleDouble type offers more than enough precision though:

It’s not until you reach a magnification level of 1031, that you need to switch to QuadDouble. The sample application allows you to enter a magnification level of up to 1038:

I know, these images don’t show a particularly interesting part of the Mandelbrot set. This is because I chose center coordinates that aren’t near the “cardioid” or other bulbs of the set. This results in fewer calculations, so the sample app can focus more on differences in precision without having to wait forever for the rendering to complete.

This sample application is build with the FireMonkey framework, so it works on all platforms natively supported by Delphi.

More Information

There is more to Neslib.MultiPrecision than described above. For more details you can look at the well-documented Neslib.MultiPrecision.pas source file. Additional usage samples can be found in the UnitTests subdirectory and the Mandelbrot sample application.

If you are interested in the technical details and algorithms used for these types, you can take a look at the QD.pdf file in the C subdirectory.

If you need arbitrary precision, or support for big integers, decimals and rationals, then I can highly recommend Rudy Velthuis’s DelphiBigNumbers libary.

17 thoughts on “High-Precision Floating-Point Types for Delphi”

  2. Hi Erik,
    In the Neslib.MultiPrecision.pas unit, the routine
    Initialize
    in the initialization part of the unit causes an exception ‘floating point stack check at 0x00432f0c” in the debugger.

    Specifically, the problem occurs within the

    class operator QuadDouble.Multiply(const A: QuadDouble; const B: Double): QuadDouble;
    begin
    _qd_mul_qd_d(A, @B, Result);
    end;

    in the line _qd_mul_qd_d(A, @B, Result);

    I am using Delphi XR5 Professional, Windows 10: Target platform: Win 32 bit.

    What am I doing wrong? Or is there a bug in the library?

    Thanks, Kind Regards, Andreas
    Germany

    Like

    Reply

    1. This may be specific to the Delphi version you are using, since I don’t have this issue with Delphi 10.4.2.
      I just pushed an update that should fix it. Let me know if this works for you.

      Like

      Reply

  3. Thanks Erik,
    after the modification from you there is no exception anymore. However, I can’t print my simple test result: it always comes out an empty string with exponent like this:
    ‘ , E+00′

    Here is a small excerpt from the console program:
    …
    VAR
    FPU_Status: UInt32;
    A : QuadDouble;
    Exp: String;
    S : String;

    Begin
    FPU_Status:= MultiPrecisionInit;
    A := Abs(QuadDouble.Pi);
    S := QuadDouble.ToString(A);
    Exp:=’3.14159265358979323846264338327950288419716939937510582097494459’;

    WriteLn(‘Expected = ‘, Exp);
    WriteLn(‘Actual = ‘, S);

    // Result: —> Actual = ‘ , E+00’

    ReadLn;
    …
    MultiPrecisionReset(FPU_Status);
    …
    End.

    Thanks, Kind Regards, Andreas

    Like

    Reply

    1. Hi Andreas,
      Must have something to do with the Delphi version.
      I have also tried it on Delphi 10.2 (which is the oldest Delphi version I have installed), and it works as expected (after making a few changes because 10.2 doesn’t support inline variables).
      This makes is hard for me to debug this.

      I would advise to get an update subscription and always use the latest Delphi version anyway, but I know that is not the answer you are looking for 😉

      If I find some time, I may install an older Delphi version and try it out. You said you are using XE5, right?

      Like

      Reply

    1. Hi Erik,
      I have been able to isolate the possible error somewhat. The routine that seems to be causing problems is:

      procedure QuadDouble.ToDigits(const S: TArray; var Exponent: Integer;
      const Precision: Integer);
      . . .
      else if (E > 300) then
      begin
      R := Ldexp(R, -53);
      R := R / IntPower(Ten, E);
      R := Ldexp(R, 53);
      end
      else
      R := R / IntPower(Ten, E);

      { Fix exponent if we are off by one }
      . . .

      Specifically, the line: R := R / IntPower(Ten, E);

      Here the record R still contains the correct digits of Pi:

      R = (3.14159265358979, 1.22464679914735e-16, -2.99476980971834e-33, 1.11245422086337e-49)

      Ten = (10, 0, 0)
      E = 0

      After division, actually by 10^0, that is by 1, all fields of R are set to +NAN:
      R = (+NAN, +NAN, +NAN, +NAN).

      I’m sure it’s because 0^0 is not mathematically defined and IntPower(Ten, E) throws an exception. Maybe the special case 0^0 would have to be caught before and the power would have to be set to 1 “by hand”. Maybe it helps you?

      Thanks, Kind Regards
      Andreas

      Like

      Reply

      1. Thank you Andreas. Now I know what to look for when I install an older Delphi version.

        Like

      2. I’m sorry Andreas, but I didn’t find a solution.

        I can reproduce the issues with XE5, but it seems to have to do with the way Delphi XE5 links C object files. Tried creating those object files in different ways, with different calling conventions but couldn’t get it to work. On Win32, it would generate invalid results, and on Win64 it would result in Access Violations.

        I don’t know in which Delphi version this was changed to the current (working) behavior. But I can’t spend too much time on trying to support an 8 year old Delphi version.

        Again, I really recommend upgrading to a more recent (especially most recent) Delphi version. But if you cannot do that, then you will have to experiment with debugging this yourself. Instructions on how to build the object files can be found in the readme.txt file in the C subdirectory.

        I’m sorry I cannot be more helpful.

        Like

  5. Hi Erik,
    first of all, many thanks for your help and toil!

    But because your “High-Precision Floating-Point Types for Delphi” would be very important also for the numerous users of older Delphi versions, I created a topic in my German “home forum” Delphi-Praxis. That is:

    „High-Precision Floating-Point Types for Delphi – Hilfe beim Linken von C object files“

    https://www.delphipraxis.net/207850-high-precision-floating-point-types-delphi-hilfe-beim-linken-von-c-object-files.html#post1488969

    There, I have asked my colleagues with special expert knowledge for help, so that this valuable treasure can be made probably available also to earlier Delphi version’s users. I hope very much that helpful suggestions for a proper solution will come.
    If something positive happens, I will be keeping you informed.

    Thank you again and kind regards
    Andreas

    Like

    Reply

  6. Hi Erik,
    sorry, I have to bother you again…
    After some discussion with my colleagues in the two forums mentioned above, I installed the German version of [B]Delphi 10.2 CE[/B] at my place to be able to do the tests with your library myself without any detours. In addition, I have downloaded today (20. May 2021) the current version of [B]Neslib.MultiPrecision-master.zip[/B] from your homepage https://github.com/neslib/Neslib.MultiPrecision.
    downloaded.

    I did my first tests with your original files, without any modification to your source code. I am sorry to tell you that [DELPHI]Neslib.MultiPrecision.pas[/DELPHI] must still contain several bugs.
    Attached is a small console test program which uses some of your tests from [DELPHI]MultiPrecision.QuadDouble.Tests[/DELPHI]:

    [DELPHI]
    Program NesLib_MP_Test_1;

    {$APPTYPE CONSOLE}

    {$R *.res}

    Uses
    System.SysUtils
    , Neslib.MultiPrecision
    // UnitTest,
    // MultiPrecision.DoubleDouble.Tests,
    , MultiPrecision.QuadDouble.Tests
    ;

    VAR
    FPU_Status: UInt32;
    T : TTestQuadDouble;

    Begin
    FPU_Status:= MultiPrecisionInit;

    Try
    T:= TTestQuadDouble.Create;

    Try

    T.TestInit; // Crash!
    T.TestToString; // Crash!
    T.TestArcTan2; // Crash!

    ReadLn;

    Finally
    T.Free;
    MultiPrecisionReset(FPU_Status);
    End;
    Except
    On E: Exception Do
    Writeln(E.ClassName, ‘: ‘, E.Message);
    End;
    End.
    [/DELPHI]

    There is always an exception in the line

    4515: [DELPHI]FromString := StrToFloat(SB.ToString, FormatSettings);[/DELPHI]

    From routine:
    [DELPHI]
    function QuadDouble.ToString(const FormatSettings: TFormatSettings;
    const Format: TMPFloatFormat; const Precision: Integer): String;
    var
    SB: TStringBuilder;
    FromString: Double;
    Off, D, DWithExtra, I, E: Integer;
    T: TArray;
    begin
    E := 0;
    SB := TStringBuilder.Create;
    . . .

    // hiere in line 4515:
    FromString := StrToFloat(SB.ToString, FormatSettings);
    . . .
    [/DELPHI]

    The problem is probably related to the DecimalSeparator. In the German version of Delphi 10.3 the following is the default:
    ThousandSeparator = ‘.’, and DecimalSeparator = ‘,’.

    The first time you call
    [DELPHI]
    A.Init(1.5);
    CheckEquals(1.5, A.X[0], 0);
    CheckEquals(0, A.X[1], 0);
    CheckEquals(0, A.X[2], 0);
    CheckEquals(0, A.X[3], 0);
    CheckEquals(‘1.50000000000000000000000000000000000000000000000000000000000000’, A);
    [/DELPHI]
    in your [DELPHI]MultiPrecision.QuadDouble.Tests.pas[/DELPHI] the content of the variable SB (of type TStringBuilder) is correct, with DecimalSeparator = ‘,’:
    SB = (‘1,50000000000000000000000000000000000000000000000000000000000000’#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0#0, 64, 2147483647)

    But in the next call
    [DELPHI]
    A.Init(DoubleDouble.Log2);
    CheckEquals(6.931471805599452862e-01, A.X[0], 0);
    CheckEquals(2.319046813846299558e-17, A.X[1], 0);
    CheckEquals(0, A.X[2], 0);
    CheckEquals(0, A.X[3], 0);
    CheckEquals(‘0.69314718055994530941723212145817599730465629273908450073995924’, A);
    [/DELPHI]
    the content of the variable SB is FALSE:
    SB = (‘0.69314718055994530941723212145817599730465629273908450073995924000000000000000000000000000000000000000000000000000000’#0’000000000’, 64, 2147483647).
    Therefore, your [DELPHI]MultiPrecision.QuadDouble.Tests.pas[/DELPHI] crashes with the following error message:
    Exception of class EConvertError:
    ‘0.69314718055994530941723212145817599730465629273908450073995924’
    is not a valid floating point value.

    [B] Next problem:[/B]
    Mandelbrot.dpr can be compiled, but the program does not show any Mandelbrot graphics.
    And if I start the program directly in the debugger and immediately click on the “Update” button, the same exception occurs:

    Exception of the class EConvertError:
    ‘-0.00677652295833245729642263781984627256356509565412970431582937’
    is not a valid floating point value.

    Erik, could you please take a look and correct this error, because it is not due to the Delphi version.

    Thanks a lot!

    Kind regards
    Andreas

    Like

    Reply

    1. Thanks for reporting this Andreas. I usually try to take international format settings into account since I am Dutch and in The Netherlands we also use a comma as a decimal separator.

      There are a couple of issues here.

      First, the unit tests assumed they were ran using format settings that use a period as a decimal separator. I modified the unit tests so they use a specific TFormatSettings record that has a period as decimal separator. Most String conversion methods accept an optional TFormatSettings parameter that you should use to specify the format of the source string. I forgot to do that in the unit tests.

      Next, when assigning a string using the implicit operator (as in DD := ‘3.14’), you had to make sure that the string was in the current system format settings. The documentation even said so:

      { Implicitly converts a string to a DoubleDouble.
      The string *must* use the current system-wide format settings. }
      class operator Implicit(const S: String): DoubleDouble; inline; static;

      However, I understand that implicit conversions are a bit different from other string conversions. So I modified these to *always* use a period as decimal separator (because you cannot use a TFormatSettings record with implicit operators). I think this makes more sense since a string literal in this case is more like a regular floating-point literal in Delphi (which also always uses a period as decimal separator).

      I pushed an update with these changes and fixes. Let me know if this solves the issue for you.

      Like

      Reply

  7. Many-many thanks, Erik!

    My first tests with your improved version of the library are very positive: It works not only with Delphi 10.3 but if I rewrite the inline declarations in the two routines
    – function MultiPrecisionInit: UInt32;
    – procedure MultiPrecisionReset(const AState: UInt32);
    also with Delphi XE5 flawlessly!

    High praise for you!

    My suggestion:
    If you would renounce the new-fashioned inline declarations in a future release, also the users of older Delphi versions could use your excellent library directly without any modifications.

    Since in my technical-scientific calculations the exactness of the calculation results has the highest priority, I will still test all your mathematical functions extensively against my up to now used multi-precision routines. I will definitely report you about my tests.

    Until then, thank you very much for your excellent work!
    You are a super programmer!

    Kind regards
    Andreas

    Like

    Reply

  8. Hi Erik,
    it’s me again…
    There is still something wrong with QuadDouble. It does not accept assignments like
    Z_Qd := ‘0.5’;
    and Z_Qd still remains NAN and therefore all operations with QuadDouble numbers return NAN as result.

    Here is an example print-out:
    ArcTan2
    ———–
    X = -0.5
    Y = 0.5
    MPA_Float = 2.356194490192344928846982537459627163147877049531329365731208444230862
    DoubleDouble = 2.3561944901923449288469825374596
    QuadDouble = NAN

    *******************
    QuadDouble.ToString(X) = NAN
    QuadDouble.ToString(‘0.5’) = NAN
    *******************
    Could you please check again?

    Thank you!
    Andreas

    Like

    Reply

    1. Hi Andreas,

      I cannot reproduce this issue. When I perform the tests above, I get the expected results.
      What Delphi version are you testing this with?

      Also, how are you able to use the code with XE5? When I call into the underlying C APIs on XE5 I get access violations. Were you able to build the object files to work with XE5?

      Finally, for administrative purposes, would you mind posting future issues on GitHub for easier tracking?

      Thanks,
      Erik

      Like

      Reply

Leave a comment