Speed vs. floating point overflow

Is it possible to enable larger floating point range for calculation?
finser1.m (524 Bytes)

I don’t know what you mean by “larger floating point range”.

Octave (like the majority of numerical math libraries) implements support for single and double precision floating point numbers according to IEEE 754.

Tks. I mean the exponent size and maybe also accuracy.
(the case 4 in the attached macro file, that seems fastest, can not perform for n>5).

Well, like I wrote before: Octave uses single and double precision floating point numbers according to IEEE 754…

As @mmuetzel said, Octave and all other standard numeric packages use the IEEE 754 standard to represent numbers. This necessarily has a maximum value which can be represented. See the function realmax.

octave:1> realmax
ans = 1.797693134862316e+308

If you want to calculate something faster you will need to take advantage of symmetries or other special aspects of your problem. In this case, the code that you want to speed up is


This is taking the logarithm of the factorial function and there is an appropriate special function gammaln which can be used. The function returns the logarithm of the gamma function and can be combined with this identity

gamma (x + 1) = factorial (x)

Putting it all together I made a fifth case

b=n*h-gammaln (h+1);

Timings for n=12 (the most difficult in the file) are

n = 12
Elapsed time is 1.17888 seconds.
Elapsed time is 0.582969 seconds.
Elapsed time is 0.004076 seconds.
Elapsed time is 8.10623e-05 seconds.
Elapsed time is 7.86781e-06 seconds.

Method 5 is approximately 500X faster than the next best method (Method 3) and gets the correct answer.

I should checked about the gammaln(), that seems better answer directly the log(factorial) without getting into large numbers.
But case could be similar but with small difference like log(factorial(h)+1), etc and
there might not be special function to server each case.

Also I try find the pi, but this case also it gives only about 50 decimals:
piser1.m (469 Bytes)

If I wanted more accuracy to see how this works, how to?

Ordinary IEEE 754 variables have ~16 decimal digits of precision. If you want more, you need to use a different representation. I think the best choice here would be to use the Symbolic package available at Octave Forge. This package has functions for variable precision arithmetic so you can decide just how many digits you need on a per variable basis. See the function vpa at Function Reference: vpa.