integral competition HP50g vs. DM42

[peacecalc]

Hello folks,

and the winner is....

ahh, I not like to spoiler the result. What was the job of the competitors?

a) Of course numerical, because the hp50g lives in both worlds, the dm42 is only working with numbers.
b) The Fresnel-Integrals have to be calculated.
c) as pair on the stack y-stack has the S(x) result and x-stack contains the result for C(x)
d) hp50g is set in approx. mode and fix 4 (is necessary for accuracy).
dm42 gets 1E-5 for the ACC-variable for controlling accuracy...

e) Definitions of the two Integrals: \[ C(x) = \int_0^x \cos\left(\frac{\pi}{2} \cdot t^2\right) dt \qquad\qquad \hbox{and} \qquad\qquad S(x) = \int_0^x \sin\left(\frac{\pi}{2} \cdot t^2\right) dt \qquad\qquad \qquad\qquad \qquad\qquad \]

f) Maybe a further information, these definite integrals can only calculated numerically. Except for x -> +- infinite.
S(3.9) = 0.4752 and C(3.9) = 0.4223;

DM42 beats the Hp50g in being three times faster by equal accuracy!

[Albert Chan]

(08-19-2020 07:28 PM)peacecalc Wrote:  S(3.9) = 0.4752 and C(3.9) = 0.4223;

Using complex numbers, we can do both C(x), S(x) at the same time

XCas> F := int(exp(i*pi/2*x*x),x)     → erf(sqrt(pi)/2 * x * (1-i)) / (1-i)

Since erf(0)=0, we have F(x) = C(x) + i * S(x)

XCas> subst(F, x=3.9)                     → 0.422332710261 + 0.475202402351*i

[peacecalc]

Hello Albert,

thank you for tip, I'll try it with my hp50g, maybe it could speed up.

The DM42 is only able to integrate real numbers as far I know, although it is able to work with complex numbers.

Quote:...we can do both C(x), S(x) at the same time

. Is that true?
I always believe our calcs are only able to that one after the other, as long as they have no multi kernel processors and no parallelising operation system.

[Stevetuc]

(08-20-2020 05:25 AM)peacecalc Wrote:  Hello Albert,

thank you for tip, I'll try it with my hp50g, maybe it could speed up.

The DM42 is only able to integrate real numbers as far I know, although it is able to work with complex numbers.

Quote:...we can do both C(x), S(x) at the same time

. Is that true?
I always believe our calcs are only able to that one after the other, as long as they have no multi kernel processors and no parallelising operation system.

By processing as C(x),S(x) as a single complex number C(x) + i * S(x)
Both results are obtained with a single operation:
https://en.m.wikipedia.org/wiki/Complex_number

[Ángel Martin]

(08-20-2020 10:02 AM)Stevetuc Wrote:  By processing as C(x),S(x) as a single complex number C(x) + i * S(x)
Both results are obtained with a single operation:
https://en.m.wikipedia.org/wiki/Complex_number

Indeed, and as example here's the 41Z routine listing for the calculation.
Input: x in X
Output: S(x) in Y, C(x) in X

Code:

1    LBL "ZCSX"
2    1
3    X<>Y
4    STO 00
5    X^2
6    PI
7    ,5
8    STO 01
9    STO 02
10    *
11    *
12    ISG 02
13    0
14    ZGHF
15    RCL 00
16    ST* Z
17    *
18    END

Execution times: 1.5 s on the 41CL

[Thomas Okken]

(08-20-2020 05:25 AM)peacecalc Wrote:  The DM42 is only able to integrate real numbers as far I know, although it is able to work with complex numbers.

That is correct. INTEG only works with real-valued functions.

[Didier Lachieze]

Hi Thomas, I noticed a small difference between Free42 and the HP-42S: when trying to integrate a complex function on the HP 42S I get an "Invalid Type" error message, however on Free42 2.5.19 I don’t get any error message and the integrate function returns 0 instead which is a bit misleading.

Here is the function I’ve used for testing with LLIM=0, ULIM=5, ACC=0.01:

Code:

00 { 18-Byte Prgm }
01▸LBL "Fp"
02 MVAR "U"
03 0
04 RCL "U"
05 X↑2
06 COMPLEX
07 E↑X
08 END

[Thomas Okken]

(08-21-2020 06:30 AM)Didier Lachieze Wrote:  I noticed a small difference between Free42 and the HP-42S: when trying to integrate a complex function on the HP 42S I get an "Invalid Type" error message, however on Free42 2.5.19 I don’t get any error message and the integrate function returns 0 instead which is a bit misleading.

I agree, that is misleading. I'll fix it in 2.5.20.

[peacecalc]

Hello Albert,

I tried your calculus on the HP50g, but I (or if I've didn't make a mistake by program that, the machine) failed. Hmm, on HP50g is only a subset of capabilities from XCAS installed. It seems to me the HP50g has the same restricting as DM42 towards integrating complex functions.

