Code covered by the BSD License  

Highlights from
HPF - a big decimal class

HPF - a big decimal class

by

John D'Errico (view profile)

 

04 May 2012 (Updated )

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

Demo script for my High Precision Floating point decimal

Demo script for my High Precision Floating point decimal

format. The user can specify any number of digits to be carried, while doing a variety of different numerical computations on these numbers. Not all MATLAB operators are defined, as my main goal here was merely to build a general tool for this purpose while learning to use MATLAB's OOP facilities. As well, it seems a useful tool to learn how one MIGHT accomplish the goal of a "big" float format.

Considerably more information is provided in the document HPF.pdf, where I provide many details of the computational methods employed in HPF.

Author: John D'Errico e-mail: woodchips@rochester.rr.com

Contents

Create a hpf number at the command line

Note that you don't need to specify the number of digits in advance. HPF has a default of 66 total digits, 2 of which are shadowed, carried as guard digits, but not reported.

DefaultNumberOfDigits 64 2
F = hpf(12)
F =
12

There are 66 digits carried in this format by default, two of which are generally hidden from view.

F.NumberOfDigits
ans =
    64     6

Note that the number 12 is an integer in MATLAB, so it is exactly convertable in the hpf format. Next, convert a true floating point number into hpf form.

F = hpf(1.2)
F =
1.199999999999999955591079014993738383054733276367187500000000000

As we see here, the number was not truly exactly 1.2. This happens because I let MATLAB create the value 1.2, which is then passed into hpf and then converted into the decimal form for that number from the hex form carried internally. The digits here are the exact representation of that number as MATLAB has stored it in the ieee form. In fact, the last 40 digits are essentially random floating point trash. Also see that this number has a few trailing zeros.

mantissa(F)
ans =
  Columns 1 through 13
     1     1     9     9     9     9     9     9     9     9     9     9     9
  Columns 14 through 26
     9     9     9     9     5     5     5     9     1     0     7     9     0
  Columns 27 through 39
     1     4     9     9     3     7     3     8     3     8     3     0     5
  Columns 40 through 52
     4     7     3     3     2     7     6     3     6     7     1     8     7
  Columns 53 through 65
     5     0     0     0     0     0     0     0     0     0     0     0     0
  Columns 66 through 70
     0     0     0     0     0

Had we specified F in a different way, we could have generated the exact decimal form. Thus

F = hpf('1.2')
F =
1.2

This number is the exact representation of the desired decimal value. A nice feature of HPF is that it will parse anumber with a long string of decimal digits correctly. Thus we can provide a number with many digits and expect that hpf will store all of them exactly as we desire.

F = hpf('123.45678909876543210123456789098765432101234567890987654321012345',[100 0])
F =
123.45678909876543210123456789098765432101234567890987654321012345

since we requested exactly 100 digits in the result, there are zeros appended to the end of the digit string.

mantissa(F)
ans =
  Columns 1 through 13
     1     2     3     4     5     6     7     8     9     0     9     8     7
  Columns 14 through 26
     6     5     4     3     2     1     0     1     2     3     4     5     6
  Columns 27 through 39
     7     8     9     0     9     8     7     6     5     4     3     2     1
  Columns 40 through 52
     0     1     2     3     4     5     6     7     8     9     0     9     8
  Columns 53 through 65
     7     6     5     4     3     2     1     0     1     2     3     4     5
  Columns 66 through 78
     0     0     0     0     0     0     0     0     0     0     0     0     0
  Columns 79 through 91
     0     0     0     0     0     0     0     0     0     0     0     0     0
  Columns 92 through 100
     0     0     0     0     0     0     0     0     0

Be careful though, since if you provide many digits, but then indicate that only a few digits be stored, then I'll do as you tell me, truncating the extra digits from that number.

F = hpf('123.45678909876543210123456789098765432101234567890987654321012345',[10 0])
mantissa(F)
F =
123.4567890
ans =
     1     2     3     4     5     6     7     8     9     0

Note that the digits of F are stored as Migits, with an initial default set as 4-migits

DefaultDecimalBase
ans =
     5

We can see the migits themselves

F.Migits
ans =
       12345       67890

Numbers are stored in a floating point form, with a separate sign and exponent from the Mantissa. The sign is stored as either -1, 0 or 1. We can extract the sign directly.

F.Sign
ans =
     1

Extracting the exponent is also easily done.

F.Exponent
ans =
                    3

That exponent is stored as an int64 integer, so numbers as large in magnitude as the largest int64 number can be formed and manipulated.

intmax('int64')
ans =
  9223372036854775807

HPF does include both Inf and NaN as permissable values.

hpf(inf)
hpf(nan)
ans =
  Inf
ans =
  NaN

In some cases a computation can have a problem. So just as 0/0 produces a NaN in matlab, it does so with HPF numbers too. The latter two examples here create HPF NaNs.

0/0
hpf(0)/0
hpf(0)/hpf(0)
ans =
   NaN
ans =
  NaN
ans =
  NaN

Inf is as easily generated, even by some functions. As large as the exponent in any floating point form can be, we can always overflow that value if we try hard enough.

