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

03202014, 06:46 AM
(This post was last modified: 03202014 12:41 PM by Ángel Martin.)
Post: #21




RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"?
(03192014 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!) (03192014 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 doubleprecision source code paper. It appears the link has been redirected to another subject, so I cannot find it online. 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 reverseengineer my own MCODE Edited: I found the document  send me your email address if you're interested. "To live or die by your own sword one must first learn to wield it aptly." 

03202014, 12:44 PM
Post: #22




RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"?
(03202014 06:46 AM)Ángel Martin Wrote: On V41 using the default settings: qf(0.01) = 2,326522873, in 1.5 seconds. That's fast ...but the true value is 2,326347874. That's also what my FOCAL program returns in 9,5 seconds. ?!? (03202014 06:46 AM)Ángel Martin Wrote: cf. the paper by Mike Giles posted at: 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 ninedigit constants just for the coefficients. (03202014 06:46 AM)Ángel Martin Wrote: Edited: I checked my files looking for the doubleprecision source code paper. It appears the link has been redirected to another subject, so I cannot find it online. 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 reverseengineer 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 

03202014, 02:21 PM
Post: #23




RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"?
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 "To live or die by your own sword one must first learn to wield it aptly." 

03202014, 09:07 PM
(This post was last modified: 03202014 09:14 PM by Dieter.)
Post: #24




RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"?
(03202014 02:21 PM)Ángel Martin Wrote: Believe it or not I typed the wrong value before. The result from ICPF is as follows: 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. (03202014 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. ;) Multidigit 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 

03212014, 05:54 AM
(This post was last modified: 03212014 10:48 AM by Ángel Martin.)
Post: #25




RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"?
(03202014 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 precompiled. 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 doubleprecision 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 oldreliable XP machine, 9years old so likely slower than yours ;) Cheers, 'AM "To live or die by your own sword one must first learn to wield it aptly." 

03212014, 01:37 PM
(This post was last modified: 03212014 06:33 PM by Dieter.)
Post: #26




RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"?
(03212014 05:54 AM)Ángel Martin Wrote: Using MCODE the advantage is obvious, as the coefficients are sort of precompiled. 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 customdesigned 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 usercode program. That's why I suggested the idea of a very simple "first guess" rational approximation with subsequent refinement. (03212014 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(1p). Easy as that. ;) That's why the results for p close to 1 are as accurate as those for 1p. 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 10digit 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  1E99 the algorithm would handle this just as well. ;) So the idea is: IF p < 0,5 THEN return QF(p) ELSE return QF(1p) (03212014 05:54 AM)Ángel Martin Wrote: Unfortunately I only used 10 digits to implement the doubleprecision expression (ran out of ROM space), but even then it returns decent results. 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 1112 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. (03212014 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 oldreliable XP machine, 9years 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 

03282014, 07:56 PM
Post: #27




RE: Accurate Bernoulli numbers on the 41C, or "how close can you get"?
(03212014 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(1p). Easy as that. ;) 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 

« Next Oldest  Next Newest »

User(s) browsing this thread: