(phixonline)-->
with javascript_semantics
function P(integer n,k)
return factorial(n)/factorial(n-k)
end function
function C(integer n,k)
return P(n,k)/factorial(k)
end function
function lstirling(atom n)
if n<10 then
return lstirling(n+1)-log(n+1)
end if
return 0.5*log(2*PI*n) + n*log(n/E + 1/(12*E*n))
end function
function P_approx(integer n, k)
return lstirling(n)-lstirling(n-k)
end function
function C_approx(integer n, k)
return lstirling(n)-lstirling(n-k)-lstirling(k)
end function
function to_s(atom v)
integer e = floor(v/log(10))
return sprintf("%.9ge%d",{power(E,v-e*log(10)),e})
end function
-- Test code
printf(1,"=> Exact results:\n")
for n=1 to 12 do
integer p = floor(n/3)
printf(1,"P(%d,%d) = %d\n",{n,p,P(n,p)})
end for
for n=10 to 60 by 10 do
integer p = floor(n/3)
printf(1,"C(%d,%d) = %d\n",{n,p,C(n,p)})
end for
printf(1,"=> Floating point approximations:\n")
constant tests = {5, 50, 500, 1000, 5000, 15000}
for i=1 to length(tests) do
integer n=tests[i], p = floor(n/3)
printf(1,"P(%d,%d) = %s\n",{n,p,to_s(P_approx(n,p))})
end for
for n=100 to 1000 by 100 do
integer p = floor(n/3)
printf(1,"C(%d,%d) = %s\n",{n,p,to_s(C_approx(n,p))})
end for