<< problem 306 - Paper-strip Game | An amazing Prime-generating Automaton - problem 308 >> |
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 toecho "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 an Intel® Core™ i7-2600K CPU @ 3.40GHz.
(compiled for x86_64 / Linux, GCC flags: -O3 -march=native -fno-exceptions -fno-rtti -std=gnu++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)
Code in various languages:
Python github.com/Meng-Gen/ProjectEuler/blob/master/307.py (written by Meng-Gen Tsai)
C++ github.com/roosephu/project-euler/blob/master/307.cpp (written by Yuping Luo)
Java github.com/thrap/project-euler/blob/master/src/Java/Problem307.java (written by Magnus Solheim Thrap)
Mathematica github.com/frrad/project-euler/blob/master/mathematica/Problem307.m (written by Frederick Robinson)
Mathematica github.com/steve98654/ProjectEuler/blob/master/307.nb
Perl github.com/shlomif/project-euler/blob/master/project-euler/307/euler-307-v1.pl (written by Shlomi Fish)
Those links are just an unordered selection of source code I found with a semi-automatic search script on Google/Bing/GitHub/whatever.
You will probably stumble upon better solutions when searching on your own. Maybe not all linked resources produce the correct result and/or exceed time/memory limits.
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 | |
black | problems are solved but access to the solution is blocked for a few days until the next problem is published | |
[new] | the flashing problem is the one I solved most recently |
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 |
576 | 577 | 578 | 579 | 580 | 581 | 582 | 583 | 584 | 585 | 586 | 587 | 588 | 589 | 590 | 591 | 592 | 593 | 594 | 595 | 596 | 597 | 598 | 599 | 600 |
601 | 602 | 603 | 604 | 605 | 606 | 607 | 608 | 609 | 610 | 611 | 612 | 613 | 614 | 615 | 616 | 617 | 618 | 619 | 620 | 621 | 622 | 623 | 624 | 625 |
626 | 627 | 628 | 629 |
I scored 13526 points (out of 15700 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 | An amazing Prime-generating Automaton - problem 308 >> |