Accurate Bernoulli numbers on the 41C, or "how close can you get"?

[Ángel Martin]

(03-19-2014 08:31 PM)Dieter Wrote:  "Fast and accurate" always sounds very promising. ;-) Do you have any figures here? How fast is, for instance, qf(0,01) on an actual HP41?

On V41 using the default settings: qf(0.01) = -2,326522873, in 1.5 seconds.
On a 41CL using TURBO=zero the same calculation takes 2.35 seconds.

and just 0.11002 seconds in TURBO50 ( I couldn't resist!)

(03-19-2014 08:31 PM)Dieter Wrote:  You say you used the inverse error function as implemented in the CUDA library. Do you have a link to this algorithm? I found some information on the library, but not on the actual algorithms.

cf. the paper by Mike Giles posted at:
http://people.maths.ox.ac.uk/gilesm/file...erfinv.pdf

Edited: I checked my files looking for the double-precision source code paper. It appears the link has been re-directed to another subject, so I cannot find it on-line. Unfortunately I didn't keep an electronic copy of the document so would need to go to the actual "paper vault" , I hope to find it or else would have to reverse-engineer my own MCODE

Edited: I found the document - send me your email address if you're interested.

[Dieter]

(03-20-2014 06:46 AM)Ángel Martin Wrote:  On V41 using the default settings: qf(0.01) = -2,326522873, in 1.5 seconds.
On a 41CL using TURBO=zero the same calculation takes 2.35 seconds.

That's fast ...but the true value is 2,326347874. That's also what my FOCAL program returns in 9,5 seconds. ?!?

(03-20-2014 06:46 AM)Ángel Martin Wrote:  cf. the paper by Mike Giles posted at:
http://people.maths.ox.ac.uk/gilesm/file...erfinv.pdf

Ah, I happened to come across that paper a few days ago. ;-) Table 5 shows the CUDA code for single precision (approx. 7 digits), which is a set of two 8th degree polynomials, requiring 18 nine-digit constants just for the coefficients.

(03-20-2014 06:46 AM)Ángel Martin Wrote:  Edited: I checked my files looking for the double-precision source code paper. It appears the link has been re-directed to another subject, so I cannot find it on-line. Unfortunately I didn't keep an electronic copy of the document so would need to go to the actual "paper vault" , I hope to find it or else would have to reverse-engineer my own MCODE

Never mind. I just would like to know the basic idea. Is it really a series expansion, or the usual polynomial/rational approximation, or something else?

Dieter

[Ángel Martin]

Believe it or not I typed the wrong value before. The result from ICPF is as follows:

ICPF(0.01) = -2,326347873, in 2.35 seconds.

You can check it out yourself on V41 - both the SandMath and the Library#4 modules will be required.

Yes it is a polynomial approximation, two to be precise depending on the region.

Cheers,
ÁM

[Dieter]

(03-20-2014 02:21 PM)Ángel Martin Wrote:  Believe it or not I typed the wrong value before. The result from ICPF is as follows:
ICPF(0.01) = -2,326347873, in 2.35 seconds.

That looks much better. ;-) Of course I did not assume that the Sandmath function would return a result with an error as large as in the previous value. BTW my current FOCAL program on a standard V41 runs in about 4,5 seconds.

(03-20-2014 02:21 PM)Ángel Martin Wrote:  Yes it is a polynomial approximation, two to be precise depending on the region.

Polynomials are very nice in compiled computer languages where it does not matter whether a constant has 2 or 15 digits, but on a calculator every single digit requires at least one byte of memory. I remember that some time ago I set up two rational approximations for the HP67/97 where the coefficients were kept in the primary and secondary registers which may be stored on a data card. ;-)

Multi-digit numeric constants in user code tend to slow down program execution, even on a HP41 where every value fits in a single line. On the other hand the method is very straightforward and fast by design. I just found an earlier HP41 program for the Normal quantile that uses two rational approximations. It runs in about 6 seconds (in FOCAL) for any input value on a hardware HP41.

BTW, if you are interested in rational approximations for the Normal quantile (or its evaluation in general) I can recommend the papers by W. T. Shaw. Several PDFs are available.

Dieter

[Ángel Martin]

(03-20-2014 09:07 PM)Dieter Wrote:  Polynomials are very nice in compiled computer languages where it does not matter whether a constant has 2 or 15 digits, but on a calculator every single digit requires at least one byte of memory.

Using MCODE the advantage is obvious, as the coefficients are sort of pre-compiled. Rational approximations allow simpler operations (basically sum and multiplication), which avoids slower routines even in MCODE. It's also nice to have the same execution time for all arguments.

