HP-15C (& LE!) 11-11-11 Mini-Challenge

[Valentin Albillo]

Hi all:

The recent release of the HP-15C LE is a momentous event which truly deserves a great celebration on my part, being as I am a life-long sworn fan of this wonderful model (most especially today which is the one-and-only 11-11-11 date you'll get to see for the next 100 years or so).

However, as I'm currently not in the mood for any great efforts and I'm feeling a little lazy today, methinks you'll have to make do with this two-bit


---

HP-15C (& LE!) 11-11-11 Mini-Challenge


---


We all know that numerical integration is all the rage these days in the forum with plenty of good algorithms being discussed and offered for consideration to the WP34S team and such. The best of those integration algorithms are powerful and complex, a many-byte affair.

However, on the other hand, the ugly duckling would be numerical derivation which, unlike integration, can usually be done by numerically evaluating a symbolically-obtained derived function y'(x) at any specified X ordinate. Regrettably, that requires some advanced symbolic manipulation which is simply unavailable in the HP-15C (&LE!) and thus numerical approximations is the only way to go.

So the challenge is:


---

You're asked to write a program for the HP-15C (& LE!) which must be as short as possible, and which can evaluate the derivative at a given ordinate X of a user-supplied function y(x), subject to the following mandatory requirements.


---

The exact requirements your program must meet are:

1. first and foremost, it must be as short as possible, anything more than a dozen steps or two means back to the drawing board for you. The program must begin with LBL A and end with RTN.

2. it must accept a given real ordinate X in the display and produce the value of the derivative of y(x) at that ordinate, where y(x) is a user-supplied function which must be accessible under LBL E.

3. it must not depend on or require any particular setting or any particular memory allocation to be previously stablished.

4. it must not require any pre-stored constants in the registers or stack contents except for the ordinate X being initially in the display at runtime.

5. and last but not least, it must strive to produce full 10-digit accuracy (except for the usual cases of particular troublesome points or ranges, etc, of course).
   
Sample of use:
1. given y(x) = e^x/Sqrt(sin(x)³+cos(x)³), find y'(1.5)
   

      	define y(x) -> LBL E, e^X, LAST X, LAST X, SIN, 3, y^x, X<>Y, COS, 3, y^x, +, SQR, /, RTN
   
   	1.5, GOSUB A ->  4.053427894, which is the correct result to all 10 digits
   
   
   
   

2. given y(x) = 7*e^(2.1*X), find y'(15)
   

   	define y(x) -> LBL E, 2.1, *, e^X, 7, *, RTN
   
   	15, GOSUB A -> 7.040338081 E14, which is the correct result to all 12 digits
   
   
   
   

3. given y(x) = sin(x)/e^x, find y'(1)
   

   	define y(x) -> LBL E, SIN, LAST X, e^X, /, RTN
   
   	1, GOSUB A -> -0.1107937653, whis is correct to all 10 digits
   
   
   
   

Within a few days I'll post my solution, which is just 11 steps long, meets all five conditions above, and produces the sample results above correct to all digits shown. Relaxing a condition or two could shorten it even further, to just 6 steps or even just 5, but let's first keep all five requirements above and see what YOU can produce ... :)

Enjoy !

Best regards from V.

.

[Gerson W. Barbosa]

Hello Valentin,

I think the following meets all your requirements:

001- f LBL A

002-   9

003-   CHS

004-   10^x

005-   STO 0

006-   Re<>Im

007-   +

008-   GSB E

009-   Re<>Im

010-   RCL/ 0

011- g RTN

Best regards,

Gerson.

[peacecalc]

Hello Gerson,

I'm sorry but I think your code doesn't function in the way you want.

I miss the "-" operation, which is nessecary for calculating: "f(x+h)/h - f(x)/h"

But your idea is genial, parallel computing with the hp 15c!

sincerely
peacecalc

[Thomas Radtke]

Looks like Gerson ingeniously implemented this.

[Gerson W. Barbosa]

My source, actually. That's what Wikipedia and Google are for :-)
Technically we should get rid of the imaginary part in the end, either by setting the complex mode off or by clearing the register X before the second Re<>Im instruction. This takes one extra step, however.

Here is an HP-42S solution:

00 { 27-Byte Prgm }

01 LBL "DRV"

02 1E-9

03 RCL* ST Y

04 STO 00

05 COMPLEX

06 XEQ "FN"

07 COMPLEX

08 RCL/ 00

09 END

Example:

00 { 8-Byte Program }

01 LBL "FN"

02 LN

03 RTN

1E-490 XEQ DRV -> 1.E490

       XEQ DRV -> 1.E-490

You might want to insert X<>Y and Rv between the steps 08 and 09. The answer will be valid only if FN returns a real result for the given X. A few more steps might handle this.

