Posts: 1,755
Threads: 112
Joined: Jan 2005
Hi all:
The recent release of the HP15C LE is a momentous event which truly deserves a great celebration on my part, being as I am a lifelong sworn fan of this wonderful model (most especially today which is the oneandonly 111111 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 twobit
HP15C (& LE!) 111111 MiniChallenge
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 manybyte affair.
However, on the other hand, the ugly duckling would be numerical derivation which, unlike integration, can usually be done by numerically evaluating a symbolicallyobtained derived function y'(x) at any specified X ordinate. Regrettably, that requires some advanced symbolic manipulation which is simply unavailable in the HP15C (&LE!) and thus numerical approximations is the only way to go.
So the challenge is:
You're asked to write a program for the HP15C (& LE!) which must be as short as possible, and which can evaluate the derivative at a given ordinate X of a usersupplied function y(x), subject to the following mandatory requirements.
The exact requirements your program must meet are:
 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.
 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 usersupplied function which must be accessible under LBL E.
 it must not depend on or require any particular setting or any particular memory allocation to be previously stablished.
 it must not require any prestored constants in the registers or stack contents except for the ordinate X being initially in the display at runtime.
 and last but not least, it must strive to produce full 10digit accuracy (except for the usual cases of particular troublesome points or ranges, etc, of course).
Sample of use:
 given y(x) = e^{x}/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
 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
 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.
.
Posts: 2,761
Threads: 100
Joined: Jul 2005
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.
Posts: 97
Threads: 9
Joined: Nov 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
Posts: 1,089
Threads: 32
Joined: Dec 2005
Looks like Gerson ingeniously implemented this.
Posts: 2,761
Threads: 100
Joined: Jul 2005
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 HP42S solution:
00 { 27Byte Prgm }
01 LBL "DRV"
02 1E9
03 RCL* ST Y
04 STO 00
05 COMPLEX
06 XEQ "FN"
07 COMPLEX
08 RCL/ 00
09 END
Example:
00 { 8Byte Program }
01 LBL "FN"
02 LN
03 RTN
1E490 XEQ DRV > 1.E490
XEQ DRV > 1.E490
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.
Posts: 2,761
Threads: 100
Joined: Jul 2005
It does work. That was Valentin's idea, though. I've only figured out which it was, I think.
Regards,
Gerson.
Posts: 97
Threads: 9
Joined: Nov 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
Posts: 2,761
Threads: 100
Joined: Jul 2005
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 ComplexStep Derivative Approximation
Regards,
Gerson.
Posts: 163
Threads: 7
Joined: Jul 2007
Or,
00 { 25Byte Prgm }
01 LBL"DER"
02 1E9
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
Posts: 3,283
Threads: 104
Joined: Jul 2005
This may help to shorten the f' function in WP 34S which is implemented in user code. Pauli?
Posts: 2,761
Threads: 100
Joined: Jul 2005
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 { 27Byte Prgm }
01 LBL "DER"
02 1E9
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.
Posts: 1,089
Threads: 32
Joined: Dec 2005
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.
Posts: 653
Threads: 26
Joined: Aug 2010
What about x = 0? This currently means that h = 0 as well, throwing a division by zero error. Setting h = 1E9 in this case is not really an option. ;)
Dieter
Posts: 2,761
Threads: 100
Joined: Jul 2005
Ooops!
00 { 28Byte Prgm }
01 LBL "DER"
02 1E9
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.
Posts: 653
Threads: 26
Joined: Aug 2010
Are negative values no practical cases ?)
What about x = 1 000 000? This would set h = 1E9 as well.
Dieter
Posts: 2,761
Threads: 100
Joined: Jul 2005
Back from a nap... which might explain it (don't program when you're sleepy :)

Back to the drawing board:
00 { 29Byte Prgm }
01 LBL "DER"
02 1E9
03 RCL* ST Y
04 X=0?
05 LASTX
06 STO 00
...
Or on the HP15C:
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.00000090
GSB B > 1.000000 90
Any idea to save a step or two?
Edited: 12 Nov 2011, 7:57 p.m.
Posts: 3,229
Threads: 42
Joined: Jul 2006
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
Posts: 735
Threads: 34
Joined: May 2007
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
Posts: 735
Threads: 34
Joined: May 2007
It appears that Valentin took the first example from this paper:
Kind regards
Thomas
Posts: 653
Threads: 26
Joined: Aug 2010
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 = 1E9 will do. Remember, h is not added to x, it's just used for its imaginary part.
Example:
x = 1234567890
h = 1E9
ln (1234567890 + 1E9 i)
= 2,9339868592 + 8,10000007371E19 i
So
f' = 8,10000007371E19 / 1E9
= 8,10000007371E10
= 1/x
It even works down to h = x * 1E99 resp. x * 1E499.
Or, in other words: your very first solution was fine already. ;)
Dieter
Edited: 13 Nov 2011, 8:09 a.m.
Posts: 2,761
Threads: 100
Joined: Jul 2005
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.
Posts: 653
Threads: 26
Joined: Aug 2010
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
Posts: 163
Threads: 7
Joined: Jul 2007
I thought so too, at first, but:
take f(x)= LN(x) and x = 1e10 then (with h=1e9):
f(1e10,1e9) = (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=1e499. But then,
Taking h=1e499 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
Posts: 2,761
Threads: 100
Joined: Jul 2005
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.
Posts: 2,761
Threads: 100
Joined: Jul 2005
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.
Posts: 735
Threads: 34
Joined: May 2007
From The ComplexStep 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 1e10 we should be happy with h = 1e5 x. It turns out this isn't exactly enough but it seems h = 1e6 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 = 1e9 in most of the cases.
Cheers
Thomas
Posts: 412
Threads: 40
Joined: Mar 2006
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
1E9
x<>y
XEQ E
x<>y
1E9
/
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...
Posts: 2,761
Threads: 100
Joined: Jul 2005
The first problem is possible as well:
A01 LBL A
A02 1E9
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 HP28S symbolic differentiator instead :)
Regards,
Gerson.
Posts: 1,755
Threads: 112
Joined: Jan 2005
Hi, all:
Thanks for your interest in my HP15C (&LE!) Minichallenge, 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 minichallenges 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=1E6 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., prestoring 1E6 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 R0R1 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 Minichallenge, not a fullfledged article on numerical derivation complete with industrialstrength code ! ... :D
When I was busy concocting this minichallenge, I first implemented it for the HP71B, where the code is as simple as
IMPT(FNF((X,H)))/H
where FNF is the userdefined function f(x) and H is a suitably small constat (1E7, say). I conducted some experiments to see which H would be best, for instance:
10 DEF FNF(X)=5*COS(10*X)+((X2)*X6)*X+10
20 DEF FND(X)=50*SIN(10*X)+3*X*X4*X6
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=1E7 first nails the exact derivative
8 20.2010555444 20.2010555444
9 20.2010555444 20.2010555444
and in virtually all cases H = 1E7 was optimal (except for X near 0, of course, where a smaller H is needed). For the HP15C, which only carries 10 (13) digits instead of 12 (15), H=1E6 does the trick.
That's all. I've got (hopefully) interesting ideas for some dozen challenges and minichallenges as well as for another dozen fullfledged 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.
