From - Wed Oct 01 07:22:42 1997 Path: news.mitre.org!blanket.mitre.org!agate!nntp.info.ucla.edu!208.134.241.18!newsfeed.internetmci.com!199.60.229.5!newsfeed.direct.ca!Supernews60!supernews.com!uunet!in5.uu.net!news-feed.newscorp.com!news.delphi.com!news From: alan Newsgroups: comp.os.vms Subject: Re: VAX Fortran...limits.... Date: Wed, 1 Oct 97 00:03:13 -0500 Organization: Delphi (info@delphi.com email, 800-695-4005 voice) Lines: 52 Message-ID: References: <009BAFB3.1FA61393.17@earth.oscs.montana.edu> <342F72CF.41C6@spammenot.wisconsin.cern.ch> <30SEP199702351329@gerg.tamu.edu> NNTP-Posting-Host: 199.93.4.4 X-To: Carl Perkins Carl Perkins writes: >}You might want to look into Stirling's approximation for the factorial >}(actually, for the Gamma-function). My trusted little old "Particle >}Properties >}Data Booklet" from 1980(*) says: >} >}Gamma(x) == (x-1)! = (approx.) sqrt(2 pi) * exp(-x) * x**(x-1/2) * ( >}1+1/(12x)+...) >} >}It also claims the given approximation is good to better than +-0.1% for >}all x >= 1/2, so a few more terms in the expansion should get you to >}your >}required precision. Double precision arithmetic might then be sufficient -- >Uh, no. It would be quicker but the problem is that the factorial would still >overflow the exponent of any of the floating point formats available nativly >(and going to a mere double precision would make that overflow happen sooner). > >What he needs is to use unlimited precision math. Off hand, the only such >stuff that he'd already have is the STR$MUL function that comes in the STR$ >RTL. -- YES!, you can use Stirling's formula, but you need to do some algebra first! Take the log of both sides. This yields (these are NOT Fortan statements, generally!!!) (You'd at least need to change a few of the 10's to 10.'s) -- log n! = log(e) [(n+1/2)(ln(n)) + ln(sqrt(2pi)) - n + 1/12n - ...] -- Then do a 10**(log n!) separately for the integer and fractional parts. -- Example: n = 100 What is 100!? -- log (100!) = 157.9696... ! using the log n! formula above (well, without the n/12 part, include it for more accuracy!) -- n! = 10**(157.9696...) = 10**(.9696) * 10**157 -- ==> 100! ~= 9.325 E+157 -- My calculator does 100! and comes up with 9.333 E+157. Thus, the formula *in this case* is accurate to better than 0.1% without the n/12 term! -- There is no need for raising n**n. However, of course, you still may end up with exponents too large to put in a Fortran variable for very large n. -- Disclaimer: JMHO. -- Alan E. Feldman alan48 delphi.com &-)