<< problem 306 - Paper-strip Game | Integer Ladders - problem 309 >> |

# Problem 307: Chip Defects

(see projecteuler.net/problem=307)

k defects are randomly distributed amongst n integrated-circuit chips produced by a factory

(any number of defects may be found on a chip and each defect is independent of the other defects).

Let p(k,n) represent the probability that there is a chip with at least 3 defects.

For instance p(3,7) approx 0.0204081633.

Find p(20 000, 1 000 000) and give your answer rounded to 10 decimal places in the form 0.abcdefghij

# My Algorithm

It took quite some effort to finally solve this problem. First I had a minor flaw in my formulas - but most time was spent

working around precision issues. The current code does *not* work with Visual C++ because I need extended floating-point precision (`long double`

).

The main idea is that instead of counting chips with 3+ defect I count the opposite: chips with 0, 1 or 2 defects.

The probability of having exactly x chips with 2 defects each is

dfrac{k!}{(k-x)!} * dfrac{p! / ((p-2x)! x!)}{2^x} * k^{-p}

Summing these values for all possible x (→ number of chips with 2 defects) will give the percentage of working chips, so

p(k,n) = 1 - sum probabilities(x,k,n)

The involved numbers are huge and exceed the range of C++'s floating point data types.

Therefore I perform all calculations in "log space":

a = e^{ln(a)}

a * b = e^{ln(a) + ln(b)}

a ^ b = e^{ln(a) * b}

I wrote a function `logFactorial(n)`

that returns ln{n!}.

A similar function `logFactorial(n, onlyTopValues)`

returns ln{dfrac{n!}{(n - onlyTopValues)!}} which is needed twice

and much faster and more precise to calculate than calling `logFactorial(n)`

twice.

You will find my old code of a Monte-Carlo simulation as well (see `monteCarlo`

).

It helped me fixing a major bug in my mathematical formula, which was way off.

As always, it's impossible to achieve the required accuracy (10 digits) with a basic Monte-Carlo simulations.

## Note

As mentionend before, a simple `double`

has some precision issues and the last digit was off by 3.

GCC supports `long double`

(which has 80 instead of 64 bits) and finds the correct result.

# Interactive test

You can submit your own input to my program and it will be instantly processed at my server:

This is equivalent to`echo "3 7" | ./307`

Output:

*Note:* the original problem's input `20000 1000000`

__cannot__ be entered

because just copying results is a soft skill reserved for idiots.

*(this interactive test is still under development, computations will be aborted after one second)*

# My code

… was written in C++11 and can be compiled with G++, Clang++. You can download it, too.

#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.

# Benchmark

The correct solution to the original Project Euler problem was found in 0.24 seconds on a Intel® Core™ i7-2600K CPU @ 3.40GHz.

(compiled for x86_64 / Linux, GCC flags: `-O3 -march=native -fno-exceptions -fno-rtti -std=c++11 -DORIGINAL`

)

See here for a comparison of all solutions.

**Note:** interactive tests run on a weaker (=slower) computer. Some interactive tests are compiled without `-DORIGINAL`

.

# Changelog

August 12, 2017 submitted solution

August 12, 2017 added comments

# Difficulty

Project Euler ranks this problem at **35%** (out of 100%).

# Links

projecteuler.net/thread=307 - **the** best forum on the subject (*note:* you have to submit the correct solution first)

# Heatmap

*Please click on a problem's number to open my solution to that problem:*

green | solutions solve the original Project Euler problem and have a perfect score of 100% at Hackerrank, too | |

yellow | solutions score less than 100% at Hackerrank (but still solve the original problem easily) | |

gray | problems are already solved but I haven't published my solution yet | |

blue | solutions are relevant for Project Euler only: there wasn't a Hackerrank version of it (at the time I solved it) or it differed too much | |

orange | problems are solved but exceed the time limit of one minute or the memory limit of 256 MByte | |

red | problems are not solved yet but I wrote a simulation to approximate the result or verified at least the given example - usually I sketched a few ideas, too |

1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | 19 | 20 | 21 | 22 | 23 | 24 | 25 |

26 | 27 | 28 | 29 | 30 | 31 | 32 | 33 | 34 | 35 | 36 | 37 | 38 | 39 | 40 | 41 | 42 | 43 | 44 | 45 | 46 | 47 | 48 | 49 | 50 |

51 | 52 | 53 | 54 | 55 | 56 | 57 | 58 | 59 | 60 | 61 | 62 | 63 | 64 | 65 | 66 | 67 | 68 | 69 | 70 | 71 | 72 | 73 | 74 | 75 |

76 | 77 | 78 | 79 | 80 | 81 | 82 | 83 | 84 | 85 | 86 | 87 | 88 | 89 | 90 | 91 | 92 | 93 | 94 | 95 | 96 | 97 | 98 | 99 | 100 |

101 | 102 | 103 | 104 | 105 | 106 | 107 | 108 | 109 | 110 | 111 | 112 | 113 | 114 | 115 | 116 | 117 | 118 | 119 | 120 | 121 | 122 | 123 | 124 | 125 |