Gerson.

[Werner]

Or,

00 { 25-Byte Prgm }

01 LBL"DER"

02 1E-9

03 COMPLEX

04 ASTO L

05 XEQ IND ST L

06 COMPLEX

07 1E9

08 *

09 END

It does not need a predefined alpha label for the function

It does not use registers, numbered or other. Just the alpha register to pass the function name

Cheers,
Werner

[Gerson W. Barbosa]

Thanks for the tip, Werner. It doesn't give correct results when |x| >= 1e4 for my example (LN), unless a slight modification is made:

00 { 27-Byte Prgm }

01 LBL "DER"

02 1E-9

03 RCL* ST Y

04 STO 00

05 COMPLEX

06 ASTO ST L

07 XEQ IND ST L

08 COMPLEX

09 RCL/ 00

10 END

Cheers,

Gerson.

Edited: 12 Nov 2011, 10:48 a.m.

[Dieter]

What about x = 0? This currently means that h = 0 as well, throwing a division by zero error. Setting h = 1E-9 in this case is not really an option. ;-)

Dieter

[Gerson W. Barbosa]

Ooops!

00 { 28-Byte Prgm }

01 LBL "DER"

02 1E-9

03 X<Y?

04 RCL* ST Y

...

This should handle practical cases. Of course the result is not valid for LN(x) because LN(0) is not defined.

[Dieter]

Are negative values no practical cases ?-)

What about x = -1 000 000? This would set h = 1E-9 as well.

Dieter

[Gerson W. Barbosa]

Back from a nap... which might explain it (don't program when you're sleepy :-)

----

Back to the drawing board:

00 { 29-Byte Prgm }

01 LBL "DER"

02 1E-9

03 RCL* ST Y

04 X=0?

05 LASTX

06 STO 00

...

Or on the HP-15C:

001- f LBL B

002-   ENTER

003-   ENTER

004-   EEX

005-   CHS

006-   9

007-   *

008- g x=0 

009- g LSTx

010-   STO 0

011- f Re<>Im

012-   +

013-   GSB E

014- f Re<>Im

015-   RCL/ 0

016- g CF 8       ; leave complex mode

017- g RTN

018- f LBL E

019- g LN

020- g RTN

 keystrokes       display

EEX 90 GSB B -> 1.000000-90

       GSB B -> 1.000000 90

Any idea to save a step or two?

Edited: 12 Nov 2011, 7:57 p.m.

[Dieter]

Gerson,

I think we are discussing a problem that actually does not exist. ;-) There is no need to make sure that h is a certain fraction of x. And x=0 does not require any special handling either. The method should work for essentially any h, regardless of the magnitude of x. Consider your ln(x) example. It doesn't matter if x is 1 or 50000 or 10^9 - in all cases h = 1E-9 will do. Remember, h is not added to x, it's just used for its imaginary part.

Example:

   x = 1234567890

   h = 1E-9

   ln (1234567890 + 1E-9 i)

     = 2,9339868592 + 8,10000007371E-19 i

So

  f' = 8,10000007371E-19 / 1E-9

     = 8,10000007371E-10

     = 1/x

It even works down to h = x * 1E-99 resp. x * 1E-499.

Or, in other words: your very first solution was fine already. ;-)

Dieter

Edited: 13 Nov 2011, 8:09 a.m.

[Gerson W. Barbosa]

Dieter,

I am just trying to avoid the results below. The second program gets them all right.

   keystrokes           display

.1 GSB A              10.00000000

.01 GSB A             100.0000000

.001 GSB A            1000.000000

.00001 GSB A          10000.00000

EEX CHS 5 GSB A       99999.99967

EEX CHS 6 GSB A       999999.9967

EEX CHS 7 GSB A       9999999.987

EEX CHS 8 GSB A       99668652.49

EEX CHS 9 GSB A       785398163.4

EEX CHS 10 GSB A      1471127674.

EEX CHS 20 GSB A      1570796327.

EEX CHS 30 GSB A      1570796327.

EEX CHS 40 GSB A      1570796327.

...

Regards,

Gerson.

[Dieter]

I see, so the problem arises for values close to zero. On the other hand you previously said "...it doesn't give correct results when |x| >= 1e4 for my example (LN)", so I took a closer look at large values. #-)

This leads to the question if there is a clever way to determine a suitable h that's small enough to ensure sufficient precision as well as large enough to prevent underflow. But I fear this cannot be done in 11 steps or less. ;-)

Dieter

[Gerson W. Barbosa]

Quote:

---

But I fear this cannot be done in 11 steps or less. ;-)

---

