Code covered by the BSD License

### Highlights fromHPF - a big decimal class

5.0

5.0 | 6 ratings Rate this file 41 Downloads (last 30 days) File Size: 1.39 MB File ID: #36534

# HPF - a big decimal class

04 May 2012 (Updated 10 Dec 2012)

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

File Information
Description

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.

Acknowledgements
Required Products MATLAB
MATLAB release MATLAB 7.14 (R2012a)
Other requirements HPF should work in any MATLAB release that employs the classdef capability for user defined classes.
Tags for This File
Everyone's Tags
Tags I've Applied
10 Sep 2012

This 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.

% (a plain text file) and then import to SSHData.mat
% Derek O'Connor 10 Sep 2012

[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

07 Jun 2012

This 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!

15 May 2012
12 May 2012
11 May 2012

Hi 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

11 May 2012

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

04 May 2012

as always, excellent job from John

04 May 2012
07 May 2012

Many repairs & updates - moved hpf.m into @hpf to enable new methods to be added easily.

22 May 2012

Added static methods, ones, zeros, eye, ten. A few minor doc changes.

11 Sep 2012

New .pdf file provided by Derek O'Connor.

12 Sep 2012

Further updates to the HPF pdf file

10 Dec 2012