version 1.15 (1.39 MB) by
John D'Errico

High precision floating point arithmetic, a new class written in MATLAB

Very often I see people asking for a tool that offers more than 16 digits or so of accuracy. MATLAB itself only allows you to use doubles or singles in standard arithmetic, so that is the normal limit. The fact is, most of the time, if you can't do it with a double, you are doing something wrong. Good practices of numerical analysis are worth far more than any high precision tool. Even so, there are times when you will have a use for a bit of extra precision. And some of you will just want to play in the huge number sandbox. While some of you may use tools like that written by Ben Barrowes, HPF is written purely in MATLAB, so no compiles are needed. For all of you, whatever your reasons, I offer HPF, a High Precision Floating point tool.

In fact, the reason I wrote HPF was for my own purposes. I wanted to learn to use the classdef tools in MATLAB that were released a few years ago. As well, I wanted to try building such a tool as a natural extension of the VPI tools I wrote some time ago. And I wanted to learn some techniques for working in a high number of digits. The result is HPF.

There are a few ideas I've introduced for how HPF interacts with the user. For example, HPF can work in any number of decimal digits, as chosen by the user. You can set the number of digits as a default. Thus, if you want to always work in 30 decimal digits, with 2 guard digits on all computations, then type this at the command prompt:

DefaultNumberOfDigits 30 2

From now on, for you HPF will always work in a total of 32 decimal digits of precision, and report the top 30 digits, thus two guard digits will be used internally. For example,

pie = hpf('pi')

pie =

3.14159265358979323846264338328

exp(pie - 3)

ans =

1.15210724618790693014572854771

HPF will recall this state the next time you start MATLAB. You can override the default state by specifying a different number of digits though.

hpf('e',12)

ans =

2.71828182846

I've stored values as part of HPF for e and pi that are accurate to 500,000 digits. In fact, those numbers were generated by HPF itself.

Finally, for speed and efficiency, HPF stores all numbers in the form of Migits, which are bundles of decimal digits. This yields a huge bonus in the speed of multiplies, since conv is employed for that purpose. We can see them here:

pie.Migits

ans =

[3141 5926 5358 9793 2384 6264 3383 2795]

The nice thing is that the use of Migits will be transparent to most users. But if you want a bit more speed in your multiples, then you can get a boost by typing this:

DefaultDecimalBase 6

From now on, HPF will employ base 1000000 migits internally, what I call 6-migits. The only problem is, you will be restricted from using numbers with more than 36000 decimal digits. Speed has a price.

Another nice use of HPF is to extract the exact decimal form that MATLAB uses to store its own numbers. For example, what number does MATLAB REALLY store internally when you type in something like 1.23?

hpf(1.23,55)

ans =

1.229999999999999982236431605997495353221893310546875000

Is HPF complete as it stands? Of course not. HPF currently represents nearly 7000 lines of MATLAB code, in the form of dozens of methods available for the class. As it is, you will find many hundreds of hours of work on my part, over the course of several years. But I've not yet written a huge number of things that might be useful to some people. For example: roots, eig, chol, det, rank, backslash, gamma, etc. And HPF offers no support for complex numbers. Even so, I hope that some will find this useful, if only to learn some of the tricks I've employed in the building thereof. Some of those tricks are described in HPF.pdf.

For example, multiplies are best done in MATLAB by conv. But divides take more work, so here I use a Newton scheme that employs only adds and multiplies, and is quadratically convergent. A similar trick is available for square roots.

Or, look into how my exponential function works. Here I've used a few tricks to enhance speed of convergence of the exponential series. Of course, there are obvious range reduction tricks, but I've gone an extra step there. I also employ a different way of summing the series for exponentials (as well as the sine and cosine series) that minimizes divides.

A lot of thought and research has gone into the methods of HPF. Those thoughts are captured in the HPFMod.pdf file, as enhanced by Derek O'Connor. Many thanks there. HPFMod.pdf is sort of a manual too, for those who want to truly understand the tool.

HPF will probably never be in what I consider to be in a final form, as I am sure there are a few bugs still unfound. Even so, the tool is working quite well on the thousands of tests I have performed. For those of you who try HPF out and do find a bug, please send me an e-mail and I will repair it immediately.

