Posts: 1,089
Threads: 32
Joined: Dec 2005
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.
Posts: 1,368
Threads: 212
Joined: Dec 2006
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.
Posts: 1,089
Threads: 32
Joined: Dec 2005
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
Posts: 1,368
Threads: 212
Joined: Dec 2006
My routines start with series and continued fraction computations of the incomplete gamma function, straight out of Numerical Recipesindeed 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 leftsided incomplete gamma by series expansion, the rightsided 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 uppertail) and the chisquare probability distribution (again, both cumulative and upper tail).
When I get a chance I will post thisI 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.54e12), 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 seriesmaybe 45 seconds topsbut I should note that the result is off by 1 ULP1.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
Posts: 1,368
Threads: 212
Joined: Dec 2006
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 12digit result should be 9.99999999998e1. This should come as no surpriserounding 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 wellerf(1) is computed swiftly as 0.84270079295, which is completely accurate. erf(0.1) is even faster, giving 1.12462916019e1 (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.
Posts: 1,368
Threads: 212
Joined: Dec 2006
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
1e250
+
STO C
RCL G
RCL/ E
RCL+ B
1e250
+
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.67773498102e3 (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.57299207048e1 (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 1erf(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=1e250 as required to avoid division by zero issues.
Les
Edited: 1 Aug 2007, 10:08 a.m. after one or more responses were posted
Posts: 1,089
Threads: 32
Joined: Dec 2005
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.
Posts: 1,368
Threads: 212
Joined: Dec 2006
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
Posts: 1,368
Threads: 212
Joined: Dec 2006
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.99999984583e1). 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 nonalternating 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.
