Align two strings by dynamic programming

salign.c

#include <stdio.h>
#include <string.h>
#include <stdlib.h>

#define xnew(x, n) \
  if ((x = calloc(n, sizeof(*(x)))) == NULL) { \
    fprintf(stderr, "no memory for %s x %u\n", #x, (unsigned) (n)); \
    exit(1); }

/* return the best match of two strings */
static int stralign(const char *s, const char *t,
    int match, int mis, int gap)
{
  int i, j, n, m, len;
  int d, sc[3];
  struct { int di, sc; } **a,  *a_;
  char *sm, *tm;

  n = strlen(s);
  m = strlen(t);
  xnew(a, n+1);
  xnew(a_, (n+1)*(m+1));
  for (i = 0; i <= n; i++)
    a[i] = a_ + (m+1)*i;

  /* dynamic programming to fill to matching table */
  for (i = 0; i <= n; i++) {
    a[i][0].di = 1;
    a[i][0].sc = gap*i;
  }
  for (j = 0; j <= m; j++) {
    a[0][j].sc = 2;
    a[0][j].sc = gap*j;
  }
  for (i = 1; i <= n; i++)
    for (j = 1; j <= m; j++) {
      sc[0] = a[i-1][j-1].sc + ((s[i-1] == t[j-1]) ? match : mis);
      sc[1] = a[i-1][j].sc + gap;
      sc[2] = a[i][j-1].sc + gap;
      if (sc[1] > sc[0]) {
	d = (sc[1] > sc[2]) ? 1 : 2;
      } else {
	d = (sc[0] > sc[2]) ? 0 : 2;
      }
      a[i][j].di = d;
      a[i][j].sc = sc[d];
    }

  /* print the matching table */
  printf("  ");
  for (j = 0; j <= n; j++)
    printf("%5c", (j == 0) ? ' ' : t[j-1]);
  printf("\n");
  for (i = 0; i <= n; i++) {
    printf("%c ", (i == 0) ? ' ' : s[i-1]);
    for (j = 0; j <= m; j++)
      printf("%5d", a[i][j].sc);
    printf("\n");
  }
  printf("score %d\n", a[n][m].sc);

  /* back trace */
  xnew(sm, n+m+2);
  xnew(tm, n+m+2);
  for (len = n+m, i = n, j = m; i > 0 || j > 0; len--) {
    d = a[i][j].di;
    if (d == 0) {
      sm[len] = s[--i];
      tm[len] = t[--j];
    } else if (d == 1) {
      sm[len] = s[--i];
      tm[len] = ' ';
    } else {
      sm[len] = ' ';
      tm[len] = t[--j];
    }
  }
  printf("%s\n%s\n", sm+len+1, tm+len+1);
  i = a[n][m].sc;
  free(sm);
  free(tm);
  free(a);
  free(a_);
  return i;
}

int main(void)
{
  stralign("exponential", "polynomial", 0, -1, -1);
  return 0;
}