jurjen@cwi.nl (Jurjen NE Bos) (02/25/91)
I am sorry for the guy who just posted his complex factorial program, but my
program was almost ready at that point, so I post it anyway.
The end of this posting contains a library.
- load the string in your HP48
- convert it using ASC->
- store it in a port by typing <port> STO (for example, 0 STO)
- turn your calculator off and on (you'll lose your stack!)
- enter 1617 ATTACH in the directory where you intend to use it
Now you have remapped the FACT function to a special version that also takes
complex arguments.
Features:
- FACT of a name or algebraic gives ! instead. To make a function that works
with both algebraics and complex numbers, use the user function:
FAC @warning: typed in by hand, download with caution
%%HP: T(3);
\<< \-> z
\<<
IF { 6 7 9 } z TYPE POS
THEN z 1 \->LIST 'FAC' APPLY
ELSE z FACT
END
\>>
\>>
This program should be called FAC.
This should be downloaded or typed in AFTER you installed the library.
Sorry for the inconvenience; as far as I know, it is impossible to
incorporate library functions into an algebraic.
- The library is only 495.5 bytes.
- For real arguments, the program is identical to !.
- This program is extremely accurate, because it uses Long Real and Long
Complex internally. The accuracy is more than 11 digits (most of the time
all digits are correct), even around the poles.
- The program is fast, because it uses internal functions only. Of course,
a complex factorial is hard to compute, so it still takes 540 milliseconds
to compute (5 3) FACT. As far as I know, the maximum computation time is
875 milliseconds. (For real arguments, it can go up to 2180 milliseconds!)
- The algorithm used is:
For real numbers, ! is called.
For numbers with real part >-6: the same algorithm as HP uses. I will
not tell you what it is, because it is probably patented. I will only
tell you that I did not see it before, and that it is fast and
accurate. Of course, the algorithm had to be adapted for complex
numbers.
For number with real part <-6, I use the formula:
pi * x
FACT(X) = --------------------
sin(pi*x) * FACT(-X)
The sine is compute using degrees for the real part of x to preserve
accuracy around the poles.
- Note that -250.000000001 FACT is slower (and less accurate) than
(-250.000000001 0) FACT.
- Weird errors or results can occur if an error occurs during the
computation. I have not see this happen, and I only expect it to happen for
extremely large arguments (absolute value > 10000).
Have fun!
Here it is:
%%HP: T(1)A(R)F(,);
"04B20AD300C034F6D607C656870264143445C01564100000000B700000000E4A
20C6000000000000000000640000000000000000000000000000000000000000
00000000000000000000000210004064143445000D0000E4A204100061000551
007F1008156000D9D20ECE819FF304C0B230040D9D20BC915A8C50D9D2088130
55920000000000000000000069F18A23991629E20156100FF1B22C230F49A299
FA288130B08A27F81629E20156200CA130D9C1529E20156200322300B5A2AEC8
1592300E4A2D9C1588130592302A170E0E304D226EE170D9D20592300E4A2D9C
15592302C2302ED1559230B21305E170592308523031F15B2130CBD5070A15B2
1303504057155B2130D9D208813019B1529E201562002C230854A2BAD1588130
CBD50B9F06BBF0655920200000000000000000810A99A28813071CA2EF1167CD
A2A99A23223086CA25923059DA2A99A2A8C50592302ED1531F15B2130D9D200E
4A2D9C1588130CA1302ED15CD2B2BAD1588130003B2D9C15F13B232230B3F153
43B2D9C152C230E3C15C63B232230B3F15093B2D9C152C230E3C159B3B232230
B3F15DD3B2D9C15E3C1531F152C23001D15014B2D9C152C230CBD5051225C7DA
2A8C5059230C75A283D152ED15E3C15CBD5032230C1BA22C23087CA22C230DEF
2672CA259230A99A2A8C50B2130FCB8F124"