Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]

Groups > comp.lang.forth > #43289 > unrolled thread

Rosetta code binomial coefficient

Back to article view | Back to comp.lang.forth

Contents

  Rosetta code binomial coefficient krishna.myneni@ccreweb.org - 2015-11-06 22:05 -0800
    Re: Rosetta code binomial coefficient mhx@iae.nl - 2015-11-06 23:24 -0800
      Re: Rosetta code binomial coefficient Albert van der Horst <albert@spenarnc.xs4all.nl> - 2015-11-07 10:58 +0000
    Re: Rosetta code binomial coefficient "WJ" <w_a_x_man@yahoo.com> - 2015-11-07 08:14 +0000
      Re: Rosetta code binomial coefficient "WJ" <w_a_x_man@yahoo.com> - 2015-11-07 14:05 +0000
        Re: Rosetta code binomial coefficient "WJ" <w_a_x_man@yahoo.com> - 2015-11-15 06:18 +0000
    Re: Rosetta code binomial coefficient Paul Rubin <no.email@nospam.invalid> - 2015-11-07 00:23 -0800
      Re: Rosetta code binomial coefficient Bernd Paysan <bernd.paysan@gmx.de> - 2015-11-07 15:11 +0100
        Re: Rosetta code binomial coefficient krishna.myneni@ccreweb.org - 2015-11-07 07:49 -0800
        Re: Rosetta code binomial coefficient krishna.myneni@ccreweb.org - 2015-11-07 08:10 -0800
        Re: Rosetta code binomial coefficient Paul Rubin <no.email@nospam.invalid> - 2015-11-07 08:50 -0800
          Re: Rosetta code binomial coefficient Bernd Paysan <bernd.paysan@gmx.de> - 2015-11-07 20:20 +0100
            Re: Rosetta code binomial coefficient krishna.myneni@ccreweb.org - 2015-11-07 17:45 -0800
            Re: Rosetta code binomial coefficient Paul Rubin <no.email@nospam.invalid> - 2015-11-07 19:50 -0800
              Re: Rosetta code binomial coefficient "WJ" <w_a_x_man@yahoo.com> - 2015-11-08 15:35 +0000
                Re: Rosetta code binomial coefficient "WJ" <w_a_x_man@yahoo.com> - 2015-11-25 20:29 +0000
                  Re: Rosetta code binomial coefficient Paul Rubin <no.email@nospam.invalid> - 2015-11-25 12:43 -0800
              Re: Rosetta code binomial coefficient Bernd Paysan <bernd.paysan@gmx.de> - 2015-11-08 18:51 +0100
                Re: Rosetta code binomial coefficient Albert van der Horst <albert@spenarnc.xs4all.nl> - 2015-11-08 18:51 +0000
                Re: Rosetta code binomial coefficient Paul Rubin <no.email@nospam.invalid> - 2015-11-08 19:55 -0800
    Re: Rosetta code binomial coefficient krishna.myneni@ccreweb.org - 2015-11-08 07:03 -0800
    Re: Rosetta code binomial coefficient krishna.myneni@ccreweb.org - 2015-11-09 06:57 -0800
      Re: Rosetta code binomial coefficient krishna.myneni@ccreweb.org - 2015-11-10 03:49 -0800
        Re: Rosetta code binomial coefficient Ron Aaron <rambamist@gmail.com> - 2015-11-10 14:01 +0200
          Re: Rosetta code binomial coefficient krishna.myneni@ccreweb.org - 2015-11-10 19:13 -0800
            Re: Rosetta code binomial coefficient Ron Aaron <rambamist@gmail.com> - 2015-11-11 07:38 +0200
    Re: Rosetta code binomial coefficient glidedog@gmail.com - 2015-11-14 05:26 -0800
      Re: Rosetta code binomial coefficient krishna.myneni@ccreweb.org - 2015-11-14 06:19 -0800
        Re: Rosetta code binomial coefficient krishna.myneni@ccreweb.org - 2015-11-14 06:22 -0800
        Re: Rosetta code binomial coefficient glidedog@gmail.com - 2015-11-14 08:16 -0800
          Re: Rosetta code binomial coefficient krishna.myneni@ccreweb.org - 2015-11-14 11:16 -0800
    Re: Rosetta code binomial coefficient "WJ" <w_a_x_man@yahoo.com> - 2016-06-28 00:02 +0000

