RE: HP48-HP50G : Error functions ERF(z), FRESNEL integrals, PROBIT, etc.

[Gil]

Last version 16.09
at last post (page 3)with comments or here

Wikipedia:
"In statistics, for non-negative 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.

[John Keith]

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:


@ ERF
\<< DUP SQ 2. * 1. SWAP UTPC 1. - SWAP 0. >
  { NEG } IFT
\>>

Code:


@ ERFC
\<< DUP SQ 2. * 1. SWAP UTPC SWAP 0. <
  { 2. - NEG } IFT
\>>

Now if you happen to have a program for the inverse error function, I would be very interested.

[Albert Chan]

It may be more accurate if Erfc derived from UTPN

Erfc(x) = « 0 0.5 ROT UTPN 2 * »
Erf(x) = « 1 SWAP ERFC - »

[Gil]

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 built-in 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

[Gil]

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.

[John Keith]

(11-10-2023 08:15 PM)Gil Wrote:  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 shorter and faster, but less accurate for x<0.

Very nice! Have you tried Albert's suggestion of deriving ERF from ERFC?

[Gil]

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...inary-part

Observation:
The Nspire CX II-T CAS from TI seems to give 14 correct digits for the corresponding integral, versus "only" 13 digits with the built-in erf/erfc functions of the Prime calculator from HP.

For Erf(z a complex) & erfc(z), see below.

[John Keith]

(11-11-2023 09:26 AM)Gil Wrote:  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...inary-part

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.

[Gil]

I already tried with X as a local variable and solver (ROOT), but with no success.

Regards,
Gil

[DavidM]

(11-11-2023 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
\<<
  0 .5 ROT UTPN -2 * 1 +
\>>

iERF
\<<
  1 - -2 / NOVAL \-> p X
  \<<
    \<< p 0 .5 X UTPN - \>>
    'X'
    2
    # 2F12Ch SYSEVAL   @ direct ROOT code
    "x of Erf" \->TAG
  \>>
\>>

ERFc
\<<
  0 .5 ROT UTPN 2 *
\>>

iERFc
\<<
  2 / NOVAL \-> p X
  \<<
    \<< p 0 .5 X UTPN - \>>
    'X'
    2
    # 2F12Ch SYSEVAL   @ direct ROOT code
    "x of Erfc" \->TAG
  \>>
\>>

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.

[Gil]

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!

[Gil]

It worked, John.

Many, many thanks.

By the way, where can we get the hidden code related to SYSEVAL?

[John Keith]

Er, that was David's program, not mine.

If you want to see the code at that address, use nosy..

[Gil]

Sorry for the confusion, David,

and my thanks to you and to all the experts for their contributions.

[Gil]

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?

[John Keith]

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.

[Gil]

I did not look at the open sources, but finally used instead https://math.stackexchange.com/questions...inary-part

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)*(1-COS(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*(1-COS(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.

[Thomas Klemm]

(11-11-2023 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 (HP-15C, 15C LE - and others)

[Albert Chan]

(11-11-2023 11:34 PM)Gil Wrote:  The erf function of the HP Prime, however, gives all the digits correct.
Does somebody know how (which algorithm) the HP Prime proceeds to obtain such a performance?

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/thread-15...#pid135589

[Gil]

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.