14-byte factorials

[Dirk Laurie <dirk.laurie@...>]

The recent discussion using factorials as an example prompted me,
just as a curiosity, to implement exact factorials up to 31! using the
fact that floating point gives up to 53 high-order bits and integer
arithmetic gives up to 64 low order bits. The total of 117 (=8*14.625)
suggests that up to 14-byte factorials may be achievable. 32! needs
118 bits, but up to 31!, which needs 113 bits turns out to be attainable.

Since Lua has no 16-byte integers, the only way to represent
the result is as a string. The code is a nice illustration of several
features newly available in Lua 5.3.

~~~ file lfac.lua
return function(n)
-- fac(n) returns n! as a Lua integer if n<=20, as a string containing
--    a big-endian unsigned 16-byte integer if 20<n<32, and as a Lua float
--    if n>=32 or if n is not integer-valued.
   local ifac,dfac = 1, 1.
   local imax = 20 -- largest n such that signed 64-bit n! does not overflow
   local dmin = 32 -- smallest n such that the paste trick does not work
   local m=n
   while n>0 do
      ifac = ifac*n
      dfac = dfac*n
      n=n-1
   end
   if n==0 and m<=imax then return ifac end
   if n>0 or m>=dmin then return dfac end
-- paste together correct bits from the float and integer versions
   local M,E = math.frexp(dfac)
   return string.pack(">I8I8",math.floor(M*(1<<(E-64))),ifac)
end
~~~

For those to whom it matters: the above code is hereby placed in the
public domain.

Example:
> lfac = require"lfac"
> f={lfac(31):byte(1,-1)}
> for k=1,16 do io.write(("%02x"):format(f[k])) end
0001956ad0aae33a4560c5cd2c000000