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 theDoubletype.QuadDouble: this is a 256-bit type that has four times the precision of theDoubletype.
Both types have the same range as the Double type (about ± 10308) but much higher precision:
| Type | Exponent Bits | Mantissa Bits | Precision (decimal digits) |
|---|---|---|---|
| Single | 9 | 24 | 7 |
| Double | 12 | 53 | 16 |
| DoubleDouble | 12 | 106 | 32 |
| QuadDouble | 12 | 212 | 64 |
The names
DoubleDoubleandQuadDoublecome from the underlying QD library. I considered naming theseFloat128andFloat256instead, but there are already official IEEE specifications for such (hardware) types. And since these are not compatible withDoubleDoubleandQuadDouble, 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
Initoverloads. For example,DD.Init(1.23);to initialize aDoubleDoublefrom aDouble. - 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.
grijjy blog 
excellent! thanks …
LikeLike
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
LikeLike
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.
LikeLike
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
LikeLike
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?
LikeLike
Yes it is XE5. That was a typo before. Thanks for your help Erik!
LikeLike
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
LikeLike
Thank you Andreas. Now I know what to look for when I install an older Delphi version.
LikeLike
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.
LikeLike
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
LikeLike
Hi Erik,
According to the motto: “two is better than one”, I have also published the topic in an other renowned German programming forum Entwickler-Ecke. It is the issue “Hilfe beim Linken von C object files in in Delphi XE5?”
https://entwickler-ecke.de/topic_Hilfe+beim+Linken+von+C+object+files+in++in+Delphi+XE5_118708.html&sid=6ec5e7d91d4224819de324c038d33e60
Kind regards
Andreas
LikeLiked by 1 person
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
LikeLike
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.
LikeLike
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
LikeLike
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
LikeLike
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
LikeLike
To be continued at:
https://github.com/neslib/Neslib.MultiPrecision
LikeLike