[Albert Chan]

(08-21-2020 12:19 PM)Thomas Okken Wrote:  I agree, that is misleading. I'll fix it in 2.5.20.

Instead of fixing error message, is it better to extend Free42 to support complex integrals ?

[rprosperi]

(08-22-2020 03:24 PM)Albert Chan Wrote:  

(08-21-2020 12:19 PM)Thomas Okken Wrote:  I agree, that is misleading. I'll fix it in 2.5.20.

Instead of fixing error message, is it better to extend Free42 to support complex integrals ?

Ignoring for a moment the fact that the 42S doesn't do that, and therefore Free42 shouldn't either, are you volunteering to do that?

[peacecalc]

Hello all,

after trying to persuade my hp50g to integrate complex functions in general (that machine has the tools for that: first working with complex numbers and the two mighty MATCH commands). I'm a little bit disapointed about my own skills, because if the integration path isn't a closed curve, the resulting value depends on the curve you choose. Here in our example it is easy because the "curve" is all along the real x-axis (tribute to jimi hendrix, because when iI wrote "all along" my mind tried to set "the watchtower").

If somebody tries to implement a general command for integrating complex functions, that command needs a second input which describes the curve.
Of course it is possible as first steps to implement a) "curves" along the real axis, or b) closed curves, because in that two cases it is easier. Nevertheless that is a nice challenge.

@Albert Chan, do you know how this is implemented in XCAS? So, I don't have to invent the wheel a second time and maybe it is possible for me to translate it in code for the HP50g.

@stevetuc, I don't understand your code, although I'm familiar with RPN, for example what does "ZGHF" mean.

[Didier Lachieze]

(08-23-2020 08:56 AM)peacecalc Wrote:  @stevetuc, I don't understand your code, although I'm familiar with RPN, for example what does "ZGHF" mean.

The 41C code above is not from stevetuc but from Angel Martin, ZHGF is a function from Angel's 41Z Deluxe module, an amazing Complex Number module for the HP-41:
ZHGF - Complex Hypergeometric function - Author: Jean-Marc Baillard

The 41Z Deluxe module image can be loaded on the 41CL, on the DM41X, or on a Clonix module to be used with a standard 41C/CV/CX. It requires the Library4 module to work.

You can find the 41Z Deluxe module manual here.

The Complex Hypergeometric function is described by it's author here, and it's application to calculate the Fresnel Integrals is provided on section 3.b with the 41C program example.

I've tested it on my DM41X and it works fine but as mentioned by Jean-Marc Baillard, the input parameter must be small, compared with your example :
S(3.9) = 0.4752 and C(3.9) = 0.4223
it returns S(3.9) = 0.4752 and C(3.9) = 0.4282
and for x>=4 the result is meaningless.

[peacecalc]

Hello Didier Lachieze,

thank for answering my last question. And very much thanks about showing the sources. Now i know the meaning of 41Z.

@Stevetuc
It is as I supposed: of course it is not possible to integrate a complex value in one step.
The Solution that Jean-Marc Baillard used to calculate the fresnel integrals is a function which is expressed as a finite serie (for |x| < 1) and a approximation |x| < 2. It isn't an numerical integration.

So all calculations with complex numbers are implented as a kind of vector calculus. You have a dual vector space with a special "cross product" rule:

\[ \left(\array{e \\ f}\right) = \left(\array{a \\ b}\right) \cdot \left(\array{c \\ d}\right) = \left(\array{ac - bd \\ bc + ad}\right)\]

So every addition, subtraction and so on doubles the computational costs.
And if you have one processor with one kernel, you have to do one step after the other.

@Thomas Okken: for the estimation of accuracy of complex integrals, why there are problems? It would be enough to show that the |R(z)| < epsilon or is there something I not see? R(z) could be a kind of circle which can small enough?

[Thomas Okken]

(08-23-2020 02:07 PM)peacecalc Wrote:  for the estimation of accuracy of complex integrals, why there are problems? It would be enough to show that the |R(z)| < epsilon or is there something I not see? R(z) could be a kind of circle which can small enough?

I don't know. Do numerical integration schemes for real-valued functions generalize to complex-valued functions just like that? This is well outside my mathematical comfort zone.

I'm also not sure just how useful this would be. We're not talking about true complex integrals here, X would still be on a real interval, so all you're accomplishing is not having to integrate Re(f(x)) and Im(f(x)) separately.

[peacecalc]

Hello Thomas,

thank you for your reply. I'm not shure that we are talking of the same thing, but in usual it is possible to convert every complex integral into the RxR space:

My textbook about complex analysis tells me:

