My code

… was written in C++11 and can be compiled with G++, Clang++. You can download it, too. Or just jump to my GitHub repository.

       
#include <iostream>
       
#include <iomanip>
       
#include <cmath>
       
 
       
// ---------- Monte-Carlo simulation ----------
       
// => not used anymore, I ran it to have a crude approximation of the result

       
#include <vector>
       
 
       
// a simple pseudo-random number generator
       
// (produces the same result no matter what compiler you have - unlike rand() from math.h)

       
unsigned int myrand()
       
{
       
  static unsigned long long seed = 0;
       
  seed = 6364136223846793005ULL * seed + 1;
       
  return (unsigned int)(seed >> 30);
       
}
       
 
       
// simulate k defects on n chips, return percentage of 3+ defects on at least one chip

       
double monteCarlo(unsigned int iterations, unsigned int k, unsigned int n)
       
{
       
  // a chip is defect if 3+ defects at once
       
  const unsigned char threshold = 3;
       
 
       
  // at least one chip isn't working (3+ defects)
       
  unsigned int bad = 0;
       
  std::vector<unsigned char> defects(n);
       
  for (unsigned int i = 0; i < iterations; i++)
       
  {
       
    // reset array
       
    for (auto& x : defects)
       
      x = 0;
       
 
       
    // spread defects randomly
       
    for (unsigned int j = 0; j < k; j++)
       
    {
       
      auto id = myrand() % n;
       
      defects[id]++;
       
      // found one more iteration with at least one chips that's out of order
       
      if (defects[id] == threshold)
       
      {
       
        bad++;
       
        break;
       
      }
       
    }
       
  }
       
 
       
  // percentage of chips with 3+ defects
       
  return bad / double(iterations);
       
}
       
 
       
// ---------- explicit computation ----------

       
 
       
typedef long double Number;
       
 
       
// return log(n!)

       
Number logFactorial(unsigned int n)
       
{
       
  Number result = 0; // = log(1);
       
  for (unsigned int i = 2; i <= n; i++)
       
    result += log(Number(i));
       
  return result;
       
}
       
 
       
// return log(n! / (n - onlyTopValues)!)

       
Number logFactorial(unsigned int n, unsigned int onlyTopValues)
       
{
       
  Number result = 0;
       
  for (auto i = n - onlyTopValues + 1; i <= n; i++)
       
    result += log(Number(i));
       
  return result;
       
}
       
 
       
int main()
       
{
       
  unsigned int chips   = 1000000;
       
  unsigned int defects =   20000;
       
  std::cin >> defects >> chips;
       
 
       
  // 10^-13
       
  Number precisionThreshold = 0.0000000000001;
       
 
       
  // total number of combinations:
       
  // chips ^ defects
       
  auto combinations = log(Number(chips)) * defects;
       
 
       
  // add probabilities of not having a chip with 3+ defects
       
  Number sum = 0;
       
 
       
  // process all possible number of chips with exactly 2 defects
       
  for (unsigned int numTwoDefects = 0; numTwoDefects <= defects / 2; numTwoDefects++)
       
  {
       
    // chips with one or two defects
       
    auto affectedChips = defects - numTwoDefects;
       
    auto permutations  = logFactorial(chips, affectedChips);
       
 
       
    // count combinations of chips with one defect
       
    auto defectsFoundOnChipsWithTwo = 2 * numTwoDefects;
       
    auto count = logFactorial(defects, defectsFoundOnChipsWithTwo);
       
 
       
    // divide by numTwoDefects!
       
    count -= logFactorial(numTwoDefects);
       
 
       
    // divide by 2^numTwoDefects
       
    auto countTwoDefects = numTwoDefects * log(Number(2));
       
    count -= countTwoDefects;
       
 
       
    // multiply both
       
    auto noDefects = permutations + count;
       
 
       
    // percentage of working chips
       
    auto ratio = noDefects - combinations;
       
 
       
    // convert from "log space" to normal numbers
       
    ratio = exp(ratio);
       
    sum  += ratio;
       
 
       
    // abort early if sufficient precision
       
    if (sum > 0.01 && ratio < precisionThreshold)
       
      break;
       
  }
       
 
       
  // percentage of 3+ defects is opposite of percentage of 2- chips
       
  auto result = 1 - sum;
       
  std::cout << std::fixed << std::setprecision(10)
       
            << result << std::endl;
       
 
       
  // run repeatedly one million Monte-Carlo simulations
       
  // the first two digits will be okay
       
  //while (true)
       
  //  std::cout << monteCarlo(100000, 20000, 1000000) << std::endl;
       
 
       
  return 0;
       
}

This solution contains 25 empty lines, 32 comments and 4 preprocessor commands.