Agreed! The steps 11 and 12 in the program above can be replaced by [f] [I] per Thomas Klemm's suggestion. Now it's only 16 steps long and decounting :-)

Gerson.

[Thomas Klemm]

From The Complex-Step Derivative Approximation:

Quote:

---

Taking the imaginary parts of both sides of this Taylor series expansion (7) and dividing it by h yields:

Similarly, for the truncation error of the derivative estimate to vanish, we require that:

---

For this example we have:

The third derivative is:

This leads us to:

or

So with a relative working precision of 1e-10 we should be happy with h = 1e-5 |x|. It turns out this isn't exactly enough but it seems h = 1e-6 |x| is.

Quote:

---

This leads to the question if there is a clever way to determine a suitable h that's small enough to ensure sufficient precision as well as large enough to prevent underflow.

---

I doubt that this can be achieved in general. OTOH we probably run only into problems in the neighborhood of singularities but Valentin stated in requirement 5:

Quote:

---

(except for the usual cases of particular troublesome points or ranges, etc, of course)

---

Therefore I guess we can live with h = 1e-9 in most of the cases.

Cheers

Thomas

[Werner]

I thought so too, at first, but:
take f(x)= LN(x) and x = 1e-10 then (with h=1e-9):

 f(1e-10,1e-9) = (-20.7182906715,1.4711276743)

and the calculated derivative is 1471127674.3 instead of 1e10.

Making h smaller will work, so you might think, let's take the smallest possible h=1e-499. But then,

Taking h=1e-499 returns f'(1e499)=0, as underflow occurs in the calculation of LN(x,h).

Scaling h by x, as Gerson did, works for the whole range.

Cheers, Werner

[Thomas Radtke]

Quote:

---

That's what Wikipedia and Google are for :-)

---

I still smile about it, too, but while I think wikipedia will serve us well now and in the long run, I'm waiting for the day Google renders itself useless for anything off mainstream. It already act at its own authority too much.

[Marcus von Cube, Germany]

This may help to shorten the f' function in WP 34S which is implemented in user code. Pauli?

[Paul Dale]

This method isn't really suitable for the 34s. It requires a complex evaluation of the user supplied function which isn't seamless unfortunately.

The 34s differentiation routine uses 4, 6 and 10 point estimates and returns the 10 point one if all evaluations succeed.

I've not attempted this challenge on the 34S but the code is nice and short:

    LBL A

    f'(x) D

    RTN

Using label D instead of E. Shorter again to not bother coding this and just enter f'(x) D :-)

- Pauli

[Gerson W. Barbosa]

It does work. That was Valentin's idea, though. I've only figured out which it was, I think.

Regards,

Gerson.

[peacecalc]

Hello Gerson,

I beg your pardon Gerson, my lessons in complex functional analysis are 30 years left. I forgot them.

Now, it's a good reason for me repeating that stuff, because it's new for me to use result from this for numerical purposes.

sincerely
peacecalc

[Gerson W. Barbosa]

Hello,

Same here. As I said, I've only implemented the formula in the Wikipedia article above. The following paper provides a detailed explanation of the method:

The Complex-Step Derivative Approximation

Regards,

Gerson.

[Thomas Klemm]

It appears that Valentin took the first example from this paper:

Kind regards

Thomas

[Thomas Klemm]

Saved one step:

001 - 42,21,11   LBL A

002 -       26   EEX

003 -        9   9

004 -       16   CHS

005 -    44  0   STO 0

006 -   42  25   I

007 -   32  15   GSB E

008 -   42  30   Re<>Im

009 - 45,10, 0   RCL ÷ 0

010 -   43  32   RTN

Nice solution!

Cheers

Thomas

[Gerson W. Barbosa]

Quote:

---

Saved one step

---

Great!

This should be closer to Valentin's original solution. I guess his solution includes

10- 53, 5, 8  CF 8

This is required to get rid of the imaginary part. Try, for instance, squaring f'(pi/6) when f(x) = sin(x). The correct result is 3/4 not 1/2 + sqrt(3)/2*i. When I saw your listing, I imagined it would not work if the complex mode had not been previously set. I had tried a similar program, except that I included SF 8 in the beginning (I didn't remember [f] [I] activates the complex mode, like Re<>Im does).

Cheers,

Gerson.

Edited: 13 Nov 2011, 9:32 a.m.

[J-F Garnier]

I learnt something today, thanks to Valentin.

For the fun, here is a solution of the last two problems for the HP32SII. At least I found a use of the HP32SII complex mode!

LBL A

1E-9

x<>y

XEQ E

x<>y

1E-9

/

RTN

LBL E  ( 2nd problem )

0

2.1

CMPLX*

CMPLXe^x

0

7

CMPLX*

RTN

LBL E  ( 3rd problem )

