RE: HP48HP50G : Error functions ERF(z), FRESNEL integrals, PROBIT, etc.

11092023, 10:32 AM
(This post was last modified: 01172024 10:02 PM by Gil.)
Post: #1




RE: HP48HP50G : Error functions ERF(z), FRESNEL integrals, PROBIT, etc.
Last version 16.09
at last post (page 3)with comments or here Wikipedia: "In statistics, for nonnegative values of x, the error function has the following interpretation: for a random variable Y that is normally distributed with mean 0 and standard deviation 1/√2, erf x is the probability that Y falls in the range [−x, x]." Therefore, with mean 0, variance (1/√2)² = 1/2 = 0.5 and UTPN command, the following short steps enable to easily and accurately calculate Erf(x) = « 0 .5 ROT UTPN 2 * 1 +» and Erfc(x) = « 1 SWAP ERF » or « 0 0.5 ROT UTPN 2 * ». For complex argument, see below. For inverse Erf & Erf, see also below. 

11092023, 12:24 PM
Post: #2




RE: HP48HP50G : Error functions Erf(x) and Erfc(x)
That's a nice short program. I have been using similar programs that were posted years ago on comp.sys.hp48 by John H. Myers.
Code:
Code:
Now if you happen to have a program for the inverse error function, I would be very interested. 

11092023, 12:46 PM
Post: #3




RE: HP48HP50G : Error functions Erf(x) and Erfc(x)
It may be more accurate if Erfc derived from UTPN
Erfc(x) = « 0 0.5 ROT UTPN 2 * » Erf(x) = « 1 SWAP ERFC  » 

11092023, 01:48 PM
(This post was last modified: 11092023 05:18 PM by Gil.)
Post: #4




HP48HP50G : Error functions Erf(x), Erfc(x) & Erf/Erfc —>x
To find the corresponding x value
of the inverse of Erf and Erfc functions, I just reversed the order of the operations and used the builtin solver. iERF « 1  2 / —> p « « p 0 .5 X UTPN  » 'X' 2 ROOT "x of Erf" —>TAG » » iERFc « 2 / —> p « « p 0 .5 X UTPN  » 'X' 2 ROOT "x of Erfc" —>TAG » » Regards, Gil 

11102023, 08:15 PM
(This post was last modified: 11102023 08:52 PM by Gil.)
Post: #5




RE: HP48HP50G : Error functions Erf(x), Erfc(x) & Erf/Erfc —>x
To shorten Joe Keith's code, instead of using the condition IF:
Modified J. Keith's Erf (x) : \<< DUP SIGN NEG 1 ROT SQ 2 * UTPC 1  * \>> (40 bytes instead of 50 bytes —> space gain of 20%). Besides, evaluation time reduced by about 10%. "My" previously given version is still sortner and faster, but less accurate for x<0. To get the very same accuracy for x<0, I will have to work on the absolute value of x and multiply the result by sign(x): "My" modified Erf(x) would be: « DUP SIGN 0 .5 4 ROLL ABS UTPN 2 * 1 + *», which is then longer and slower than Joe Keith's modified version. As usual, Joe Keith always wins. 

11112023, 12:32 AM
Post: #6




RE: HP48HP50G : Error functions Erf(x), Erfc(x) & Erf/Erfc —>x
(11102023 08:15 PM)Gil Wrote: Modified J. Keith's Erf (x) : Very nice! Have you tried Albert's suggestion of deriving ERF from ERFC? 

11112023, 09:26 AM
(This post was last modified: 11132023 03:55 PM by Gil.)
Post: #7




