How can I find an algorithm to solve this problem using C++: given an integer z<=10^100, find the smallest row of Pascal's triangle that contains the number z.
1
1 1
1 2 1
1 3 3 1
1 4 6 4 1
For example if z=6 => result is on the 4th row.
Another way to describe the problem: given integer z<=10^100, find the smallest integer n: exist integer k so that C(k,n) = z.
C(k,n) is combination of n things taken k at a time without repetition
EDIT This solution needs Logarithmic time, it's O(Log z)
. Or maybe O( (Log z)^2 )
.
Say you are looking for n,k
where Binomial(n,k)==z
for a given z.
Each row has its largest value in the middle, so starting from n=0
you increase the row number, n
, as long as the middle value is smaller than the given number. Actually, 10^100 isn't that big, so before row 340 you find a position n0,k0=n0/2
where the value from the triangle is larger than or equal to z
: Binomial(n0,k0)>=z
You walk to the left, i.e. you decrease the column number k
, until eventually you find a value smaller than z
. If there was a matching value in that row you would have hit it by now. k
is not very large, less than 170, so this step won't be executed more often than that and does not present a performance problem.
From here you walk down, increasing n
. Here you will find a steadily increasing value of Binomial[n,k]
. Continue with 3 until the value gets bigger than or equal to z
, then goto 2.
EDIT: This step 3 can loop for a very long time when the row number n
is large, so instead of checking each n
linearly you can do a binary search for n
with Binomial(n,k) >= z > Binomial(n-1,k)
, then it only needs Log(n)
time.
A Python implementation looks like this, C++ is similar but somewhat more cumbersome because you need to use an additional library for arbitrary precision integers:
# Calculate (n-k+1)* ... *n
def getnk( n, k ):
a = n
for u in range( n-k+1, n ):
a = a * u
return a
# Find n such that Binomial(n,k) >= z and Binomial(n-1,k) < z
def find_n( z, k, n0 ):
kfactorial = k
for u in range(2, k):
kfactorial *= u
xk = z * kfactorial
nk0 = getnk( n0, k )
n1=n0*2
nk1 = getnk( n1, k )
# duplicate n while the value is too small
while nk1 < xk:
nk0=nk1
n0=n1
n1*=2
nk1 = getnk( n1, k )
# do a binary search
while n1 > n0 + 1:
n2 = (n0+n1) // 2
nk2 = getnk( n2, k )
if nk2 < xk:
n0 = n2
nk0 = nk2
else:
n1 = n2
nk1 = nk2
return n1, nk1 // kfactorial
def find_pos( z ):
n=0
k=0
nk=1
# start by finding a row where the middle value is bigger than z
while nk < z:
# increase n
n = n + 1
nk = nk * n // (n-k)
if nk >= z:
break
# increase both n and k
n = n + 1
k = k + 1
nk = nk * n // k
# check all subsequent rows for a matching value
while nk != z:
if nk > z:
# decrease k
k = k - 1
nk = nk * (k+1) // (n-k)
else:
# increase n
# either linearly
# n = n + 1
# nk = nk * n // (n-k)
# or using binary search:
n, nk = find_n( z, k, n )
return n, k
z = 56476362530291763837811509925185051642180136064700011445902684545741089307844616509330834616
print( find_pos(z) )
It should print
(5864079763474581, 6)
Stirling estimation for n! can be used to find first row in triangle with binomial coefficient bigger or equal to a given x. Using this estimation we can derive lower and upper bound for
and then by observation that this is the maximum coefficient in row that expands 2n:
P( 2n, 0), P( 2n, 1), P( 2n, 2), ..., P( 2n, 2n -1), P( 2n, 2n)
we can find first row with maximum binomial coefficient bigger or equal to a given x. This is the first row in which x can be looking for, this is not possible that x can be found in the row smaller than this. Note: this may be right hint and give an answer immediately in some cases. At the moment I cannot see other way than to start a brute force search from this row.
template <class T>
T binomial_coefficient(unsigned long n, unsigned long k) {
unsigned long i;
T b;
if (0 == k || n == k) {
return 1;
}
if (k > n) {
return 0;
}
if (k > (n - k)) {
k = n - k;
}
if (1 == k) {
return n;
}
b = 1;
for (i = 1; i <= k; ++i) {
b *= (n - (k - i));
if (b < 0) return -1; /* Overflow */
b /= i;
}
return b;
}
Stirling:
double stirling_lower_bound( int n) {
double n_ = n / 2.0;
double res = pow( 2.0, 2 * n_);
res /= sqrt( n_ * M_PI);
return res * exp( ( -1.0) / ( 6 * n_));
}
double stirling_upper_bound( int n) {
double n_ = n / 2.0;
double res = pow( 2.0, 2 * n_) ;
res /= sqrt( n_ * M_PI);
return res * exp( 1.0 / ( 24 * n_));
}
int stirling_estimate( double x) {
int n = 1;
while ( stirling_lower_bound( n) <= x) {
if ( stirling_upper_bound( n) > x) return n;
++n;
}
return n;
}
usage:
long int search_coefficient( unsigned long int &n, unsigned long int x) {
unsigned long int k = n / 2;
long long middle_coefficient = binomial_coefficient<long long>( n, k);
if( middle_coefficient == x) return k;
unsigned long int right = binomial_coefficient<unsigned long>( n, ++k);
while ( x != right) {
while( x < right || x < ( right * ( n + 1) / ( k + 1))) {
right = right * ( n + 1) / ( ++k) - right;
}
if ( right == x) return k;
right = right * ( ++n) / ( ++k);
if( right > x) return -1;
}
return k;
}
/*
*
*/
int main(int argc, char** argv) {
long long x2 = 1365;
unsigned long int n = stirling_estimate( x2);
long int k = search_coefficient( n, x2);
std::cout << "row:" << n <<", column: " << k;
return 0;
}
output:
row:15, column: 11
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With