A Tiny (WP34S) program does a lot

[BarryMead]

While looking at a graph of the "Factorial" Gamma(X+1) function I noticed that between the points 0 and 1 where the value of the function is 1 (one), the function dips below one. I was curious how to find the minimum value between these two points. Using SLV with the f'(x) function made it quite easy. This tiny program illustrates how powerful a small WP-34S program can be. It took only 15 program steps on the WP-34S and over 40 on the HP-15C.

Code:


LBL A                //Entry point for the program
# 000                //Lower search limit 0
# 001                //Upper search limit 1
SLV 00               //Find the zero point of the derivative
RTN                  //Return if SLV succeeds
ERR 20               //If SLV fails show "No root Found"
LBL 00               //Function to find the first derivative
f'(x) 01
RTN
LBL 01               //Function to calculate the Factorial of value
x!
RTN
LBL'[DELTA]X'        //Function to set the Delta X of derivative
# 001
SDR 005              //Set Delta X to .00001
END

It takes a while to run on the real WP-34S (about 46 seconds), but is almost instantaneous on the emulator.
Obviously when a function reaches a "Minimum" value, its first derivative (slope) will be zero.
Since I don't know what the first derivative of the Factorial Gamma(X+1) looks like, I used the f'(x) function
to approximate it. Note: The Optional '[Delta]X' routine helps to improve the accuracy of the result.
And the answer is: 0.4616321449683623 where the factorial reaches it's minimum value of: 0.8856031944108887
Thanks to Paul Dale and Dieter for the excellent enhancements to the program.

[Paul Dale]

The derivative of the gamma function is the digamma function (ψ), well is trivial to get the derivative from this. This function was implemented both in C and in XROM but it was omitted from the final firmware. It remains a compile time option.

Oh, SLV is a condition, it skips a step if it can't find a solution. Better would be two RTN statements after the SLV not one.


- Pauli

[Dieter]

(11-28-2014 03:16 PM)BarryMead Wrote:  While looking at a graph of the "Factorial" function

The x! function actually evaluates Γ(x+1), just as the 34C and 15C did. For non-negative integer arguments this agrees with the factorial function.

(11-28-2014 03:16 PM)BarryMead Wrote:  And the answer is: .4616321449682268 where the factorial reaches it's minimum value of: .8856031944108887

For the record: the exact 16-digit value is x = 0.4616321449683623 with a minimum of 0.8856031944108887.

Dieter

[BarryMead]

(11-28-2014 10:29 PM)Dieter Wrote:  For the record: the exact 16-digit value is x = 0.4616321449683623 with a minimum of 0.8856031944108887.

What method/software did you use to obtain that result. I tried various values for '[Delta]X' and was not able to get any more accuracy with smaller values.

I did notice that the value you quoted for the factorial was the same as my findings. Perhaps the result was
in error because it could not be discriminated from the correct value in the SLV function which probably stops
searching when it sees no change in the value returned from the derivative function.

[BarryMead]

(11-28-2014 10:05 PM)Paul Dale Wrote:  Oh, SLV is a condition, it skips a step if it can't find a solution. Better would be two RTN statements after the SLV not one.

Thanks, I forgot about the conditional skip on return from the SLV function. I corrected the program above to include this extra RTN instruction.

[BarryMead]

I am VERY IMPRESSED with how much "behind the scenes" computation power is going on with these
simple looking WP-34S functions. Essentially the calculator sets up a double nested loop submitting
intelligent guesses to the outer function, which calls an inner loop to compute the first derivative of the guesses.

The result is empowering, even though I did not know what the derivative (DiGamma X+1) of the Factorial function (Gamma of X+1) was, the calculator was able to figure it out with it's intelligent internal looping. Even though I chose
a bad example for a function with no algebraically defined derivative, there must be functions with no known derivative to which this concept could be applied.

Thanks for your informative replies. I learn something new every day.

[BarryMead]

(11-28-2014 10:29 PM)Dieter Wrote:  For the record: the exact 16-digit value is x = 0.4616321449683623 with a minimum of 0.8856031944108887.

When I put the calculator into Double Precision mode, I did in fact get the 16-digit answer you reported above.

I corrected the original posting to show the correct value you reported.

[Gerson W. Barbosa]

(11-28-2014 10:05 PM)Paul Dale Wrote:  The derivative of the gamma function is the digamma function (ψ), well is trivial to get the derivative from this. This function was implemented both in C and in XROM but it was omitted from the final firmware. It remains a compile time option.

Oh, SLV is a condition, it skips a step if it can't find a solution. Better would be two RTN statements after the SLV not one.


- Pauli

On the HP 50g, which has digamma built-in, I remember I just solved 'Psi(X)=0.' for X:

X = 1.46163214497 and GAMMA(X), or (X - 1)!, = 0.885603194411

If only there were a closed form solution...

8*EulerGamma/(10 + 1/(2*(163 + √2)))

Gerson.

[Gerson W. Barbosa]

The above approximation can easily be evaluated on the HP 50g:

white-shift MTH [NXT] SPECI 1 Psi [+/-] 8 * 2 [sqrt] 163 + 2 * [1/X] 10 + /

--> 0.461632144989

[Dieter]

(11-29-2014 12:34 AM)BarryMead Wrote:  What method/software did you use to obtain that result.

Wolfram Alpha. Simply enter "Minimum of Gamma(x+1) near 1.5". ;-)

(11-29-2014 12:34 AM)BarryMead Wrote:  I did notice that the value you quoted for the factorial was the same as my findings. Perhaps the result was in error because it could not be discriminated from the correct value in the SLV function which probably stops searching when it sees no change in the value returned from the derivative function.

It will stop as soon as the solved function becomes zero. Depending on the second derivative, there may a more or less wide range of 16-digit values that return the same function result near or equal to zero.