Page 2 of 2 — ← Prev page 1 [2]

#43366

On Saturday, November 7, 2015 at 12:05:04 AM UTC-6, krishna...@ccreweb.org wrote:
> The Forth entry on Rosetta Code for computing the binomial coefficient is beautiful. Compare its compactness of implementation with those in other languages. This is a very nice code example in Forth. Although it meets the requirements of the entry, single integers for 32 bit word size quickly lead to failure, e.g. "60 30 choose". The implementation is better suited for 64-bit word size.
> 
> : choose ( n k -- nCk ) 1 swap 0 ?do over i - i 1+ */ loop nip ;
> 
...

Here is the code from FSL Algorithm 59, permcomb.x. It allows computation of the binomial coefficient for up to n = 512, on a 32-bit system, by using a double length sum and an intermediate triple length product, with standard Forth.

512 256 combs d.
9191530535263022392  ok

Krishna

=========================================

VARIABLE temp  \ private to combs

: combs ( u1 u2--ud)  1 S>D 2SWAP
                       2DUP - MIN
                       SWAP temp !
                       1+ 1 ?DO  temp @  I M*/
                                 -1 temp +!
                            LOOP ;
\
\ return nCr, where u1=n u2=r. All arguments are unsigned, result is double.
\
\ This is an iterative version of the recurrence;
\      r=0 --> nCr = 1
\      r>0 --> nCr = nC(r-1)(n-r+1)/r
\
\ This recurrance was chosen in favour of the more common
\      nCr = n!/(n-r)! r!
\ to avoid excessively large intermediate results. Use of integer maths
\ necessitates that the multiplication be done before the division, to avoid
\ truncation errors, hence the use of M*/, which has a triple length
\ intermediate result. Advantage is taken of the symmetry of the function
\ to minimise the number of iterations.

[toc] | [prev] | [next] | [standalone]

#43430

On Saturday, November 7, 2015 at 12:05:04 AM UTC-6, krishna...@ccreweb.org wrote:
> The Forth entry on Rosetta Code for computing the binomial coefficient is beautiful. Compare its compactness of implementation with those in other languages. This is a very nice code example in Forth. Although it meets the requirements of the entry, single integers for 32 bit word size quickly lead to failure, e.g. "60 30 choose". The implementation is better suited for 64-bit word size.
> 
> : choose ( n k -- nCk ) 1 swap 0 ?do over i - i 1+ */ loop nip ;
> 

The FSL module permcomb.x, FSL Algorithm #59 by Gordon Charlton, uses a double length number and M*/ for a triple length intermediate product, to compute the correct double length result for n up to 66 (all k) on a 32-bit machine. Pascal's triangle formed by addition should allow at least n = 67 when double length integers are used on a 32-bit system.

Pascal's triangle suggests a strategy for computing binomial coefficients for larger n:

1) Use single length integers to store a Pascal's triangle array of binomial coefficients, up to the maximum unsigned integer representation for 32/64 bit.

2) Store an extended portion of Pascal's triangle computed with double length integers, until the maximum double length representation is reached.

3) Use floating point calculations (with addition only) to obtain binomial coefficients up to some maximum n value, N. No coefficients are stored.

4) Use logarithmic forms for calculating binomial coefficients for n > N, up to the maximum IEEE double prec. fp value.

5) Provide words which can return single integer, double length integer, and floating point binomial coefficients: two variants, ones without a check on the upper limit of n for each type of output, and ones with a bounds check on n.

Krishna

[toc] | [prev] | [next] | [standalone]

#43445

On Monday, November 9, 2015 at 8:57:43 AM UTC-6, krishna...@ccreweb.org wrote:
> On Saturday, November 7, 2015 at 12:05:04 AM UTC-6, krishna...@ccreweb.org wrote:
> > The Forth entry on Rosetta Code for computing the binomial coefficient is beautiful. Compare its compactness of implementation with those in other languages. This is a very nice code example in Forth. Although it meets the requirements of the entry, single integers for 32 bit word size quickly lead to failure, e.g. "60 30 choose". The implementation is better suited for 64-bit word size.
> > 
> > : choose ( n k -- nCk ) 1 swap 0 ?do over i - i 1+ */ loop nip ;
> > 
> 
> The FSL module permcomb.x, FSL Algorithm #59 by Gordon Charlton, uses a double length number and M*/ for a triple length intermediate product, to compute the correct double length result for n up to 66 (all k) on a 32-bit machine. Pascal's triangle formed by addition should allow at least n = 67 when double length integers are used on a 32-bit system.
> 
> Pascal's triangle suggests a strategy for computing binomial coefficients for larger n:
> 
> 1) Use single length integers to store a Pascal's triangle array of binomial coefficients, up to the maximum unsigned integer representation for 32/64 bit.
> 
> 2) Store an extended portion of Pascal's triangle computed with double length integers, until the maximum double length representation is reached.
> 
> 3) Use floating point calculations (with addition only) to obtain binomial coefficients up to some maximum n value, N. No coefficients are stored.
> 
> 4) Use logarithmic forms for calculating binomial coefficients for n > N, up to the maximum IEEE double prec. fp value.
> 
> 5) Provide words which can return single integer, double length integer, and floating point binomial coefficients: two variants, ones without a check on the upper limit of n for each type of output, and ones with a bounds check on n.
> 
> Krishna

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

[toc] | [prev] | [next] | [standalone]

#43446

On 11/10/2015 13:49, krishna.myneni@ccreweb.org wrote:
> =======
> binomial(66, 33);
> (%o5)                         7219428434016265740
>
> binomial(67, 33);
> (%o6)                        14226520737620288370
>
> binomial(130,65);
> (%o3)               95067625827960698145584333020095113100

In 8th no need for special number words, just implement the algo:

: perm \ p k -- P(p,k)
   over swap \ p p k
   n:- n:1+   \ p p-k+1
   over n:1-
   ' n:* -rot loop ;

: comb \ p k -- C(p,k) == P(p,k) / fact(k)
   tuck perm  \ k P(p,k)
   swap 2 swap
   ' n:/ -rot loop ;

ok> 130 65 comb . cr
95067625827960698145584333020095113100
ok> 200 65 comb . cr
355410124430498052252434369613867238863579093582797880

[toc] | [prev] | [next] | [standalone]

#43462

On Tuesday, November 10, 2015 at 6:01:20 AM UTC-6, Ron Aaron wrote:
...
> In 8th no need for special number words, just implement the algo:
> 
> : perm \ p k -- P(p,k)
>    over swap \ p p k
>    n:- n:1+   \ p p-k+1
>    over n:1-
>    ' n:* -rot loop ;
> 
> : comb \ p k -- C(p,k) == P(p,k) / fact(k)
>    tuck perm  \ k P(p,k)
>    swap 2 swap
>    ' n:/ -rot loop ;
...

Congrats on your 8th language. I haven't looked at it, but from your examples above, it appears to be a Forth like language that supports a higher level of abstraction, and includes arbitrary precision arithmetic. I use a variety of low and high level languages/environments for my computing needs, from assembly language to Forth to scientific computing environments such as R and Maxima. At present, I'm not looking to learn another environment. For reasons of convenience, familiarity, and ability to write code over a range of levels, from bare metal to modular programming, Forth suits quite a few of my computing needs, and I occasionally want to port functions/library routines from higher level systems to Forth. Sometimes this is just for a check on the output of higher level systems where I can't see what's going on under the hood very easily.

Cheers,
Krishna

[toc] | [prev] | [next] | [standalone]

#43464

On 11/11/2015 05:13, krishna.myneni@ccreweb.org wrote:

>
> Congrats on your 8th language. I haven't looked at it, but from your examples above, it appears to be a Forth like language that supports a higher level of abstraction, and includes arbitrary precision arithmetic. I use a variety of low and high level languages/environments for my computing needs, from assembly language to Forth to scientific computing environments such as R and Maxima. At present, I'm not looking to learn another environment. For reasons of convenience, familiarity, and ability to write code over a range of levels, from bare metal to modular programming, Forth suits quite a few of my computing needs, and I occasionally want to port functions/library routines from higher level systems to Forth. Sometimes this is just for a check on the output of higher level systems where I can't see what's going on under the hood very easily.

Hi, Krishna -

Indeed, it's exactly as you've described it.  The tack we took with 8th 
regarding numbers is that the end-user usually just wants to get the 
work done, and doesn't really care about 'integers' vs 'floats' vs 'big 
numbers' etc.  So 8th takes care of the low-level details and as far as 
the user is concerned there is (usually) no reason to care.

8th is well suited to applications development (and the reason it was 
written to begin with); however, it's not suitable for bare-metal stuff 
though there are accessor words for things like Bluetooth etc.

Best regards,
Ron

[toc] | [prev] | [next] | [standalone]

#43532

On 11/7/15 1:05 AM, krishna...@ccreweb.org wrote:
> The Forth entry on Rosetta Code for computing the binomial
> coefficient is beautiful. Compare its compactness of implementation
> with those in other languages. This is a very nice code example in
> Forth.

> : choose ( n k -- nCk ) 1 swap 0 ?do over i - i 1+ */ loop nip ;

An Oforth solution can be pretty good I think:

: choose // n k -- nCk
  | i | 1 swap loop: i [ over i - 1 + * i / ] nip ;
200 100 choose .
90548514656103281165404177077484163874504589675413336841320

For those not familiar with Oforth:

1) | i | is an uninitialized local variable declaration, not Forth's i

2) " ( limit ) loop: y [ ... ]" is a looping construct
    - loop: removes one stack item which becomes the loop limit
    - y must be a declared local variable
    - y steps from 1 through loop limit inclusive
    - the code in the [ ... ] block is executed loop limit times

The same algorithm in standard Forth would be:
: choose \ n k -- nCk
  1 swap 1+ 1 do over i - 1 + * i / loop nip ;

Note use of the same stack manipulation words as in Oforth
and exactly the same code inside the loop.

-Doug

[toc] | [prev] | [next] | [standalone]

#43533

On Saturday, November 14, 2015 at 7:26:53 AM UTC-6, glid...@gmail.com wrote:
> On 11/7/15 1:05 AM, krishna...@ccreweb.org wrote:
> > The Forth entry on Rosetta Code for computing the binomial
> > coefficient is beautiful. Compare its compactness of implementation
> > with those in other languages. This is a very nice code example in
> > Forth.
> 
> > : choose ( n k -- nCk ) 1 swap 0 ?do over i - i 1+ */ loop nip ;
> 
> An Oforth solution can be pretty good I think:
> 
> : choose // n k -- nCk
>   | i | 1 swap loop: i [ over i - 1 + * i / ] nip ;
> 200 100 choose .
> 90548514656103281165404177077484163874504589675413336841320
> 
...
> The same algorithm in standard Forth would be:
> : choose \ n k -- nCk
>   1 swap 1+ 1 do over i - 1 + * i / loop nip ;
> 
...

Interesting example of Oforth. However, please note that your standard Forth version, above, does not use */ for an intermediate double length, as does the Rosetta code example, so it overflows earlier, e.g.

30 15 choose u.

gives the correct answer in the Rosetta code example, but not in your standard Forth version.

Regards,
Krishna

[toc] | [prev] | [next] | [standalone]

#43534

On Saturday, November 14, 2015 at 8:19:42 AM UTC-6, krishna...@ccreweb.org wrote:
> On Saturday, November 14, 2015 at 7:26:53 AM UTC-6, glid...@gmail.com wrote:
...
> > The same algorithm in standard Forth would be:
> > : choose \ n k -- nCk
> >   1 swap 1+ 1 do over i - 1 + * i / loop nip ;
> > 
> ...
> 
> Interesting example of Oforth. However, please note that your standard Forth version, above, does not use */ for an intermediate double length, as does the Rosetta code example, so it overflows earlier, e.g.
> 
> 30 15 choose u.
> 
> gives the correct answer in the Rosetta code example, but not in your standard Forth version.
> 