John D'Errico (2021). HPF - a big decimal class (https://www.mathworks.com/matlabcentral/fileexchange/36534-hpf-a-big-decimal-class), MATLAB Central File Exchange. Retrieved .

Created with
R2012a

Compatible with any release

**Inspired by:**
Multiple Precision Toolbox for MATLAB, num2strexact (exact version of num2str), Variable Precision Integer Arithmetic

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!Create scripts with code, output, and formatted text in a single executable document.

Alexey NenashevHi John! It's a really good job!

Toder noticed that hpf('-1',[16 0]) + hpf('-1e-20',[16 0]) gives ans = 1. I suppose it's because of

A.Sign = 1;

B.Sign = 1;

at lines 4495-4496 in file hpf.m

Ralph ColemanI think there is missing code in the nthroot method: it returns no result

Ameer HamzaRichardThank you for this, John!

Bug report -- I'm running into errors for mod(hpf(-1),1), ceil(hpf(1)), and floor(hpf(-1))

ToderHi John. I really like hpf (and vpi). I've noticed that

hpf('-1',[16 0]) + hpf('-1e-20',[16 0])

gives ans = 1. Can you help me understand why the sign changes?

Cemil KözkurtNitinMojgan RostaminiaSrighakollapu Vali MHi, I want to know whether I can pass hpf value to other functions and can those functions return hpf values?

Luis MendoJohn D'ErricoNote that as Fabio found, if you try to form

hpf(1/3)

MATLAB first creates 1/3 as a double. Then it passes the result to HPF. This fails because 1/3 is already corrupted with garbage past the 16th digit or so. The fixis as Fabio found,

fabio sacchiright, i find the way .... (:-)

F = hpf('1',100)/3

F =

0.3333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333

fabio sacchiops, the result of this ratio :

hpf(1/3)

is this :

ans = 0.3333333333333333148296162562473909929394721984863281250000000000

mmhhh surely there is an explanation of this but .... IT'S WRONG !!!

Can i use hpf in other way (preferably more exactly)

ErichThe power function doesn't handle negative integer exponents. Rather than create a new function, as David Verrelli did, I modified the power function in hpf.m (starting around line 4807) (All code referencing the variable "yneg" is added by me):

% y is small and integer

if isa(y,'hpf')

% y is small enough to safely convert to a flint

y = double(y);

end

yneg = false;

if y < 0

y = -y;

yneg = true;

end

% y is a numeric integer, but not zero or 1 or an inf or nan

% so convert y to binary, forming the powers of x by repeated squaring

ybin = dec2bin(y) == '1';

if ybin(end)

result = hpf(x,NDig);

else

result = hpf('1',NDig);

end

xsq = x.*x;

for i = (numel(ybin)-1):-1:1

if ybin(i)

result = result.*xsq;

end

if i > 1

xsq = xsq.*xsq;

end

end

if yneg

result = 1 ./ result;

end

end

end % function power

Joseph KirkHi John, this submission is fantastic! Great work!

I may have found a bug on Line 2997:

k = find(D ~= 0,'1','last');

should be

k = find(D ~= 0,1,'last');

Libor SedaDavid VerrelliAn additional function, saved as saveHPF.m to save an array of HPF data into a tab-delimited text file. There is certainly scope to improve this, such as including options for comma-separated variables (CSV), etc..

~~~~~~~~~~

% Output HPF array data to ASCII file

% Equivalent of

% save(outFile, varName, '-ascii','-double','-tabs')

% Input arguments:

% outFile output filename (and, optionally, path)

% arrayHPF array of HPF variables (no more than two dimensional)

% Output variable:

% res Boolean flag of function success

%

% David I. Verrelli (MQ), 2016-02-21

function res = saveHPF(outFile, arrayHPF)

res = false;

fid = fopen(outFile, 'w');

numRows = size(arrayHPF, 1);

for row = 1 : numRows

numCols = size(arrayHPF, 2);

for col = 1 : numCols

if col ~= numCols,

fprintf( fid , '%s\t' , num2str(arrayHPF(row, col)) );

else

fprintf( fid , '%s\r\n' , num2str(arrayHPF(row, col)) );

end;

end;

end;

fclose(fid);

res = true;

end

~~~~~~~~~~

[licence: CC BY]

David VerrelliAn additional function, saved as powerHPF.m to perform element-by-element exponentiation for HPF arguments with negative exponents allowed. There may well be scope to improve this further.

~~~~~~~~~~

% Compute elementwise power for HPF arguments in which the exponent y may be negative

% David I. Verrelli (MQ), 2016-02-21

function outcome = powerHPF(x,y)

numX = length(x(:));

numY = length(y(:));

numEls = max( numX , numY );

%outcome = nan * hpf(zeros(numEls, 1));

outcome = hpf( nan(numEls, 1) );

for i = 1 : numEls

if numX == 1,

elsX = 1;

else

elsX = i;

end;

if numY == 1,

elsY = 1;

else

elsY = i;

end;

xVal = x(elsX);

yVal = y(elsY);

if yVal < 0,

outcome(i) = 1 / (xVal ^ -yVal);

else

outcome(i) = xVal ^ yVal;

end;

end;

if numX > 1,

outcome = reshape(outcome, size(x));

elseif numY > 1,

outcome = reshape(outcome, size(y));

end;

end

~~~~~~~~~~

[licence: CC BY]

David VerrelliThere is apparently a small bug in the code for the eps function within hpf.m

On about line 1587, I believe the command

D(I) = eps(X(i));

should be amended to

D(i) = eps(X(i));

David VerrelliThis is a very helpful and extensive contribution. It was much easier to use than I had anticipated, due to the fact that it was implemented as a class.

As with any major project, there are opportunities to add more features; there may also be a couple of minor bugs that would not affect most users. I will expand on these points in separate comments here.

—DIV

Eneru YI've problem with fprintf, norm(), ecc.: "Function is not defined for 'hpf' inputs."

John D'ErricoShaun,

I'll admit that extending HPF to the complex plane would be a nice idea. I suppose that had I written it that way from scratch, a complex HPF would be my choice today. For a few seconds the other day, I toyed with the idea that I could simply extend it with a new class on top. However, there are too many places where the code would need to be updated, and then carefully tested, etc. SQRT or LOG, for example, both of which explicitly trap for negative inputs.

Sadly, I'll admit I spent perhaps at least a half man year writing this suite of code (then several more man-months all over several years re-writing it, because I wrote the entire toolbox three times, with three subtly different implementations to decide on which I felt was the best implementation.) My fear is that I would be forced to invest another large fraction of a man year, choosing algorithms that are robust and stable in the complex plane, then testing them all. And Bessel functions have their own intriguing set of quirks.

So, I am honestly sorry to say that while I hugely enjoyed writing HPF, that I could not do it again. I'll defer to the next person who wishes to write a version of a tool like HPF. Anyone who wishes to do so is welcome to use any and all of HPF, any ideas used in my implementation.

ShaunNice John

Been meaning to play with this for a while

I also vote for ...

bessel/complex methods

But, if you did nothing else, 5 Stars

John D'ErricoSorry, but not at this time. I've not implemented bessel functions in HPF, and, unless you only had real positive Z, besseli will result in complex results. Since I've not implemented a complex version of HPF, that too will be a problem.

Eric DiazIs there anyway to get this to work with the special besseli functions?

John D'ErricoI fixed that problem and posted the fix for that, I am sure. Well, I know that I fixed it.

However, it is conceivable I made a mistake, and did not update the file properly, or you might be using an older release. As well, there was a recent problem with the FEX where a file I uploaded did not get properly posted, so perhaps this might be their fault. :)

So I've just now uploaded the current release. As I just tested, your example SHOULD give the proper result.

hpf(-1, [20 0]) + hpf(-1e-21, [20 0])

ans =

-1

Denis GrebenkovDear John,

Thanks for this nice code! I have the same problem of sign changing as pointed out by Stephen Lucas: the command

> hpf(-1, [20 0]) + hpf(-1e-21, [20 0])

returns +1, instead of -1.

Is there anything I can do?

Best regards,

Denis

PS: An extension to complex-valued operations would be great...

John D'ErricoHi Eric,

it should have been:

