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

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.

Answer 2

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 3

I guess the answer to this is to use BigInt:

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

println(fact(BigInt(100))) 

Which gives result:

93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000

Tested on: http://forio.com/julia/repl/

As stated in some other answers, Python implicitly converts int(s) that exceed the maximum size to bigint(s) for you and so you get the result you expect rather than failing silently.

Julia on the other hand seems to be more explicit about this and favors performance over "expected behavior". Julia is a Dynamic Language with OPtional Type Annotations and Inference.

Answer 4

Python automatically uses a BigInt that can hold arbitrarily large numbers. In Julia, you have to do it yourself. I would think it is corrected like this

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

println(fact(BigInt(100))) 

Answer 5

One of the reasons Julia is fast is that they avoid features that would compromise performance. This is one of them. In Python, the interpreted constantly checks to see whether it should automatically switch to the BigInt library. That constant checking comes at a cost.

Here is a function that does what you want:

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

println( fact(100) )
println( fact(0) )

I took the liberty of correcting a bug in your program: Zero factorial is defined, and is 1. Incidentally, you could write your function like this:

function !(n)
    if n == 0
        1
    else
        big(n) * !(n-1)
    end
end

println( !(100) )
println( !(0) )

I would not personally do this, because "foo!" functions are conventionally used for functions that modify the arguments. But I wanted to offer the option. Lastly, I can't resist offering the one-line alternative:

fact(n) = n == 0 ? 1 : big(n) * !(n-1)

println( fact(100) )
println( fact(0) )

Answer 6

By the way, I think an idiomatic way of doing this in Julia would be to make use of the type system to compile different versions for different types:

fact(n) = n <= zero(n) ? one(n) : n*fact(n-one(n)) 
# one(n) gives you a one, as it were, of the same type as n

Then, a different version of that function is compiled and called depending on the type of the input, and the user would have to decide which type to use, and thus which version of the function to call:

julia> fact(10)
3628800

julia> fact(100)
0

julia> fact(BigInt(100))
93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000

The overhead of BigInt vs machine (Int64) arithmetic, well described by @abarnert, can be seen when we look at the (LLVM) compiled versions of the Int64 vs the BigInt fact():

julia> code_llvm(fact,(Int64,))

define i64 @"julia_fact;23421"(i64) {
top:
  %1 = icmp sgt i64 %0, 0, !dbg !10800
  br i1 %1, label %L, label %if, !dbg !10800

if:                                               ; preds = %top
  ret i64 1, !dbg !10800

L:                                                ; preds = %top
  %2 = add i64 %0, -1, !dbg !10800
  %3 = call i64 @"julia_fact;23398"(i64 %2), !dbg !10800
  %4 = mul i64 %3, %0, !dbg !10800
  ret i64 %4, !dbg !10800
}



julia> code_llvm(fact,(BigInt,))

define %jl_value_t* @"julia_fact;23422"(%jl_value_t*, %jl_value_t**, i32) {
top:
  %3 = alloca [6 x %jl_value_t*], align 8
  %.sub = getelementptr inbounds [6 x %jl_value_t*]* %3, i64 0, i64 0
  %4 = getelementptr [6 x %jl_value_t*]* %3, i64 0, i64 2, !dbg !10803
  store %jl_value_t* inttoptr (i64 8 to %jl_value_t*), %jl_value_t** %.sub, align 8
  %5 = load %jl_value_t*** @jl_pgcstack, align 8, !dbg !10803
  %6 = getelementptr [6 x %jl_value_t*]* %3, i64 0, i64 1, !dbg !10803
  %.c = bitcast %jl_value_t** %5 to %jl_value_t*, !dbg !10803
  store %jl_value_t* %.c, %jl_value_t** %6, align 8, !dbg !10803
  store %jl_value_t** %.sub, %jl_value_t*** @jl_pgcstack, align 8, !dbg !10803
  store %jl_value_t* null, %jl_value_t** %4, align 8, !dbg !10803
  %7 = getelementptr [6 x %jl_value_t*]* %3, i64 0, i64 3
  store %jl_value_t* null, %jl_value_t** %7, align 8
  %8 = getelementptr [6 x %jl_value_t*]* %3, i64 0, i64 4
  store %jl_value_t* null, %jl_value_t** %8, align 8
  %9 = getelementptr [6 x %jl_value_t*]* %3, i64 0, i64 5
  store %jl_value_t* null, %jl_value_t** %9, align 8
  %10 = load %jl_value_t** %1, align 8, !dbg !10803
  %11 = call %jl_value_t* @julia_BigInt2(i64 0), !dbg !10804
  store %jl_value_t* %11, %jl_value_t** %4, align 8, !dbg !10804
  %12 = getelementptr inbounds %jl_value_t* %10, i64 1, i32 0, !dbg !10804
  %13 = getelementptr inbounds %jl_value_t* %11, i64 1, i32 0, !dbg !10804
  %14 = call i32 inttoptr (i64 4535902144 to i32 (%jl_value_t**, %jl_value_t**)*)(%jl_value_t** %12, %jl_value_t** %13), !dbg !10804
  %15 = icmp sgt i32 %14, 0, !dbg !10804
  br i1 %15, label %L, label %if, !dbg !10804

if:                                               ; preds = %top
  %16 = call %jl_value_t* @julia_BigInt2(i64 1), !dbg !10804
  %17 = load %jl_value_t** %6, align 8, !dbg !10804
  %18 = getelementptr inbounds %jl_value_t* %17, i64 0, i32 0, !dbg !10804
  store %jl_value_t** %18, %jl_value_t*** @jl_pgcstack, align 8, !dbg !10804
  ret %jl_value_t* %16, !dbg !10804

L:                                                ; preds = %top
  store %jl_value_t* %10, %jl_value_t** %7, align 8, !dbg !10804
  store %jl_value_t* %10, %jl_value_t** %8, align 8, !dbg !10804
  %19 = call %jl_value_t* @julia_BigInt2(i64 1), !dbg !10804
  store %jl_value_t* %19, %jl_value_t** %9, align 8, !dbg !10804
  %20 = call %jl_value_t* @"julia_-;23402"(%jl_value_t* inttoptr (i64 140544121125120 to %jl_value_t*), %jl_value_t** %8, i32 2), !dbg !10804
  store %jl_value_t* %20, %jl_value_t** %8, align 8, !dbg !10804
  %21 = call %jl_value_t* @"julia_fact;23400"(%jl_value_t* inttoptr (i64 140544559367232 to %jl_value_t*), %jl_value_t** %8, i32 1), !dbg !10804
  store %jl_value_t* %21, %jl_value_t** %8, align 8, !dbg !10804
  %22 = call %jl_value_t* @"julia_*;23401"(%jl_value_t* inttoptr (i64 140544121124768 to %jl_value_t*), %jl_value_t** %7, i32 2), !dbg !10804
  %23 = load %jl_value_t** %6, align 8, !dbg !10804
  %24 = getelementptr inbounds %jl_value_t* %23, i64 0, i32 0, !dbg !10804
  store %jl_value_t** %24, %jl_value_t*** @jl_pgcstack, align 8, !dbg !10804
  ret %jl_value_t* %22, !dbg !10804
}