I mean on a 32-bit system, for the particular example I gave. The same problem would exist on a 64-bit system, but for a larger value of n.

[toc] | [prev] | [next] | [standalone]

#43537

On Saturday, November 14, 2015 at 9:19:42 AM UTC-5, krishna...@ccreweb.org wrote:
> On Saturday, November 14, 2015 at 7:26:53 AM UTC-6, glid...@gmail.com wrote:
> > On 11/7/15 1:05 AM, krishna...@ccreweb.org wrote:
> > > The Forth entry on Rosetta Code for computing the binomial
> > > coefficient is beautiful. Compare its compactness of implementation
> > > with those in other languages. This is a very nice code example in
> > > Forth.
> > 
> > > : choose ( n k -- nCk ) 1 swap 0 ?do over i - i 1+ */ loop nip ;
> > 
> > An Oforth solution can be pretty good I think:
> > 
> > : choose // n k -- nCk
> >   | i | 1 swap loop: i [ over i - 1 + * i / ] nip ;
> > 200 100 choose .
> > 90548514656103281165404177077484163874504589675413336841320
> > 
> ...
> > The same algorithm in standard Forth would be:
> > : choose \ n k -- nCk
> >   1 swap 1+ 1 do over i - 1 + * i / loop nip ;
> > 
> ...
> 
> Interesting example of Oforth. However, please note that your standard Forth version, above, does not use */ for an intermediate double length, as does the Rosetta code example, so it overflows earlier, e.g.

Yes. I knew that.  The idea was to illustrate how Forth-like Oforth is.
The Oforth solution has little, if any, issues with overflow as it will automatically use bignums as required.  Nothing like */ required.  Not saying Oforth is better than standard Forth necessarily.  
But it can do this kind of calculation in a straightforward manner.

-Doug

[toc] | [prev] | [next] | [standalone]

#43545

On Saturday, November 14, 2015 at 10:16:28 AM UTC-6, glid...@gmail.com wrote:
> ... The idea was to illustrate how Forth-like Oforth is.
> The Oforth solution has little, if any, issues with overflow as it will automatically use bignums as required.  Nothing like */ required.  Not saying Oforth is better than standard Forth necessarily.  
> But it can do this kind of calculation in a straightforward manner.
> 

Yes. Oforth appears to be a higher level language than standard Forth, and frees the user from having to worry about the mundane details of arithmetic precision. The application space for Oforth is likely a bit different than the one for standard Forth, or for other lower level languages. For applications where efficiency is a big concern, the coder must bound the problem and use the least computationally intensive solution, or at least have full control of the computational complexity, rather than the system deciding the complexity. On the other hand, where the necessary bounds are not known in advance, it is very useful to have a system that handles such details.

Cheers,
Krishna

[toc] | [prev] | [next] | [standalone]

#48590

krishna.myneni@ccreweb.org wrote:

> The Forth entry on Rosetta Code for computing the binomial coefficient is
> beautiful. Compare its
> compactness of implementation with those in other languages
> . This is a very nice code example in Forth. Although
> it meets the requirements of the entry
> , single integers for 32 bit word size quickly lead to failure, e.
> g. "60 30 choose". The implementation is better suited for 64-bit word size.
> 
> : choose ( n k -- nCk ) 1 swap 0 ?do over i - i 1+ */ loop nip ;

OCaml:

#load "nums.cma";;
open Num;;
 
let rec choose n = function
    0 -> Int 1
  | k -> Int (n-k+1) */ (choose n (k-1)) // Int k ;;

string_of_num (choose 200 100);;
 ===>
"90548514656103281165404177077484163874504589675413336841320"

-- 
If confirmed, Garland would be the fourth Jewish justice on the nation's
highest court, which is comprised entirely of Jews and Catholics. The three
current Jewish members of the Supreme Court are Ruth Bader Ginsburg, Elana
Kagan, and Stephen Breyer.
www.jta.org/2016/03/16/news-opinion/united-states/obama-to-name-jewish-federal-judge-to-supreme-court

[toc] | [prev] | [standalone]

Page 2 of 2 — ← Prev page 1 [2]

Back to top | Article view | comp.lang.forth

csiph-web