HP Forums

Full Version: HP 50g Double factorial
You're currently viewing a stripped down version of our content. View the full version with proper formatting.
Pages: 1 2
Hi,

currently I extend my 50g user's guide with additional information. Among others, I would like to add a User RPL program to calculate the double factorial (https://en.wikipedia.org/wiki/Double_factorial).
My version:
Code:

@program DF
@in: number > 0
@out: double factorial of number
\<<
  \-> N
  \<<
     IF 1 > THEN @valid number?
       N 2 MOD @even or odd?
       IF 0 == THEN @even number
          2 N 2 / ^ N 2 / ! *
       ELSE @odd number
          N 1 + ! 2 N 1 + 2 / ^ N 1 + 2 / ! * /          
       END
     ELSE @not valid
       1
     END 
  \>>
\>>
It based on the fourmulars n!!=2^(n/2)*(n/2)! for even numbers and n!!=(n+1)!/2^(n+1/2)*(n+1/2)! for odd numbers and it works, but can you tell me where you could optimize it or how you would solve this task?
But no additional libraries should be used and it should be a User RPL program.

Many thanks and kind regards
Joerg
Hi Joerg , another way is :

Code:
 « -> n 'IFTE(n MOD 2,(n+1)!/(2^((n+1)/2)*((n+1)/2)!),2^(n/2)*(n/2)!)'»

returns 1 for zero, and an error for n<0
And works in symbolic mode with a symbolic variable
Hi Gilles,

(05-01-2019 01:40 PM)Gilles Wrote: [ -> ]
Code:
 « -> n 'IFTE(n MOD 2,(n+1)!/(2^((n+1)/2)*((n+1)/2)!),2^(n/2)*(n/2)!)'»
returns 1 for zero, and an error for n<0
And works in symbolic mode with a symbolic variable

thanks for the answer with the short and sweet version.
I like that a question mark is displayed for negative numbers. However, at -1 and -2 not. This is certainly related to the formula?

Kind regards
Joerg
(05-01-2019 03:57 PM)joeres Wrote: [ -> ]thanks for the answer with the short and sweet version.
I like that a question mark is displayed for negative numbers. However, at -1 and -2 not. This is certainly related to the formula?

Kind regards
Joerg

Yes it is. We can do :

Code:
«  -> n 
 «
  IF 'n<0 OR FP(n)' THEN  #203h DOERR ELSE
   'IFTE(n MOD 2,(n+1)!/(2^((n+1)/2)*((n+1)/2)!),2^(n/2)*(n/2)!)' ->NUM
  END 
 »
»

A negative or not integer input will result an error "Bad argument value".
Gilles, I'm curious. What are these rectangles in your post supposed to be?

Ces petits rectangles sont censés représenter quoi au juste (voir image attachée) ?

[attachment=7191]
Here's another solution:

Considering that  n!!  is equivalent to  n·(n−2)·(n−4)· ... ·1  for every positive integer  n ≠ 0  (odd or even), all it needs is a simple loop which calculates the next factor and multiplies it with the running product.

1             @ Initialize result.
SWAP 2 FOR f  @ Loop down from n to 2.
  f *         @ Multiply with current factor.
-2 STEP       @ Decrement factor by 2 ; repeat until factor < 2.

Taking into account the special case where  n = 0 , we finally have:

Code:

@ NAME  : DF (Double factorial)
@ STACK : ( n --> n!! )
\<<
  1                 @ Initialize result.
  IF SWAP DUP 0 ==  @ If n = 0,
  THEN DROP         @   then n!! = 1 (just drop n, 1 already on stack).
  ELSE 2 FOR f      @ Otherwise: Loop down from n to 2.
    f *             @   Multiply with current factor.
  -2 STEP END       @   Decrement factor by 2 ; repeat until factor < 2.
\>>
Hi joeres, the program you posted contains a small bug:

After having created and initialized the local variable N, you need to put it on the stack again for comparison, because \-> consumes its argument.

So, instead of:

\<<
  \-> N
  \<<
     IF 1 > THEN

(...)

you have to do:

\<<
  \-> N
  \<<
     IF N 1 > THEN

(...)
(05-01-2019 04:14 PM)Gilles Wrote: [ -> ]Yes it is. We can do :

Code:
«  -> n 
 «
  IF 'n<0 OR FP(n)' THEN  #203h DOERR ELSE
   'IFTE(n MOD 2,(n+1)!/(2^((n+1)/2)*((n+1)/2)!),2^(n/2)*(n/2)!)' ->NUM
  END 
 »
»

