Error function 35s vs. 32SII
#1

FYI, I implemented a simple program to calculate the error function on both the 35s and the 32SII.

Results (Argument 5, full precision)

32SII: 14.5 s
35s: 22.3 s
OTOH, using the equation editor and integrate it (Argument 5, four digits precision) results in:
32SII: 5.87 s
35s: 4.9
So, while RPL programs seem to run a lot faster on the old pioneer, at least one of the internal functions of the 35s is able to beat its grandfather.
#2

Thomas, what algorithm do you use to calculate erf and erfc?

I have both a series and continued fraction implementation that are perceptibly slower on the 35s vs. the 33s, but a typical calculation doesn't seem to take any more than a couple or three seconds tops, even faster on the 33s. Haven't tested on 32sii but in my experience it runs about as fast as the 33s for most things.

22s, even 14s, seems a really slow way to compute this function.

Les

Edited: 31 July 2007, 3:04 a.m.

#3

Dear Les,

I'm using the taylor series, argument 5 is completely unusual for most applications and takes 51 iterations. E.g., only 5 are needed for 0.1

Of course, the program can be improved (there's not even an error handling)

Thomas

LBL A
STO X
0
STO S
STO N
1
STO F
STO P
STO V
100
STO L

(A012)

RCL N
x=0?
GTO A016
STO * F

(A016)

2
*
1
+
STO T
RCL V
RCL X
RCL T
y^x
*
RCL F
/
RCL T
/
RCL P
/
STO B
2
STO * P
-1
STO * V
RCL B
STO + S
RCL S
RCL L
x=y?
GTO A048
RCL S
STO L
1
STO + N
GTO A012

(A048)

RCL S
2
PI
*
SQRT
/
RTN

#4

My routines start with series and continued fraction computations of the incomplete gamma function, straight out of Numerical Recipes--indeed when I program continued fraction calculations I use the modified Lentz method in preference to the better known and very old Wallis method, for the reasons NR gives. I compute the left-sided incomplete gamma by series expansion, the right-sided by continued fraction. The error functions are of course special cases of the incomplete gamma functions, as are the normal distribution probability function (both cumulative and upper-tail) and the chisquare probability distribution (again, both cumulative and upper tail).

When I get a chance I will post this--I have code for everything except the chisquare stuff. All told just over 100 steps tops, but there is no error handling, and the user needs to decide which function (series or continued fraction?) is preferable based on his input and his knowledge of the properties of these calculations.

I agree that it is unusual to compute erf(5) by series expansion. It makes more sense to compute erfc(5) (about 1.54e-12), and compute erf(5) = 1 - erf(5). That said, my series routines don't take anywhere near over 20 seconds to compute erf(5) directly by series--maybe 4-5 seconds tops--but I should note that the result is off by 1 ULP--1.00000000001.

It also looks like you are using the alternating series expansion. The other series expansion I know is not alternating and has an exponential term out front. I don't know if it converges any faster.

Les

#5

Thomas, here is my offering. The reference is formula 7.1.6 of Abramowitz and Stegun, whereas you seem to use 7.1.5. The idea of each term being a recursion on the prior term I get from the Numerical Recipes code for the series computation of the incomplete gamma function.

(E001) LBL E
STO T
STO S
x^2
STO Z
1
STO D
RCL S
(E009) STO O
2
STO+ D
RCL Z
STO* T
2
STO* T
RCL D
STO/ T
RCL T
STO+ S
RCL O
RCL S
x#y?
GTO E009
RCL Z
+/-
e^x
*
2
*
PI
SQRT
/
(E033) RTN

5 XEQ E ENTER returns 1.00000000002 in a few seconds (certainly no more than 10), whereas the actual 12-digit result should be 9.99999999998e-1. This should come as no surprise--rounding error is bound to crop up. For input greater than about 1.8, one should compute erfc by the continued fraction and get erf as the complement (erf(z) = 1 - erfc(z))

But for smaller arguments, the above does very well--erf(1) is computed swiftly as 0.84270079295, which is completely accurate. erf(0.1) is even faster, giving 1.12462916019e-1 (the last digit should be an 8).

I guess what caught my attention is that I knew one could compute these series faster on the 33S and 32sii than what you were reporting. I must confess that the comparative sluggishness of the 35S (though still much faster than the 41 series, the 11C, the 15C, the Spice series, and even the 42S) is a disappointment. I love the speed of the 33S, despite its foibles, and really hoped this would port to the 35s intact. Alas, that hasn't happened.

I may port my continued fraction code for the incomplete gamma function, modified for the special case of erfc, and post it later if anyone is interested.

Les

Edited: 31 July 2007, 12:49 p.m.

#6

Here is that promised continued fraction computation for erfc:

(C001) LBL C
STO X
x^2
0.5
+
STO B
1/x
STO C
STO D
1e250
STO E
CLx
STO H
(C014) 1
STO+ H
STO+ B
STO+ B
0.5
RCL- H
RCL* H
STO G
RCL* C
RCL+ B
1e-250
+
STO C
RCL G
RCL/ E
RCL+ B
1e-250
+
STO E
RCL C
1/x
x<> C
/
STO* D
1
x#y?
(C040) GTO C014
RCL X
ENTER
x^2
+/-
e^x
*
RCL* D
PI
SQRT
/
(C051) RTN

This is best employed when input is greater than about 1.8. For example, 2 XEQ C ENTER returns 4.67773498102e-3 (the last digit should be 5) in a few short seconds. The routine gets faster as the input gets larger. Nonetheless, 1 XEQ C chugs along for a few seconds to yield 1.57299207048e-1 (the last digits should be 50 so it is off 2 ULP). Below that, if not sooner, the continued fraction really should give way to the series computation of erf and erfc(x) computed as 1-erf(z).

My routine started out as my port of the Numerical Recipes continued fraction computation of the incomplete gamma function, modified for the specific case of the complementary error function. The modified Lentz algorithm is used with TINY=1e-250 as required to avoid division by zero issues.

Les

Edited: 1 Aug 2007, 10:08 a.m. after one or more responses were posted

#7

Thanks for your code! This makes me wonder if there's a repository for such pearls somewhere? I fear, the work shown here quickly gets lost.

#8

Thomas, my plan is eventually submit the incomplete gamma routine on which these are based, along with subroutines to compute the special cases of the error function, chisquare distribution, and normal distribution, to Dave for inclusion in the software library on this site. I also hope to port a somewhat long 42S program I have that computes the incomplete beta function, along with its special cases (t distribution and F distribution).

I mentioned in an earlier thread that I hoped the expanded flexibility of the 35s programming environment will make it easier to write meaningful programs. I find the line number addressing heaven sent. I have three programs in my 33S that use up almost all the labels right there. The 35s is now free of that limitation.

Les

#9

Here, for completeness, is my take on the the Taylor series (formula 7.1.5 in Abramowitz and Stegun):

(T001) LBL T
STO T
STO S
x^2
+/-
STO Z
0
STO D
1
STO E
RCL S
(T012) STO O
1
STO+ D
STO+ E
STO+ E
RCL Z
STO* T
RCL D
STO/ T
RCL T
RCL/ E
STO+ S
RCL O
RCL S
x#y?
(T027) GTO T012
2
*
PI
SQRT
/
(T033) RTN

This converges rapidly and accurately for small values of input but once you get much beyond 3 results become not only inaccurate, but pretty meaningless.

For example, 4 XEQ T ENTER returns 1.00000031185 (whereas the actual result is actually about 9.99999984583e-1). I assume that with the alternating series, rounding errors accumulated, or, in the alternative, differences between terms begin to cancel out so the loop terminates too soon.

For this reason, I really think the non-alternating series is better in generally, though there is no doubt this series is very fast for smaller values of input.

Les

Edited: 1 Aug 2007, 3:19 p.m.



Possibly Related Threads…
Thread Author Replies Views Last Post
  HP50g: Writing a function that returns a function Chris de Castro 2 2,063 12-10-2013, 06:49 PM
Last Post: Han
  MS advert shows spreadsheet with obvious error BruceH 3 2,146 11-14-2013, 09:50 AM
Last Post: Bill (Smithville, NJ)
  HP Prime: Rounding error in determinant Stephan Matthys 3 1,607 10-25-2013, 09:29 PM
Last Post: Walter B
  Prime Error or Mine? toml_12953 12 3,531 10-22-2013, 10:35 AM
Last Post: toml_12953
  Explaination on How to Reset Caculator in Users guie: error Harold A Climer 5 2,512 10-15-2013, 02:11 AM
Last Post: cyrille de Brébisson
  Repair of HP-34C - Error 0 and Error 9 Jeff Kearns 3 1,644 10-11-2013, 12:29 PM
Last Post: Randy
  Do You Think a Listing Error Was Made? Jim Johnson 13 3,775 09-04-2013, 09:23 AM
Last Post: Eddie W. Shore
  HP Prime: Error reports : Here Patrice 111 22,317 07-24-2013, 05:52 PM
Last Post: Thomas Klemm
  PARTFRAC error. Romeu Medeiros 2 1,295 06-23-2013, 10:27 PM
Last Post: Romeu Medeiros
  WP-34s Documentation error Marcel Samek 9 2,707 06-13-2013, 10:25 AM
Last Post: Miguel Toro

Forum Jump: