12

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?
15
  • 1
    This might be a wild guess but 100! is a crazy large number, so you might overflow your Julia program, and it might return 0 as a result of your program crashing. Have you tried with a smaller number, like 5 and see what happens?
    – Dyrborg
    Jan 15, 2014 at 2:01
  • It works with a smaller number, I'll clarify that.
    – Mike Vella
    Jan 15, 2014 at 2:04
  • 5
    See JMW's post here.
    – DSM
    Jan 15, 2014 at 2:08
  • 1
    Based on the julia documentation it handles numbers up to a maximum of 128 bits (docs.julialang.org/en/latest/manual/…) - a number MUCH smaller than 100!. Try and see if 30! works but 35! doesn't. The overflow limits lies between there as far as I can see.
    – Dyrborg
    Jan 15, 2014 at 2:10
  • 2
    @JamesMills: Are you kidding? 100! is about 9E157, which takes 526 bits.
    – abarnert
    Jan 15, 2014 at 2:15

6 Answers 6

20

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.

4
  • There's one thing I still don't understand - why wasn't Julia designed to throw an error in this situation? Is it for performance optimization or an oversight?
    – Mike Vella
    Jan 15, 2014 at 2:19
  • @MikeVella: Probably because it's relying on the underlying C representation. Let me edit again to explain.
    – abarnert
    Jan 15, 2014 at 2:20
  • 7
    It's quite intentional. Partly for performance and partly because it's sometimes quite useful: docs.julialang.org/en/latest/manual/faq/…. Jan 15, 2014 at 2:24
  • @StefanKarpinski: Thanks, that link covers everything I wanted to cover and better (which makes sense, since it's written as part of the Julia docs, rather than a quickly answer from someone who actually prefers Python…). I'll edit the answer.
    – abarnert
    Jan 15, 2014 at 2:44
18

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.

3
  • Thanks for this, Can I ask how you work this out? It seems like quite a complicated thing to figure out.
    – Mike Vella
    Jan 15, 2014 at 23:03
  • 5
    I updated with code. The understanding of the problem is from taking math classes, and reading about cryptography and modular arithmetic.
    – ivarne
    Jan 16, 2014 at 8:06
  • 4
    A much easier way to compute the factorization of 100! in Julia is to ask for it directly: factor(factorial(big(100))). Aug 4, 2014 at 4:21
10

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.

6
  • Modern Python only has bigints. (Although in the PyPy 2.1 implementation of Python 3.2, they do get optimized to int8, int32, or int64 in some cases, that still isn't visible to the user.)
    – abarnert
    Jan 15, 2014 at 2:33
  • Python version(s) what onwards use bigints everywhere? Jan 15, 2014 at 2:35
  • 3.0. Although IIRC the implementation was drastically changed in 3.1, and as usually 3.0 should be ignored. Meanwhile, 2.6 is where they made it harder to distinguish the two types (without making it impossible—the repr and type are obviously different, and isinstance doesn't have any magic), and I forget which version (2.3?) added the automatic conversion.
    – abarnert
    Jan 15, 2014 at 2:36
  • Right - confirmed. So Python 3+ uses a single type <type 'int'> which I presume covers any size of integer thrown at it. Jan 15, 2014 at 2:37
  • Well, you can't exactly throw Graham's number at it; eventually you'll get a MemoryError. But otherwise, yeah, that's the idea.
    – abarnert
    Jan 15, 2014 at 2:54
2

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))) 
1
  • syntax: malformed function argument (call BigInt n) Jan 15, 2014 at 2:29
0

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) )
1
  • I think you forgot to correct the function name in the last sample :)
    – xzczd
    Apr 29, 2014 at 7:49
0

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
}

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.