---

poly.c

Go to the documentation of this file.

#include <stdio.h>
#include <math.h>
#include "xfig.h"


#define QUAD(x) ((x)*(x)*(x)*(x))
#define CUBE(x) ((x)*(x)*(x))
#define SQ(x) ((x)*(x))

#define N 10
#define XMIN 0
#define XMAX 3
#define YMIN 0
#define YMAX 3

#define OK(i,j) (!(i == 0 && j == 0) && !(i == 1 && j == 0) && !(i == 0 && j == 1) && \
                 !(i == 3 && j == 0) && !(i == 3 && j == 1) && !(i == 3 && j == 3))


void ludcmp(double a[N+1][N+1], int indx[], double *d)
{
  int i,imax,j,k;
  double big,dum,sum,temp;
  double vv[N+1];
  
  *d=1.0;
  for (i=1;i<=N;i++) {
    big=0.0;
    for (j=1;j<=N;j++)
      if ((temp=fabs(a[i][j])) > big) big=temp;
    if (big == 0.0) {
      fprintf(stderr, "Singular matrix in routine ludcmp\n");
      exit(1);
    }
    vv[i]=1.0/big;
  }
  for (j=1;j<=N;j++) {
    for (i=1;i<j;i++) {
      sum=a[i][j];
      for (k=1;k<i;k++) sum -= a[i][k]*a[k][j];
      a[i][j]=sum;
    }
    big=0.0;
    for (i=j;i<=N;i++) {
      sum=a[i][j];
      for (k=1;k<j;k++)
        sum -= a[i][k]*a[k][j];
      a[i][j]=sum;
      if ( (dum=vv[i]*fabs(sum)) >= big) {
        big=dum;
        imax=i;
      }
    }
    if (j != imax) {
      for (k=1;k<=N;k++) {
        dum=a[imax][k];
        a[imax][k]=a[j][k];
        a[j][k]=dum;
      }
      *d = -(*d);
      vv[imax]=vv[j];
    }
    indx[j]=imax;
    if (a[j][j] == 0.0) a[j][j]=1e-30;
    if (j != N) {
      dum=1.0/(a[j][j]);
      for (i=j+1;i<=N;i++) a[i][j] *= dum;
    }
  }
}


void lubksb(double a[N+1][N+1], int indx[], double b[])
{
  int i,ii=0,ip,j;
  double sum;
  
  for (i=1;i<=N;i++) {
    ip=indx[i];
    sum=b[ip];
    b[ip]=b[i];
    if (ii)
      for (j=ii;j<=i-1;j++) sum -= a[i][j]*b[j];
    else if (sum) ii=i;
    b[i]=sum;
  }
  for (i=N;i>=1;i--) {
    sum=b[i];
    for (j=i+1;j<=N;j++) sum -= a[i][j]*b[j];
    b[i]=sum/a[i][i];
  }
}


main(int argc, char *argv[])
{
  int i, j, k = 1, l;
  double a[N+1][N+1], y[N+1][N+1], col[N+1], p;
  int indx[N+1];
  FILE *fptr;
  char s[256];

  for (i = XMIN; i <= XMAX; i++)
    for (j = YMIN; j <= YMAX; j++)
      if (OK(i,j)) {
        l = 1;
/*
        a[k][l++] = QUAD(i);
        a[k][l++] = QUAD(j);
        a[k][l++] = CUBE(i)*j;
        a[k][l++] = i*CUBE(j);
        a[k][l++] = SQ(i)*SQ(j);
*/
        a[k][l++] = CUBE(i);
        a[k][l++] = CUBE(j);
        a[k][l++] = SQ(i)*j;
        a[k][l++] = i*SQ(j);
        a[k][l++] = SQ(i);
        a[k][l++] = SQ(j);
        a[k][l++] = i*j;
        a[k][l++] = i;
        a[k][l++] = j;
        a[k][l++] = 1.0;
        k++;
      }
  for (i = 1; i <= N; i++) {
    for (j = 1; j <= N; j++)
      printf("%g\t", a[i][j]);
    printf("\n");
  }

  ludcmp(a, indx, &p);
  for (j = 1; j <= N; j++) {
    for (i = 1; i <= N; i++) col[i] = 0.0;
    col[j] = 1.0;
    lubksb(a, indx, col);
    for (i = 1; i <= N; i++) y[i][j] = col[i];
  }

  printf("\n");
  for (i = 1; i <= N; i++) {
    for (j = 1; j <= N; j++)
      printf("%g\t", y[i][j]);
    printf("\n");
  }

  printf("\nu[i][j] = ");
  k = 1;
  for (i = XMIN; i <= XMAX; i++)
    for (j = YMIN; j <= YMAX; j++)
      if (OK(i,j)) {
        if (fabs(y[N][k]) > 1e-10) {
          printf("%g*u", y[N][k]);
          if (i == 0)
            printf("[i]");
          else if (i == 1)
            printf("[i+imin]");
          else
            printf("[i+%d*imin]", i);

          if (j == 0)
            printf("[j] ");
          else if (j == 1)
            printf("[j+jmin] ");
          else
            printf("[j+%d*jmin] ", j);
        }
        k++;
      }
  printf("\n");

#define XFIG(x) (100*(x) + 100)
#define YFIG(y) (100*(y) + 100)

  fptr = fopen("poly.fig", "wt");
  fprintf(fptr, "#FIG 2.1\n80 1\n");
  XfigStartCompound(fptr, XFIG((double)XMIN), YFIG((double)YMIN), 
                    XFIG((double)XMAX + 1), YFIG((double)YMAX + 1));
  for (i = XMIN; i <= XMAX + 1; i++)
    XfigLine(fptr, XFIG((double)i), YFIG((double)YMIN),
             XFIG((double)i), YFIG((double)YMAX + 1),
             XFIG_SOLID, 1, XFIG_BLACK);
  for (j = YMIN; j <= YMAX + 1; j++)
    XfigLine(fptr, XFIG((double)XMIN), YFIG((double)j),
             XFIG((double)XMAX + 1), YFIG((double)j),
             XFIG_SOLID, 1, XFIG_BLACK);

  XfigText(fptr, XFIG((double)0.5), YFIG((double)0.6),
           XFIG_CENTER, XFIG_TIMES_ROMAN, 24, 
           XFIG_BLACK, 0.0, "x");  
  k = 1;
  for (i = XMIN; i <= XMAX; i++)
    for (j = YMIN; j <= YMAX; j++)
      if (OK(i,j)) {
        if (fabs(y[N][k]) > 1e-10) {
          sprintf(s, "%g", y[N][k]);
          XfigText(fptr, XFIG((double)i + 0.5), YFIG((double)j + 0.6),
                   XFIG_CENTER, XFIG_TIMES_ROMAN, 24, 
                   XFIG_BLACK, 0.0, s);
        }
        k++;
      }
  XfigEndCompound(fptr);
  fclose(fptr);
}

---