RE: HP48HP50G : Error functions Erf(x), Erfc(x) & Erf/Erfc —>x
About Albert's version:
It's the one I use from the very beginning. The published version is the derived one (c —> complementary, ie "1..."). For clarity, here are my non modified versions for argument a real number: ERF « 0 .5 ROT UTPN 2 * 1 + » iERF (reversed operations) « 1  2 / —> p « « p 0 .5 X UTPN  » 'X' 2 ROOT "x of Erf" —>TAG » » ERFc « 0 .5 ROT UTPN 2 * » iERFc (reversed operations) « 2 / —> p « « p 0 .5 X UTPN  » 'X' 2 ROOT "x of Erfc" —>TAG » » Question 1: Is there an elegant way to solve for the introduced variable X without creating that new variable (that I can of course delete at the end of the programs iERF and iERFc by 'X' PURGE)? Question 2: How could I calculate erf(i*x) for x=0.5 in '2/sqrt(pi)*S(0.,.5*i,e^t^2.,t)' ? A way should be to implement the explanations given in https://math.stackexchange.com/questions...inarypart Observation: The Nspire CX IIT CAS from TI seems to give 14 correct digits for the corresponding integral, versus "only" 13 digits with the builtin erf/erfc functions of the Prime calculator from HP. For Erf(z a complex) & erfc(z), see below. 

11112023, 01:24 PM
Post: #8




RE: HP48HP50G : Error functions Erf(x), Erfc(x) & Erf/Erfc —>x
(11112023 09:26 AM)Gil Wrote: Question 1: Could you make x a local variable? I've never tried using local variables with ROOT but it should work. The formula in the stackexchange link seems straightforward but the resulting program would be big and slow. Unfortunately, neither UTPC nor UTPN allow complex arguments. The vast majority of calculations using ERF involve small, positive real numbers. If I ever need to calculate complex ERF I will probably just use Mathematica or some such. 

11112023, 03:28 PM
Post: #9




RE: HP48HP50G : Error functions Erf(x), Erfc(x) & Erf/Erfc —>x
I already tried with X as a local variable and solver (ROOT), but with no success.
Regards, Gil 

11112023, 03:43 PM
Post: #10




RE: HP48HP50G : Error functions Erf(x), Erfc(x) & Erf/Erfc —>x
(11112023 01:24 PM)John Keith Wrote: Could you make x a local variable? I've never tried using local variables with ROOT but it should work. This can be done, but it's a bit more convoluted than it should be. It appears that the ROOT command has overloaded argument signatures that all expect global IDs only (as opposed to local). This was probably done to save space, as there are already 8 dispatches for ROOT, and allowing locals would have doubled that number. The underlying code, though, uses SysRPL SAFE@, which works fine with either globals or locals. So getting this to work properly involves a little trickery. Basically you have to instantiate the local, then bypass the dispatcher with a SYSEVAL call instead of the UserRPL ROOT command. The address to pass to SYSEVAL depends on the platform you are targeting. I've used Gil's code from above as the basis here, and modified it for a 50g. As usual, be careful with SYSEVAL commands. One small typo or using a code for the wrong platform will likely crash your calculator. Code: ERF I chose to use NOVAL when instantiating the local X, just to make it easier to see. Any value could have been used there, since whatever is there will be overwritten by the ROOT solver. The 50g needs # 2F12Ch as the argument to SYSEVAL. A rev. R 48gx needs # 32F77h instead. 

11112023, 05:18 PM
Post: #11




RE: HP48HP50G : Error functions Erf(x), Erfc(x) & Erf/Erfc —>x
Great, John!
I never understood nor used SYSRPL, but will try it on my HP50G real calculator that I just bought (I use up to now EMU48). Regarding NIVAL Very often in my programs I used an initial, fictive value equal to 0 for my local variable, instead of the more appropriate NOVAL. Never late to learn and improve oneself. Thanks again! 

11112023, 07:41 PM
Post: #12




RE: HP48HP50G : Error functions Erf(x), Erfc(x) & Erf/Erfc —>x
It worked, John.
Many, many thanks. By the way, where can we get the hidden code related to SYSEVAL? 

11112023, 08:29 PM
Post: #13




