ANOVA Using Raw C#

Recently, over a span of several days, I had been reviewing how to implement the incomplete beta function, Ix(a,b), from scratch using the C# language. Once I had Ix(a,b) working (which was no small task), I figured I might as well use my Ix(a,b) implementation to do an end-to-end demo of Analysis of Variance (ANOVA). The key idea is that in an ANOVA problem you compute the right tail of the F-distribution using the Ix(a,b) function.

The goal of an ANOVA problem is to analyze three of more samples of numeric data and determine if all three samples have the same source population means. For example, if you have sample IQ test scores from White, Black, and Asian students you can use ANOVA to estimate the likelihood that all three groups have the same mean IQ.

Briefly, starting with grouped data, you do several intermediate calculations and end up with an F-statistic and two degrees of freedom (df1, df2) values. Then the area under F from negative infinity to the F-stat (the left tail) is given by:

x = df2 / (df2 + df1 * fstat)
a = df2 / 2
b = df1 / 2
left_tail = IncBeta(x, a, b)

The area from the F-stat to positive infinity (the right tail) is 1 – left_tail and is the estimate that all group population means are equal. In my demo, the p-value is 0.0004 which means it is highly unlikely that all three groups have the same population mean. (It’s pretty obvious that the mean of Group 1 is smaller than that of Groups 2 and 3).

Like all classical statistics techniques, ANOVA makes several assumptions (variances of group data are equal, samples are Gaussian distributed, etc.) which can’t be completely verified, and furthermore, interpretation of the ANOVA results is very subtle and tricky.

But still, ANOVA can be useful in some practical scenarios.



ANOVA is based on the variance/variability of numbers. Artist Edwin Georgi (1896-1964) used color variability to good effect. Here are three examples of his work. Notice that almost all of the colors are unnatural, such as the green, orange, and purple skin tones.


Demo code. Replace “lt”, “gt”, “lte”, “gte” with Boolean operator symbols.

using System;
namespace AnovaCSharp
{
  internal class AnovaProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin ANOVA using C# demo \n");
      Console.WriteLine("The goal is 3 or more population");
      Console.WriteLine("means are all equal. \n");

      Console.WriteLine("The sample data is: \n");

      double[][] data = new double[3][]; // 3 groups
      data[0] = new double[] { 3, 4, 6, 5 };
      data[1] = new double[] { 8, 12, 9, 11, 10, 8 };
      data[2] = new double[] { 13, 9, 11, 8, 12 };
      string[] colNames = new string[] { "Group1", "Group2",
        "Group3" };
      ShowData(data, colNames); 

      Console.WriteLine("\nCalculating F-statistic");
      int[] df = null; // degrees of freedom
      double fstat = Fstat(data, out df);
      Console.WriteLine("F stat = " + fstat.ToString("F3"));

      Console.Write("The degrees of freedom are ");
      Console.WriteLine(df[0] + ", " + df[1]);

      Console.WriteLine("\nCalculating p-value from F stat");
      double pValue = FDist(fstat, df[0], df[1]);

      Console.Write("p-value = ");
      Console.WriteLine(pValue.ToString("F8"));

      Console.WriteLine("\nThe p-value is prob all ");
      Console.WriteLine("group means are equal.");
      Console.WriteLine("Interpreting the p-value is subtle.");

      Console.WriteLine("\nEnd ANOVA demo ");
      Console.ReadLine();
    } // Main()

    static double Fstat(double[][] data, out int[] df)
    {
      // calculate F statistic and degrees freedom num, denom
      // assumes data has specific structure:
      // data[0] : 3, 4, etc (group 1)
      // data[1] : 8, 12, etc. (group 2)
      // etc.

      int K = data.Length; // number groups
      int[] n = new int[K]; // number items each group
      int N = 0; // total number data points
      for (int i = 0; i "lt" K; ++i)
      {
        n[i] = data[i].Length;
        N += data[i].Length;
      }

      // 1. group means and grand mean
      double[] means = new double[K];
      double gMean = 0.0;
      for (int i = 0; i "lt" K; ++i)
      {
        for (int j = 0; j "lt" data[i].Length; ++j)
        {
          means[i] += data[i][j];
          gMean += data[i][j];
        }
        means[i] /= n[i];

        Console.WriteLine("Group [" + i + "] mean = " + 
          means[i].ToString("F2"));
      }
      gMean /= N;
      Console.WriteLine("Overall mean = " +
        gMean.ToString("F2"));

      // 2. SSb and MSb
      double SSb = 0.0;
      for (int i = 0; i "lt" K; ++i)
        SSb += n[i] * (means[i] - gMean) * (means[i] - gMean);
      double MSb = SSb / (K - 1);

      // 3. SSw and MSw
      double SSw = 0.0;
      for (int i = 0; i "lt" K; ++i)
        for (int j = 0; j "lt" data[i].Length; ++j)
          SSw += (data[i][j] - means[i]) * (data[i][j] -
            means[i]);
      double MSw = SSw / (N - K);

      // for demo only
      Console.WriteLine("");
      Console.WriteLine("Calculated SSb = " + SSb.ToString("F4"));
      Console.WriteLine("Calculated MSb = " + MSb.ToString("F4"));
      Console.WriteLine("Calculated SSw = " + SSw.ToString("F4"));
      Console.WriteLine("Calculated MSw = " + MSw.ToString("F4"));
      Console.WriteLine("");
      Console.Write("F stat = MSb / MSw");
      Console.WriteLine(" = " + MSb.ToString("F4") + " / " +
        MSw.ToString("F4"));
      //Console.WriteLine("");

      df = new int[2]; // store df values
      df[0] = K - 1;
      df[1] = N - K;

      double F = MSb / MSw;
      return F;
    } // Fstat

    static double FDist(double fstat, double df1, double df2)
    {
      // right tail of F-dist past fstat
      // https://www.itl.nist.gov/div898/handbook/..
      //   eda/section3/eda3665.htm
      double k = df2 / (df2 + df1 * fstat);
      double leftTail = 1 - IncBeta(k, df2 / 2, df1 / 2);
      return 1 - leftTail;  // the right tail
    }

    static double ContFraction(double x, double a, double b,
      int maxTerms = 200)
    {
      // 1. pre-compute 200 d values
      double[] d = new double[maxTerms];  // d[0] not used

      int end = maxTerms / 2;  // 100
      for (int m = 0; m "lt" end; ++m)  // [0,99]
      {
        int i = 2 * m;  // even di
        int j = i + 1;  // odd di
        d[i] = (m * (b - m) * x) / ((a + 2 * m - 1) *
          (a + 2 * m));
        d[j] = -1 * ((a + m) * (a + b + m) * x) / ((a + 2 * m) *
          (a + 2 * m + 1));
      }

      // 2. work backwards
      double[] t = new double[maxTerms];  // t[0] not used
      //t[4] = 1 + d[4] / 1;
      //t[3] = 1 + d[3] / t[4];
      //t[2] = 1 + d[2] / t[3];
      //t[1] = 1 + d[1] / t[2];
      //double cf = 1 / t[1];

      t[maxTerms - 1] = 1 + d[maxTerms - 1];
      for (int j = maxTerms - 2; j "gte" 1; --j)
        t[j] = 1 + d[j] / t[j + 1];

      // ShowVector(t);
      // Console.ReadLine();

      double cf = 1 / t[1];
      return cf;
    }

    static double IncompleteBeta(double x, double a, double b)
    {
      // double top = Math.Pow(x, a) * Math.Pow((1 - x), b);
      // double bot = a * Beta(a, b);
      double cf = ContFraction(x, a, b);

      double logTop = (a * Math.Log(x)) +
        (b * Math.Log(1 - x));
      double logBot = Math.Log(a) + LogBeta(a, b);
      double logLeft = logTop - logBot;
      double left = Math.Exp(logLeft);

      return left * cf;  // should be [0.0, 1.0]
    }

    static double IncBeta(double x, double a, double b)
    {
      // pick the form that converges best
      if (x "lt" (a + 1.0) / (a + b + 2.0))
        return IncompleteBeta(x, a, b);
      else
        return 1.0 - IncompleteBeta(1 - x, b, a);
    }

    static double LogGamma(double z)
    {
      // plain Gamma() can easily overflow
      // Lanczos approximation g=5, n=7
      double[] coef = new double[7] { 1.000000000190015,
        76.18009172947146, -86.50532032941677,
        24.01409824083091, -1.231739572450155,
        0.1208650973866179e-2, -0.5395239384953e-5 };

      double LogSqrtTwoPi = 0.91893853320467274178;

      if (z "lt" 0.5) // Gamma(z) = Pi / (Sin(Pi*z))* Gamma(1-z))
        return Math.Log(Math.PI / Math.Sin(Math.PI * z)) -
          LogGamma(1.0 - z);

      double zz = z - 1.0;
      double b = zz + 5.5; // g + 0.5
      double sum = coef[0];

      for (int i = 1; i "lt" coef.Length; ++i)
        sum += coef[i] / (zz + i);

      return (LogSqrtTwoPi + Math.Log(sum) - b) +
        (Math.Log(b) * (zz + 0.5));
    }

    static double LogBeta(double a, double b)
    {
      // Beta(a,b) = (Gamma(a) * Gamma(b)) / (Gamma(a+b))
      // but plain Gamma() can easily overflow, so:
      // Log( Beta(a,b) ) 
      //   = Log( (Gamma(a) * Gamma(b)) / (Gamma(a+b) ) 
      //   = LogGamma(a) + LogGamma(b) - LogGamma(a+b)
      double logbeta = LogGamma(a) + LogGamma(b) -
        LogGamma(a + b);
      return logbeta;
    }

    //static double Beta(double a, double b)
    //{
    //  return Math.Exp(LogBeta(a, b));
    //}

    static void ShowData(double[][] data, string[] colNames)
    {
      for (int i = 0; i "lt" data.Length; ++i)
      {
        Console.Write(colNames[i] + ": ");
        for (int j = 0; j "lt" data[i].Length; ++j)
        {
          Console.Write(data[i][j].ToString("F1").PadLeft(7));
        }
        Console.WriteLine("");
      }
    }

    //static void ShowVector(double[] v)
    //{
    //  for (int i = 0; i "lt" v.Length; ++i)
    //    Console.Write(v[i].ToString("F4") + " ");
    //  Console.WriteLine("");
    //}

  } // Program

} // ns
Posted in Miscellaneous | Leave a comment

Why PyTorch Layer Definition Order Sometimes Matters

Does the order in which you define layers of a PyTorch neural network matter? Answer: sometimes.

This idea was pointed out to me in a recent workshop I taught. One of the workshop examples was a multi-class classification problem where the goal was to predict the job type of an employee (0 = mgmt, 1 = supp, 2 = tech) from sex (-1, +1), age, city (100 = anaheim, 010 = boulder, 001 = concord), and income. The example defined a 6-(10-10)-3 network like so:

class Net(T.nn.Module):
  def __init__(self):
    super(Net, self).__init__()
    self.hid1 = T.nn.Linear(6, 10)  # 6-(10-10)-3
    self.hid2 = T.nn.Linear(10, 10)
    self.oupt = T.nn.Linear(10, 3)
   
    T.nn.init.xavier_uniform_(self.hid1.weight)
    T.nn.init.zeros_(self.hid1.bias)
    T.nn.init.xavier_uniform_(self.hid2.weight)
    T.nn.init.zeros_(self.hid2.bias)
    T.nn.init.xavier_uniform_(self.oupt.weight)
    T.nn.init.zeros_(self.oupt.bias)

  def forward(self, x):
    z = T.tanh(self.hid1(x))
    z = T.tanh(self.hid2(z))
    z = self.oupt(z)  # no softmax: CrossEntropyLoss() 
    return z

I pointed out that the order in which the hid1, hid2, oupt layers are defined int __init__() doesn’t matter — it’s the order in the forward() method that matters. For example, defining layers in this order has no effect:

  super(Net, self).__init__()
  self.hid2 = T.nn.Linear(10, 10)
  self.oupt = T.nn.Linear(10, 3)
  self.hid1 = T.nn.Linear(6, 10) 

But I was only partially correct. It’s true that the order in which the layers are defined doesn’t matter but only because I use explicit weight and bias initialization.

Suppose instead of explicit weight and bias initialization, the network uses default initialization:

class Net(T.nn.Module):
  def __init__(self):
    super(Net, self).__init__()
    self.hid1 = T.nn.Linear(6, 10) 
    self.hid2 = T.nn.Linear(10, 10)
    self.oupt = T.nn.Linear(10, 3)
   
  def forward(self, x):
    z = T.tanh(self.hid1(x))
    z = T.tanh(self.hid2(z))
    z = self.oupt(z) 
    return z

In this scenario the weights and biases are initialized in the order in which they’re listed in the __init__() method. Therefore, the following definition would yield slightly different results:

class Net(T.nn.Module):
  def __init__(self):
    super(Net, self).__init__()
    self.hid2 = T.nn.Linear(10, 10)
    self.oupt = T.nn.Linear(10, 3)
    self.hid1 = T.nn.Linear(6, 10)
   
  def forward(self, x):
    z = T.tanh(self.hid1(x))
    z = T.tanh(self.hid2(z))
    z = self.oupt(z) 
    return z

This somewhat subtle idea was pointed out to me by one of the workshop attendees (“Vi”, from Denmark).