A negative or not integer input will result an error "Bad argument value".

Thanks Gilles, it works fine.

Kind regards
Joerg
(05-01-2019 05:05 PM)Giuseppe Donnini Wrote: [ -> ]
1          @ Initialize result.
SWAP 2     @ Loop down from n to 2.
FOR f f *  @ Multiply with current factor.
-2 STEP    @ Decrement factor by 2 ; repeat until factor < 2.

Taking into account the special case where  n = 0 , we finally have ...

I am not familiar with RPL code, but is there a loop construct where you test the condition first ?
That way, the 0 special case need not be added.

BTW, I like your code. Doing odd n!! as n! / (n-1)!! might overflow the result.
(05-01-2019 05:05 PM)Giuseppe Donnini Wrote: [ -> ]Here's another solution:

Considering that  n!!  is equivalent to  n·(n-2)·(n-4)· ... ·1  for every positive integer  n ≠ 0  (odd or even), all it needs is a simple loop which calculates the next factor and multiplies it with the running product. Such a solution could do entirely without local variables and algebraic objects, requiring only stack manipulations, and very few in number:

1          @ Initialize result.
SWAP 2     @ Loop down from n to 2.
FOR f f *  @ Multiply with current factor.
-2 STEP    @ Decrement factor by 2 ; repeat until factor < 2.

Taking into account the special case where  n = 0 , we finally have:

Code:

@ NAME  : DF (Double factorial)
@ STACK : ( n --> n!! )
\<<
  1                 @ Initialize result.
  IF SWAP DUP 0 ==  @ If n = 0,
  THEN DROP         @   then n!! = 1 (just drop n, 1 already on stack).
  ELSE 2            @ Otherwise: Loop down from n to 2.
  FOR f f *         @   Multiply with current factor.
  -2 STEP           @   Decrement factor by 2 ; repeat until factor < 2.
\>>

Hi Giuseppe,

Thank you, very interesting and nice solution. I tested it (END added). Works fine, but if I use very large numbers, the program runs in an infinite loop. The other programs output the message "Integer too large". But I can easily query this.

Kind regards
Joerg
(05-01-2019 05:50 PM)Giuseppe Donnini Wrote: [ -> ]Hi joeres, the program you posted contains a small bug:

After having created and initialized the local variable N, you need to put it on the stack again for comparison, because \-> consumes its argument.

So, instead of:

\<<
  \-> N
  \<<
     IF 1 > THEN

(...)

you have to do:

\<<
  \-> N
  \<<
     IF N 1 > THEN

(...)

Many thanks for the hint. I have corrected it.

Kind regards
Joerg
(05-01-2019 04:14 PM)Gilles Wrote: [ -> ]A negative or not integer input will result an error "Bad argument value".

Not necessarily. Smile See this thread.

Here is a version for the HP-48 or HP 50:

Code:

\<< \pi EVAL DUP2 * COS 1. - .25 * SWAP .5 * SWAP ^ SWAP .5 * DUP 2. SWAP ^ SWAP 1. + GAMMA * *
\>>

It will handle non-integer arguments, and complex numbers if complex mode is enabled. Negative even integers return results near negative infinity.

See the graphs here.
(05-01-2019 06:26 PM)joeres Wrote: [ -> ]If I use very large numbers, the program runs in an infinite loop.
I cannot confirm this behavior. If n!! is too large, I get MAXR, i.e. 9.99999999999E499.
You must be using ZINTs, then. As a matter of fact, I wrote the code with an HP-48 in mind, so it is primarily intended for floats.

(05-01-2019 06:26 PM)joeres Wrote: [ -> ]The other programs output the message "Integer too large".
Again, this must be a ZINT-specific message.

P.S. Thanks for reminding me of the missing "END", I just added it.
(05-01-2019 04:59 PM)grsbanks Wrote: [ -> ]Gilles, I'm curious. What are these rectangles in your post supposed to be?
Ces petits rectangles sont censés représenter quoi au juste (voir image attachée) ?

I think it's because I use HPUserEdit and drag and drop from it and then, I change the special chars by hand :/

[Image: fske.png]
(05-01-2019 06:17 PM)Albert Chan Wrote: [ -> ]I am not familiar with RPL code, but is there a loop construct where you test the condition first ? That way, the 0 special case need not be added.

Thanks, Albert, for your valuable comments. Sure there is! Here we go:

\<<
  1 SWAP          @ Initialize result.
  WHILE
    DUP 2 \>=     @ Repeat while factor >= 2.
  REPEAT
    DUP ROT *     @ Multiply product with current factor.
    SWAP 2 -      @ Decrement factor by 2.
  END DROP        @ Drop remaining factor (by now invalid, i.e. < 2).
