// ////////////////////////////////////////////////////////
// # Title
// Exploring Pascal's pyramid
//
// # URL
// https://projecteuler.net/problem=154
// http://euler.stephan-brumme.com/154/
//
// # Problem
// A triangular pyramid is constructed using spherical balls so that each ball rests on exactly three balls of the next lower level.
//
// ![Pascal's pyramid](p154_pyramid.gif)
//
// Then, we calculate the number of paths leading from the apex to each position:
//
// A path starts at the apex and progresses downwards to any of the three spheres directly below the current position.
//
// Consequently, the number of paths to reach a certain position is the sum of the numbers immediately above it
// (depending on the position, there are up to three numbers above it).
//
// The result is Pascal's pyramid and the numbers at each level n are the coefficients of the trinomial expansion `(x + y + z)^n`.
//
// How many coefficients in the expansion of `(x + y + z)^200000` are multiples of `10^12`?
//
// # Solved by
// Stephan Brumme
// July 2017
//
// # Algorithm
// Wikipedia told me that a value of Pascal's pyramid can be derived from Pascal's triangle ( https://en.wikipedia.org/wiki/Pascal%27s_pyramid )
// `P(level, i, j) = T(level, i) * T(i, j)`
//
// And a value of Pascal's triangle is computed with the binomial coefficient:
// `T(n, k) = dfrac{n!}{(n-k)! * k!}`
//
// Those values get huge - but I don't need to know their true value.
// All I need to know is whether `P(level, i, j)` can be factorized into at least 12x 2 and 12x 5 because `2^12 * 5^12 = 10^12`.
//
// My program stores for each number between 1 and 200000 (= ''layer'') how often it contains 2 and 5 (see ''mulPrime1'' and ''mulPrime2'').
// Repeatedly dividing the current number by 2 and 5 produces the desired result.
//
// Then it stores for each factorial between 1 and 200000 how often it contains 2 and 5 (see ''sum1'' and ''sum2'').
// Since `s^x * s^y = s^{x+y}` it boils down to summing the exponents up to the current point.
//
// The same idea make ''choose'' very simple:
// ''choose(T(n, k)) = sums[n] - (sums[n-k] + sums[k])''
// When I call ''choose'' with ''sums = sum1'' then I know how many 2s T(n,k) contains.
// And the same happens with ''sums = sum2'' for the number of 5s.
//
// Whenever I have at least 12x 2s and 12x 5s then I found a multiple of `10^12`.
//
// There a two simple optimizations:
// - if `T(level, i)` is a multiple of `10^12` then any `T(level, i) * T(i, j)` is as well
// - `T(n, k) = T(n, n - k)` ==> the pyramid is symmetric and I only need to compute its left half
//
// # Note
// There are more symmetries in the pyramid that I didn't exploit: I can see six "areas" with identical values.
// Proper use would make my program probably 3 times faster.
//
// # Hackerrank
// Hackerrank wants you to count all numbers that are multiples of `p_1^{a_1} * p_2^{a_2}` where `p_1` and `p_2` are distinct prime numbers.
// `10^12 = 2^12 * 5^12`, therefore my program finds the correct solution for the original problem when entering ''200000 2 12 5 12''.
#include
#include
typedef std::vector Exponents;
// return exponent of a prime factor of C(n,k)
// looks a bit like the logarithm:
// C(n,k) = n! / ((n-k)! * k!)
unsigned int choose(const Exponents& sums, unsigned int n, unsigned int k)
{
return sums[n] - (sums[n - k] + sums[k]);
}
int main()
{
unsigned int layer = 200000;
// 10^12 = 2^12 * 5^12
unsigned int prime1 = 2;
unsigned int exponent1 = 12;
unsigned int prime2 = 5;
unsigned int exponent2 = 12;
std::cin >> layer >> prime1 >> exponent1 >> prime2 >> exponent2;
// analyze for each number between 0 and layer how often they contain prime1 and prime2
Exponents mulPrime1 = { 0 };
Exponents mulPrime2 = { 0 };
for (unsigned int x = 1; x <= layer; x++)
{
auto current = x;
unsigned int count = 0;
// extract first prime (=2) as often as possible
while (current % prime1 == 0)
{
current /= prime1;
count++;
}
mulPrime1.push_back(count);
count = 0;
// extract second prime (=5) as often as possible
while (current % prime2 == 0)
{
current /= prime2;
count++;
}
mulPrime2.push_back(count);
}
// sum1[x] = sum of mulPrime1[0 ... x]
Exponents sum1;
unsigned int count = 0;
for (auto x : mulPrime1)
{
count += x;
sum1.push_back(count);
}
// the same stuff for the other prime
Exponents sum2;
count = 0;
for (auto x : mulPrime2)
{
count += x;
sum2.push_back(count);
}
unsigned long long result = 0;
for (unsigned int i = 0; i <= layer; i++)
{
// how often is each prime used by C(layer, i) ?
auto found1 = choose(sum1, layer, i);
auto found2 = choose(sum2, layer, i);
// already enough ?
if (found1 >= exponent1 && found2 >= exponent2)
{
// no need to enter the inner-most loop, each iteration would succeed
result += i + 1;
continue;
}
// note: abort early because of mirrored values
for (unsigned int j = 0; j <= (i+1) / 2; j++)
{
if (found1 + choose(sum1, i, j) >= exponent1 &&
found2 + choose(sum2, i, j) >= exponent2)
{
// found a match
result++;
// left and right side are identical
if (j < i / 2)
result++;
}
}
}
// and we're done !
std::cout << result << std::endl;
return 0;
}