The moral of the story is that it’s better to use explicit weight and bias initialization when possible.



Out of order PyTorch layer definitions are one thing. Out of order (in the sense of broken) robots in science fiction movies is another thing.

Left: In “Red Planet” (2000), a team is sent to Mars. Their dog-like robot “AMEE” (Autonomous Mapping Exploration and Evasion) is damaged by a gamma-ray burst and attacks the crew.

Center: In “Saturn 3” (1980), the two crew members on an outpost on moon of Saturn are menaced by an out of control robot named Hector that doesn’t seem to have a head.

Right: In “Virus” (1999), the crew of a salvage ship finds an apparently deserted Russian ship with lots of robots that have been infected by an alien intelligence.


Posted in PyTorch | Leave a comment

Recap of the 2022 Machine Learning Week Conference

I presented a workshop titled “Deep Learning in Practice: A Hands-On Introduction” at the 2022 Machine Learning Week conference. The event ran June 19-24 and was in Las Vegas.

Machine Learning Week is co-located with the Predictive Analytics World conference. The events share the same conference space but have different areas of emphasis. PAW focuses on industry applications of all kinds of predictive systems. MLW focuses on more technical topics, mostly related to neural technologies and advanced classical ML techniques.


Left: The workshop room before attendees arrived. Right: The event had a small Expo event with about 20 companies represented.

In addition to my hands-on workshop at MLW, later in the week of the events, I presented a talk on advanced neural architectures at PAW.

The workshop had 24 attendees. The attendees came from all over the world, and from many different companies: a guy from the U.S. FDA, a couple of people from the Brazilian Army, a guy from an Israeli ML company, people from Denmark, Sweden, and one other Scandinavian country I can’t recall, a guy from IBM, and so on.

In the first part of the workshop, we walked through installing Anaconda Python (3.7.6), PyTorch (1.10.0), and Keras (2.8.0), with an emphasis on all the things that can go wrong. My co-presenters Leo B and Sam L did a good job of helping everyone succeed in getting a working environment up and running — no small task.

The bulk of the workshop centered around an example of multi-class classification. We examined data loaders, minus-one-plus-one vs zero-one encoding for binary predictors, Dictionary vs tuple batches, numeric predictor normalization details, neural architecture, implicit vs. explicit weight initialization, NLLLoss vs. CrossEntropyLoss, tanh vs. relu hidden activation, saving training checkpoints with exact reproducibility, global accuracy vs class-accuracy, dealing with unbalanced data, and many other topics.

The main challenge for me in presenting the workshop was that the attendees ranged from raw beginners through very advanced data scientists. But all things considered, I think the workshop went quite well.

Posted in Conferences, PAW | Leave a comment

Researchers Generate Realistic Images from Text on the Pure AI Web Site

I contributed to an article titled “Researchers Generate Realistic Images from Text” on the Pure AI web site. See https://pureai.com/articles/2022/06/01/images-from-text.aspx.

The article describes how researchers at Google have demonstrated a new technique that generates photo-realistic images from arbitrary text. The technique is implemented in a system called Imagen (“image generation”) that combines image diffusion generation with text transformer architecture language understanding.

Some of the examples are quite astonishing:


The system accepts natural language text and creates an image pixel-by-pixel.

The article notes that image generation systems like Imagen and inevitable successors raise interesting ethical questions because there are clearly ways that realistic fake images can be misused. And Imagen relies on TA text encoders trained on raw uncurated web-scale data, and so the system inherits the biases and limitations common to all large language models.


Image generation using diffusion. A low-resolution image source (left) is reduced to noise. The diffusion system learns to reconstruct the image from noise. The result is a high-quality reproduction (right).

I contributed a few opinions: The Pure AI editors spoke to Dr. James McCaffrey from Microsoft Research. McCaffrey has extensive experience working with image generation systems and natural language systems. He commented, “The images produced by the Imagen system are very impressive, and in my opinion represent a significant step forward in image generation state of the art.”

McCaffrey also noted, “The ability to create fake images that can fool humans likely means it will become increasingly important to the develop machine learning security systems that can distinguish real images from fake images.”



Photoshop guru James Fridman accepts requests from his fans to correct their photos, with clever results.


Posted in Machine Learning | Leave a comment

Computing the Incomplete Beta Function from Scratch in Excel

I made the mistake of wondering how difficult it would be to compute the regularized incomplete Beta function from scratch. After several hours of thrashing around, I knew the answer: extremely difficult.

I remember reading several references which stated that computing incomplete Beta was one of the most difficult challenges in all of numerical programming. Based on my experience, I think that statement is correct.

My initial exploration was to use Excel to calculate an example. I could compare my Excel result with the result of the betainc() function in the Python SciPy library. The point was to make sure I understood how the incomplete Beta function works before trying to implement the function in code.

The incomplete Beta function is usually written as Ix(a,b) or I(x; a, b). It would take pages to explain what Ix(a,b) is so for the purposes of this blog post, Ix(a,b) is just a function that accepts three floating point values, x, a, and b, and returns a floating point value between 0.0 and 1.0.


The SciPy library has a built-in inc_beta() function but I wanted to explore how difficult it’d be to implement from scrarch.

For example, Ix(a,b) with x = 0.5, a = 5.0, b = 3.0 returns 0.226562 — as computed via the SciPy library built-in inc_beta() function.

There are dozens of versions of equations for Ix(a,b) — all of them nutso complicated. I used one (eq. 26.5.8) called a “continued fraction” from the famous A and S Handbook of Mathematics (1964):


The incomplete beta function in the A and S Handbook.

The left part of the equation has a B(a,b) term which is the plain Beta() function. The right part between the curly braces is a continued fraction. It’s shortcut notation for:


Understanding the continued fraction syntax.

The d(i) values are computed differently for even and odd values of i.

I arbitrarily decided to stop at d4 for my experiment. After the d1 to d4 values were computed, I worked backwards and called the terms t4 to t1. The final continued fraction value was 2.7602.

To compute beta(5,3) I used the built-in Excel GAMMA() function where

beta(a,b) = [ gamma(a) * gamma(b) ] / gamma(a+b)

Putting it all together, my from-scratch version of Ix(a,b) for x=0.5, a=5.0, b=3.0 gave 0.226428 which was close to the true value (from the SciPy library) of 0.226562. I could get a better result by using more values of d beyond just d4.


Computing incomplete beta in Excel.

Anyway, after the experiments I felt confident that I could implement an end-to-end Ix(a,b) function completely from scratch. I also know it would take a lot of effort (likely a full day at the very least — perhaps as long as a week). Of course I can always use SciPy but it’s nice to know I could implement from scratch if needed.

And now that the idea of implementing incomplete Beta in code is in my head, I know that it will be nagging at me until I actually do it.

