Factorial function works in Python, returns 0 for Julia

Question

I define a factorial function as follows in Python:

def fact(n):
    if n == 1:
        return n
    else:
        return n * fact(n-1)

print(fact(100))

and as follows in Julia:

function fact(n)
    if n == 1
        n
    else
        n * fact(n-1)
    end
end

println(fact(100))

The python program returns a very large number for the evaluation of 100 (as expected). Julia returns 0. With a smaller number (like 10) they both work.

I have two questions:

1. Why does Python handle this OK and Julia not.
2. Why doesn't Julia throw an error and just print 0 instead?

Answer 1

Julia has separate fixed-size integer types, plus a BigInt type. The default type is Int64, which is of course 64 bits.

Since 100! takes about 526 bits, it obviously overflows an Int64.

You can solve this problem by just doing fact(BigInt(100)) (assuming you've required it), or of course you can do the conversion in the fact function.

---

Python used to be the same, once upon a time. It had separate types int, which was 16 or 32 or 64 bits depending on your machine, and long, which was arbitrary-length. If you ran your program on Python 1.5, it would either wrap around just like Julia, or raise an exception. The solution would be to call fact(100L), or to do the conversion to long inside the fact function.

However, at some point in the 2.x series, Python tied the two types together, so any int that overflows automatically becomes a long. And then, in 3.0, it merged the two types entirely, so there is no separate long anymore.

---

So, why does Julia just overflow instead of raising an error?

The FAQ actually explains Why does Julia use native machine integer arithmetic. Which includes the wraparound behavior on overflow.

---

By "native machine arithmetic", people generally mean "what C does on almost all 2s-complement machines". Especially in languages like Julia and Python that were originally built on top of C, and stuck pretty close to the metal. In the case of Julia, this is not just a "default", but an intentional choice.

In C (at least as it was at the time), it's actually up to the implementation what happens if you overflow a signed integer type like int64… but on almost any platform that natively uses 2's complement arithmetic (which is almost any platform you'll see today), the exact same thing happens: it just truncates everything above the top 64 bits, meaning you wrap around from positive to negative. In fact, unsigned integer types are required to work this way in C. (C, meanwhile, works this way because that's how most CPUs work.)

In C (unlike most CPUs' machine languages), there is no way to detect that you've gotten an overflow after the fact. So, if you want to raise an OverflowError, you have to write some logic that detects that the multiplication will overflow before doing it. And you have to run that logic on every single multiplication. You may be able to optimize this for some platforms by writing inline assembly code. Or you can cast to a larger type, but (a) that tends to make your code slower, and (b) it doesn't work if you're already using the largest type (which int64 is on many platforms today).

In Python, making each multiplication up to 4x slower (usually less, but it can be that high) is no big deal, because Python spends more time fetching the bytecode and unboxing the integer objects than multiplying anyway. But Julia is meant to be faster than that.

As John Myles White explains in Computers are Machines:

In many ways, Julia sets itself apart from other new languages by its attempt to recover some of the power that was lost in the transition from C to languages like Python. But the transition comes with a substantial learning curve.

---

But there's another reason for this: overflowing signed arithmetic is actual useful in many cases. Not nearly as many as overflowing unsigned arithmetic (which is why C has defined unsigned arithmetic to work that way since before the first ANSI spec), but there are use cases.

And, even though you probably want type conversions more often than you want rollover, it is a lot easier to do the type conversions manually than the rollover. If you've ever done it in Python, picking the operand for % and getting the signs right is certainly easy to get wrong; casting to BigInt is pretty hard to screw up.

---

And finally, in a strongly-typed language, like both Python and Julia, type stability is important. One of the reasons Python 3 exists was that the old str type magically converting to unicode caused problems. It's far less common for your int type magically converting to long to cause problems, but it can happen (e.g., when you're grabbing a value off the wire, or via a C API, and expect to write the result out in the same format). Python's dev team argued over this when doing the int/long unification, quoting "practicality beats purity" and various other bits of the Zen, and ultimately deciding that the old behavior caused more problems than the new behavior would. Julia's designed made the opposite decision.

Answer 2

Nobody answers why the result in Julia is 0.

Julia does not check integer multiplication for overflow and thus the multiplication for 64 bit integers is performed mod 2^63. See this FAQ entry

when you write out the multiplication for factorial you get

1*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16*17*18*19*20*21*22*23*24*25*26*27*28*29*30*31*32*33*34*35*36*37*38*39*40*41*42*43*44*45*46*47*48*49*50*51*52*53*54*55*56*57*58*59*60*61*62*63*64*65*66*67*68*69*70*71*72*73*74*75*76*77*78*79*80*81*82*83*84*85*86*87*88*89*90*91*92*93*94*95*96*97*98*99*100

This can also be written as prime factors

2^97 * 3^48 * 5^24 * 7^16 * 11^9 * 13^7 * 17^5 * 19^5 * 23^4 * 29^3 * 31^3 * 37^2 * 41^2 * 43^2 * 47^2 * 53^1 * 59^1 * 61^1 * 67^1 * 71^1 * 73^1 * 79^1 * 83^1 * 89^1 * 97^1

If you look at the exponent of 2 you get 97. Modular arithmetic gives that you can do the mod function at any step of the calculation, and it will not affect the result. 2^97 mod 2^63 == 0 which multiplied with the rest of the chain is also 0.

UPDATE: I was of course too lazy to do this calculation on paper.

d = Dict{Int,Int}()
for i=1:100
   for (k,v) in factor(i)
       d[k] = get(d,k,0) + v
   end
end
for k in sort(collect(keys(d)))
    print("$k^$(d[k])*")
end

Julia has a very convenient factor() function in its standard library.