previous Back to the Homepage of Factorial Algorithms. You might also whish to visit our FactorialCalculator page and see how today an ultrafast calculator works. next


hp41 RETRO-FACTORIAL

Once upon a time Hewlett-Packard build RPN-Calculators, like the legendary HP-41C.

Even today you can put an HP-41 on a Windows desktop and use it like in the good old days. The (virtual) Hardware by Warren Furlow and Software. And here you can find many other calculator emulators.

In 1977 Peter Henrici decscribed in Computational Analysis with the HP-25 Pocket Calculator a simple routine for calculating the Gamma function Gamma(x) = (x-1)!, x > 0, x real. Henrici used the Stirling formula tuned by the QD-algorithm (a Padé-approximation). With minor modifications this becomes the following code for the HP-41C.

1STO 02+*1 / xEXP
21GTO 061 / x12
3STO 01LBL 13RCL 02+PI
47STO 0241 / x*
5RCL 026*CHSRCL 02
6LBL 06*+RCL 02/
7x > y1 / xRCL 02LNSQRT
8GTO 13RCL 02*+*
9STO * 013RCL 02RCL 01
1015**/

On my page comparing different approximation formulas for the factorial function I opted for Stieltjes' formula. So let us compare! The règles de jeu are: no more then 50 `HP-steps´ and only single-digit numbers. The formula used is discussed here. The HP-41C version is this sequence of keystrokes (the blue keystrokes are common to both routines):

1STO 02+*LNRCL 02
21GTO 061 / x1/
3STO 01LBL 13RCL 02SQRT
47STO 026RCL 02*
5RCL 024**RCL 01
6LBL 06*++/
7x > y1 / x2EXP 
8GTO 13RCL 02*2 
9STO * 01+1 / xPI 
1015RCL 02* 

Maybe a creative reader joins this friendly competition? (Please send me your solution.) Alternatively, you can use a very modern approach which combines object-oriented and functional programming (note the use of lambda expressions below) and use the cute little program below (written in C#), which simulates this evaluation of the factorial on the HP-41 calculator.

   1:  using System;
   2:  using System.Linq;
   3:  using System.Collections.Generic;
   4:   
   5:  namespace RPNCalculator
   6:  {
   7:      public delegate void VoidFunc<T>(T arg);
   8:   
   9:      public static class FPSwitch
  10:      {
  11:          public static void Switch<T>
  12:          (
  13:              this IEnumerable<T> source,
  14:              Func<T, int> expr,
  15:              params VoidFunc<T>[] funcArray
  16:          )
  17:          {
  18:              foreach (T item in source)
  19:              {
  20:                  funcArray[expr(item)](item);
  21:              }
  22:          }
  23:      }
  24:   
  25:      public enum TokenType
  26:      {
  27:          Operand, Enter, Exchange,
  28:          Add, Subtract, Multiply, Divide, Invert,
  29:          Exp, Ln, Sqrt, PI
  30:      };
  31:   
  32:      public class Token
  33:      {
  34:          public TokenType Type;
  35:          public double Operand;
  36:   
  37:          public Token(TokenType type)
  38:          {
  39:              this.Type = type;
  40:          }
  41:   
  42:          public Token(double op)
  43:          {
  44:              this.Type = TokenType.Operand;
  45:              this.Operand = op;
  46:          }
  47:   
  48:          // The RPN Token Processor
  49:          public static double Processor(Token[] prog)
  50:          {
  51:              Stack<double> st = new Stack<double>();
  52:   
  53:              prog.Switch
  54:              (
  55:                  i => (int)i.Type,
  56:                  s => st.Push(s.Operand),
  57:                  s => st.Push(st.Peek()),
  58:                  s =>
  59:                  {
  60:                      double a = st.Pop(), b = st.Pop();
  61:                      st.Push(a); st.Push(b);
  62:                  },
  63:                  s => st.Push(st.Pop() + st.Pop()),
  64:                  s => st.Push(-st.Pop() + st.Pop()),
  65:                  s => st.Push(st.Pop() * st.Pop()),
  66:                  s => st.Push(1 / st.Pop() * st.Pop()),
  67:                  s => st.Push(1 / st.Pop()),
  68:                  s => st.Push(Math.Exp(st.Pop())),
  69:                  s => st.Push(Math.Log(st.Pop())),
  70:                  s => st.Push(Math.Sqrt(st.Pop())),
  71:                  s => st.Push(Math.PI)
  72:              );
  73:   
  74:              return st.Pop();
  75:          }
  76:      }

That's all! Now we can start writing programs for this HP-41 simulator. Here is the RETRO-FACTORIAL program:

  77:      class Factorial
  78:      {
  79:          static double factorial(double x)
  80:          {
  81:              Token[] program = 
  82:              {
  83:                  new Token(x),
  84:                  new Token(1),
  85:                  new Token(TokenType.Add),
  86:                  new Token(TokenType.Enter),
  87:                  new Token(TokenType.Enter),
  88:                  new Token(TokenType.Enter),
  89:                  new Token(TokenType.Enter),
  90:                  new Token(4),
  91:                  new Token(TokenType.Multiply),
  92:                  new Token(TokenType.Invert),
  93:                  new Token(TokenType.Add),
  94:                  new Token(5),
  95:                  new Token(TokenType.Multiply),
  96:                  new Token(TokenType.Invert),
  97:                  new Token(TokenType.Exchange),
  98:                  new Token(6),
  99:                  new Token(TokenType.Multiply),
 100:                  new Token(TokenType.Add),
 101:                  new Token(2),
 102:                  new Token(TokenType.Multiply),
 103:                  new Token(TokenType.Invert),
 104:                  new Token(TokenType.Exchange),
 105:                  new Token(TokenType.Enter),
 106:                  new Token(TokenType.Ln),
 107:                  new Token(1),
 108:                  new Token(TokenType.Subtract),
 109:                  new Token(TokenType.Multiply),
 110:                  new Token(TokenType.Add),
 111:                  new Token(TokenType.Exp),
 112:                  new Token(TokenType.Exchange),
 113:                  new Token(TokenType.PI),
 114:                  new Token(2),
 115:                  new Token(TokenType.Multiply),
 116:                  new Token(TokenType.Exchange),
 117:                  new Token(TokenType.Divide),
 118:                  new Token(TokenType.Sqrt),
 119:                  new Token(TokenType.Multiply)
 120:              };
 121:   
 122:              return Token.Processor(program);
 123:          }
 124:   
 125:          static void Main()
 126:          {
 127:              for (int i = 0; i < 12; i++)
 128:              {
 129:                  Console.WriteLine(i + " -> " + factorial(i));
 130:              }
 131:          }
 132:      }
 133:  }

previous Back to the Homepage of Factorial Algorithms. You might also whish to visit our FactorialCalculator page and see how today an ultrafast calculator works. next