Update: After I wrote this blog post, I did in fact implement Ix(a,b) from scratch. See jamesmccaffrey.wordpress.com/2022/06/29/computing-the-regularized-incomplete-beta-function-from-scratch-using-csharp/ .



There are pros and cons to do-it-yourself implementation of complex software functions like the incomplete beta function from scratch. And the same applies to do-it-yourself car repair.


Posted in Miscellaneous | Leave a comment

Runs Testing Using C# Simulation in Visual Studio Magazine

I wrote an article titled “Runs Testing Using C# Simulation” in the June 2022 edition of Microsoft Visual Studio Magazine. See https://visualstudiomagazine.com/articles/2022/06/17/runs-testing.aspx.

Suppose you observe a sequence of people entering a store and you record the color of shirt worn by each person. You see this sequence of 24 shirts:

0, 0, 3, 3, 2, 1, 1, 2, 0, 0, 3, 3, 2, 0, 0, 1, 1, 2, 2, 2, 2, 0, 1, 1

where 0 = white shirt, 1 = yellow, 2 = red, 3 = black. You’re interested in knowing if the pattern is random or not. The pattern doesn’t look very random, but how can you quantify that observation? This is an example of a runs test.

A run is a set of consecutive data points that are the same. The sequence above has 13 runs: (0 0), (3 3), (2), (1 1), (2), (0 0), (3 3), (2), (0 0), (1 1), (2 2 2 2), (0), (1 1).

I present a demo program that estimates the probability of seeing fewer than 13 runs in a random sequence in two different ways. The first technique uses the fact that over many iterations the number of runs is approximately Normal (bell-shaped curve) distributed. The statistical probability estimate is 0.0019 — only about 2 out of 1,000 — possible, but very unlikely.

The second technique uses a raw counting technique. In the 1,000,000 iterations there were 13 or fewer runs in a random sequence only 4,888 + 1,370 313 + 59 + 13 + 1 = 6,644 times which is an estimated probability of 0.0066 — only about 7 out of 1,000. Again, possible, but very unlikely.

It’s important to remember that the results of a runs test on a sequence can only suggest that a sequence is probably not generated randomly. Therefore, conclusions should be conservative and similar to, “It is unlikely the sequence was generated by a random process, not “The pattern isn’t random”.



Fractal artwork combines randomness with patterns. Here are two wonderful illustrations by artist Jorge Abalo (“batjorge”) that represent alien landscapes and life. Left: “Deep Bonds”. Right: “Tryptamine Trance”.


Posted in Miscellaneous | Leave a comment

Computing the Regularized Incomplete Beta Function From Scratch Using C#

The regularized incomplete beta function is very strange. One form of the function looks like I(x=0.5, a=5.0, b=3.0) = 0.2265. Another form is Ix(a,b). The x value is between 0 and 1, and a and b are positive. The result is a value between 0.0 and 1.0.

The Ix(a,b) function has no obvious intuitive meaning and isn’t used directly. But Ix(a,b) is used as a helper function in many numerical functions. For example, you can use Ix(a,b) to find the area under the F-distribution for an ANOVA problem, where a and b are derived from df values and x is derived from the computed F statistic.

One weekend I sat down and decided to implement Ix(a,b) from scratch using the C# language. This was primarily for mental exercise and entertainment. (Yes, my idea of entertainment is probably like yours if you’re reading this blog post).


The continued fractions definition of Ix(a,b)

Anyway, after about six hours I had a working demo up and running.

It was crazy complicated both conceptually and from an engineering perspective. It would take many pages to explain fully, so instead, here are a few key ideas of my from-scratch implementation:

* I used the continued fractions definition (eq. 26.5.8) of Ix(a,b) from the A_and_S Handbook.
* I computed the continued fraction forward with a default of 200 terms rather than use Lenz’s algorithm.
* I used Log(Gamma(a,b)) to compute Beta(a,b) to avoid overflow.
* I implemented a LogGamma() function using the Lanczos g=5, n=7 approximation.
* I didn’t add any significant error-checking: Ix(a,b) can blow up in many ways.

I tested my from-scratch version by comparing its results with results from the Python language betainc() function in the SciPy library.

Even though the Ix(a,b) function is extraordinarily tricky, the implementation was surprisingly short. (But adding error checks would at least quadruple the code – probably more like 10x more code).

All in all it was a really fun and interesting exploration and I learned a lot.



I’ve always been fascinated and entertained by coin operated arcade games. Here are three views of the Rock-Ola World Series game, manufactured in 1934. It’s completely mechanical and quite amazing. It’s like a pinball machine. The game tracks balls and strikes, number of outs, runners advance, and a score is kept. In some ways, the device is a forerunner of the first mechanical computers.


Demo code. Replace “lt”, “gt”, “lte”, “gte” with Boolean operator symbols.

using System;

namespace IncBeta
{
  internal class IncBetaProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin IncBeta() from scratch demo \n");

      double a = 5.0; double b = 3.0; double x = 0.5;
      Console.WriteLine("\nI(x=0.5, a=5.0, b=3.0) via SciPy =" +
        " 0.22656250");
      double ib = IncBeta(x, a, b);
      Console.WriteLine("I(x=0.5, a=5.0, b=3.0) scratch  =  " +
        ib.ToString("F8"));

      a = 24.0; b = 36.0; x = 0.2;
      Console.WriteLine("\nI(x=0.2, a=24.0, b=36.0) via SciPy =" +
        " 0.00022272");
      ib = IncBeta(x, a, b);
      Console.WriteLine("I(x=0.2, a=24.0, b=36.0) scratch  =  " +
        ib.ToString("F8"));

      a = 60.0; b = 60.0; x = 0.7;
      Console.WriteLine("\nI(x=0.7, a=60.0, b=60.0) via SciPy =" +
        " 0.999997499205322");
      ib = IncBeta(x, a, b);
      Console.WriteLine("I(x=0.7, a=60.0, b=60.0) scratch  =  " +
        ib);