How does your method hold up in the vicinity of one? That's where the CUDA approach really surpasses any other approximation I tried in the past. Unfortunately I only used 10 digits to implement the double-precision expression (ran out of ROM space), but even then it returns decent results.

For instance,

ICPF(0.9999) = 3,719016392
ICPF(0.99999) = 4,264890427

accurate to the 8th. and 7th decimal digits respectively, according to the results in Wolfram Alpha:
https://www.wolframalpha.com/input/?i=sq...9999+-1%29

BTW, V41 is not good for benchmarking execution speed - even using default settings it's still a function of the PC's CPU. I'm using an old-reliable XP machine, 9-years old so likely slower than yours ;-)

Cheers,
'AM

[Dieter]

(03-21-2014 05:54 AM)Ángel Martin Wrote:  Using MCODE the advantage is obvious, as the coefficients are sort of pre-compiled. Rational approximations allow simpler operations (basically sum and multiplication), which avoids slower routines even in MCODE. It's also nice to have the same execution time for all arguments.

Yes, that's why I like these methods in Pascal and VB. Over the years I have been using custom-designed rational approximations for various applications. BTW this does not require Maple or Mathematica, it can also be done with standard software like Excel. I really like polynomial/rational approximations - but less so in a user-code program. That's why I suggested the idea of a very simple "first guess" rational approximation with subsequent refinement.

(03-21-2014 05:54 AM)Ángel Martin Wrote:  How does your method hold up in the vicinity of one? That's where the CUDA approach really surpasses any other approximation I tried in the past.

There generally is no problem near one if you use the obvious solution: In all the methods I use for the Normal quantile, only the lower half (p<0,5) is considered and the other half is simply evaluated as -QF(1-p). Easy as that. ;-)

That's why the results for p close to 1 are as accurate as those for 1-p.

QF(0,9) = -1,281551566
QF(0,99) = -2,236347874
QF(0,999) = -3,090232306
QF(0,9999) = -3,719016485
QF(0,99999) = -4,264890794
QF(0,999999) = -4,753424309
QF(0,9999999) = -5,199337582
QF(0,99999999) = -5,612001244
QF(0,999999999) = -5,997807015
QF(0,9999999999) = -6,361340902

All these are exact. Which of course does not mean that the 10-digit FOCAL routine might not have errors in the last place here and there.

In general, it simply does not matter whether QF(0,99999) or QF(0,00001) is calculated. If it was possible to enter 1 - 1E-99 the algorithm would handle this just as well. ;-)

So the idea is:

IF p < 0,5
THEN return QF(p)
ELSE return -QF(1-p)

(03-21-2014 05:54 AM)Ángel Martin Wrote:  Unfortunately I only used 10 digits to implement the double-precision expression (ran out of ROM space), but even then it returns decent results.
...
ICPF(0.9999) = 3,719016392
ICPF(0.99999) = 4,264890427
...
accurate to the 8th. and 7th decimal digits respectively

The FOCAL program I mentioned (initial guess + 1 correction step) is mostly within 1 or 2 ULP, even with just 10 digits working precision.

The only exception I noticed are results between 0,7 and 1 which may be off by 5 or 6 ULP. That's because at this point the CDF that's used in the correction step would require 11-12 digits for an exact quantile.

The FOCAL program using two rational approximations should be within 2 ULP over the whole range. But I did not do any thorough tests.

(03-21-2014 05:54 AM)Ángel Martin Wrote:  BTW, V41 is not good for benchmarking execution speed - even using default settings it's still a function of the PC's CPU. I'm using an old-reliable XP machine, 9-years old so likely slower than yours ;-)

I wouldn't bet on it. The device I'm writing this on shows a BIOS date in late 2002 when powering up. ;-) With standard settings, V41 is about twice as fast as my "real" HP41. Due to the TIME function in V41 the FOCAL program using a rational approximation runs in 2,67 seconds.

Dieter

[Dieter]

(03-21-2014 01:37 PM)Dieter Wrote:  There generally is no problem near one if you use the obvious solution: In all the methods I use for the Normal quantile, only the lower half (p<0,5) is considered and the other half is simply evaluated as -QF(1-p). Easy as that. ;-)
...
So the idea is:
IF p < 0,5
THEN return QF(p)
ELSE return -QF(1-p)

BTW, Ángel: accurate results near 0 and 1 are easy to get if Sandmath would offer the complementary error function erfc(x) and its inverse - instead of erf(x). So what about ERFC and IERFC in a future release?

Quote:The device I'm writing this on shows a BIOS date in late 2002 when powering up.

I now remember I did a BIOS update about 10 years ago. So the hardware must be even older. ;-)

Dieter