import mir.bignum.fp: Fp;
import mir.math.common: approxEqual;
import mir.math.numeric: prod;
import mir.ndslice.topology: iota;

static assert(is(typeof(factorial(33)) == Fp!128));
static assert(is(typeof(factorial!256(33)) == Fp!256));
static immutable double f33 = 8.68331761881188649551819440128e+36;
static assert(approxEqual(cast(double) factorial(33), f33));

assert(cast(double) factorial(0) == 1);
assert(cast(double) factorial(0, 100) == 1);
assert(cast(double) factorial(1, 100) == 100);
assert(approxEqual(cast(double) factorial(100, 1000), iota([100], 1000).prod!double));