      Console.WriteLine("\nEnd demo ");
      Console.ReadLine();
    } // Main()

    static double ContFraction(double x, double a, double b,
      int maxTerms = 200)
    {
      // 1. pre-compute d values
      double[] d = new double[maxTerms];  // d[0] not used
      
      int end = maxTerms / 2;  // 10
      for (int m = 0; m "lt" end; ++m)  // [0,9]
      {
        int i = 2 * m;  // [0,18]  // even
        int j = i + 1;  // [1,19]  // odd
        d[i] = (m * (b - m) * x) / ((a + 2 * m - 1) *
          (a + 2 * m));
        d[j] = -1 * ((a + m) * (a + b + m) * x) / ((a + 2 * m) *
          (a + 2 * m + 1));
      }

      // 2. work backwards
      double[] t = new double[maxTerms];  // t[0] not used
      //t[4] = 1 + d[4] / 1;
      //t[3] = 1 + d[3] / t[4];
      //t[2] = 1 + d[2] / t[3];
      //t[1] = 1 + d[1] / t[2];
      //double cf = 1 / t[1];

      t[maxTerms - 1] = 1 + d[maxTerms - 1];
      for (int j = maxTerms - 2; j "gte" 1; --j)
        t[j] = 1 + d[j] / t[j + 1];

      // ShowVector(t);
      // Console.ReadLine();

      double cf = 1 / t[1];
      return cf;
    }

    static double IncBeta(double x, double a, double b)
    {
      // double top = Math.Pow(x, a) * Math.Pow((1 - x), b);
      // double bot = a * Beta(a, b);
      double cf = ContFraction(x, a, b);

      double logTop = (a * Math.Log(x)) + (b * Math.Log(1 - x));
      double logBot = Math.Log(a) + LogBeta(a, b);
      double logLeft = logTop - logBot;
      double left = Math.Exp(logLeft);

      return left * cf;  // should be [0.0, 1.0]
    }

    static double LogGamma(double z)
    {
      // plain Gamma() can easily overflow
      // Lanczos approximation g=5, n=7
      double[] coef = new double[7] { 1.000000000190015,
        76.18009172947146, -86.50532032941677,
        24.01409824083091, -1.231739572450155,
        0.1208650973866179e-2, -0.5395239384953e-5 };

      double LogSqrtTwoPi = 0.91893853320467274178;

      if (z "lt" 0.5) // Gamma(z) = Pi / (Sin(Pi*z))* Gamma(1-z))
        return Math.Log(Math.PI / Math.Sin(Math.PI * z)) -
          LogGamma(1.0 - z);

      double zz = z - 1.0;
      double b = zz + 5.5; // g + 0.5
      double sum = coef[0];

      for (int i = 1; i "lt" coef.Length; ++i)
        sum += coef[i] / (zz + i);

      return (LogSqrtTwoPi + Math.Log(sum) - b) +
        (Math.Log(b) * (zz + 0.5));
    }

    static double LogBeta(double a, double b)
    {
      // plain Beta(a,b) = (Gamma(a) * Gamma(b)) / (Gamma(a+b))
      // but plain Gamma() can easily overflow, so:
      // Log( Beta(a,b) ) 
      //   = Log( (Gamma(a) * Gamma(b)) / (Gamma(a+b) ) 
      //   = LogGamma(a) + LogGamma(b) - LogGamma(a+b)
      double logbeta = LogGamma(a) + LogGamma(b) -
        LogGamma(a + b);
      return logbeta;
    }

    //static double Beta(double a, double b)
    //{
    //  return Math.Exp(LogBeta(a, b));
    //}

    static void ShowVector(double[] v)
    {
      for (int i = 0; i "lt" v.Length; ++i)
        Console.Write(v[i].ToString("F4") + " ");
      Console.WriteLine("");
    }

  } // Program

} // ns
Posted in Miscellaneous | Leave a comment

Generating the mth Lexicographical Element of a Combination Using the Combinadic

Suppose you are working with mathematical (n=5, k=3) combinations. There are Choose(5,3) = 5! / (3! * (5-3)!) = 120 / 12 = 10 different combination elements. In lexicographical order they are:

[0] 0 1 2
[1] 0 1 3
[2] 0 1 4
[3] 0 2 3
[4] 0 2 4
[5] 0 3 4
[6] 1 2 3
[7] 1 2 4
[8] 1 3 4
[9] 2 3 4

Note: I use 0-based indexing.

Several years ago I learned about a neat technique from Laci Lovasz of Microsoft Research. He suggested that the combinadic of a number can be used to generate the mth (n,k) combination. He did not use the term “combinadic” so I coined the term because the concept is similar to the factoradic of a number.

I published an article in the old Microsoft MSDN Magazine but MSDN has been vaporized into Internet ether. I summarize that old article in this blog post.

If you want just element m = [4] you could start with (0 1 2) and assuming you have a Successor() function, loop 4 times to get to (0 2 4). It turns out that writing a Successor() function for combinations is a bit tricky but feasible. However, because the number of combinations for order n and subset size k become astronomically large for even moderate sized values of n and k, the iteration approach won’t work except for small values of m. For example, if you have (n=200, k=10) there are Choose(200, 10) = 22,451,004,309,013,280 possible combination elements

It turns out that there’s a clever way to directly generate combination element [m] using an old math idea I call the combinadic. The essence of the algorithm is to take m, compute its “combinadic”, and then use that to compute the combination.

Combinadics

You can think of a combinadic as an alternate representation of an integer. Consider the integer 859. It can be represented as

859 = (8 * 100) + (5 * 10) + (9 * 1)

Or another way of looking at it is as based on a fixed radix (base) of powers of 10:

859 = (8 * 10^2) + (5 * 10^1) + (9 * 10^0)

In short, any number can be uniquely represented as a linear sum of powers of 10. The combinadic of an integer is its representation based on a variable base corresponding to the values of the Choose(n,k) function. For example if (n=7, k=4) then the integer 27 can be represented as

27 = Choose(6,4) + Choose(5,3) + Choose(2,2) + Choose(1,1)
   = 15 + 10 + 1 + 1

With (n=7, k=4), any number m between 0 and 34 (the total number of combination elements for n and k) can be uniquely represented as:

m = Choose(c1,4) + Choose(c2,3) + Choose(c3,2) + Choose(c4,1)

where n less-than c1 less-than c2 less-than c3 less-than c4. Notice that n is analogous to the base because all combinadic digits are between 0 and n-1 (just like all digits in ordinary base 10 are between 0 and 9). The k value determines the number of terms in the combinadic. The combinadic of a number can be calculated fairly quickly.

Here’s an example of how a combinadic is calculated. Suppose you are working with (n=7, k=4) combinations, and m = 8. You want the combinadic of 8 because, as it turns out, the combinadic can be converted to combination element [8].

The combinadic of 8 will have the form:

8 = Choose(c1,4) + Choose(c2,3) + Choose(c3,2) + Choose(c4,1)

The first step is to determine the value of c1. We try c1 = 6 (the largest number less than n = 7) and get Choose(6,4) = 15, which is too large because we’re over 8. Next, we try c1 = 5 and get Choose(5,4) = 5, which is less than 8, so bingo, c1 = 5.

At this point we have used up 5 of the original number m=8 so we have 3 left to account for. To determine the value of c2, we try 4 (the largest number less than the 5 we got for c1), but get Choose(4,3) = 4, which is barely too large. Working down we get to c2 = 3 and Choose(3,3) = 1, so c2 = 3.