0

ENTER

CMPLX+  

CMPLXSIN

Rv

Rv

CMPLXe^x

CMPLX/

RTN

I didn't succeed to write it on the HP35s...

[Gerson W. Barbosa]

The first problem is possible as well:

A01 LBL A

A02 1E-9

A03 STO B

A04 x<>y

A05 XEQ E

A06 x<>y

A07 RCL/ B

A08 RTN

E01 LBL E

E02 STO A

E03 RCL B

E04 x<>y

E05 CMPLXCOS

E06 RCL B

E07 RCL A

E08 CMPLXSIN

E09 CMPLX+

E10 0

E11 3

E12 CMPLX*

E13 3

E14 RCL* B

E15 3

E16 RCL* A

E17 CMPLXCOS

E18 CMPLX+

E19 3

E20 RCL* B

E21 3

E22 RCL* A

E23 CMPLXSIN

E24 CMPLX-

E25 0

E26 4

E27 CMPLX/

E28 2

E29 RCL* B

E30 2

E31 RCL* A

E32 CMPLXe^x

E33 CMPLX/

E34 CMPLX1/x

E35 0

E36 ENTER

E37 0.5

E38 CMPLXy^x

E39 RTN

1.5 XEQ A -> 4.05342789391

Better use the original expression and save the intermediate results to registers. Or use the HP-28S symbolic differentiator instead :-)

Regards,

Gerson.

[Valentin Albillo]

Hi, all:

Thanks for your interest in my HP-15C (&LE!) Mini-challenge, I'm glad many of you enjoyed it and produced a number of correct solutions, some of which virtually duplicated my original one, as well as producing solutions for other models. That's the spirit of these challenges and mini-challenges and I thank you for joining in.

As for my original solution, this is it:

  01  LBL A

  02  EEX

  03   6

  04  CHS

  05  STO I

  06   I

  07  GSB E

  08  Re<>Im

  09  CF 8

  10  RCL/ I

  11  RTN

which is 11 steps and meets all 5 requirements. I prefer to use h=1E-6 instead of smaller values and I think that CF 8 must be included lest you'll leave the machine in complex mode, which would be utterly unexpected to the user as both the input and the output are assumed and expected to be real.

Partially relaxing condition (1), i.e., supressing LBL and RTN, would save two steps, while relaxing condition (4), i.e., pre-storing 1E-6 in a register, would save three additional steps, leaving the grand total at a mere 6 steps or so. An additional step could be saved by supressing CF 8 but as stated above, I think this would seriously detract from the usability.

As for the register used to store the constant, I prefer using permanent register RI instead of R0 or R1 (also permanent) though, as the constant is less than one, the best solution is using register RAN#, which frees RI for indexing purposes and R0-R1 for general use or matrix operations. However, as there's no RCL\ RAN# instruction this means an additional step to perform the division.

Finally, I knew that for X arguments near 0 accuracy would get worse and so require extra steps to cater for that case, thus this is why I included the exception at condition (5) in order to keep it short and sweet. After all, this was intended as a Mini-challenge, not a full-fledged article on numerical derivation complete with industrial-strength code ! ... :D

When I was busy concocting this mini-challenge, I first implemented it for the HP-71B, where the code is as simple as

    IMPT(FNF((X,H)))/H

where FNF is the user-defined function f(x) and H is a suitably small constat (1E-7, say). I conducted some experiments to see which H would be best, for instance:

  10 DEF FNF(X)=5*COS(10*X)+((X-2)*X-6)*X+10

  20 DEF FND(X)=-50*SIN(10*X)+3*X*X-4*X-6

  30 FOR I=1 TO 9 @ H=10^(-I) @ DISP I;IMPT(FNF((1,H)))/H,FND(1) @ NEXT I

  >RUN

   1  24.9567129442     20.2010555444

   2  20.24631331       20.2010555444

   3  20.2015078976     20.2010555444

   4  20.201060068      20.2010555444

   5  20.2010555897     20.2010555444

   6  20.2010555449     20.2010555444

   7  20.2010555444     20.2010555444  << H=1E-7 first nails the exact derivative

   8  20.2010555444     20.2010555444

   9  20.2010555444     20.2010555444

and in virtually all cases H = 1E-7 was optimal (except for X near 0, of course, where a smaller H is needed). For the HP-15C, which only carries 10 (13) digits instead of 12 (15), H=1E-6 does the trick.

That's all. I've got (hopefully) interesting ideas for some dozen challenges and mini-challenges as well as for another dozen full-fledged articles but regrettably next to no or very little time to write them down & post them hera, and actually no suitable place to publish the long articles so it might be a while till I make another lengthy contribution ... :(

Thanks again for your interest and best regards from V.