Typeof(n*n) in factorial

The implementation of factorial in Base is

function factorial(n::Integer)
    n < 0 && throw(DomainError(n, "`n` must be nonnegative."))
    f::typeof(n*n) = 1
    for i::typeof(n*n) = 2:n
        f *= i
    end
    return f
end

Why typeof(n*n) instead of typeof(n)? This maybe a dumb question but I really don’t get it.

There are numbers who’s type is not closed under multiplication. This makes sure you get something of the right type. In practice, this will almost always just be typof(n) though.

4 Likes

Can you please elaborate a little more? Example?

One example would be values with physical units, as provided by the Unitful.jl package.

The area of a square has not the same type as the length of its sides:

julia> using Unitful

julia> length = 1u"m"
1 m

julia> typeof(length)
Quantity{Int64,𝐋,Unitful.FreeUnits{(m,),𝐋,nothing}}

julia> area = length*length
1 m^2

julia> typeof(area)
Quantity{Int64,𝐋^2,Unitful.FreeUnits{(m^2,),𝐋^2,nothing}}

julia> typeof(area) == typeof(length)
false

edit: Hmm, maybe scratch that…since it is not applicable to the factorial method.

2 Likes

It’s fine. It still gave the example I was looking for.

Well, but now I don’t understand why it is necessary in the factorial case :slight_smile:

I really hope someone will answer.

If there is some strange implementation of Integer for which n*n is not closed under multiplication, why should (n*n)*n be? And why should 1 be of the same type as n*n?

1 Like

Working with result types is actually one of the most irksome parts of Julia. Even restricting to <:Number, the interfaces are not formalized. So while it may be strange, making a <:Number type that is not closed under multiplication does not formally violate any contract.

While in this particular case, typeof(n) would be fine.

That said, I am curious where you got that code from, currently in Julia factorial uses a lookup table and looks very different.

2 Likes

If I recall that’s just the fallback, and probably used by BigInt, but Int64 etc. has a lookup.

I see that it’s actually in intfuncs.jl, but I am not sure what would call this.

Int & friends use lookup tables, BigInt falls back to using MPZ.

That function gets called for user-defined types that are <:Integer. For example:

using StaticNumbers
factorial(static(20))
1 Like

Is there a particular type(?) that benefits from typeof(n*n) instead of typeof(n)?

I think the point is that you cannot know for certain that n*n will always have the same type as n, for example in types created by users. Therefore this is an extra safeguard. You want to make sure that things will work out, instead of relying on that it “seems to work in most cases.”

1 Like

Seems like this was introduced by @jeff.bezanson in https://github.com/JuliaLang/julia/commit/b6687823d0e4a61f44890325e55aded0d5e07f29 all the way back in 2014. It does seem a bit questionable, whether the typeof(n*n) is sensible here. Maybe

f = one(n)
for i in oftype(n, 2):n
    f *= i
end

might be a better generic implementation? It won’t be type stable if n is not closed under multiplication, but assuming n and 1 can always be converted to typeof(n*n) is probably worse, as seen in the Unitful example. It is also not guaranteed that (n*n) * n can be converted to typeof(n*n), so typeof(n*n) seems rather arbitrary as output type.

1 Like

Thinking about it, oftype(n, 2) might need to be oneunit(n) + oneunit(n) because otherwise, that fails for static(true) from StaticNumbers.

It would also fail for plain old true from Base…

I thought that had a custom implementation, but apparently not.

It doesn’t seem to need one.

julia> @code_llvm factorial(true)

;  @ intfuncs.jl:889 within `factorial'
define i8 @julia_factorial_19411(i8) {
top:
;  @ intfuncs.jl:895 within `factorial'
  ret i8 1
}
2 Likes

But if that is violated, then you cannot know if n * n * n will be the same type, etc, so I am not sure that is worth much.

Frankly, I don’t understand the motivation for this — if anything, I would just do some widening, eg Int8, Int16, etc to Int.