recip = reciprocal(hpf(den,result(1).NumberOfDigits));

I'll post the fix.

v = 1:5;

d = hpf(23,17);

v./d

ans =

HPF array of size: 1 5

|1,1| 0.043478260869565217

|1,2| 0.086956521739130435

|1,3| 0.13043478260869565

|1,4| 0.17391304347826087

|1,5| 0.21739130434782609

Erik BenklerDear John,

I receive an error message (Improper index matrix reference.) on line 5061 of hpf.m, when dividing a double type Nx1 vector by a scalar hpf:

recip = reciprocal(hpf(den,num(1).NumberOfDigits)); %shouldn't this be result(1) instead of num(1), or simply reciprocal(den) ?

Best regards,

Erik

Erik BenklerStephen LucasI lovely contribution, but I have an example where I'm not getting the expected result. Consider:

a=hpf('0.1'); b=hpf('1e-400'); c=b-a, d=-b-a

In default precision (64 4), due to roundoff I would expect (and want!) the answer to be -0.1 in both cases. But Matlab is returning 0.1 for d, not -0.1.

if I use b=hpf('1e-50') I get the correct answer.

louis kovalevskyThai V. HoangThanks for your quick response. Unfortunately, I still cannot get the updated version from the official download today. Could you please send it to me somehow? Thank you!

John D'ErricoYes, that was indeed a bug, now fixed. I'm uploading a new version now (to appear today I hope.)

So with the new release, in 100 digits of precision...

exp(hpf('-2.22222222222222222222222222e-137',100))

ans =

1

And in 200 digits, we see a tiny difference from 1, as expected.

exp(hpf('-2.22222222222222222222222222e-137',200))

ans =

0.99999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999997777777777777777777777777780000000000000000000000000000000000000

This release adds one new toy: the rat function, implemented for hpf. So for a 100 digit rational fraction approximation to pi...

[N,D] = rat(hpf('pi',100),hpf('1e-101'))

N =

394372834342725903069943709807632345074473102456264

D =

125532772013612015195543173729505082616186012726141

N/D

ans =

3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068

hpf('pi',100)

ans =

3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068

Thai V. HoangHello John,

Could you please check the result of th following expression for me. I got -1 in my Linux Mint 13 Maya x64. Maybe this is a bug? Thanks

exp(hpf('-2.22222222222222222222222222e-137'))

ShaoyangAdamThank you very much.

John D'ErricoAdam,

I'd give you 10 to 1 odds (or more) that you have too old a MATLAB release to use it. I'm sorry, but it was written using 2012a, and it uses some of the capabilities from that release. The error that you got is consistent with what you would see from an older release.

So you need to upgrade your version of MATLAB.

John

AdamDear John D'Errico

I downloaded your program, but when I tell the MatLab that calculate the following line:

pie=hpf('pi')

I got Error with the next message:

File: hpf.m Line: 863 Column: 9

