Powers of Binomials in C

Expanding (a + b)^n yields the summation of the 2^n terms. This program will show you the expanded formula. The example below shows (1.6 + 0.4)^3 expanded and solved:

$ ./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;
}