1/0
hpf(1)/0
cot(hpf(0))
tan(hpf('pi',[64 1])/2)
exp(hpf('1e1000'))
ans =
   Inf
ans =
  Inf
ans =
  Inf
ans =
  Inf
ans =
  Inf

Most numeric classes are currently supported for conversion into a hpf format, with the exception of complex numbers. (Sorry, perhaps in a future version.)

Special numbers in hpf

The numbers pi and e are useful enough that I'm willing to store the true value for those digits, up to a realistic limit. Thus, there are 100000 digits available for pi, and 100000 digits stored away for e.

Here I will define a 100 digit version of the number pi. Of course, I cannot simply type hpf(pi), as that would convert the double precision version of pi that MATLAB supplies into a hpf. Instead, HPF looks for 'pi' and 'e' as keys here.

PI = hpf('pi',100)
PI =
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068

Properties of a hpf number

properties('hpf')
Properties for class hpf:
    NumberOfDigits
    DecimalBase
    Base
    Numeric
    Sign
    Exponent
    Migits

Methods currently defined for hpf numbers.

(This will grow with time of course.)
methods('hpf')
Methods for class hpf:

abs                csch               isinf              round              
acos               cubrt              isnan              roundn             
acosd              cumprod            isnumeric          sec                
acosh              cumsum             iszero             secd               
acot               det                le                 sech               
acotd              diag               linspace           sign               
acsc               disp               log                sin                
acscd              display            log10              sind               
adjustdecimalbase  double             log2               single             
asec               eps                lt                 sinh               
asecd              eq                 lu                 sort               
asin               erf                mantissa           sqrt               
asind              exp                max                sum                
asinh              factorial          min                tan                
atan               fix                minus              tand               
atand              floor              mod                tanh               
atanh              fminsearch         mpower             times              
augmentdigits      fractionalpart     mrdivide           uint16             
ceil               full               mtimes             uint32             
chol               ge                 ne                 uint64             
cos                gt                 nthroot            uint8              
cosd               hpf                num2str            uminus             
cosh               int16              plus               uplus              
cot                int2str            power              vpi                
cotd               int32              prod               
coth               int64              rat                
csc                int8               rdivide            
cscd               isfinite           reciprocal         

Static methods:

eye                ten                
ones               zeros              

Arithmetic operations are defined on a hpf number

2*PI
ans =
6.283185307179586476925286766559005768394338798750211641949889184615632812572417997256069650684234136

So if we subtract the value of 2*pi, as defined by MATLAB, we expect a result that is on the order of eps.

2*PI - 2*pi
ans =
0.0000000000000002449293598294706354452131864550002116419498891846156328125724179972560696506842341359642800000000000

Many of the standard functions in mathematics have been defined for hpf numbers

For example, I could get 100 digits in the value of e as simply

e_100 = hpf('e',100)
e_100 =
2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427

But I could have also found that value as

exp(hpf(1,100))
ans =
2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427

Changing the number of digits in an hpf number

We can specify the number of digits to be carried in an hpf number in a variety of ways.

First, I'll define 1/3 as a hpf number, with 100 decimal digits shown. in fact, there will be 104 digits stored, the last four of which are simply there to be conservative, for an extra bit of accuracy.

F = hpf('1',100)/3
F.Migits
F.NumberOfDigits
F =
0.3333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333333
ans =
  Columns 1 through 6
       33333       33333       33333       33333       33333       33333
  Columns 7 through 12
       33333       33333       33333       33333       33333       33333
  Columns 13 through 18
       33333       33333       33333       33333       33333       33333
  Columns 19 through 21
       33333       33333       33333
ans =
   100     5

Reduce the number of digits in F

F = augmentdigits(F,[50,2])
F.NumberOfDigits
F =
0.33333333333333333333333333333333333333333333333333
ans =
    50     5

Of course, we can increase the number of digits stored again, but the information in those truncated digits is lost. Future computations will now be carried out with the new precision.

F = augmentdigits(F,[250,4])
F.NumberOfDigits
F.Migits
F =
0.3333333333333333333333333333333333333333333333333333333
ans =
   250     5
ans =
  Columns 1 through 6
       33333       33333       33333       33333       33333       33333
  Columns 7 through 12
       33333       33333       33333       33333       33333           0
  Columns 13 through 18
           0           0           0           0           0           0
  Columns 19 through 24
           0           0           0           0           0           0
  Columns 25 through 30
           0           0           0           0           0           0
  Columns 31 through 36
           0           0           0           0           0           0
  Columns 37 through 42
           0           0           0           0           0           0
  Columns 43 through 48
           0           0           0           0           0           0
  Columns 49 through 51
           0           0           0

Trig functions