\>>

Saves 9 bytes, but negative arguments will return 1, too, whereas my former implementation left them unchanged.
(05-01-2019 07:00 PM)John Keith Wrote: [ -> ]Not necessarily. Smile See this thread.

Here is a version for the HP-48 or HP 50:

Code:

\<< \pi EVAL DUP2 * COS 1. - .25 * SWAP .5 * SWAP ^ SWAP .5 * DUP 2. SWAP ^ SWAP 1. + GAMMA * *
\>>

It will handle non-integer arguments, and complex numbers if complex mode is enabled. Negative even integers return results near negative infinity.

See the graphs here.

Hi John,
Thank you for this information. I needed something to understand everything reasonably. The connection between the gamma function and the double faculty I did not know. Mathematics and RPL, two really incredible things ....

many Greetings
Joerg
(05-01-2019 07:01 PM)Giuseppe Donnini Wrote: [ -> ]I cannot confirm this behavior. If n!! is too large, I get MAXR, i.e. 9.99999999999E499.
You must be using ZINTs, then. As a matter of fact, I wrote the code with an HP-48 in mind, so it is primarily intended for floats.

Again, this must be a ZINT-specific message.

Sorry Giuseppe, you are right. The answer will be given after about 21 seconds for number 9999. I was just confused, the other program versions end immediately.

Kind regards,
Joerg
This post gave me an opportunity to develop an HP-67 version. And here we go!



*LBL A:
001: 31 25 11  f LBL A
002:   33 00   STO 0
003:      84   R/S

*LBL B:
004: 31 25 12  f LBL B
005:     02    2
006:     81    ÷
007:     41    ENTER
008:   31 83  f INT
009:   32 51  g x=y?
010:   22 02    GTO 2

LBL 1:
011: 31 25 01  f LBL 1
012:   33 01    STO 1
013:     02     2
014:     71     ×
015:     01     1
016:     61     +
017:   35 81  h  N!
018:     02     2
019:   34 01    RCL 1
020:   35 63  h  yˣ
021:   35 82  h LST x
022:   35 81  h  N!
023:     71     ×
024:     81     ÷
025:   35 22  h  RTN

LBL 2:
026: 31 25 02  f LBL 2
027:     02     2
028:   35 52  h x≷y
029:   35 63  h  yˣ
030:   35 82  h LST x
031:   35 81  h  N!
032:     71     ×
033:   35 22  h  RTN

All you have to do…

Enter x Press A and then B
x!! is in the display.

You're welcome.
Hello

Just for the sake of using SEQ and ΠLIST:

Code:
<< → n
    << n 2 MOD
       { :: n :: n                 @ if n is even: 
         1 n 2 SEQ }                @ a list from 1 to n is created. 
       { :: n 1 − :: n             @ if n is odd:
         0 n 2 SEQ                  @ a list from '-1' to 'n − 1' is created.
         1 ADD                       @ adds 1 to each element of the list
         DUP                           @ duplicates the list (now from 0 to n)
         NOT                             @ new list: 0 becomes 1 ; others values become 0 
         ADD } IFTE                        @ adds the 2 lists, element by element
       ΠLIST                       @ returns the product of the elements of the list = the double factorial of n
    >>
>>

To avoid the {1} bug with ΠLIST, ΠLIST can be replaced by LPROD (from ListExt library).

:: are used to replace ' ' which can not be used in a list.

Or more conventionally, ' ' within chevrons:

<< 'n' … >> << 'n − 1' … >> IFTE
(01-27-2024 05:57 AM)FLISZT Wrote: [ -> ]...(from ListExt library).

If you've got ListExt installed, there's a fairly simple way to do this if we are only interested in non-negative integers (ie. no gamma function required):
Code:
\<<
  { 2 - } { 1 > } LWHL
  LPROD
\>>

The first line simply builds the list of multiplicands, given n as the argument on stack level 1:
0 yields { }
1 yields { }
2 yields { 2 }
3 yields { 3 }
4 yields { 4 2 }
5 yields { 5 3 }
6 yields { 6 4 2 }
7 yields { 7 5 3 }
...

By definition, LPROD provides 1 as the result when given an empty list. That works nicely here. Also, using "1 >" as the check condition, the superfluous final "1" is skipped for multiplication with odd n.

Note that this is just the core of the function; if input validation is needed, that should be added as required.
Pages: 1 2
Reference URL's