\[ \displaystyle \int_C f(z) dz = \int^{t_1}_{t_o} \left( u(x,y) \frac{dx}{dt} - v(x,y)\frac{dy}{dt}\right) dt + i\cdot \int_{t_o}^{t_1} \left( v(x,y)\frac{dx}{dt} + u(x,y) \frac{dy}{dt} \right) dt
\]

\[z = x + iy; \quad dz = dx + idy; \quad u(x,y) = \Re f(z); \quad v(x,y) = \Im f(z); \quad x= x(t); \quad y=y(t); \quad u,v,x,y,t \in IR \]
C is an arbitray curve in the complex plane and t is the parameter. When our TR is able to operate with complex numbers and even operate with regular expressions (like the CAS from HP50g), then I see the chance to make a RPL program which can integrate (not for all cases) exact complex integrals. For closed curves for instance it is possible to calculate complex residues.

And if not the build in numerical integration will give us results, yes you have to calculate two integrals, but it should works.

And that should be possible for even non-analytical functions, because they are defined that their integral over a closed curve is not zero.

[Albert Chan]

(08-23-2020 08:56 AM)peacecalc Wrote:  @Albert Chan, do you know how this is implemented in XCAS? So, I don't have to invent the wheel a second time and maybe it is possible for me to translate it in code for the HP50g.

I just downloaded giac source, usual.cc, function erf0():

Quote:if (absx<=3){
      // numerical computation of int(exp(-t^2),t=0..x)
      // by series expansion at x=0
      // x*sum( (-1)^n*(x^2)^n/n!/(2*n+1),n=0..inf)

Code:

def hyp1f1(a,b,x):
    s = t = c = 1
    while 1:
        t *= a/(b*c)*x
        if s == s+t: return s
        a, b, c, s = a+1, b+1, c+1, s+t

>>> k = 2/pi**0.5
>>> erf = lambda x: k*x * hyp1f1(1/2, 3/2, -x*x)
>>> CiS = lambda x, erf=erf: (1+1j)/2 * erf((1-1j)/k * x)
>>> CiS(3.9)
(0.42233270272270274 + 0.47520241152278453j)

CiS(3.9) called erf(z = (1-1j) * 3.4562850092657564), |z| ≈ 5 > 3
As expected, we lose many digits precision due to massive cancellation.

For big |z|, XCas use continued fraction form:

Quote:// 2*exp(z^2)*int(exp(-t^2),t=z..inf) = 1/(z+1/2/(z+1/(z+3/2/(z+...))))

int(exp(-t^2),t=z..inf) = sqrt(pi)/2 * erfc(z), we got this:

Code:

def erfc(x, n=20, res=0):   # continued fraction, bottom up
    while n>0: res = n/(x+res); n -= 0.5
    return 0.5*k * exp(-x*x) / (x+res)

>>> CiS(3.9, erf = lambda x: 1 - erfc(x))
(0.42233271026093344 + 0.47520240235068834j)

Or, we could calculate generalized continued fraction top down, without hard-coded maximum n.

Code:

def erfc2(x, eps=1/16):     # continued fraction, top down
    a1, b1, a2, b2 = 1, 0, x, 1
    f1, f2, n = 0, 1/x, 0.5
    while f2 != f2 + (f1-f2)*eps:
        a1, b1, a2, b2 = a2, b2, n*a1+x*a2, n*b1+x*b2
        f1, f2, n = f2, b2/a2, n+0.5
    return 0.5*k * exp(-x*x) * f2

>>> CiS(3.9, erf = lambda x: 1 - erfc2(x))
(0.42233271026093344 + 0.47520240235068834j)

[ttw]

In these comparisons, does the DM42 use the same algorithm as the HP50g? If not, one is comparing algorithm speed as well as hardware.

[peacecalc]

Hello ttw,

I cannot answer this, but I guess the two machines use at least a similar algorithmus.
The manuals of these calcs told us: the integration routine doesn't use the left and right limit of the integration interval. The algorithmus can be a kind of gaussian quadrature, because that one doesn't use the limits on the left and right side.

Romberg and other method use the limits of the integration aera.

But maybe I'm wrong.

[John Keith]

For numerical integration they are most likely the same algorithm, probably even the same code except for minor OS differences. The math routines on HP calculators have been mostly unchanged from the HP-71 Math ROM through the Home mode of the Prime.

If all 12 digits of the results on both calculators are the same, that would tend to confirm that the algorithms are the same.

The main reason that the DM42 is faster than the HP 50 is that the DM42 is running Free42, which is a simulator of the HP42 written in C. A more appropriate comparison would be the original HP42 vs the HP-48g as both calculators were of the same era and use the same processor and (internally) a similar OS.

Another valid comparison would be the DM42 vs the Home mode of the Prime, as the Prime Home mode uses the HP math routines rewritten in C for the ARM processor. The HP 50 is a Saturn emulator running on an ARM processor which is clearly slower than native code.