The standard trigonometric functions are all defined (although I've probably missed your favorite special function as there are so many. I seem to be constantly adding new functions as I notice that one is missing.)

The sine of 30 degrees = pi/6 radians is 1/2

sin(hpf('pi',200)/6)
sind(hpf('30',1000))
ans =
0.5
ans =
0.5

The cosine of 45 degrees = pi/4 radians is sqrt(2)/2

C = cosd(hpf('45',500))
% and of course, the square of that number should be 1/2
C.^2
C =
0.70710678118654752440084436210484903928483593768847403658833986899536623923105351942519376716382078636750692311545614851246241802792536860632206074854996791570661133296375279637789997525057639103028573505477998580298513726729843100736425870932044459930477616461524215435716072541988130181399762570399484362669827316590441482031030762917619752737287514387998086491778761016876592850567718730170424942358019344998534950240751527201389515822712391153424646845931079028923155579833435650650780928449361862
ans =
0.5

The tangent of 45 degrees = pi/4 radians is 1

tand(hpf('pi',1000)/4)
atand(hpf('1',1000))
ans =
0.01370864253439405605227298869790930363121380384841842673544536711960521807164775059427280020891113254561562378582045422164719174631109765713950037630957898860802480020923765874675922829584873864221797843978765892606013704345929094437512442369399787095125097994006840738875250901794241128544100024806068605721788387338083182793243192696944371011114307326529305276957593055314700017836809662218854752179847180761861543376311460311997654051066068212362010944818083667667870312358230187038986881975475528319227379662049032167313596293521706598712408936830751971770520631033465064783718752512525241556030042759546136161110199432687468919207329256921678814504836256955001984830490726781332113305296750494156109972970571233516623036276312379019005767225536128258770802849775452059663312639073442324922770841878513726292896357741375245258143298927264569637440177459635257380514260980739539297643828167434020129368034683562815450708160231110620759427376454042367056840328205686160894406218718726859525538377348
ans =
45.00000000000000000000000000000061126127699329378291627080405353152961002069794566320284014987471586052921522972883322289189941252744527347380985012593057152946666228033261577480426283714271065393213449833936939677437166541523596761541420694844119118451822347307618593269157362614758226121617600011016392320700428591759140538802010327262537812895257965616974942371276177308128567919439465090686362505853869732633506914223762428271613477563266694960330121597931964084357505973850173535781669912167706917407337403560066173710033526063926319316534276889372163811835099753296679926834353991235661061045441264751144652206643842483595530083773086694052909636811439835155571285930860308793060637890021090891380636492407776295192078845422380592318320307672367231862136068772502805614320021785865993228945458469309758700596486223048360859093858162813729399889969337289369183726107403360333791009804296254519712364059904368532189367831417474788343677658152703552888266673618731191211537941671369877978713250536

Basic arithmetic is provided, and much more.

As long as one term in the expression is an HPF number, the computation is done with an HPF result. Be careful though. Here I forced 5 to be an hpf number, with the current DefaultNumberOfDigits, because then I took a square root of 5. I can add an integer to that, since HPF converts 1 on the fly into HPF.

DefaultNumberOfDigits 64 4
phi = (1 + sqrt(hpf(5)))/2
phi =
1.618033988749894848204586834365634158823368255109555301291650298

Beware however, that while the following result will be in HPF form, with the desired number of digits, that the result will be incorrect past the 16th digit or so, because here sqrt(5) is computed as a double precision number, and only then converted to an HPF form.

(hpf(1) + sqrt(5))/2
ans =
1.618033988749894902525738871190696954727172851562500000000000000

If you want to be confident that the number computed is accurate, then do it a second time, but with more digits of precision, using a few more spare digits in the number. Here, I'll subtract the more accurate result from phi. All the reported digits are identical to those from our original computation.

phi - (1 + sqrt(hpf(5,[64 10])))/2
ans =
    0

A test of some identities

Working in 1000 decimal digits, with no shadow digits to hide any flaws...

DefaultNumberOfDigits 1000 0
x = sqrt(hpf(rand))
x =
0.35635209595778450416563669407896462587713144320364254245988443923232647798600668882384640419170789466741025019641200123731418935451583012300799336534546783620886726673640409688350916789668607299430414966561066515151300071462259211330092512071132659912109375

Should be 1

sin(x).^2 + cos(x).^2
ans =
0.9999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999996

Should be -2

(log(cosh(x)/exp(x) - 1/2) + log(hpf(2)))/x
ans =
-2.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000004

Computing pi

As I've shown, I already store the value of pi internally in HPF to a rather large number of digits. But, one could use HPF to compute pi yourself. Here I'll use Viete's formula to generate 64 digits of pi. http://en.wikipedia.org/wiki/Pi

a = sqrt(hpf(2,[64 3]));
b = a/2;
tol = 10*eps(b)
niter = 0;
while 1
  niter = niter + 1;
  a = sqrt(2 + a);
  b0 = b;
  b = b*a/2;

  if abs(b - b0) <= tol
    break
  end
end
piest = 2/b
tol =
1.e-69
piest =
3.141592653589793238462643383281337594640838857701917327089225014

It took niter iterations before it had converged to my specified tolerance.

niter
niter =
   114

As you can see, the approximation yields the same digits as those I have stored for pi.

hpf('pi',64)
ans =
3.141592653589793238462643383279502884197169399375105820974944592

Another choice might be Chudnovski's formula: http://en.wikipedia.org/wiki/Chudnovsky_algorithm

DefaultNumberOfDigits 5000 10
DefaultDecimalBase 5
a1 = hpf(13591409); % small integers, so exact
a2 = hpf(545140134);
a3 = hpf(640320);
a4 = a3.*a3.*a3;
num = hpf(1);
den = sqrt(a4);
tol = eps(num);
piinv = hpf(0);
k = 0;
while 1

  piterm = num.*(a1+a2.*k)./den;
  piinv = piinv + piterm;

  if abs(piterm) < tol
    break
  end

  k = k + 1;
  num = num.*prod(hpf(6*k+(-5:0)));
  den = den.*a4.*(3*k-2).*(3*k-1).*(3*k).*hpf(k).^3;

  num = -num;

end
piest = reciprocal(12.*piinv)
piest =
3.1415926535897932384626433832793909304074940369500141177724882830002667582175677534238500349122595480051550203519998032783479590370264426370833005237512528693652545212823962637546181312689885611034582341403118413868252540236632365314377422474157438814781175156133543630619768306671517119588104634503781851570210541885221432240484817565938951252388200066329475978083991623403328565758927052102881260001431839872946181084932609454905832715472143706767923321430146686786739354427444764254713849040522793379140027577820849235788652223132965057543889604386531005970894426308513794096965069183274365193023857508441010964343904364272134602548145155763129743974517318238901293156295689299600905565000531440272263479619052100042622789767910812753272533062668899788815419526581979802628287376273673678698142045936909914452642278994236613222963254859934900481793562809094399062916691578743888036362139182637860557845546984341923751244757686981701424169909602146615315722506505296285725910134635838143671352405930179121462940961619088359359590468600325771551614362378865518548404163186522733883576905762099953879157911963730870298902941589173492099444071912204144489177765166716184774592758984317402041270031592104908836762334867461228814146520436359175394701333029043327046710713444200032732837249452235255538978431241127920369378417333800372763557442051478877416277388623856271405632167925274554395608316355696022009358304975304420693258131472304989358654163436328452296342690100852300273995074398598077718487156180560635132143670735367399011224307015200431526637311459066091054241751260623673153265564632294078541298289341910532887585140926616444624022300099272373263019611956697782054220149223870079769394266750197443165839257944819272773780578514832274071290174237770846206760403197876566256199070096625830769052377406164405750177099156461767577326128757602959597959030511841184865371682939650228386518506518097130061602940217425866266085122879869172707458319561840914924659063906711082376286165737812519108891661540285386642450959950324237641449405417809720346174624898813119802847205858827410198456038767517719645353006186137127188382468686737919963423996265324691923251878373162338390351559817127347828399207543309425187844711387257596072094615388289330860099970878233050901714345524313283466391786485644602878954570058576533480090631966858900752752635228031236516534442727028113756841443064952217654097337251306459159013117235603527549845836037787650425760858515894769102957544452777212803967425130874984441926242383015248793359081735556436575765301481633588907421671058396894588567950745723814103789254902448833043327066541839321531944049982869458738741148014888623180661246910771033574907026970710902891926861466867617671769599979009963039748374783947377939139347839905000195930464802403954695904493574614859871092840160144788354952542411098239009129960292686287391201934737209913996687432245266303156498046790020640525376163096570097625967410257896612016576511082018826099081062256278089813695566018112723842779997993167778864747997367250177300117873298369830258608671526466017020214333915589868697440611043996513341542852038499158620486662004869551448060034543684427440240951754268883156304842795201361167505424916365859690784536654020214757927481021631003436746002290395089847890790961641391885451987575452914486636916100307580944858430029451919766531546724436367184145670339496115106297638694286806152677506777661580099491137400910406899139705991100980822374108175644321266728042542743932805386552999778971003216154225736759976072522523239354727020499724969116939131011448887403402819611106792307607465029801293077657657119517797598446506617681505190464563166166528279369325602825852896965436734171138406509449856585117537926773091109465944922594728770638986122204917644854708348273150851094946244927266872988786304606248600324356990231801886682835359936683191730727401704286561263068471983681253033067015540122953847416465229360047933677775259822972822356739742436798823514824268731647855855891418658109272895166718227588884676493618620758008394826844403189786302315925891292576063970434298224799431953201276403231803300382510607363842312894469194428278866574766506172827857337275331570641734648776443444117292084779474095756946476898990342210162091890562209643485408363473350031868288091502019046456363574759052437943237158508156185363726085017688799793405503418077755086797061272115300625954591858209108650082224930005218322117299749337237211122066220839135183735649512122087043444934843948741401728773139864621833724091166496561889219644661881746554681409992762982921708499365532960476373372068877008741488928358870411437428199207066854130188439425924966466269506537381395842942665206286791247225636862249816530271193165138075349384092470611668509567768338218289053715282075043667965164873317687641163488956525489210842957402488275188953189881677939806455838017138543281669343154542608448032804552270331251213596720957871699614196874063205096518823720041120523433665938331042259024003284899174

The numner of iterations required for this accuracy is

k
k =
   354

So roughly 5000/k digits per iteration

5000/k
ans =
       14.124

With an error of:

piest - hpf('pi')
ans =
-0.00000000000000000000000000000011195378967536242509170320245630930754964806864124520418479042985751997699306616128250336874588557252413959464205888437722824808502958141954225735094151335396039382692373028849913427910819210481241170234904091785545802766753105131010598554847749598106168164844978569089553942998555212706603165747243920623535896647654436125964276220293136820848760994549144174165338560928987328307114148345983124062679866106457873444175574231936529407089341644300824846575655327895966979191967060465857171950715617263309871898280174174120784310568109658632491376705587775649103115010274624972402604882917039135855788250297509359732420432401670907733482411451698065541096023577963611489631567199420691901769858505762707853830252146504307053982566939734180039271522121343236436494628211400133346801008268040321886098595701919904177188300881472909369793246121742965881957778440325940536124781369879199456309360668379416571221137060279174312032620828152069717804041091441302738154208117960592018311090696932395444292997710652179539107514051579741530111415535730394714050289203119907348374359978366086936712080565667061576765178344869785468393283686094826770077760112502926993093488284566480236747792733022091398182876729091196028140729770414654392267510365550268404426136248310146073351559345018956423695221426301878201968966831383300230576962564356764996049968427968073645875971177263435319458141895113847889965231540677735581609442076663092286152020423648450335689959618348735098973745423970310928147398838888669974076437668918031475460090263866151993342243837669235520735756442130956439222751723108195419590341470312211070567903974355862092403485539135646369753642745351613631593045471665952801810027859318908755828472116137070104514312523946389470587392561314701733200161288477642320130692273023606501352593617387071477392233161224944482143318899303294437014040843594800026218355881053449703127246174517320246049710527467296944710117145551158056608149650750019472440836719681669677567921029961005041260522330235756328606167767998070843569210201577860128949670618205425365341760970119145172428791644878103771871097094379057714522535406188917738919508541004864353913646377418004119368529834792953800025537152071719948634674996951114359159483599018582139075848240218379129775248542940267942410255578713069210943278865249475010762485435256663857121690240413365061069654168418215293669061877229171209353636615595533489938353796635088796278141037074452687797593482301122031431323782566123661407088344524114957171675811705532406723556913898870386985819839299566361915907188584630237950744601309876483160534327860718161762869872478653473712507460785731599524594814578647802902215074882667514407535017229810850410064409114345373943554827594877136434205560144745245246534414062990536993399983057694225602020129452304771955195255423737985942053271783948331476731002497947005708133662739015404460362838651893788800888204711909743088114838269240251879673336481135993206523458501412113192581693763015029556093009534723344734512908219341707603212606862680195476122217191714547564470074917104067463167744196363879640039491740156475006315283124236377737385999867806564544668639115630776133582342151674842488875756416716041219449628271865054998331014896197701062753168021000025053088807860425619868139080847301347172769573026309438902604708518888099414921416690675477813041459743665996376006595822232803105702334749132537000740155753476728159726955033628647442987094709757255608247161675247664018657165882673323989054423621725782597255789891619320829305978518315397321607675883749860518804967439696868777514672320536032436716737790528444483060840914733089750340314502052744304506174964190366787473921717520948544635472642592181947281207328024666525410725770654024420166254027376356675502512235659341865206908269004552444028629626015886633401483957344182650893565319691117174579333832776892183758475542735490644340367343085860254278103263783097329979671550552508324146685736479997378919323210044162557689780074638196519438396261745821441594912757689102257379249093808000543042596897147175426283305629543051015353021847333450642321119060730324234849827163664172507875518875791333155869813114381808148794192305589346592992592376533989472591454063996034029786193794418689038472910273411966145905198812062477841779650130979306144148521721494178825908455527348532213561521123959247461754888800300494148838380004516988568236790302465622952530609969152105545338892333940511694562114834499852533940740225156658476176609399719258079491323838655354250259206424189781582464127114794501965547692268586095880690299495029760526813924474915602217800900987795555175914031682907675810211203600121131053035979632273771466531720557944716473227961678005740164582542346958046063792330185873846132310057582421681928552541271052785300328889648120859041165441635613566739940755919645832753145294987681917418379403860244319387107090699544403379368985888866176334821440040137275028789816993003268540081766929675658408646841283612980398128053000000000000000000000

Of course, we can use the rat function to give a simple rational fraction approximation for pi also. Perhaps you have tried rat out for doubles.

[N,D] = rat(pi)
N =
   355
D =
   113

Rat now works on HPF numbers, and as you can see, here it yields a rational approximation for pi accurate to 100 digits. Its not terribly fast, since there are a lot of divides. Sorry about that, but 100 digits is no problem at all.

[N,D] = rat(hpf('pi',100),hpf('1e-101'))
hpf('pi',100)
N/D
N =
394372834342725903069943709807632345074473102456264
D =
125532772013612015195543173729505082616186012726141
ans =
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068
ans =
3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068

A timing test

Time for a multiply in HPF should be roughly inversely quadratic although 35000 digits is not enough to really stress the system

DefaultNumberOfDigits 35000 0
DefaultDecimalBase 1
tic,hpf('pi')*hpf('e');toc

DefaultDecimalBase 2
tic,hpf('pi')*hpf('e');toc

DefaultDecimalBase 3
tic,hpf('pi')*hpf('e');toc

DefaultDecimalBase 4
tic,hpf('pi')*hpf('e');toc

DefaultDecimalBase 5
tic,hpf('pi')*hpf('e');toc

DefaultDecimalBase 6
tic,hpf('pi')*hpf('e');toc
Elapsed time is 0.227864 seconds.
Elapsed time is 0.060248 seconds.
Elapsed time is 0.045643 seconds.
Elapsed time is 0.027505 seconds.
Elapsed time is 0.024214 seconds.
Elapsed time is 0.024402 seconds.

Linear algebraic tools: chol, LU, det

Start with a Hilbert matrix of order 5

N = 5;
k = repmat(hpf(1:N,100),N,1);
A = reciprocal(k + k' - 1);

Cholesky factor

L = chol(A);

Max error

max(max(abs(A - L'*L)))
ans =
5.93560228369665931268754727907941298465405429655213843281917621573723e-33

LU factorization

[L,U,P] = lu(A);
max(max(abs(P*A - L*U)))
ans =
1.e-102

Determinant. For this one, lets create a little higher order matrix. Note that the matrix is NOT truly singular, but the deterinant is vanishingly tiny. Had I used too high of a order for this matrix, we would have seen problems even in 100 digits of precision.

N = 25;
k = repmat(hpf(1:N,100),N,1);
A = reciprocal(k + k' - 1);
det(A)
ans =
-1.339885345032204390143551076857085237546956006581482557366597988101633171274685669261626122466988338e-357

Use, and "Abuse" of HPF

Of course, you can do many interesting things with the HPF toolbox. Above I used it to compute pi to many digits of precision. You can also use this tool to do computations to many digits of precision on numbers that you know inexactly, an abuse of both those numbers and this toolbox. As well, you can also use this toolbox to do computations that are better done using well thought out numerical methods. Good numerical analysis will always trump brute force computation. But, at times brute force can be too much of a temptation.

For example, from a recent Project Euler problem (318), how many leading nines are there in the fractional part of (sqrt(2) + sqrt(3)).^n, where n is a large even number?

Thus, for various values of n...

DefaultNumberOfDigits 2000 3 session
k = sqrt(hpf(2)) + sqrt(hpf(3));
for n = [2 4 8 16 32 64 128 256 512 1024]
  f = fractionalpart(k.^n)
  % How many 9's does that fractional part begin with?
  find(mantissa(f) ~= 9,1,'first') - 1
end
f =
0.89897948556635619639456814941170967512964847039917284858587868124577889314176966169444681852915598558452077409123274046978323335955485891507629014207641528247435599958732976112872434022943130369278062789231143591485827043116338218714787978831933342229962919918843434509414525024340036907233297824859619140625
ans =
     0
f =
0.9897948556635619639456814941163804342516449309419118573183352454123299988050466163636589237776443590044182357597873202954285007996812168105398763076743547752788437066471125426758210093876510928058948101223203771369078318808707841140607029180252440700356725458889473477002521351042830173738704418570312642568701146431327022179014523251651550437894320090365806570906800083455239556017911247996781013889269948962658971805273640991424536444902768343406983077134241018991362880636387847731319132159190220637493752734003157787831228238602414188897337948713307944155354872279046888651488078636475620442070066928863525390625
ans =
     1
f =
0.99989585502907246667678642326346640202059464972498721482550430414880335712101187890182202075285883691322644487657043300573409155095902252650905563503096152453636935703788403372954235965591695965212737472742324476138449976189744579252256546336897457828660294878964398774092858641853273072402005751783133111974892014984158055653411969755223144969650165868620532327912934654085384390108808596402365596697135320991339516891419131002482728919202898478128305579192392224895256953255483390166068606189395253712359404752341185739112044876121252001869560305291771395380713186100080088405104159463589480768538718620742168122660094314878776199638818517313109817372794144656809168641792997410959815236601089408069264169844943300711850408437100279821355301105298268973810824996429002899289894028566754610693268153653037882405397670473860007937835821417582992048672759572691499407539240223491762766049976398517131295520558076903022341392954276424172276819105047686994087545164685239211130524079438814691371059228180156147719133552296344061827882977608803740506253345593154749277448774450358644595732934455054385644742250301486025361591968713752654883951865441536496527827256391398060385453930677143071490721448757454936639987863600254058837890625
ans =
     3
f =
0.99999998915382503050323345208335571309999378373841032658331506526935899950368792002839666611742452268792868920428877055693521384423257770932922164404707314916213148887981860710735987034418377245911741505782467024053437923653210591539826102589892741621219766168873840580288651551123272931071445077689341228963791733578230268368862078006562791257128683016258939102373556615988394888257819992708938367176037361035868885889591155296488947728105660303443904709840254423729065718846307870252995049874614456243088804895002502481515058757235086038249880831988274259483102197345456909352053964606324990770024036150206064557090194683354120380419321856413062043025071595193640286047766572861333956731748327193628387091470622377525418265491138104971563284367938617807339485254947827801788326890862203265691403981524577989811248562240904838497879443769269082944276403151164474041208933159873341662498387548259484610829339233234086781760967941428149693400263424929942113471416962712202769399329689954817890347631864920425546535168766643860774372661845257268912778075308608262772348115368338566021236360871062239634150457075607228326272988793146458693510102392544624160956685025028286294348378091151170649929334674927320199573762009536125134541668621048815635237237905586612658249943517731306580005385957644889554690879795391374589675205743589892965501416506104073316054983386664851130485872134404403131101139424587412545284688588114401435812243317608653629847525051601803620265838994218576186042459195134535753913171533211862180736820337053690589153105125205207622050396947512728512026457000036792212913053766683918576630187607991250338894827182011222766279482137369584510478012186959529070678790771151270901360523907767876688487285360153918568124950465287658644862408732653643940477335901772988632943297884615554141818854158704544713294553199880333316530139047026075890667625854830918921719788218568716252754126772201877663887235946328044438438000718707659499045896682730401286199968195252925913924875642861210000
ans =
     7
f =
0.99999999999999887786956940402896454228809976549498323958385413592701454013587092285638730298942863800626761678151507577786227333017795171212448988159568756331004448353049755674842654961380637209339288254809956794478849551703699149917847633135430698072688400835228699864854009278912980331625388820409601439214584705167391346490688043534101643512425589843422878385689891227320746420923768488983034587169668846353692415716166694676520074833263563984472326202005657216291835481870921908800245836828584759888304607765595574894115008550324812273942693194513692497165769675352456361456646512515014613590008695330534578215263482438706489918615418601702686754718626370482199432546824085045677683709947229232852050843355417959348576936170711598956451003155110108751132515446458738158772354943624511240513616624170418813441323727092371045880775581187836014398313317421851034802895541334680473803266680636919312385737485541371770213396378042366544905399811419598139512878165971412850997381938096571082420451132180751209261963042504307510453874881078097005577834926108616494890791612528529007287389317753357735919862846022195528456319834586035125238715245536375793462333194851355568310828138219050771685803384466086927615254728149643307397300972574563746358050610527074152884514414107340190928912638853925788738999432371676534700943452950162710444358803317475675851268535116766727075954575031748106004174754506185098066780397321022632108316161474402991501884684690818230070708881080217268691037151370668477393703840573969143892433993016894874757278740466634679612774099739206371480255870439233955398381878930288256828163335627009917258098884999225640920137024989470340963317555387899536800550867489617343973344505637019451070727310870312951868575209553499074149860726293832745564212031773684915889114573475469973322716869716413355977663834025838373249876793916953482184051103054838339087503158598391911974919899761627703995255271900295843605416261972673480466488949435410021837450141462585971752379576000000000000
ans =
    14
f =
0.92255881405522260521823021687830463937216218481217210207014336947428099253623053489857751426708977364849372070278062664911918395336870098729644363516772318012127530034775618380105323664669138280489011238374709343764447762868725859553697051405573569642134164060138793783059567455057291253218595876799665911011893408424786907609811184492828941580383064605057634433835231618701033448911054332930137067007203644933010693505577463132898101408777156192565635362770058875504674922458769005235894954833497164413862765893256684547681290967340494358013015659425703203102626368603716395299924637068076734432600359982450619797043755158063365778487247840301970526118918818309986892635531639355582822706522766241120012663557866448443731828472027447333985391175033859929140455281084108037050638397693942559985585447298740663007462453425664582851249722034102069889460182518007890665513829078555866319643086737651922130537405084000496595290675777976438474320210971911017770332658687688979314952003821182546252567562858983074420930569582396729458554694195035911417453673500277953197281981582248755015809737256901097030733323859002725626462077013397547263173729701846112906834508029967077209294554608313904463353765363869292864587024332281345120355361711934282332290540975511640301154588267649591404186742717617037032831959947054959225184096129361823124769804676246547222649947654583639309976589733018318878232821417606216607473512834306884407323575965986026296936676801760801399011593073544588091423473087518717242873811388002785487299949506904473424931835023955032299980452578440071371833339372588233023320553030865279464898964806915139781157900252692024664583261490364035729676589725287514213975805844850883591788195916101810407185365444186479930524773971543460678474537961226951916089151860724517102877904289541211850969809454443780587028104085369662127707119142480490450893954674787822651121164907323423120790691861522832346577974071149874827323309731045178198146869882093144011981181560000000000000000000000000000
ans =
     1
f =
0.84838511579154507028137197233380596132950297757530524917256051841381686676363009679676643375422667846605018537262244943964796064341937267737412051290900859959698957787077296167218164804826222869958912925003989188804671248155091408499307712105915700656466436820524764802547208204386510337310883422483939058199447321451757536738300261656449836131003871243664628138916814826985408579988845178128728917410051851047601676966233790390484412112893949777813302558088383500385596430955004775223028599949439821261146467821930227149862498381624733590870318521778825011469792456922290186667943608961716221640246269220905692016164311253704924226872224765682007859832980428580437594645935044636310086290603100341494564908501255141336661858865571236013448910112040357940244022712247886194101583589963366623913446283471489657730526610807121757823859299914840076667925416886423897461898385002049454172231132767791638088685302277433882938266624379792922119104561170776874859175664762449724147709155157571055921114843885202391950503554031969904199773893990462607803740237115814729001629210846450613891658317640714771207896506870029644106368424013678114887192977401332053644236789162553353373101910653977423530121027858201636543801688635829460730538277062543034482960179449266310175734099211444383328870900985257635037323467835244892209980922984929807002731032635150959799386080027279663628717652060876508131200016377222158809533412729155393765902039880835766827154289419168672608686824656687659327434066366255757908253416213645686001225664302297859498767756249232527760869807687086064338008161013281112609761770450151077613706711042010935408691610424064284826646800229815312187453178932539749151943895112305396689481100888079610172645208543973494078120338055633823410271876368186293455148214632579238633107063211061238366540065754244449619241828669773008531054856117893342022523056981848212756301886493555831354893309747104849476375385129471413149216039881714000000000000000000000000000000000000000000000000000000000000
ans =
     0
f =
0.33116915405007761467367272029938444598558730982855983039607839028087699376726417906599974568634880813836521370564520998256241179711801999788916456174265079625198504111430165133178432379650666477687974760140255959382094051641614608444334476356903369609403817787630278680589083278538008861050757246894018046748016421178409342518818474129772766849868141739695744270700948592627682380700585964724000104377117075314806308417568995704472909152707076054456595306716307473308080134242876719882112966682704031237346985106986975502678137756542712575726645348433191242720095021524620392731445347416864903684538163108416819922607844291002446871949703657172501489677419151791386211043629796490253438262825569657151206286428695642058802257457941844746809798810150288959947427268535845807857218574262052951771396274731339118323458956771254730291440470605169664406853294397350384880999841360538445724409619707426980503679112721966963389144476253180443029548148621483645466142086953022959693903750787493658777950687666480535192908211102944449847378254230195375365579643867405394923230720573554498880424945906012032044726823016497657110042980200740726380317231864376664482230961790458874344967097787825128326938917472893101585549546667748162011928459993975520164533984111224120084767612590421487343717760094266188682524860071278974318080985536194859590686409113766051843083327623393527079011592259300452236629969343664393963081142118459046574929163148519970640906918601712280380931054446080823406306367171581465363940343628770135688446194707114166071053427354895741925219659039317273571416658151559254294693703510249617713860532314533236308404340889149699166425559411759312057163673263708976303617980274337568583434666040152886549289370207405984484705776681999415643608234986798451939805121089776647099286438847585595059730541233783764867091175178497726315403098910672029653912590269283975534060000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
ans =
     0
f =
0.99165400522392012265781253424858004682312891222003893113885843353505209976276480932842970646230664492107776431446497783458207563213632359457156785761645822746366182976510831151260419365219698554169348109789930519936372737363239452020604998935127769172070102959879635687778685549582950298078922697532013331533902871876244717714533915447452157327300264185398676975117364476198476071090783535333336836200503920565586371675465455632030737117876048499414332974916191796145403123406819012066976407198179590902045894333341314683828296433882409250530952279180907611944728632694557425241446079095964827491000362839352839000599143381170574570289069622084175172986615366924227689299520527739017425856664141023012478254841538549784641468217863094697044001538376894646594588663872549503524001637857862486604841320714027861660740885294353986386566206933842378863029253019975169445826243179682724545547286757943299149506879513901057677615689104357937183704426052977035433726081545238817505639745546244507678447968704943821454257350496480956424418695522376747345323026102651246442576832879917750300123929280446098249355096972762061499455863117687034044007772291563027954370677998162977861680620626608218691916680817740901322240422003546597238372606610204701047301136306610258494759207129702899925164612300076554128483170127050442656795728902124831229586469396583562249864794984100185348418811171606679166483654631555696658033150823366465139682282330147445962848827879141616015580021682032882638979858287190995999021901903889884036877414500988877020634820090341134649650637194096723033672221177720446256838097755653007710313335385641413666699555134570753352065511025472703506886269536504351504980980926795922761608596020673210614389496881259252065202750722729561517600000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
ans =
     2
f =
0.398846026065970500309467997100544781302796629868038337100274252839438445834490857469262493847692008809807896610796113271209019792043510620311765648608027446528641043249511593184967554346008246172630337096512579287507220961745795654605571857714471033085668398762687094365833526167609276186019421901225210021589416731908085691877619986202344590214814373444224588113088753165453514965596242934751059778953182459435178069689655917599319453334683407142938449668774548995132159240277406442402174367747339351327434300352349536422980500666182528015661510214274710416136411789409468001559519515979667926618501567550991610707608333804165032646214279368608788346128977065398902015122696377908180627412149988263638249457437411555375845308071869448312997663393542796953922913214244253667251917785863084930009912829712043716929862720691110600439710047497737726102297876820457999833330055642399380148590691119473350932280342709643789453941597383727140823074765816203540107868971001220913907630316578338158681208778560634388505979886938800391297591694061293720347283759494909072697309396705499462779184051420832644310403205239357627310882812676290280278246158347042392447516021441042466929174874373164955443748974707018283798419278631249340750799064329981518995011933576991081067655424406185689969703490487116980824957827304979932398884066267740064726901089155046876362166604962666202503401527386644073900876044716110613323339564607820997415838791747549269950461853910965937548828707450902170363569918291367273
ans =
     0

The problem is brute force will fail you in this quest, unless you find the key. That, of course, is the essence of all Project Euler problems.

Contact us