x86 32-bit machine code, 69 bytes
\x6a\x01\x58\x89\x03\x89\xc1\xc1\xe1\x02\x29\xcc\x89\xe7\x89\xde\xf3\xa4
\x89\xc5\x89\xd1\x49\x7e\x25\x31\xf6\x89\xf7\xf7\xdf\x8b\x3c\xb4\x11\x3c
\xb3\x19\xff\x46\x39\xee\x7c\xf1\xf7\xdf\x74\x07\x83\x04\xb3\x01\x46\x72
\xf9\x39\xf0\x0f\x4c\xc6\xeb\xd8\x8d\x24\xac\x4a\x7f\xc1\xc3
Try it online!
Unlike other machine code entries, I implemented arbitrary-precision multiplication from scratch.
The algorithm is the most basic one. For \$a × b\$, \$a\$ is added \$b\$ times with a loop. It doesn't scale well for really big numbers, but 125! is calculated within 0.3 seconds (user time) in TIO.
The function outputs an arbitrary-precision integer in base \$2 ^ {32}\$. The challenge doesn't require the function to output a string, so I can use any external method to convert the output to printable string apart from the challenge. The huge footer in TIO is the source code of mini-GMP that I copy-pasted for bignum-to-decimal-string conversion. TIO has GMP, but only for 64-bit, unfortunately.
The code is quite straightforward. Here is a pseudocode to illustrate what the assembly code is doing.
for x = y!:
x = 1
i from y to 1:
z = x
j from y - 1 to 1:
x += z
assembly (nasm)
; input: edx, output: ebx (bignum array), eax (number of "limbs")
; custom calling convention, everything except `ebx` and `esp` gets dirty
_fac:
push 1
pop eax
mov [ebx], eax
.L0:
mov ecx, eax
shl ecx, 2
sub esp, ecx
mov edi, esp
mov esi, ebx
rep movsb
mov ebp, eax
mov ecx, edx
.L1:
dec ecx
jle .1
xor esi, esi
mov edi, esi
.L2:
neg edi
mov edi, [esp + esi * 4]
adc [ebx + esi * 4], edi
sbb edi, edi
inc esi
cmp esi, ebp
jl .L2
neg edi
jz .0
.L3:
add dword [ebx + esi * 4], 1
inc esi
jc .L3
.0:
cmp eax, esi
cmovl eax, esi
jmp .L1
.1:
lea esp, [esp + ebp * 4]
dec edx
jg .L0
ret