We used up 1 of the remaining 3 we had to account for, so we have 2 left to consume. Using the same ideas we’ll get c3 = 2 with Choose(2,2) = 1, so we have 1 left to account for. And then we’ll find that c4 = 1 because Choose(1,1) = 1. Putting our four c values together we conclude that the combinadic of m=8 for (n=7, k=4) combinations is ( 5 3 2 1 ).

Duals

Before I explain how to use the combinadic of a number to determine the mth lexicographical element of a combination, I need to explain the idea of the dual of each lexicographic index. Suppose (n=7, k=4). There are Choose(7,4) = 35 combination elements, indexed from 0 to 34. The dual indexes are the ones on opposite ends so to speak: indexes 0 and 34 are duals, indexes 1 and 33 are duals, indexes 2 and 32 are duals, and so forth. Notice that each pair of dual indexes sum to 34, so if you know any index it is easy to find its dual.

Now, continuing the first example above for the number m=27 with (n=7, k=4), suppose you are able to find the combinadic of 27 and get ( 6 5 2 1 ). Now suppose you subtract each digit in the combinadic from n-1 = 6 and get ( 0 1 4 5 ). Amazingly, this gives you the combination element [7], the dual index of 27. Putting these ideas together you have an elegant algorithm to determine an arbitrarily specified combination element for given n and k values. To find the combination element for index m, first find its dual and call it x. Next, find the combinadic of x. Then subtract each digit of the combinadic of x from n-1 and the result is the mth lexicographic combination element.

The table below shows the relationships among m, the dual of m, combination element [m], the combinadic of m, and (n-1) – ci for (n=5, k=3).

m dual(m) Element(m) combinadic(m) (n-1) - ci
==============================================
[0]  9    { 0 1 2 }   ( 2 1 0 )     ( 2 3 4 )
[1]  8    { 0 1 3 }   ( 3 1 0 )     ( 1 3 4 )
[2]  7    { 0 1 4 }   ( 3 2 0 )     ( 1 2 4 )
[3]  6    { 0 2 3 }   ( 3 2 1 )     ( 1 2 3 )
[4]  5    { 0 2 4 }   ( 4 1 0 )     ( 0 3 4 )
[5]  4    { 0 3 4 }   ( 4 2 0 )     ( 0 2 4 )
[6]  3    { 1 2 3 }   ( 4 2 1 )     ( 0 2 3 )
[7]  2    { 1 2 4 }   ( 4 3 0 )     ( 0 1 4 )
[8]  1    { 1 3 4 }   ( 4 3 1 )     ( 0 1 3 )
[9]  0    { 2 3 4 }   ( 4 3 2 )     ( 0 1 2 )

Implementation

Here is C# code that implements the algorithm, where a combination is represented as an array of int values. The implementation uses the BigInteger type because Choose() results can be huge.

Replace “lt”, “gt”, “lte”, “gte” with Boolean operator symbols (my lame blog editor chokes on symbols).

using System;
using System.Numerics;

namespace Combinations // classic style template
{
  internal class CombinationsProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin combinadic demo ");

      BigInteger numCombs = Choose(100, 10);
      Console.WriteLine("Number (n=100, k=10) combinations: ");
      Console.WriteLine(numCombs.ToString("N0"));

      int n = 100;
      int k = 10;
      BigInteger m = BigInteger.Parse("9999999999");
      int[] c = Element(m, n, k);
      Console.WriteLine("\nCombination m=[" + m + "]" +
        " of C(n=100,k=10): ");
      ShowComb(c);

      Console.WriteLine("\nEnd combinadic demo \n");
      Console.ReadLine();
    } // Main

    static void ShowComb(int[] comb)
    {
      int n = comb.Length;
      for (int i = 0; i "lt" n; ++i)
        Console.Write(comb[i] + " ");
      Console.WriteLine("");
    }

    static int[] Element(BigInteger m, int n, int k)
    {
      // compute element [m] using the combinadic
      BigInteger maxM = Choose(n, k) - 1;

      if (m "gt" maxM)
        throw new Exception("m value too large in Element");

      int[] ans = new int[k];

      int a = n;
      int b = k;
      BigInteger x = maxM - m; // x is the "dual" of m

      for (int i = 0; i "lt" k; ++i)
      {
        ans[i] = LargestV(a, b, x); // see helper below    
        x = x - Choose(ans[i], b);
        a = ans[i];
        b = b - 1;
      }

      for (int i = 0; i "lt" k; ++i)
        ans[i] = (n - 1) - ans[i];

      return ans;
    }

    static int LargestV(int a, int b, BigInteger x)
    {
      // helper for Element()
      int v = a - 1;
      while (Choose(v, b) "gt" x)
        --v;
      return v;
    }

    static BigInteger Choose(int n, int k)
    {
      // number combinations
      if (n "lt" 0 || k "lt" 0)
        throw new Exception("Negative arg in Choose()");
      if (n "lt" k) return 0; // special
      if (n == k) return 1; // short-circuit

      int delta, iMax;

      if (k "lt" n - k) // ex: Choose(100,3)
      {
        delta = n - k; iMax = k;
      }
      else           // ex: Choose(100,97)
      {
        delta = k; iMax = n - k;
      }

      BigInteger ans = delta + 1;
      for (int i = 2; i "lte" iMax; ++i)
        ans = (ans * (delta + i)) / i;

      return ans;
    }

  } // Program

} // ns
Posted in Miscellaneous | Leave a comment

Generating the kth Lexicographical Element of a Permutation Using the Factoradic

Suppose you are working with mathematical permutations of order n = 3. There are 3! = 6 different permutation elements. In lexicographical order they are:

[0] 0 1 2
[1] 0 2 1
[2] 1 0 2
[3] 1 2 0
[4] 2 0 1
[5] 2 1 0

Note: I use 0-based indexing.

Several years ago I learned about a neat technique from Peter Cameron of Queen Mary, University of London. He suggested that the factoradic of a number can be used to generate the kth permutation of order n. He was not sure where the term “factoradic” originated so I use it too. I published an article in the old Microsoft MSDN Magazine but MSDN has been vaporized into Internet ether. I summarize that old article in this blog post.

If you want just element [4] you could start with (0 1 2) and assuming you have a Successor() function, loop 4 times to get to (2 0 1). It turns out that writing a Successor() function for permutations is a bit tricky but feasible. However, because the number of permutations for order n become astronomically large for even moderate sized values of n, the iteration approach won’t work except for small values of n.

It turns out that there’s a clever way to directly generate permutation element [k] using an old idea called the factoradic. The essence of the algorithm is to take k, compute its “factoradic”, and then use that to compute the permutation.

You can think of a factoradic as an alternate representation of an integer. Consider the integer 859. It can be represented as

859 = (8 * 100) + (5 * 10) + (9 * 1)

Or another way of looking at it is as based on a fixed radix (base) of powers of 10:

859 = (8 * 10^2) + (5 * 10^1) + (9 * 10^0)

In short, any number can be uniquely represented as a linear combination of powers of 10. The factoradic of an integer is its representation based on a variable base corresponding to the values of n factorial. It turns out that any integer k can be uniquely represented in the form

k = (a0 * 1!) + (a1 * 2!) + (a2 * 3!) + (a3 * 4!) + . . .
  = (a0 * 1) + (a1 * 2) + (a2 * 6) + (a3 * 24) + . . .

For example, the integer 859 can be represented as

 
 = (1 * 1!) + (0 * 2!) + (3 * 3!) + (0 * 4!) + (1 * 5!) + (1 * 6!)
 = (1 * 1) + (3 * 6) + (1 * 120) + (1 * 720)

So 859 in factoradic form is { 1 1 0 3 0 1 } where the right-most digit is the value of the 1!s. It will be useful to append a trailing 0 onto the right end of all factoradics so { 1 1 0 3 0 1 0 } is the final form for 859.

Furthermore, it turns out that there is a one-to-one mapping between the factoradic of an integer k and the kth permutation of order n, meaning that each factoradic uniquely determines a permutation. To illustrate this, the following table shows the values of k, the factoradic of k, and the kth permutation for order n=4.

 k     factoradic    permutation
--------------------------------
[ 0]   { 0 0 0 0 }   ( 0 1 2 3 )
[ 1]   { 0 0 1 0 }   ( 0 1 3 2 )
[ 2]   { 0 1 0 0 }   ( 0 2 1 3 )
[ 3]   { 0 1 1 0 }   ( 0 2 3 1 )
[ 4]   { 0 2 0 0 }   ( 0 3 1 2 )
[ 5]   { 0 2 1 0 }   ( 0 3 2 1 )
[ 6]   { 1 0 0 0 }   ( 1 0 2 3 )
[ 7]   { 1 0 1 0 }   ( 1 0 3 2 )
[ 8]   { 1 1 0 0 }   ( 1 2 0 3 )
[ 9]   { 1 1 1 0 }   ( 1 2 3 0 )
[10]   { 1 2 0 0 }   ( 1 3 0 2 )
[11]   { 1 2 1 0 }   ( 1 3 2 0 )
[12]   { 2 0 0 0 }   ( 2 0 1 3 )
[13]   { 2 0 1 0 }   ( 2 0 3 1 )
[14]   { 2 1 0 0 }   ( 2 1 0 3 )
[15]   { 2 1 1 0 }   ( 2 1 3 0 )
[16]   { 2 2 0 0 }   ( 2 3 0 1 )
[17]   { 2 2 1 0 }   ( 2 3 1 0 )
[18]   { 3 0 0 0 }   ( 3 0 1 2 )
[19]   { 3 0 1 0 }   ( 3 0 2 1 )
[20]   { 3 1 0 0 }   ( 3 1 0 2 )
[21]   { 3 1 1 0 }   ( 3 1 2 0 )
[22]   { 3 2 0 0 }   ( 3 2 0 1 )
[23]   { 3 2 1 0 }   ( 3 2 1 0 )

For example, for k = 5, the factoradic is { 0 2 1 0 } = (0 * 3!) + (2 * 2!) + (1 * 1!) + 0 and the [5] permutation of order 4 is ( 0 3 2 1 ).

The clever and efficient way to compute the kth permutation of order n is to first find the factoradic of k and then to generate the corresponding permutation from the factoradic.

The trickiest part of the algorithm is the computation of the permutation that corresponds to the factoradic. Here’s an example that shows how the algorithm converts the factoradic { 1 2 3 2 1 1 0 } into its corresponding permutation of order n=7 which is ( 1 3 5 4 2 6 0 ).

First create a temp[] array and copy into it the factoradic values incremented by 1:

[ 2 3 4 3 2 2 1 ]

The result[] permutation array is seeded with a 1 value in the right-most cell:

[ ? ? ? ? ? ? 1 ] 

Now starting with the second value from the right-most value of the factoradic (skipping over the right-most because it will always be 1 since it came from the padded 0 value), add it to the result[] array:

[ ? ? ? ? ? 2 1 ]

Now the algorithm scans through all the values to the right of the new value and increments by 1 all values that are greater than or equal to the new value. Continuing this process generates:

[ ? ? ? ? 2 3 1 ]
[ ? ? ? 3 2 4 1 ]
[ ? ? 4 3 2 5 1 ]
[ ? 3 5 4 2 6 1 ]
[ 2 4 6 5 3 7 1 ]

Last, the algorithm traverses the result[] array and decrements all values by 1 to put the resulting permutation in 0-based form:

( 1 3 5 4 2 6 0 )

To summarize, if you want to generate the kth permutation of order n, first compute the factoradic of k, and then use that result to compute the corresponding permutation. The example above started with k = 1,047 and then computed its factoradic = { 1 2 3 2 1 1 0 } and then computed the permutation ( 1 3 5 4 2 6 0 ). So the permutation [1047] of order 7 is ( 1 3 5 4 2 6 0 ).

Here is C# code that implements the algorithm, where a permutation is represented as an array of int values. The implementation uses the BigInteger type because Factorial() can be huge.

Replace “lt”, “gt”, “lte”, “gte” with Boolean operator symbols (my lame blog editor chokes on symbols).

using System; 
using System.Numerics;
namespace Permutations // classic style template
{
  internal class PermutationsProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin factoradic demo ");

      int n = 7;  // order
      BigInteger k = BigInteger.Parse("1047");  // k could be huge
      int[] perm = Element(k, n);
      Console.WriteLine("\nPermutation[" + k + "] of order 7: ");
      ShowPerm(perm);

      Console.WriteLine("\nEnd permutations demo ");
      Console.ReadLine();
    } // Main

    static void ShowPerm(int[] perm)
    {
      int order = perm.Length;
      for (int i = 0; i "lt" order; ++i)
        Console.Write(perm[i] + " ");
      Console.WriteLine("");
    }

    static int[] Element(BigInteger k, int order)
    {
      if (k "gte" Factorial(order))
        throw new Exception("k too large in Element");
      int[] result = new int[order];

      int[] factoradic = new int[order]; // factoradic of k
      for (int j = 1; j "lte" order; ++j)  // note: skip [0]
      {
        factoradic[order - j] = (int)(k % j);  // always an int
        k /= j;
      }

      for (int i = 0; i "lt" order; ++i)
        ++factoradic[i];

      result[order - 1] = 1; // last value set to 1
      for (int i = order - 2; i "gte" 0; --i)
      {
        result[i] = factoradic[i];  // copy factoradic
        for (int j = i + 1; j "lt" order; ++j) 
        {
          if (result[j] "gte" result[i]) 
            ++result[j];
        }
      }

      for (int i = 0; i "lt" order; ++i) // make 0-based
        --result[i];

      return result;
    }

    static BigInteger Factorial(int n)
    {
      if (n == 0 || n == 1)
        return BigInteger.One;

      BigInteger ans = BigInteger.Parse("1");  // alternative
      for (int i = 1; i "lte" n; ++i)
        ans *= i;
      return ans;
    }

  } // Program
} // ns
Posted in Miscellaneous | 1 Comment

