Well, I wouldn't have used such colorful language, but I share
Prof. Fateman's skepticism on the utility of lgamma for computing
factorials.
I tried the expression as given, (see attachment, if you want to try
yourself or point out errors) and it begins to diverge from the correct
value at n=17 (which is a bit more than 48 bits for the integer).
Interestingly, if I switch to gamma(), instead of exp(lgamma()), it
doesn't diverge until n=23, which corresponds to a 75 bit number.
(Presumably there are a lot of zeros at the bottom of a factorial.)
So, that *would* work if you were in C and restricted to 64 bit long
longs. But, as Fateman points out, with that restriction you might as
well keep a short table.
Of course, in Python we expect thousands of bits to work
transparently, so I would definitely avoid converting from lgamma to
an int.
BTW, I don't understand the "far superior" implementation of fact
later in the posting. When I try it, I just get regular linear
recursion. Presumably there is a typo in my translation or in the
original posting.
On Saturday, May 3, 2014 4:35:41 PM UTC-7, Richard Fateman wrote:
>
>
>
> On Tuesday, April 1, 2014 2:58:34 PM UTC-7, Ondřej Čertík wrote:
>>
>> On Mon, Mar 31, 2014 at 9:14 PM, Matthew Rocklin <[email protected]>
>> wrote:
>> > http://www.evanmiller.org/mathematical-hacker.html
>> >
>> > I reference that blog post pretty often. I fully intend to reference
>> it
>> > again in my talk (if it is accepted).
>>
>
> That's unfortunate. I was totally unacquainted with Evan Miller. I read
> this particular
> blog entry. I did not follow any of the links, but I think that the essay
> is a piece of crap.
> It is uninformed, shallow, and, I daresay, just wrong. I hope you
> re-think referring to it.
>
>
>
>> >
>> > The interesting thing about the Factorial / Gamma / loggamma example is
>> that
>> > to find the solution you need to find someone who knows both that n! =
>> > Gamma(n+ 1) and who knows that a loggamma routine is commonly found in
>> lower
>> > level languages. Those bits of information are usually held by
>> different
>> > experts. Ondrej said "Of course, that's obvious" when I first reposted
>> the
>> > article on G+.
>>
>
> Well, let me quote from the article..
> "
> *if it is ever necessary to compute a factorial, a programmer should be
> taught to take advantage of the system's log-gamma function, as the
> following C code does: *
>
> *
> long int fac(unsigned long int n) {
> return lround(exp(lgamma(n+1)));
> }
> *
>
> * Again, no recursion is required as long as one...*
>
> Here are a few things that are wrong. Rounding an approximate floating
> point
> number to a integer may not get you the exact integer factorial you want.
> The program will certainly not work for arbitrary unsigned long int n
> because
> the exp() will overflow. Say certainly by n=175.
> If you want a floating-point approximation to any factorial that can be
> represented as an IEEE double, you could do it by checking if N<175 and
> then index into an array of length 175. [I'd check that 175 more carefully
> if I were really doing this. And the number would change a little for
> 64-bit ints.
>
> If you want exact integer factorials on a 32 bit computer, you only need
> an array of 12, since
> 13! > 2^32.
>
> So the program is stupid. Why? Because the author is mathematically
> and computationally ignorant. (And seemingly quite full of himself.)
>
> Say you want to solve the much more interesting problem of computing
> factorials
> of any size integer. Then you might see
> http://www.luschny.de/math/factorial/FastFactorialFunctions.htm
> which tells how it might be done better.
>
> Here's a far superior factorial, in lisp, and yes it is recursive. I
> suppose you
> could write it iteratively, but why bother. A Lisp compiler will do that
> for you.
>
> (defun fact(n)(f n 1))
>
> (defun f (n m) (if (<= n 1) 1 (* n (f (- n m) m))))
>
> It is used in some computer algebra systems.
>
>
> The snide characterizations that lisp programmers know more or less
> mathematics
> or that what they know is more or less relevant is just junk.
>
> Of the 3 people quoted -- Raymond, Graham, Yegge...
> I've never heard of Yegge; I think that Graham is smart, and Raymond just
> says things
> to be, uh, smart-ass.
>
>
>
>
>
>
>>
>> That's funny, I forgot that I said that and just had the same reaction.
>>
>> However, we are still missing rational function approximation in SymPy
>> or mpmath.
>> That's what is used to implement things like log_gamma or erf (error
>> function) in
>> Fortran or C. So I have just created an issue for it:
>>
>> https://github.com/sympy/sympy/issues/7359
>>
>> and I spent time explaining exactly how it works in it and giving
>> examples,
>> including for example the implementation of erf() in gfortran.
>>
>> Aaron, you were asking about examples of numerical cancellation. I
>> worked out one in the
>> issue as well.
>>
>> Ondrej
>>
>> >
>> > You're right that this is similar to my last talk. The last one though
>> was
>> > mostly about an application (numerical linear algebra). I actually
>> want to
>> > talk a bit more about the philosophy and some of the more abstract
>> tools
>> > that people might actually use. Your first impression is a valuable
>> one
>> > though, I should go through my last talk and make sure that I'm not
>> > repeating too much that shouldn't be repeated.
>> >
>> >
>> > On Mon, Mar 31, 2014 at 6:03 PM, Aaron Meurer <[email protected]>
>> wrote:
>> >>
>> >> That's a good point. One of the nicest things about symbolics, when
>> you
>> >> can get it, is that it can make things drastically more efficient by
>> doing
>> >> mathematical simplifications. Evaluating integrals symbolically is a
>> nice
>> >> example of this (especially for SymPy, which has some pretty nice
>> algorithms
>> >> to compute definite integrals).
>> >>
>> >> I'm reminded of a popular blog post (I can't find a link right now)
>> about
>> >> how know math is important for programmers. It has the example of how
>> all
>> >> these programming languages show how they they compute factorial, and
>> how
>> >> tail recursion can make it linear or whatever, but the actual best way
>> to
>> >> compute it is to use loggamma, which gives the answer in constant
>> time.
>> >>
>> >> Aaron Meurer
>> >>
>> >>
>> >> On Mon, Mar 31, 2014 at 7:51 PM, Tim Lahey <[email protected]> wrote:
>> >>>
>> >>>
>> >>>
>> >>> On 31 Mar 2014, at 20:29, Aaron Meurer wrote:
>> >>>
>> >>>> On Mon, Mar 31, 2014 at 11:32 AM, Matthew Rocklin
>> >>>> <[email protected]>wrote:
>> >>>>
>> >>>>> I like that you emphasized the utility for numerics, I think that
>> this
>> >>>>> is
>> >>>>> likely to be a selling point for the SciPy crowd.
>> >>>>>
>> >>>>
>> >>>> Yes, this was very intentional. I may need some help gathering up
>> some
>> >>>> nice
>> >>>> motivating examples if this is accepted.
>> >>>
>> >>>
>> >>> One motivating example for me is the integration of products of
>> functions
>> >>> over areas and volumes. For finite elements, you'll get products of
>> pairs of
>> >>> trial functions (usually polynomials). It's even more useful for
>> products of
>> >>> trig functions. Performing the integration of any of theses is easy
>> enough
>> >>> with numerical integration, but it's much more efficient to calculate
>> the
>> >>> integrals symbolically and then perform the evaluation for each
>> element.
>> >>>
>> >>> Cheers,
>> >>>
>> >>> Tim.
>> >>>
>> >>>
>> >>> --
>> >>> You received this message because you are subscribed to the Google
>> Groups
>> >>> "sympy" group.
>> >>> To unsubscribe from this group and stop receiving emails from it,
>> send an
>> >>> email to [email protected].
>> >>> To post to this group, send email to [email protected].
>> >>> Visit this group at http://groups.google.com/group/sympy.
>> >>> To view this discussion on the web visit
>> >>>
>> https://groups.google.com/d/msgid/sympy/AE5C61B5-8C69-49AE-BAFF-85E5AC924F8F%40gmail.com.
>>
>>
>> >>>
>> >>> For more options, visit https://groups.google.com/d/optout.
>> >>
>> >>
>> >> --
>> >> You received this message because you are subscribed to the Google
>> Groups
>> >> "sympy" group.
>> >> To unsubscribe from this group and stop receiving emails from it, send
>> an
>> >> email to [email protected].
>> >> To post to this group, send email to [email protected].
>> >> Visit this group at http://groups.google.com/group/sympy.
>> >> To view this discussion on the web visit
>> >>
>> https://groups.google.com/d/msgid/sympy/CAKgW%3D6K8F9weBP-7yjpzrN9hQ5FR54TK_QdCPqbzBdwo20zKdg%40mail.gmail.com.
>>
>>
>> >>
>> >> For more options, visit https://groups.google.com/d/optout.
>> >
>> >
>> > --
>> > You received this message because you are subscribed to the Google
>> Groups
>> > "sympy" group.
>> > To unsubscribe from this group and stop receiving emails from it, send
>> an
>> > email to [email protected].
>> > To post to this group, send email to [email protected].
>> > Visit this group at http://groups.google.com/group/sympy.
>> > To view this discussion on the web visit
>> >
>> https://groups.google.com/d/msgid/sympy/CAJ8oX-HZABscP43PDuoAeWhcpBYW2w39G8qPx0%2BZCbKcYbJdJg%40mail.gmail.com.
>>
>>
>> >
>> > For more options, visit https://groups.google.com/d/optout.
>>
>
--
You received this message because you are subscribed to the Google Groups
"sympy" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To post to this group, send email to [email protected].
Visit this group at http://groups.google.com/group/sympy.
To view this discussion on the web visit
https://groups.google.com/d/msgid/sympy/f999c5c5-b99b-46bf-a40b-5d1eb5ff11e6%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
import math
def fact(n):
"""Factorial algorithm from R. Fateman sympy group post. Gives
correct values, but seems to result in ordinary linear
behavior."""
return f(n, 1)
def f(n, m):
"""Semantics?"""
if n<=1:
return 1
return n*f(n-m, m)
def miller_fac(n):
"""Factorial approximation, based on Evan Miller blog post."""
return int(math.exp(math.lgamma(n+1)) + 0.5)
def miller2_fac(n):
"""Slight modification - use gamma instead of lgamma."""
return int(math.gamma(n+1) + 0.5)
def test_factorial_implementations(fimpl, f_max):
for n in xrange(f_max+1):
base = math.factorial(n) # this should be exact
f2 = fimpl(n)
if base != f2:
print "Discrepancy found for", n, base, f2
print "lg of factorial", math.log(base, 2)
if __name__ == "__main__":
print "Original expression, using exp(lgamma)"
test_factorial_implementations(miller_fac, 20)
print "Using Gamma"
test_factorial_implementations(miller2_fac, 25)
