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

+- HP Forums (https://archived.hpcalc.org/museumforum)
+-- Forum: HP Museum Forums (https://archived.hpcalc.org/museumforum/forum-1.html)
+--- Forum: Old HP Forum Archives (https://archived.hpcalc.org/museumforum/forum-2.html)

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

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) = ex/Sqrt(sin(x)3+cos(x)3), 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)/ex, 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.
```.
```

Re: HP-15C (& LE!) 11-11-11 Mini-Challenge - Gerson W. Barbosa - 11-11-2011

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.

Re: HP-15C (& LE!) 11-11-11 Mini-Challenge - peacecalc - 11-12-2011

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

Re: HP-15C (& LE!) 11-11-11 Mini-Challenge - Thomas Radtke - 11-12-2011

Looks like Gerson ingeniously implemented this.

Re: HP-15C (& LE!) 11-11-11 Mini-Challenge - Gerson W. Barbosa - 11-12-2011

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.

Re: HP-15C (& LE!) 11-11-11 Mini-Challenge - Gerson W. Barbosa - 11-12-2011

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

Regards,

Gerson.

Re: HP-15C (& LE!) 11-11-11 Mini-Challenge - peacecalc - 11-12-2011

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

Re: HP-15C (& LE!) 11-11-11 Mini-Challenge - Gerson W. Barbosa - 11-12-2011

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:

Regards,

Gerson.

Re: HP-15C (& LE!) 11-11-11 Mini-Challenge - Werner - 11-12-2011

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

Re: HP-15C (& LE!) 11-11-11 Mini-Challenge - Marcus von Cube, Germany - 11-12-2011

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

Re: HP-15C (& LE!) 11-11-11 Mini-Challenge - Gerson W. Barbosa - 11-12-2011

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.

OT: Google et al. [Re: HP-15C (& LE!) 11-11-11 Mini-Challenge] - Thomas Radtke - 11-12-2011

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.

Re: HP-15C (& LE!) 11-11-11 Mini-Challenge - Dieter - 11-12-2011

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

Re: HP-15C (& LE!) 11-11-11 Mini-Challenge - Gerson W. Barbosa - 11-12-2011

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.

Re: HP-15C (& LE!) 11-11-11 Mini-Challenge - Dieter - 11-12-2011

Are negative values no practical cases ?-)

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

Dieter

Re: HP-15C (& LE!) 11-11-11 Mini-Challenge - Gerson W. Barbosa - 11-12-2011

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.

Re: HP-15C (& LE!) 11-11-11 Mini-Challenge - Paul Dale - 11-13-2011

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

Re: HP-15C (& LE!) 11-11-11 Mini-Challenge - Thomas Klemm - 11-13-2011

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

Re: HP-15C (& LE!) 11-11-11 Mini-Challenge - Thomas Klemm - 11-13-2011

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

Kind regards

Thomas

Re: HP-15C (& LE!) 11-11-11 Mini-Challenge - Dieter - 11-13-2011

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.

Re: HP-15C (& LE!) 11-11-11 Mini-Challenge - Gerson W. Barbosa - 11-13-2011

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.

Re: HP-15C (& LE!) 11-11-11 Mini-Challenge - Dieter - 11-13-2011

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

Re: HP-15C (& LE!) 11-11-11 Mini-Challenge - Werner - 11-13-2011

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

Re: HP-15C (& LE!) 11-11-11 Mini-Challenge - Gerson W. Barbosa - 11-13-2011

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.

Re: HP-15C (& LE!) 11-11-11 Mini-Challenge - Gerson W. Barbosa - 11-13-2011

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.

Re: HP-15C (& LE!) 11-11-11 Mini-Challenge - Thomas Klemm - 11-13-2011

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:





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

Re: HP-15C (& LE!) 11-11-11 Mini-Challenge - J-F Garnier - 11-13-2011

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

Re: HP-15C (& LE!) 11-11-11 Mini-Challenge - Gerson W. Barbosa - 11-13-2011

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.

HP-15C (& LE!) 11-11-11 Mini-Challenge ORIGINAL SOLUTION & COMMENTS - Valentin Albillo - 11-14-2011

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.

```
```