Lightweight Derangements Using C#

I was doing some work with permutations and I remembered “derangements”, an old topic from my college days. A derangement is a kind of permutation where all elements are out of order. For example, all 6 zero-based permutations of order 3 are:

0 1 2
0 2 1
1 0 2
1 2 0
2 0 1
2 1 0

Of these six permutations, the only derangements are:

1 2 0
2 0 1

For example, permutation (0 2 1) isn’t a derangement because element 0 is in it’s original position. And permutation (2 1 0) isn’t a derangement because element 1 is in it’s original position. And so on. Note: I use 0-based indexing (common in computer science) even though 1-based indexing is more common in mathematics.

In general, about one-third of permutations of order n are also derangements.

I implemented a demo. Given a particular derangement, to find the next derangement, I used this mini-algorithm:

loop
  get next permutation
  if it's a derangement return it
end-loop
if none found
  return null
end-if

The algorithm relies on the fact that all derangements are permutations, and the fact that it’s possible to write a function that returns the next permutation to a given permutation.

The number of permutations for order n is factorial(n), usually written as n! as in 3! = 3 * 2 * 1 = 6. The number of derangements for order n doesn’t have a name but it’s usually written as !n (the exclamation is in front). Computing !n is a fascinating topic. I used a recurrence relation.

My demo to find the next derangement is crude in the sense that it generates the next permutation and if that permutation is a derangement then the next derangement has been found, but if the next permutation isn’t a derangement, try the next permutation. This approach isn’t horrible because about one-third of permutations are derangements and so there’s not too much wasted effort. I discovered a neat algorithm to directly generate a next derangement in a 2004 paper “Constant Time Generation of Derangements” by J. Korsch and P. LaFolette. The algorithm can be used to find all derangements but doesn’t give the next derangement in lexicographical order. The algorithm is moderately complex but I’ll implement a demo when I get some free time.

Good fun.



A mathematical derangement is a type of permutation. A permutation is a type of mutation. In early science fiction movies, a common theme is giant insects mutated by radiation. Left: “Them!” (1954) featured giants ants rampaging through the desert, and later, Los Angeles. Center: “Beginning of the End” (1957) tells the story of giant locusts in the Chicago area. Right: In “Monster from Green Hell” (1957), a wasp is sent into space, irradiated, and lands back in Africa.


Demo code. Note: Some of this code is extremely tricky and I haven’t thoroughly tested it so there are likely some edge cases that are incorrect. Replace “lt”, “gt”, “lte”, “gte” with Boolean operator symbols.

using System;  // .NET 6.0
using System.Numerics;

namespace Derangements // classic style template
{
  internal class DerangementsProgram
  {
    static void Main(string[] args)
    {
      Console.WriteLine("\nBegin lightweight derangements demo ");

      int[] d = MakeDerange(4);
      Console.WriteLine("\nFirst derangement order 4 is: ");
      ShowDerange(d);

      BigInteger numDerange = NumDerange(4);
      Console.WriteLine("\nNumber of derangements order 4 is: ");
      Console.WriteLine(numDerange);

      numDerange = NumDerange(20);
      Console.WriteLine("\nNumber of derangements order 20 is: ");
      Console.WriteLine(numDerange.ToString("N0"));

      Console.WriteLine("\nAll derangements order 4: ");
      while (d != null)
      {
        ShowDerange(d);
        d = NextDerange(d);
      }

      Console.WriteLine("\nEnd derangements demo ");
      Console.ReadLine();
    } // Main

    static int[] MakePerm(int order)
    {
      int[] result = new int[order];
      for (int i = 0; i "lt" order; i++)
        result[i] = i;
      return result;
    }

    static int[] MakeDerange(int order)
    {
      int[] start = new int[order];
      for (int i = 0; i "lt" order; ++i)
        start[i] = i;

      int[] result = NextDerange(start);
      return result;
    }

    static void ShowDerange(int[] derange)
    {
      int order = derange.Length;
      for (int i = 0; i "lt" order; ++i)
        Console.Write(derange[i] + " ");
      Console.WriteLine("");
    }

    static bool IsLast(int[] perm)
    {
      // is perm(6) like [5,4,3,2,1,0] ?
      int order = perm.Length;
      if (perm[0] != order - 1) return false;
      if (perm[order - 1] != 0) return false;
      for (int i = 0; i "lt" order - 1; ++i)
      {
        if (perm[i] "lt" perm[i + 1])
          return false;
      }
      return true;
    }

    static bool IsDerange(int[] perm)
    {
      // is perm[] also a derangement?
      if (perm == null)
        return false;

      int order = perm.Length;
      for (int i = 0; i "lt" order; ++i)
        if (perm[i] == i)
          return false;
      return true;  // passed all checks
    }

    static int[] NextDerange(int[] derange)
    {
      int order = derange.Length;
      int[] result = new int[order];
      for (int i = 0; i "lt" order; ++i)  // make a copy
        result[i] = derange[i];

      while (result != null)
      {
        result = NextPerm(result);  // could be null

        if (IsDerange(result) == true)
          return result;
      }

      return null;  // didn't find a next
    }

    static int[] NextPerm(int[] perm)
    {
      int order = perm.Length;  // [0,1,2,3,4] is order 5

      // is perm [4,3,2,1,0] and so no successor?
      if (IsLast(perm) == true)
        return null;

      int[] result = new int[order];
      for (int k = 0; k "lt" order; ++k) 
        result[k] = perm[k];

      int left, right;

      left = order - 2;  // Find left value 
      while ((result[left] "gt" result[left + 1]) &&
        (left "gte" 1))
        --left;

      right = order - 1; 
      while (result[left] "gt" result[right])
        --right;

      int tmp = result[left];  // swap 
      result[left] = result[right];
      result[right] = tmp;

      int i = left + 1;  // order the tail
      int j = order - 1;
      while (i "lt" j)
      {
        tmp = result[i];
        result[i++] = result[j];
        result[j--] = tmp;
      }

      return result;
    }

    static BigInteger NumDerange(int n)
    {
      if (n == 0)
        return BigInteger.One;
      else if (n == 1)
        return BigInteger.Zero;
      else if (n == 2)
        return BigInteger.One;

      return (n - 1) * (NumDerange(n - 1) + NumDerange(n - 2));
    }
    
  } // Program

} // ns
Posted in Miscellaneous | Leave a comment