$ ./binomial 1.6 0.4 3
(1.600000 + 0.400000)^3
= (1.600000 + 0.400000)(1.600000 + 0.400000)(1.600000 + 0.400000)
= (1.600000)(1.600000)(1.600000) + (1.600000)(1.600000)(0.400000) + (1.600000)(0.400000)(1.600000) + (1.600000)(0.400000)(0.400000) + (0.400000)(1.600000)(1.600000) + (0.400000)(1.600000)(0.400000) + (0.400000)(0.400000)(1.600000) + (0.400000)(0.400000)(0.400000)
= 8.000000
Rather than use the FOIL method to compute the 2^20 terms of the expansion, it makes use of the fact that all of the terms consist of of either 4 or 2; a binary dichotomy. Instead an algorithm to convert a the value of a variable of type long to a binary representation is used.
$ time ./binomial 4.0 2.0 20 | tail -n 1
= 3656158440062976.000000
real 0m17.482s
user 0m16.650s
sys 0m0.830s
It may be slow, but at least accurately computed (4 + 2)^20 properly using a very similar algorithm to that of a person doing algebra with pencil and paper who had never heard of the Binomial Theorem.
For solving Powers of Binomials problems, using Binomial Theorem instead would have been an approach which used less running time, of course. Or the following:
double result = pow(a + b, exp);
Here is my binomial.c file:
/*
* Copyright (c) 2009, Sean A.O. Harney
*
* Calculate (a+b)^exp
* TODO: Breaks with exponents larger than 30, check for arch specific limit when parsing cmdline args
* TODO: Change it to use C99 fixed width integer types instead of long.
*/
#include <stdio.h>
#include <stdlib.h>
double binomial_expansion(double a, double b, int exp);
double calcterm(unsigned long x, int width, double a, double b);
int main(int argc, char *argv[])
{
double a, b;
int exp; //positive.
char *endptr;
if (argc < 4)
{
fprintf(stderr, "Usage:\t%s a b exp\n", argv[0]);
exit(1);
}
a = strtod(argv[1], &endptr); // C89
if (endptr == NULL)
{
perror("strtod()");
exit(1);
}
b = strtod(argv[2], &endptr);
if (endptr == NULL)
{
perror("strtod()");
exit(1);
}
exp = strtol(argv[3], &endptr, 10); // atoi() does not error check
if (endptr == NULL)
{
perror("strtol()");
exit(1);
}
if (exp <= 0)
{
fprintf(stderr, "exp must be positive integer.\n");
exit(1);
}
binomial_expansion(a, b, exp);
exit(0);
}
// return value as computed
double binomial_expansion(double a, double b, int exp)
{
long num_terms; // will be 2^exp terms
long i;
double total = 0;
printf("(%f + %f)^%d\n\t= ", a, b, exp);
for (i = 0, num_terms = 2; i < exp; i++)
{
if (i < exp - 1)
num_terms *= 2; //raise num_terms by one power for all iteration but last
printf("(%f + %f)", a, b);
}
printf("\n\t= ");
for (i = 0; i < num_terms; i++)
{
total += calcterm(i, exp, a, b);
if (i < num_terms - 1)
printf(" + ");
}
printf("\n\t= %f\n", total);
return total;
}
// modified from a convert int to binary string representation function
double calcterm(unsigned long x, int width, double a, double b)
{
double s[64]; // max exp size
int i = width - 1;
double total = 0.0;
do
{
s[i--] = (x & 1) ? b : a; // b for 1, a for 0
x >>= 1;
}
while (x > 0);
while (i >= 0)
s[i--] = a; // fill with a (0), all must be fixed width
for (i = 0; i < width; i++)
{
if (i == 0)
total = s[0];
else
total *= s[i];
printf("(%f)", s[i]);
}
return total;
}