RE: HP48HP50G : Error functions Erf(x), Erfc(x) & Erf/Erfc —>x  
11112023, 08:46 PM
Post: #14




RE: HP48HP50G : Error functions Erf(x), Erfc(x) & Erf/Erfc —>x
Sorry for the confusion, David,
and my thanks to you and to all the experts for their contributions. 

11112023, 11:34 PM
(This post was last modified: 11122023 04:32 PM by Gil.)
Post: #15




RE: HP48HP50G : Error functions Erf(x), Erfc(x) & Erf/Erfc —>x
From Wikipedia, we can write the following program, with z a complex number:
erf(z) « —> z « 'Sum(n=0,50,(1)^n/((2*n+1)*n!)*z^(2*n+1))' 2 * SQRT(PI) / —>NUM » » Note that a loop of 50, 60 or 100 does not make a real difference, because of the roundings. Most of the given digits are correct, but not all. However, even for z complex, the error function erf(z) of the HP Prime gives all the digits correct. Does somebody know how (which algorithm) the HP Prime proceeds to obtain such a performance? 

11122023, 02:19 PM
Post: #16




RE: HP48HP50G : Error functions Erf(x), Erfc(x) & Erf/Erfc —>x
The problem with that formula is the term ((2*n+1)*n! which rapidly grows into numbers with more than 12 digits. Large rounding errors are thus inevitable. I don't know offhand what algorithm the Prime uses but it probably comes from Xcas/GIAC which is open source.


11122023, 03:57 PM
(This post was last modified: 11122023 06:52 PM by Gil.)
Post: #17




RE: HP48HP50G : Error functions Erf(x), Erfc(x) & Erf/Erfc —>x
I did not look at the open sources, but finally used instead https://math.stackexchange.com/questions...inarypart
creating the following HP50G function/program, that gives most acceptable results up to the last two digits (errors due to "loop roundings"), here now a mended, faster version with summation sign (sigma) ERF(z) « C>R > x y « x ERF 'e^x^2/(2*pi*x)*(1COS(2*x*y)+i*SIN(2*x*y))' >NUM + '2/pi*e^x^2' >NUM 'Sigma(k=1,30,e^(k^2/4)/(k^2+4*x^2)*(2*x*(1COS(2*x*y)*COSH(k*y))+k*SIN(2*x*y)*SINH(k*y)+i*(2*x*SIN(2*x*y)*COSH(k*y)+k*COS(2*x*y)*SINH(k*y))))' >NUM * + >NUM » » And ERF(x) « 0 .5 ROT UTPN 2 * 1 + » The HP Prime must use a similar algorithm. 

11122023, 04:15 PM
Post: #18




RE: HP48HP50G : Error functions Erf(x), Erfc(x) & Erf/Erfc —>x
(11112023 11:34 PM)Gil Wrote: Most of the given digits are correct, but not all. You could try sum it from right to left. Cf. Programming Exercise (HP15C, 15C LE  and others) 

11122023, 04:25 PM
Post: #19




RE: HP48HP50G : Error functions Erf(x), Erfc(x) & Erf/Erfc —>x
(11112023 11:34 PM)Gil Wrote: The erf function of the HP Prime, however, gives all the digits correct. For small argument, erf(x) = 2/√(pi) * x * hyp1f1(1/2, 3/2, x*x) For large argument, erfc via continued fraction. (erf from erfc very accurate!) see https://www.hpmuseum.org/forum/thread15...#pid135589 

11122023, 04:35 PM
(This post was last modified: 11122023 04:48 PM by Gil.)
Post: #20




RE: HP48HP50G : Error functions Erf(x), Erfc(x) & Erf/Erfc —>x
I was amazed by the accuracy of the Hp Prime for erf(z) until I checked by myself, numerically on the HP50G, the "appropriate" expansion and verified numerically the nice/fast "convergence" to expected the solution.


« Next Oldest  Next Newest »

User(s) browsing this thread: 3 Guest(s)