Completely refactored since rev 0. Bonuses claimed: non-square and over 9000.
#define B i%w][j]
Calculating a number like 4000!/(2000!)^2 using multiplication and division seemed rather complicated to me. Instead I went for an approach based on iterating through the rows of Pascal's triangle, using only addition.
Here's an example using an
int32 type, just to show the algorithm clearly:
Where n+m is the number of rows of Pascal's triangle (counting the
1 at the top as row zero.) At the end of the program, a contains rows x and x-1 of Pascal's triangle, the even row in the even cells, and the odd row in the odd cells.
It's important that
w is an odd number, so that on each pass through the array it alternates between writing to the odd and even cells. Also, it is important to ensure that the 1 used to initialise the array is not overwritten on the first pass. This requires that
w/2 be even. I will explain this later.
This program will give the right answer, up to about n=m=16, after which the integer type overflows.
This works in the same way, but uses an array of
char to store the decimal digits (in big endian format.) In order to minimise the code needed to print the answer, the calculations are performed using the ASCII values (add 48 to each value.)
Note that the value of
w is enough to hold rows 9000 and 8999 the same time (a total of 9001+9000=18001 big numbers.) We actually get as far as row 18000 of Pascal's triangle, so by that stage there will not be enough space in the array to hold all 9001+9000 big numbers, and some wrapping around will occur. But this does not affect the central cell of row 9000, which is the only cell we are interested in.
Besides adding whitespace, I have explicitly spelt out the
#define and ungolfed
w/4 to the correct value of 5417. This value is chosen because 18000C9000 has 5417 digits.
for(a[w/2]=1;i<w*(n+m);i+=2) //Initialise array with a 1 in least significant bit of the middle cell of [a]. Loop through the numbers of each row of Pascal's triangle.
for(j=5418;j--;a[1+i%w][j]-=c*10) //Loop through the digits of each number. The code to subtract 10 from the digit if carry flag `c` is executed after the following line.
c=(a[1+i%w][j]=(a[i%w][j]+a[2+i%w][j]+c)%48+48)>57; //Update the cell with a digit that is the sum of the ones to the left and right, taking into account the carry, and using mod 48 arithmetic in order to store the number as its ASCII code. The return value of the assigment to a[1+i%w][j] is compared with ASCII 57 (`9`) and if it is greater the carry flag `c` is set.
My first program took about an hour to print 4000C2000. This one can handle more digits and a higher number of rows to comply with the over 9000 rule. Consequently it takes about 15 hours to print 4000C2000 and would take several days to print 18000C9000.
For 2000x2000 it prints the following (formatted to 80 columns for clarity.) This can be verified to be correct here: http://www.wolframalpha.com/input/?i=4000C2000
First 50 rows of 80 zeros (total 4000) omitted