Expression or statement is incorrect--possibly unbalanced (, {, or [.

I do not find out what could be the problem. I did the exactly what you write in the ReadMe.rtf file. Do you have an idea?

Thank you your response.

John D'ErricoRaj - You are not listening to what I have said. HPF CANNOT be used with MATLAB's version of FZERO. It will NOT work. I would need to re-write FZERO for that to happen, and I have never done so. I said this before.

Why would you expect to find a demo for code that does not exist?

This does not say it is impossible to re-write FZERO, but that I have not chosen to do so.

Raj RajHi John,

Thanks for taking time in replying.

However, in your file "hpf_demo.m", I couldnt find any demo on the use of FZERO with HPF.

Are you referring to this file?

John D'ErricoHi Raj,

The problem is in FZERO, which explicitly tests its argument to see if it is a double, rather than a floating point numeric class like HPF. In fact, FZERO even fails to run on single input.

f = @sin;

fzero(f,hpf(3,50))

Error using fzero (line 128)

FZERO only accepts inputs of data type double.

To make it work, I would need to write an overloaded FZERO, or TMW would need to allow FZERO to accept other than double input for a starting value.

(Note that the way you are trying to call fzero will not work if an overloaded FZERO were supplied by me, as MATLAB would always use the base version then. In order to call an overloaded operator, you need to pass in a variable of the proper class. See the example case where I tried to use FZERO for HPF input.)

Raj RajMay be here is the application I encountered:

f1=@(x) ((5*x+7))+hpf(9*10^15+6*10^6)

fzero(f1, [0.00])

Raj RajHi John,

Can we use this function with fzero. Here is the test case:

f1=@(x) hpf(5*x+7)

f2=@(x) (5*x+7)

fzero(f2, [0.00])

fzero(f1, [0.00])

Appreciate your response.

Thanks!

John D'ErricoHI Yamasani,

It looks like you have a few misunderstandings. What you need to appreciate is how MATLAB works with numbers as an interpreted language. Also how assignments work when there are different classes involved.

First of all, when you do something like this...

V = fun(1 + 2.3);

Here MATLAB performs the operation 1 + 2.3, as a double. The result will be 3.3 (approximately), but it WILL be a double. Then that result is passed into the function fun. Here fun could be ANY function, including a call to HPF, to create an HPF number.

The same holds true for something like...

V = hpf(1.4);

MATLAB FIRST takes the number 1.4. It stores that as a double. Note that 1.4 is NOT exactly storable as a double. In fact, the representation of 1.4 in MATLAB is...

X = hpf(1.4)

X =

1.399999999999999911182158029987476766109466552734375000000000000

Remember that numbers are stored in binary form in MATLAB. So if you really wanted 1.4, you needed to do it like this

X = hpf('1.4')

X =

1.4

Here X is EXACTLY 1.4, stored in a decimal form inside HPF.

The same thing applied when you do an operation like

1.4 + hpf(1)

Here 1 is an integer, so HPF is capable of storing that exact value as an HPF version of the integer 1. But then adding 1.4 to that number, remember that MATLAB is an interpreted language. So again, MATLAB has a problem, because MATLAB sees the 1.4 as a DOUBLE number, only then does it realize that it must add that number to an HPF number.

1.4 + hpf(1)

ans =

2.399999999999999911182158029987476766109466552734375000000000000

So this does work properly:

hpf('1.4') + 1

ans =

2.4

2.4 is indeed the exact representation internally.

Next, you have another issue because IF the variable E exists already as a double. See what happens here:

E = zeros(1,2);

E(1) = hpf('1.4')

E =

1.4 0

whos E

Name Size Bytes Class Attributes

E 1x2 16 double

E is a double. It is NOT an HPF number. In fact, E is not represented as exactly 1.4. So if I look at the value stored for E, again, it is the wrong number, not exactly the true decimal 1.4. Since E existed already, and you chose to insert the HPF value into the first element of E, MATLAB decided to convert it BACK into a double.

hpf(E(1))

ans =

HPF element: (1,1)

1.399999999999999911182158029987476766109466552734375000000000000

In your last example, I don't know what you did. My guess is you were confused, and do not really have what you showed me.

Again, you must be careful about what MATLAB does to the numbers you are working with. And these are things totally beyond my control in the design of HPF. Remember that MATLAB is a language like any other programming language (but unlike what we see in the movies.) It has strict rules about how it will operate on what you pass it, and it cannot know that what you wanted it to do is not what you told it to do.

Yamasani RahulHello John

I am stuck up with aserious problem.HPF is working as expected under normal conditions,but when trying to store my number ina an array it is rounding off the number

Ex:E=1.4+hpf(0.002233112233445566)

>>E=1.402233112233445566

but if

E(1)=1.4+hpf(0.002233112233445566)

>>E()=1.4

Please help me

John D'ErricoI've made the tool as fast as possible, in fact, writing the entire tool using 3 different fundamental schemes for the best compromise in speed and capability. This was the best of the alternatives.

When you are working with large enough numbers, then setting the DefaultDecimalBase to 6 can help some, although there is little difference on numbers as small as only 100 digits or so.

Yamasani RahulThis is extraordinary

But processing of this tool is too slow

Derek O'ConnorThis is really excellent. I wish I'd had HPF when I was teaching numerical algorithms.

Here are some functions I use to test all such packages:

function cond1SSH = TestSSH;

% Testing John D'Errico's HPF on Sea Surface Heights problem.

% For more information on this see: http://www.scribd.com/doc/26135665/

% Download data from http://www.derekroconnor.net/NA/Notes/etaana.dat

% (a plain text file) and then import to SSHData.mat

% Derek O'Connor 10 Sep 2012

load('SSHData.mat');

[m,n] = size(SSHData); n = max(m,n);

for d = 20:30

shpc = sum(hpf(SSHData,d));

disp([shpc.NumberOfDigits(1), shpc.Migits, shpc.Exponent]);

end

cond1SSH = n*sum(abs(SSHData))/abs(sum(SSHData));

disp(['Condition Number of Data:']); disp(cond1SSH);

function z = Rump(d);

% Testing John D'Errico's HPF

% Derek O'Connor 10 Sep 2012

x = hpf(77617,d); % d = Number of digits to use

y = hpf(33096,d);

z = (33375/100)*y^6 + x^2*(11*x^2*y^2-y^6-121*y^4-2) + (55/10)*y^8 + x/(2*y);

% Correct answer: z = -54767/66192

function z = Judd(d)

% Testing John D'Errico's HPF

% Derek O'Connor 10 Sep 2012

x=hpf(192119201,d); % d = Number of digits to use

y=hpf(35675640,d);

z = (1682*x*y^4 + 3*x^3 + 29*x*y^2 - 2*x^5 + 832)/107751;

% Correct answer: z = 1783

JonathanThis is excellent! I had code that needed higher precision calculation. So, I downloaded your hpf class. Without modifying my code at all, just modifying the class of the input parameters, my regression tests worked as expected. That is how a class like this is supposed to work.

Thank you!

JeffDavid WilsonJohn D'ErricoHi Michael,

This is not really intended as a replacement for VPI, as the VPI tool concentrates entirely on large integers. In HPF you will see I've not given you tools to factor integers, or to test for primality, modular inverses, modular roots, or a powermod function.

Here I've concentrated on numerical methods, so I've provided computations for trig functions, exponentials, logs, erf, etc. HPF might be used to get around the dynamic range limitation of a double, or for someone who desperately wants to work in a given number of digits. It would seem to be a good tool for a student to learn about numerical precision problems. As well, I've found HPF to be quite useful in showing what MATLAB stores when we represent a number as an IEEE 754 double. (My own version of a tool like num2strexact that need not be compiled.) And it is a nice tool to test an algorithm to determine if a problem is due to a lack of numerical precision. I've even heard of a version of chol that will run with HPF numbers.

In general, I'd expect that most users of HPF will work with numbers in the range of 30 to perhaps 50 or so digits at most, except for those individuals who love to play in the huge digit sandbox.

As far as VPI goes, I plan on offering a new version to give a speed bump. (When I do, I'd like to play with some quadratic sieve factoring schemes to update factor.) I'd been debating in my own mind if I should leave the old VPI up there or not, for those users who have older releases of MATLAB.

By the way, that speed bump for VPI can come from one of two sources. For example, I could re-write VPI as a wrapper class for an existing Big Integer tool. I looked at such a question when I wrote HPF. (In fact, I wrote HPF in three separate versions, choosing the migit version as it ran the fastest.) Alternatively, I could in theory write the new VPI version to use migits, much like HPF uses now. The problem with migits is, they work well for numbers where the length of the mantissa does not change.

The issue is there are limits on the size of a number that can be operated on with migits, and integer arithmetic can cause integers to grow arbitrarily large. As such I would be forced to limit the size of the migits to be rather small. Thus with 3-migits, one can generate integers up to a few hundred millions of decimal digits. That would surely be adequate for the conceivable future, but perhaps one day it might be a limitation. (Somebody always wants to push the limits of any computational tool.)

Finally, a migital version of VPI would be simpler (and more efficient) to write than my HPF implementation, since it would not be forced to work in a fixed umber of digits. And there is no need for guard digits when you work with pure integers. (Is migital a word? I guess it is now.)

I'd be willing to take any advice people have to offer on these matters. Vote now.

John

MichaelJohn can you comment to existing VPI users on whether this is an official successor to VPI and if/when we should make the switch?

peteras always, excellent job from John

Christophe Lauwerys