RE: HP48-HP50G : Error functions ERF(x or z), ERFC (x or z) & iERF/iERFc —>x +- HP Forums (https://www.hpmuseum.org/forum) +-- Forum: HP Software Libraries (/forum-10.html) +--- Forum: General Software Library (/forum-13.html) +--- Thread: RE: HP48-HP50G : Error functions ERF(x or z), ERFC (x or z) & iERF/iERFc —>x (/thread-20807.html) |
RE: HP48-HP50G : Error functions ERF(x or z), ERFC (x or z) & iERF/iERFc —>x - Gil - 11-09-2023 04:32 AM 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. HP48-HP50G : Error functions Erf(x), Erfc(x) & Erf/Erfc —>x - Gil - 11-09-2023 07:48 AM 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 RE: HP48-HP50G : Error functions Erf(x), Erfc(x) & Erf/Erfc —>x - Gil - 11-11-2023 03:26 AM 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/712434/erfaib-error-function-separate-into-real-and-imaginary-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. RE: HP48-HP50G : Error functions Erf(x), Erfc(x) & Erf/Erfc —>x - Gil - 11-12-2023 09:57 AM I did not look at the open sources, but finally used instead https://math.stackexchange.com/questions/712434/erfaib-error-function-separate-into-real-and-imaginary-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. RE: HP48-HP50G : Error functions ERF(x or z), ERFC (x or z) & iERF/iERFc —>x - Gil - 11-12-2023 04:16 PM Summary Now, both ERF and ERFc functions allow as single argument a real x (the returned result will be a real "probability") or a complex number in the form (a b) (the returned result will be a complex number); Their inverse functions iERF and iERFc, however, only admit here as single argument a real number/"probability" (the returned result will be the real number x). Note For real argument/numbers, the results calculated here for ERF and Erfc may differ by 1 unit in the last digit in comparison with the built-in functions erf(x) and erfc(x) of the HP PRIME. For complex arguments, the results calculated here for ERF and ERFc may slightly differ in their last two last digits (the HP PRIME is then more accurate). RE: RE: HP48-HP50G : Error functions ERF(x or z), ERFC (x or z) & iERF/iERFc —>x - Gil - 11-14-2023 07:35 AM Check the following: ERF(+5.1) :.999999999999 ERF (+5.2'):1 (with the HP50G roundings; in fact,correct is exactly 0.99999999999980751) 'ERFc (up to now) : 0 (1-"1" above) (in fact should be 1.9249061099972305*10^-13) But of course iERFc(0) does not give the expected value of 5.2. We have to try and get iERFc(ERFc(x)) = x (or a near value). That was the point of Albert Chan to which I did not pay attention. Finally, I changed the programs, starting by calculating ERFc, then setting ERF=1-ERFc. RE: RE: HP48-HP50G : Error functions ERF(x or z), ERFC (x or z) & iERF/iERFc —>x - Gil - 11-14-2023 02:37 PM Again, a lesson... Yes, I saw and appreciate the huge improvement for x very small. Sorry for my ignorance, but I did not understand the code to try and adapt it to HP50G. Unfortunately. I try here y = ierfc(1-x) -- assume good ierfc With x=10^(-9) —> y=:x of Erfc: 8.86019642276E-10 local h = x - (x-1+1) —> should be zero? return h==0 and y or y + sqrt(pi)/2 * exp(y*y) * h But y small, y² very small & exp(y²) = 1 y=y+sqrt (pi)/2 *0 * 0 = y+0 = y end function erf_taylor(x) local y = 1 - erfc(x) -- assume good erfc —>y=.000000001128 on HP50G return h==0 and y or y + 2/sqrt(pi) * exp(-x*x) * h h=10^(-9)-8.86019642276E-10 But X small, x² very small, exp x² = 1 & exp (-x²) =1 —> y=y+2/sqrt(pi)*0*0=y+0=y RE: RE: HP48-HP50G : Error functions ERF(x or z), ERFC (x or z) & iERF/iERFc —>x - Gil - 11-15-2023 08:13 AM 1) I corrected ERF and Erfc programs, in particular regarding complex numbers. For ERFc(z), with z=0+bi, i.e. (0 b) in the program, I cheated and calculated ERFc(0.0000000001 b), and replaced at the end the real part (about 0.9999...) by 1. 2) I checked the results with Wolfram Alpha and the HP PRIME. - On the HP PRIME, approximate(erfc(5*i)), ie with no real part in the complex number input (argument), returns -82982738 79.76*i (a complex number with no real part). - On my HP50G program ERFc(0 5) returns (1.,-8298273880.6 1 ), with 1 as the real part of that complex number. - Whereas Wolfram Alpha returns 1 - 8.2982738806 7680 ×10^9 i. Fazit 1) As no real part of the HP PRIME result is shown, that part has to be considered being equal to 0, and the PRIME result as wrong here (a bug?). 2) Imaginary part of the HP PRIME result is "worse" than my HP50G program (algorithm of the HP PRIME should be improved?). 3) Regarding Albert Chan observations Still to be taken into account... RE: RE: HP48-HP50G : Error functions ERF(x or z), ERFC (x or z) & iERF/iERFc —>x - Gil - 11-16-2023 05:26 PM I changed the way of calculating erf, now with the corresponding integral (much longer, of course), instead of using the UTPN command. Reason : Try 0. 5 1E-9 UTPN —> .499999999436 × 2 —> .999999998872 erf=1-.999999998872 =.000000001128. Some digits got logically lost by this procedure. Now, with integral from -x to x gives the expected answer in our case : 1.12837916709E-9. When executing .000000001 ERF¦c We get two results: "ERFc (.000000001)" .999999998872 "ERF (.000000001)" 1.12837916709E-9 No more necessary to use ERFc, but of course you should not delete it, as it is used in ERF¦C. The input can be a real or a COMPLEX number. Note that iERFc will, logically, not give back the initial value x RE: RE: HP48-HP50G : Error functions ERF(x or z), ERFC (x or z) & iERF/iERFc —>x - Gil - 11-19-2023 10:16 AM HP50: ERF(x), ERFc(x), iERFc(x) and iERFc(x), with x a real number Now significand or mantissa with 11 to 12 accurate digits (even for the inverse functions iERF and iERFc). Note that using UTPN does not give always the full significant digits (or best solution). The critical point of x lies then at about 0.114294. Depending on the values of x, an integral was used accordingly (instead of UTPN). —> For computational time, EMU48 is then highly recommended. For ERF(z) and ERFc(z), with z a complex number (a b): - the imaginary part of the resulting complex number should also give the correct 11-12 digits of the mantissa; - the real part of the resulting complex number, however, may not show all the mantissa digits (only some incomplete digits), depending on the case. This issue still should be fixed up. RE: RE: HP48-HP50G : Error functions ERF(x or z), ERFC (x or z) & iERF/iERFc —>x - Gil - 11-21-2023 06:40 PM New version 6. Better, but still not perfect. 1) Some input checking. 2) Improvements on final output and in the calculations. 3) Inverse functions, named now IERF and IERFc with capital I (and no more iERF and iERFc, as i has another meaning). 4) In last version 5, I inadvertently overwrote the original inverse_ERF. With my apologises. RE: RE: HP48-HP50G : Error functions ERF(x or z), ERFC (x or z) & iERF/iERFc —>x - Gil - 11-23-2023 06:57 AM New version 7 - with some improvements - and inclusion of other related functions. 1) Same input checking as in version 6. 2) Improvements on output "outlook". 3) Inverse functions, named now IERF and IERFc with capital I (and no more iERF and iERFc , as small i has another meaning —> see 4) below). 4) Inclusion of 3 functions related to the error function: - iERFc - ERFi - w.FADDEEVA-KRAMP For definitions, see Wikipedia/French and Wolfram Alpha, or 10) for the codes. 5) Inverse_ERF(p) Work with positive values of p only IF p-1+1=p Then Inverse_ERF(p) = ERFc(1-p) (we want to be sure that we don't lose digits by calculating 1-p; see, for instance, the results for p=5.0000000005E-4: by Wolfram Alpha —> 4.4311349177240266E-4 by integral as p-1+1≠p —> 4.43113491773E-4 by built-in UTPN (in HP48-50G), but not using Inverse_ERFc(1-p) —> 4.43113492614E-4 by Inverse_ERFc(1-p) —> 4.43113491729E-4) Else use Solver with initial value of p with always the Integral calculation, that gives better results. END IF initially p was negative Then result =: -result Else result = last, positive result END 6) Inverse_ERFc(p) IF p > 1 Then p:= 2-p (we calculate initially for resulting positive x values) Else p= initial value p END IF p now > 0.9 (it's an approximate value) Then Solver with always the Integral calculation, that gives better results (example results for x =.9999: by Wolfram Alpha: 8.8622692777289E-5 by integral: 8.86226927774E-5 by built-in UTPN in HP48-50G: 8.86226927139E-5) Else use the built-in UTPN function in HP48-50G. END IF you proceded initially to the change of p into 2-p Then final result = -result Else final result = same, positive result END 7) ERF¦c, together with ERF (a help program used internally in ERF¦c, but not to be used directly) Both programs still should be improved (some final digits induly may not appear in the outputs ERF¦c 8) All functions take one single argument. 9) All functions admit complex arguments on the form (a b), less the inverse functions IERF and IERFc (that take only real arguments). 10) Codes listing ERF¦c 1 input: real x or a complex (a b) 2 outputs : the two corresponding "1-p" and "p" values ERFc and ERF \\<< \"1 Arg .real x or .complex (a b)\" DROP DTAG EVAL DUP DUP DUP TYPE 1 == IF THEN \"ERF \" SWAP + ELSE \"ERF (\" SWAP + \")\" + END DUP \"F\" \"Fc\" SREPL DROP 4 ROLLD SWAP DUP DUP IM SWAP RE \\-> im re \\<< CASE re DUP 0 \\>= 1 -1 IFTE SWAP ABS '\\oo' EVAL == THEN 1 * DUP 1 SWAP - ROT 5 ROLL DROP2 UNROT END DROP im 0 == re ABS .01 \\<= AND IF THEN DROP '\\.S(-re,re,e^-X^2,X)/\\v/\\pi' \\->NUM ELSE im 0 == re 0 \\<= AND IF THEN -1 ELSE 1 END SWAP im 0 == IF THEN ABS END ERFc 1 SWAP - * END im 0 == IF THEN ROT ERFc UNROT ELSE ROT DROP DUP 1 SWAP - UNROT END END \\>> 'IERR' PURGE \\>> IERF 1 real input 1 real output \\<< \"1 Arg .a real a, |a| \\<= 1 (not a complex !)\" DROP DUP DUP IM 0 \\=/ SWAP RE ABS 1 > OR IF THEN DROP \"Not OK as: .|Real| > 1 or .Complex with IM\\=/0 or |RE| > 1\" DOERR END DUP ABS 1 > IF THEN \"|a| must be \\<= 1\" DOERR END DUP \"InvERF (\" SWAP + \")\" + SWAP DUP 0 \\>= 1 -1 IFTE SWAP ABS \\-> p \\<< CASE p 1 - 1 + p == THEN 1 p - IERFc SWAP DROP END p '\\.S(-x,x,e^-X^2.,X)/\\v/\\pi' - 'x' 2 ROOT { x IERR } PURGE END * \\>> \\>> IERFc 1 complex or real input 1 complex output \\<< \"1 Arg .a real a in [0 2] (not a complex !)\" DROP \\-> a \\<< \"InvERFc \" a TYPE 1 \\=/ IF THEN \"(\" + a + \")\" ELSE a END + a IM 0 \\=/ a RE 0 < OR a RE 2 > OR IF THEN DROP \"Not OK as: .Real <0, >2 or .Complex with IM\\=/0 or RE <0, >2 \" DOERR ELSE PUSH -105 CF a 1 \\<= IF THEN 1 ELSE 2 a - 'a' STO -1 END CASE a 0 == THEN '\\oo' EVAL END a 1 == THEN 0 END a .9 \\<= THEN \\<< a 0 .5 x UTPN 2 * - \\>> 'x' 2. ROOT END 1 a - RE '\\.S(-x,x,e^-X^2.,X)/\\v/\\pi' - 'x' 2 ROOT END * { x IERR } PURGE END POP \\>> \\>> iERFc 1 complex or real input 1 complex output \\<< DUP \\-> z \\<< \"iERFc \" z TYPE 1 == IF THEN z ELSE \"(\" + z + \")\" END + SWAP ERF\\166c 4 ROLL 3 DROPN z * z SQ NEG EXP \\pi \\v/ / \\->NUM - DUP NEG \"Wiki_French\" \\->TAG SWAP \"Wolfram\" \\->TAG \\>> \\>> ERFi 1 input complex or real input 1 complex output \\<< DUP \\-> z \\<< \"ERFi \" z TYPE 1 == IF THEN z ELSE \"(\" + z + \")\" END + SWAP i * \\->NUM ERF\\166c 4 ROLLD 3 DROPN i / \\->NUM \\>> \\>> W.FADD 1 complex or real input 1 complex output \\<< DUP \\-> z \\<< \"wFADDEEVA-KRAMP \" z TYPE 1 == IF THEN z ELSE \"(\" + z + \")\" END + SWAP i NEG * \\->NUM ERF\\166c 4 ROLL 3 DROPN z SQ NEG EXP * \\->NUM \\>> \\>> ERFc Not for the user Program used internally by ERF¦c \\<< \"1 Arg .real x or .complex (a b)\" DROP DUP IM 0 == IF THEN RE 0 .5 ROT UTPN 2 * ELSE C\\->R SWAP DUP 0 == IF THEN .00000000001 SWAP ELSE 1 END \\-> y x z \\<< 0 .5 x UTPN 2 * 'e^-x^2./(2.*\\pi*x)*(1.-COS(2.*x*y)+i*SIN(2.*x*y))' \\->NUM '2./\\pi*e^-x^2.' \\->NUM '\\GS(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 * + - z 0 == IF THEN C\\->R SWAP DROP 1 SWAP R\\->C END \\>> END \\>> 11) Request for Inverse_ERF (IERF) and Inverse_ERFc (IERFc) For both above mentioned functions and without having to use an integral (which is quite time consuming if using the real calculator), I would duly appreciate to have corresponding, possible algorithms and codes allowing to get — in all possible situations — the same or better results. RE: RE: HP48-HP50G : Error functions ERF(x or z), ERFC (x or z) & iERF/iERFc —>x - Gil - 12-03-2023 08:58 AM New version 8 - with improvements - and inclusion of new other related functions. 1) Same input checking as in versions 6/7 + other input checkings,insuring that we are well in presence of not to large arguments (for real, max (IxI) OK if <= 9; for complex z=a+ib, max (IaI) and max(IbI) OK if both <= 9); 1b) The above max number (argument) 9 can be changed in the help program called max.arg. 1c) The help programs are labelled/written with small letters. See also 1g. 1d) The names of the 3 help programs are: complex (to deal with complex numbers), max.arg (above mentioned) and euler (to deal with complex numbers, to try and get the most exact result, when exact result should be ± (1 0) or ± (0 1)). 1e) The help programs are internal programs. They are not to be used manually, or changed, with the exception or the help program max.arg (see 1b). 1f) Instances of results with help program euler, in exact or approximate mode, with a dot after the 2 in the first two example: 'e^(2.*i*pi)' euler —> 1., and not (1.,2.54352307471E-12) 'e^(1.5*i*pi)' euler (0.,-1.), and not (-3.08985769397E-12,-1.) but 'e^(2.2*i*pi)' euler —> (.809016994372,.587785252296). 1g) The name of all the user fonctions are written with (some) capital letters. See also 1c. 2) Improvements on output "outlook". 3) Inverse functions, named now INV¦ and INV¦c (and no more iERF and iERFc , as small i has another meaning —> see 4) below). 4) Besides the 3 functions related to the error function: - iERFc - ERFi - w.FADDEEVA-KRAMP, inclusion now of the following 3 new functions: - PROBIT - ERFcx - F_DAWSON For definitions, see Wikipedia/French and Wolfram Alpha, or 10) for the codes. 5) Inverse_ERF(p), labelled now INV¦ Initial program will convert p (p>=0 or p<0) into IpI. IF p-1+1=p Then Inverse_ERF(p) = ERFc(1-p) (we want to be sure that we don't lose digits by calculating 1-p; see, for instance, the results for p=5.0000000005E-4: by Wolfram Alpha —> 4.4311349177240266E-4 by integral as p-1+1≠p —> 4.43113491773E-4 by built-in UTPN (in HP48-50G), but not using Inverse_ERFc(1-p) —> 4.43113492614E-4 by Inverse_ERFc(1-p) —> 4.43113491729E-4) Else use Solver with initial value of p with always the Integral calculation, that gives better results. END IF initially p was negative Then result =: -result Else result = last, positive result END 6) Inverse_ERFc(p), called now INV¦c IF p > 1 Then p:= 2-p (we calculate initially for resulting positive x values) Else p= initial value p END IF p now > 0.9065 (it's an approximate value) Then Solver with always the Integral calculation, that gives better results (example results for x =.9999: by Wolfram Alpha: 8.8622692777289E-5 by integral: 8.86226927774E-5 by built-in UTPN in HP48-50G: 8.86226927139E-5) Else use the built-in UTPN function in HP48-50G. END IF you proceded initially to the change of p into 2-p Then final result = -result Else final result = same, positive result END 6b) ERF and ERFc The programs take into account the fact that the integral calculation gives more accurate result when IxI<= 0.01 then making use of the built-in function UTPN. 7) Initially, all programs allowed ± infinity as argument, supposing that the imaginary part of the argument did not exist. Now, it's not any more the case, but the "infinity" tests still appear, though they are not any more used in the present program execution. That part regarding "infinity" argument should be carefully removed in a cleaner version (or, again, treated separately, taking into account the fact that the argument may be a complex number (a+ib) with a=±infinity and b≠ infinity, or a≠infinity and b= ± infinity, or both a=±infinity and b=±infinity). 8) As in previous versions, all functions take one single argument. 9) Now all functions admit as argument, besides a real number, a complex number in the form: - (a b) - 'a+ib' - 'r*e^(i*theta)' - 'r*EXP(i*theta)', even INV¦ and INV¦c functions (programs). RE: RE: HP48-HP50G : Error functions ERF(x or z), ERFC (x or z) & iERF/iERFc —>x - Gil - 12-10-2023 04:52 PM New version 9 MAIN FUNCTIONS: all labelled with some capitals. All of them accept real or complex format. Automatic limitation of inputs: to insure good results (up to 9-10 digits). Many help programs/sub routines: all labelled with small letters. Arg/input = ± infinity, ±i × infinity: ERF and ERFi functions OK for these extreme cases. Reverse functions of ERF and ERFc only accept, as argument, real x or z=a+ib with b=0. All other functions accept, as argument, real x or z=a+ib with a, b free real. Extreme cases for arg/input = ± infinity, ±i × infinity for functions iERFc ERFcx FADDEEVA DAWSON: still to be handled. RE: RE: HP48-HP50G : Error functions ERF(x or z), ERFC (x or z) & iERF/iERFc —>x - Gil - 12-12-2023 06:55 AM Version 9b Now inverse functions of the error function ERF and ERFc, called respectively INV¦ & INV¦c, admit as argument p a complex number p = a+ib, b≠0. For the argument p of INV¦(p): IpI < 1. For the argument p of INV¦c(p): I1-pI < 1. On that purpose, a new subroutine called 'inv¦z' was used. Better here to keep all the 50 terms of the numerator and of the denominator. To save memory, however, you can use (approximate) real numbers, instead of the (exact) integers, in the 2 {} lists of the above mentioned subroutine 'inv¦z'. Regarding PROBIT function A slight change in the program, but the argument of PROBIT function was kept to remain a real number p in [0 1]. RE: RE: HP48-HP50G : Error functions ERF(x or z), ERFC (x or z) & iERF/iERFc —>x - Gil - 12-13-2023 12:16 PM Version 9d Much more compact subroutine 'inv¦z' for inverse calculation of ERF & ERFc, when argument z is a complex with abs(z) < 1 : now, to save memory place, the coefficients in the formulae development used in 'inv¦z' derive also from a formulae. Strongly recommended to use the EMU48 Android application, if you want to get the results within a couple of seconds. RE: RE: HP48-HP50G : Error functions ERF(x or z), ERFC (x or z) & iERF/iERFc —>x - Gil - 12-14-2023 05:22 AM Version 9f In the previous versions 9, the simplification of the input argument was incorrect when in presence of a complex arguments z=a+ib, with b≠0 (but everything was correctly calculated). Now the 'arg.num' subroutine was duly fixed up. Please note for INVERSE_ERF, called 'INV¦', when using complex argument z=a+ib, b≠0, that the results are quite accurate for abs(argument z) < 0.91. Please note that this last uploaded version 9f has a "print" error at variable VERSIO. Please change it from Version 9e 2023.12.14 campart @hotmail.com into Version 9f 2023.12.14 campart @hotmail.com RE: RE: HP48-HP50G : Error functions ERF(x or z), ERFC (x or z) & iERF/iERFc —>x - Gil - 12-16-2023 07:10 AM Version 10d (correction of versions 10b & 10c, which were deleted; see also note at the end of this post) (version 10 non "published") Thanks Albert Chan's insight, here is a simplified version of inverse_ERF (as usual called INV¦) & inverse_ERFc (as usual called INV¦c), using Newton's method in the subroutine 'inv¦z' when the argument p is a complex number z=a+ib, with b≠0. However, for argument abs(p) a real number in [0 1], in case of inverse_ERF (called INV¦), & for argument (p) a real number in [0 2], in case of inverse_ERFc (called INV¦c), the usual formulae were left. Please note that, for precision purpose, erf(z=a+ib, with b≠0) is imited to the following argument condition: abs(a)<= 9 and abs (b)<= 9. As the calculation of 'INV¦' use erf(z=a+ib), INV¦(z=a+ib, with b≠0) should be theorically limited to the following argument condition: abs(a)<= 9 and abs (b)<= 9. But, in order to avoid an intermediary result > 9 with Newton's method (in subroutine 'inv¦z'), the present version 10c sets for INV¦(z=a+ib, with b≠0) the following more restrictive argument condition: abs(a)<= 8.8 and abs (b)<= 8.8. The same applies to INV¦(z=a+ib, with b≠0), whose argument condition is now limited, in this version 10d, to: abs(a)<= 8 and abs (b)<= 8. RE: RE: HP48-HP50G : Error functions ERF(x or z), ERFC (x or z) & iERF/iERFc —>x - Gil - 12-20-2023 06:00 PM New version 11b All the explanations relative to the argument were somewhat changed. RE: RE: HP48-HP50G : Error functions ERF(x or z), ERFC (x or z) & iERF/iERFc —>x - Gil - 12-21-2023 06:50 AM New version 11c Slight changes (uniformization) in the DOERR messages. Subroutine inv¦z was completed (in bold) for special cases. The rest remains mostly the same for the functions ERF & ERFc, their inverse INV¦ & INV¦c & PROBIT ERFi iERFc ERFcx FADDEEVA DAWSON according to Wolfram Alpha or Wikipedia definitions. As previously, all of them take a single argument input and accept, as argument input, a real or complex number, according to the explanation given at the beginning of each function. RE: RE: HP48-HP50G : Error functions ERF(z), inv, PROBIT, etc. - Gil - 12-23-2023 10:36 AM New version: 12c The 7 main changes: 1) As all other functions, PROBIT function accepts now, as argument, a complex number. Reason: PROBIT function is defined as k*Inverse_ERF, and the function Inverse_ERF does accept — theorically — a complex number as argument. 2) As PROBIT, the functions Inverse_ERF (called 'INV¦') and Inverse_ERFc (called 'INV¦c') use Newton's recursive algorithm when in presence of a complex number argument. Therefore, in order to prevent endless loops, a counter (with a maximum loop by default set at 90) in the subroutine 'inv¦z' has been introduced in the DO... UNTIL. 3) The second criterium to stop the DO... UNTIL in the subroutine 'inv¦z' is relative to finding a converging, acceptable solution. Now, this criterium was "facilitated", with the critical difference (to interrupt the DO... UNTIL routine) between two consecutive complex numbers found (a_n, b_n) and (a_n+1, b_n+1) passing from abs(a_n - a_n+1) & abs(b_n - b_n+1) = 10^-10 to 10^-9. 4) To avoid name confusion, the previously named function iERFc was "more correctly" renamed i¹ERFc, with a special mention at the beginning of the program, as the initial letter i does not stand here for the inverse function (inverse_ERFc). See Wikipedia or Wolfram Alpha for definition and signs, Wolfram definition reversing the terms. Note: in the menu, you will see only the first characters i¹ERF (without the final letter c) of its full, real name i¹ERFc. Reminder: the inverse_ERFc function is called here 'INV¦c' (and called INV¦ for the inverse_ERF function). 5) Similarly to the function i¹ERFc, a new function i²ERFc was created. See Wikipedia or Wolfram Alpha for definition and signs, Wolfram definition reversing the terms. Note 1: in the menu, you will see only the first characters i²ERFc (without the final letter c) of its full, real name i²ERFc. Note 2: the output values of i²ERFc function could not be checked with another program, so be careful when using that function. 6) Some parameters might be exceptionally changed (such as <—max.x, <—max.z <—loop.max or critic). Therefore, their default values are specially reminded, in the beginning of each program, in form of a string, should you change some of the corresponding values and want afterwards to get back their original values. 7) Special attention was given to have the same number of POP and PUSH commands in the different programs, in order to always finish/interrupt a program with 'ENVSTACK' containing, at the very end, an empty list { }. PS Two "cosmetical" corrections in version 12c vs precious (now deleted) version 12b. RE: RE: HP48-HP50G : Error functions ERF(z), inv, PROBIT, etc. - Gil - 12-23-2023 04:02 PM New version 13 The error function directory includes now also Owen's T function. RE: RE: HP48-HP50G : Error functions ERF(z), inv, PROBIT, etc. - Gil - 12-30-2023 10:53 AM New version 14c The directory includes now, besides the following functions ERF(z, z a complex number) & ERFc(z) & their inverses INV¦(z) & INV¦c(z) PROBIT(z) i¹ERFc(z) i²ERFc(z) ERFi(z) ERFcx(z) FADDEEVA(z) DAWSON(z) OWEN(h & a, but both arguments only real numbers) the Fresnel integrals FresnelS & FresnelC with the following restriction: 1 Arg .Real x, |x| less or equal 10 . Or complex (x,0) with |x| less or equal 10 . Or complex z=a+ib, b not equal 0 with |a| & |b| less or equal 5 Arg here normalized: 't^2*pi/2' instead of 't^2' —> lim x—>±infinity = 1/2 instead of 1/2*sqrt (pi) / sqrt(2) Complex form accepted for all function (with the exception of OWEN's T function) : - (a b) - 'a+i*b' - 'r*e^(i*pi)' - 'r*EXP(i*pi)' In this version 14c, OWEN_T function was also improved regarding its 2 arguments (calculation). Most functions are restricted in the argument range (a DOERR message with the corresponding explanation will appear for non valid argument), but all of them accept the special argument ± infinity. Names/labels of subroutines called with small letters & functions names mostly with capitals. Program specially designed for EMU48. RE: RE: HP48-HP50G : Error functions ERF(z), FRESNEL integrals, PROBIT, etc. - Gil - 12-31-2023 11:30 AM New version 14 e Now 2 ways of calculating the Fresnel integrals Fresnel_S & Fresnel_C: a) as in previous version, the normalzed form, with always 1 single input: - the function argument x (real) or z (complex) in level 1; b) new, since this version 14e, the Not normalized version, with then 2 inputs: - the letter N in stack level 2 (with N for Not normalized) - the function argument x (real) or z (complex) in l stack evel 1. RE: RE: HP48-HP50G : Error functions ERF(z), FRESNEL integrals, PROBIT, etc. - Gil - 01-02-2024 04:33 PM New version 15 b Now normalized Fresnel Integrals for real argument x or z=(x, zero), with no limit for x. RE: RE: HP48-HP50G : Error functions ERF(z), FRESNEL integrals, PROBIT, etc. - Gil - 01-03-2024 10:08 AM New version 16 The subroutine with the calculation of FRESNEL Integrals for very large real numbers ('fresnelFASTxLARGE') contained severe errors that were duly fixed up in this version. My apologises for the inconvenience. Now the Fresnel results given by the HP50G are very similar with the ones computed by Wolfram. Besides, in this new version, you can always compute both the normalized and Not normalized Fresnel Integrais F_C & F_S (just adding, when calculating the Non normalized Fresnel Integrals, the letter N [N for Non normalized] in stack level 2, before the argument x or z always in stack level 1). Regarding always the Fresnel integrals: when argument x is a real number : any value from -infinity to +infinity; when argument z is a complex number a+ib, b≠0: abs(a) & abs(b) less or equal to 5 RE: RE: HP48-HP50G : Error functions ERF(z), FRESNEL integrals, PROBIT, etc. - Gil - 01-04-2024 02:28 PM New version 16e Minor changes, mostly in the (result) presentation. In particular, the subroutine 'fresnelFASTxLARGE' for the calculation of FRESNEL Integrals when dealing with very large real numbers was renamed 'FRESNEL.xLARGE.fast', with capital letters at the beginning, as it can be used as an independent function (but only for real numbers x, with abs(x) in [4 infinity). However, to avoid confusion, this function or subroutine 'FRESNEL.xLARGE.fast' was put at the last page of the ERROR directory, as the main function FRESNEL.z to be used normally (at page 2 of the error directory) will automatically access the above mentioned subroutine as soon as x > 10. RE: RE: HP48-HP50G : Error functions ERF(z), FRESNEL integrals, PROBIT, etc. - Gil - 01-05-2024 08:34 AM New version 16.08 3 not at all important changes: a) clearer explanations at the beginning of - FRESNEL.z main program - & subroutine function FRESNEL...; b) change of the name of the above subroutine into: - 'fresnelFASTxLARGE' (at the end, preference was given to small initial letters, indicating clearly that fresnel is a subroutine); c) of course, FRESNEL.z refers now consequently to new name fresnelFASTxLARGE. Though meaningless or useless, here are some results found by the HP50G programs that might be compared with Wolfram Alpha: {"ERFc (33)" 1.93206244519E-475 "ERF (33)" 1. " " "InvERFc (1.93206244519E-475)" 33. " " "Probit (1)" '+infinity' " " "i¹ERFc (5)" :Wiki (Wolframopp): 1.4813429336E-13 " " "i¹ERFc (+infinity )" "Undef by Wolfram as '0/0' (or —>0: try i¹ERFc(x) with x=10 x=20 x=33)" " " "i¹ERFc (-infinity )" :Wiki (Wolfram—>opp): '+infinity ' " " "i²ERFc (25)" :Wiki (Wolfram—>opp): 3.3068803E-277 " " "normaliz.FRESNELS (10)" .468169978589 " normaliz.FRESNELC (10)" .499898694222 " " "fresnelFASTxLARGE" "normaliz.FRESNELS (5000*e^(pi*i) =-5000)" -.499936338031 " normaliz.FRESNELC (5000*e^(pi*i) =-5000)" -.500000002602 " " "fresnelFASTxLARGE" "normaliz.FRESNELS (500000000000)" .5 " normaliz.FRESNELC (500000000000)" .500000000001 " " "normaliz.FRESNELS (-infinity )" '-1/2' " normaliz.FRESNELC (-infinity )" '-1/2' " " "not_norm.FRESNELS (+infinity )" 'sqrt(2*pi)/4' " not_norm.FRESNELC (+infinity )" 'sqrt(2*pi)/4' " " "ERFi ((4.,2.))" (-20442.123626,3999.26730293) " " "ERFi (4+2*i =(4.,2.))" (-20442.123626,3999.26730293) " " "ERFcx ((5.,0.) =5)" .110704637733 " " "w_Faaddeva-Kramp (5)" (1.3887943865E-11,.11524596183) " " "F_Dawson (5)" .102134074424 " " "P(X>h & 0<Y<aX): Owen_T (h=5, a=3)" 1.4332578594E-7 Enjoy the use of the differents function with your HP50G phone application! |