Add nChoosek function to compute binomial coefficients.

README.md

@@ -80,10 +80,13 @@ Some basic functions are implemented for analyzing data
 
 The normal distribution is implemented.
 
-The Chi^2 distribution is in development.
+The Chi^2 distribution is implemented.
 
 ##  Functions
 These are special functions not included with Ada's elementary functions
 that are useful or needed by other routines here.
 * gamma2n - Divides the positive integer argument by two and computes the gamma function.  Used to support chi^2.
 * lngamma2n - Natural log of gamma2n.  Allows larger values of n without overflow.
+* factorial - Computes the factorial of n.
+* lnfact - Computes the natual log of the factorial of n allowing larger values of n without overflow.
+* nChoosek - Computes binomial coefficent based on n and k.

src/BBS-Numerical-functions_real.adb

@@ -62,7 +62,7 @@ package body BBS.Numerical.functions_real is
    --  around n = 35.
    --
    function factorial(n : Natural) return f'Base is
-      base : f'Base := 1.0;
+      base : f'Base := 1.0;  --  0! is defined as 1
    begin
       for i in 1 .. n loop
          base := base*f'Base(i);
@@ -74,12 +74,43 @@ package body BBS.Numerical.functions_real is
    --  values of n before overflowing.
    --
    function lnfact(n : Natural) return f'Base is
-      base : f'Base := 0.0;
+      base : f'Base := 0.0;  --  0! is defined as 1
    begin
       for i in 1 .. n loop
          base := base + elem.Log(f'Base(i));
       end loop;
       return base;
    end;
-
+   --
+   --  Compute the binomial coefficient - n choose k.  Note that the result is
+   --  an integer value, but f'Base is used to allow greater range.
+   --
+   --   (n)     n!
+   --   ( )  --------
+   --   (k)  k!(n-k)!
+   --
+   --  Expanding the factorials and cancelling terms leads to the following
+   --  formula.  Note that since it is symmetrical with respect to k, the limit
+   --  for the product can be the lesser of k or (n-k).
+   --
+   --   (n)    k   n + 1 - i
+   --   ( ) = PROD ---------
+   --   (k)   i=1      i
+   --
+   function nChoosek(n, k : Natural) return f'Base is
+      base : f'Base := 1.0;
+      top  : constant f'Base := f'Base(n + 1);
+      limit : Integer;
+   begin
+      if k < (n-k) then
+         limit := k;
+      else
+         limit := n - k;
+      end if;
+      for i in 1 .. limit loop
+         base := base*(top - f'Base(i))/f'Base(i);
+      end loop;
+      return base;
+   end;
+   --
 end BBS.Numerical.functions_real;

src/BBS-Numerical-functions_real.ads

@@ -36,4 +36,24 @@ package BBS.Numerical.functions_real is
    --
    function lnfact(n : Natural) return f'Base;
    --
+   --  Compute the binomial coefficient - n choose k.  Note that the result is
+   --  an integer value, but f'Base is used to allow greater range.
+   --
+   --   (n)      n!
+   --   ( ) = --------
+   --   (k)   k!(n-k)!
+   --
+   --  Expanding the factorials and cancelling terms leads to the following
+   --  formula.  Note that since it is symmetrical with respect to k, the limit
+   --  for the product can be the lesser of k or (n-k).
+   --
+   --   (n)    k   n + 1 - i
+   --   ( ) = PROD ---------
+   --   (k)   i=1      i
+   --
+   --  This also allows larger values of n and k without overflowing.
+   --
+   function nChoosek(n, k : Natural) return f'Base
+     with pre => (n >= k);
+   --
 end BBS.Numerical.functions_real;

test/test_stats.adb

@@ -53,4 +53,13 @@ begin
       Ada.Text_IO.New_Line;
    end loop;
    Ada.Text_IO.New_Line;
+   for n in 0 .. 10 loop
+      Ada.Text_IO.put("  ");
+      float_io.put(real(n), 2, 0, 0);
+      for k in 0 .. n loop
+         Ada.Text_IO.Put("  ");
+         float_io.Put(funct.nchoosek(n, k), 4, 0, 0);
+      end loop;
+      Ada.Text_IO.New_Line;
+   end loop;
 end test_stats;