126 | 127 | 128 | 129 | 130 | 131 | 132 | 133 | 134 | 135 | 136 | 137 | 138 | 139 | 140 | 141 | 142 | 143 | 144 | 145 | 146 | 147 | 148 | 149 | 150 |

151 | 152 | 153 | 154 | 155 | 156 | 157 | 158 | 159 | 160 | 161 | 162 | 163 | 164 | 165 | 166 | 167 | 168 | 169 | 170 | 171 | 172 | 173 | 174 | 175 |

176 | 177 | 178 | 179 | 180 | 181 | 182 | 183 | 184 | 185 | 186 | 187 | 188 | 189 | 190 | 191 | 192 | 193 | 194 | 195 | 196 | 197 | 198 | 199 | 200 |

201 | 202 | 203 | 204 | 205 | 206 | 207 | 208 | 209 | 210 | 211 | 212 | 213 | 214 | 215 | 216 | 217 | 218 | 219 | 220 | 221 | 222 | 223 | 224 | 225 |

226 | 227 | 228 | 229 | 230 | 231 | 232 | 233 | 234 | 235 | 236 | 237 | 238 | 239 | 240 | 241 | 242 | 243 | 244 | 245 | 246 | 247 | 248 | 249 | 250 |

251 | 252 | 253 | 254 | 255 | 256 | 257 | 258 | 259 | 260 | 261 | 262 | 263 | 264 | 265 | 266 | 267 | 268 | 269 | 270 | 271 | 272 | 273 | 274 | 275 |

276 | 277 | 278 | 279 | 280 | 281 | 282 | 283 | 284 | 285 | 286 | 287 | 288 | 289 | 290 | 291 | 292 | 293 | 294 | 295 | 296 | 297 | 298 | 299 | 300 |

301 | 302 | 303 | 304 | 305 | 306 | 307 | 308 | 309 | 310 | 311 | 312 | 313 | 314 | 315 | 316 | 317 | 318 | 319 | 320 | 321 | 322 | 323 | 324 | 325 |

326 | 327 | 328 | 329 | 330 | 331 | 332 | 333 | 334 | 335 | 336 | 337 | 338 | 339 | 340 | 341 | 342 | 343 | 344 | 345 | 346 | 347 | 348 | 349 | 350 |

351 | 352 | 353 | 354 | 355 | 356 | 357 | 358 | 359 | 360 | 361 | 362 | 363 | 364 | 365 | 366 | 367 | 368 | 369 | 370 | 371 | 372 | 373 | 374 | 375 |

376 | 377 | 378 | 379 | 380 | 381 | 382 | 383 | 384 | 385 | 386 | 387 | 388 | 389 | 390 | 391 | 392 | 393 | 394 | 395 | 396 | 397 | 398 | 399 | 400 |

401 | 402 | 403 | 404 | 405 | 406 | 407 | 408 | 409 | 410 | 411 | 412 | 413 | 414 | 415 | 416 | 417 | 418 | 419 | 420 | 421 | 422 | 423 | 424 | 425 |

426 | 427 | 428 | 429 | 430 | 431 | 432 | 433 | 434 | 435 | 436 | 437 | 438 | 439 | 440 | 441 | 442 | 443 | 444 | 445 | 446 | 447 | 448 | 449 | 450 |

451 | 452 | 453 | 454 | 455 | 456 | 457 | 458 | 459 | 460 | 461 | 462 | 463 | 464 | 465 | 466 | 467 | 468 | 469 | 470 | 471 | 472 | 473 | 474 | 475 |

476 | 477 | 478 | 479 | 480 | 481 | 482 | 483 | 484 | 485 | 486 | 487 | 488 | 489 | 490 | 491 | 492 | 493 | 494 | 495 | 496 | 497 | 498 | 499 | 500 |

501 | 502 | 503 | 504 | 505 | 506 | 507 | 508 | 509 | 510 | 511 | 512 | 513 | 514 | 515 | 516 | 517 | 518 | 519 | 520 | 521 | 522 | 523 | 524 | 525 |

526 | 527 | 528 | 529 | 530 | 531 | 532 | 533 | 534 | 535 | 536 | 537 | 538 | 539 | 540 | 541 | 542 | 543 | 544 | 545 | 546 | 547 | 548 | 549 | 550 |

551 | 552 | 553 | 554 | 555 | 556 | 557 | 558 | 559 | 560 | 561 | 562 | 563 | 564 | 565 | 566 | 567 | 568 | 569 | 570 | 571 | 572 | 573 | 574 | 575 |

I scored 13,386 points (out of 15600 possible points, top rank was 17 out of ≈60000 in August 2017) at Hackerrank's Project Euler+.

My username at Project Euler is

**stephanbrumme**while it's stbrumme at Hackerrank.

Look at my progress and performance pages to get more details.

# Copyright

I hope you enjoy my code and learn something - or give me feedback how I can improve my solutions.

All of my solutions can be used for any purpose and I am in no way liable for any damages caused.

You can even remove my name and claim it's yours. But then you shall burn in hell.

The problems and most of the problems' images were created by Project Euler.

Thanks for all their endless effort !!!

<< problem 306 - Paper-strip Game | Integer Ladders - problem 309 >> |