Below is the double length integer version of Bernd Paysan's code for computing Pascal's triangle. On 32-bit systems, it provides valid binomial coefficients up to n = 67, and on 64-bit systems, up to n = 130.
Now I need to flesh out the floating point code to deal with n > 130:
: fchoose ( k n -- r ) \ or ( n k -- ) ( F: -- r )
dup dmaxchoose# < IF
dout-of D>F
ELSE
\ Compute the combinations for n >= dmaxchoose#
0e
THEN ;
Krishna
============================================
kForth (32-bit system):
========================
33 66 dout-of d.
7219428434016265740 ok
33 67 dout-of d. \ kForth needs ud.
-4220223336089263246 ok
gforth (64-bit system):
=======================
65 130 dout-of ud. 95067625827960698145584333020095113100 ok
For comparison,
Maxima:
=======
binomial(66, 33);
(%o5) 7219428434016265740
binomial(67, 33);
(%o6) 14226520737620288370
binomial(130,65);
(%o3) 95067625827960698145584333020095113100
The double integer version of B. Paysan's code:
\ Double length Pascal's triangle
cell 4 = [IF] 67 [ELSE] 130 [THEN] constant dmaxchoose#
dmaxchoose# 1+ tri cells 2* buffer: dpascal-triangle
1 s>d dpascal-triangle 2!
: dtri' ( k n -- addr ) tri + cells 2* dpascal-triangle + ;
: dout-of ( k n -- dresult )
2dup u> IF 2drop 0 s>d ELSE dtri' 2@ THEN ;
: dout-of+ ( k n -- dresult )
1- over 1- over dout-of 2>r dout-of 2r> d+ ;
: dfill-line ( n -- )
dup 1+ 0 DO I over 2dup dout-of+ 2over dtri' 2! 2drop LOOP drop ;
: dfill-triangle ( n -- )
1+ 1 DO I dfill-line LOOP ;
dmaxchoose# dfill-triangle