However, this is not quite the case here – the second derivative is roughly 1 at this point, so changing x by 1 ULP = 10^–16 will change the first derivative by a similar amount. Actually the exact first derivative (i.e. the digamma function) of your original result is not quite zero:
digamma(1.4616321449682268) = –1.31 · 10^–13

The sign actually changes between
digamma(1.4616321449683623) = –3,99 · 10^–17
and
digamma(1.4616321449683624) = +5,68 · 10^–17

The difference probably is due to the fact that the 34s f'(x) function provides just an approximation of the first derivative.

Since you posted a more precise result: Wolfram Alpha returns
x = 0.46163 21449 68362 34126 26595 42325 72132 84682
y = 0.88560 31944 10888 70027 88159 00582 58873 32080

Dieter

[Dieter]

(11-29-2014 12:39 AM)BarryMead Wrote:  

(11-28-2014 10:05 PM)Paul Dale Wrote:  Oh, SLV is a condition, it skips a step if it can't find a solution. Better would be two RTN statements after the SLV not one.

Thanks, I forgot about the conditional skip on return from the SLV function. I corrected the program above to include this extra RTN instruction.

In this case the user will not notice if SLV was not able to find a solution.
This should work better:

Code:

  LBL A
  #000    // lower guess = 0
  #001    // upper guess = 1
  SLV 00  // solve for derivative = 0
  RTN     // quit if solution found
  ERR 20  // otherwise  show error message "no root found"

Dieter

[BarryMead]

(11-29-2014 06:42 AM)Dieter Wrote:  This should work better:

Code:

  LBL A
  #000    // lower guess = 0
  #001    // upper guess = 1
  SLV 00  // solve for derivative = 0
  RTN     // quit if solution found
  ERR 20  // otherwise  show error message "no root found"

Thank you Dieter. I am learning a lot. (I never knew there were so many canned error codes in the
calculator, and the opcode count went back down to 17 with the # inputs) Great Suggestions.

[Thomas Klemm]

You could as well use the complex factorial to calculate the derivative:

Code:

LBL'DFX'
RCL 00
x<>y
CPX x!
x<>y
RCL/ 00
RTN

1e-8 STO 00
0 ENTER
1
SLV'DFX'

Kind regards
Thomas

[Paul Dale]

(11-29-2014 10:06 AM)BarryMead Wrote:  and the opcode count went back down to 17 with the # inputs)

The delta X routine can save a step and run faster too:

Code:

# 001
SDR 005

And another step can be saved by leave out the RTN before the END. END does the RTN for you.


- Pauli

[BarryMead]

(11-29-2014 11:03 AM)Paul Dale Wrote:  The delta X routine can save a step and run faster too:

Code:

# 001
SDR 005

And another step can be saved by leave out the RTN before the END. END does the RTN for you.

- Pauli

The "Focal" language is so rich with hidden features and cool shortcuts. A WHOLE BOOK could be
written on the "Focal" language and it's hidden treasures.

[BarryMead]

(11-28-2014 10:05 PM)Paul Dale Wrote:  The derivative of the gamma function is the digamma function (ψ), well is trivial to get the derivative from this.

Paul: How did you get the greek letter PSI into this posting?
Do you have a Unicode table posted on your wall?
What is the ISO spec for the Unicode version used by this fourm?

I found a "Unicode Test" page in test forum with a full Unicode listing made by Katie Wasserman. I guess one could simply
cut and paste from that page to get the greek unicode characters.

Is there an easier way?

[Dieter]

(11-29-2014 07:44 PM)BarryMead Wrote:  Paul: How did you get the greek letter PSI into this posting?

The same way I got the Γ into mine, I assume. ;-)

(11-29-2014 07:44 PM)BarryMead Wrote:  Do you have an extended ASCII table posted on your wall?
What is the ISO spec for this table so that I can access greek characters in these postings?

This forum supports Unicode, so inserting special characters is really simple: If you are using some kind of Windows, simply use "charmap" for an overview of the available characters. Set a standard font like Arial, select "Unicode" as character set and browse the table. Or search for "Psi". This will display all characters with "psi" in their description (so also ePSIlon gets displayed). Select, copy and paste the desired character, and you're done: Ψ.

There also is a comprehensive overview of unicode characters on Wikipedia.

Dieter

[walter b]

(11-29-2014 07:44 PM)BarryMead Wrote:  

(11-28-2014 10:05 PM)Paul Dale Wrote:  The derivative of the gamma function is the digamma function (ψ), well is trivial to get the derivative from this.

Paul: How did you get the greek letter PSI into this posting?

I use to copy the respective characters from a file I know (e.g. a manual ). All the special characters on 1 or 2 pages ...

Just my 20 m€.
d:-)

[Dieter]

(11-29-2014 08:18 PM)walter b Wrote:  I use to copy the respective characters from a file I know (e.g. a manual ).

Wow! So you copied your complete signature...

Quote:Μάλα γὰρ φιλοσόφου τοῦτο τὸ πάθος, τὸ θαυμάζειν: οὐ γὰρ ἄλλη ἀρχὴ φιλοσοφίας ἢ αὕτη.

...from the 34s manual ?-)

Dieter

[walter b]

(11-29-2014 08:41 PM)Dieter Wrote:  Wow! So you copied your complete signature...

Quote:Μάλα γὰρ φιλοσόφου τοῦτο τὸ πάθος, τὸ θαυμάζειν: οὐ γὰρ ἄλλη ἀρχὴ φιλοσοφίας ἢ αὕτη.

...from the 34s manual ?-)

Yes, in principle. Almost. I certainly used a Word file. Accented Greek letters aren't contained in the WP 34S but some are in the 43S manual.

d:-)