#include    <stdio.h>
#include    <stdlib.h>
#define     MxVAR   8
#define     MxDG2   8
#define     MxDG3   5
#define     MxSR1   5
#define     MxSR2   5
#define     MxSR3   5
#define     NELM2   45      /* nelm2[MxDG2-1] */
#define     NELM3   55      /* nelm3[MxDG3-1] */

static double a[MxVAR+1][MxVAR+1+MxSR1];
static double b[NELM2][NELM2+MxSR2], c[NELM3][NELM3+MxSR3];
static double uu[MxDG2*2+1], vv[MxDG2*2+1];
static double xx[MxDG3*2+1], yy[MxDG3*2+1], zz[MxDG3*2+1];
static int kf1=-1, kf2=-1, kf3=-1;
static int nv1, nd2, nd3, ns1, ns2, ns3, ne2, ne3;

int sm1opn(int nvar, int nsrc)
{   int i1, i2;
    if (kf1 > 0)
        { fprintf(stderr, "sm1opn: illegal call sequence\n"); return(-1); }
    if ((nvar<1) || (nsrc<1))
        { fprintf(stderr, "sm1opn: erroneous argument\n"); return(-2); }
    if ((nvar>MxVAR) || (nsrc>MxSR1))
        { fprintf(stderr, "sm1opn: too many variables\n"); return(-2); }
    nv1 = nvar+1; ns1 = nsrc;
    for (i1=0; i1<(MxVAR+1); i1++) {
        for (i2=0; i2<(MxVAR+1+MxSR1); i2++) a[i1][i2] = 0.;
    }
    kf1 = 1; return(0);
}

int sm1ex(double var[], double src[])
{   int i1, i2;
    if (kf1 < 1)
        { fprintf(stderr, "sm1ex: illegal call sequence\n"); return(-1); }
    a[0][0] += 1.;
    for (i2=1; i2<nv1; i2++) a[0][i2] += var[i2-1];
    for (i2=0; i2<ns1; i2++) a[0][nv1+i2] += src[i2];
    for (i1=1; i1<nv1; i1++) {
        for (i2=i1; i2<nv1; i2++) a[i1][i2] += var[i1-1]*var[i2-1];
        for (i2=0; i2<ns1; i2++) a[i1][nv1+i2] += var[i1-1]*src[i2];
    }
    kf1 = 2; return(0);
}

int sm1cls(double *coef)
{   int i1, i2, i3;
    if (kf1 < 1)
        { fprintf(stderr, "sm1cls: illegal call sequence\n"); return(-1); }
    if (a[0][0] < (double)nv1)
        { fprintf(stderr, "sm1cls: not enough data\n"); return(-2); }
    for (i1=1; i1<nv1; i1++) for (i2=0; i2<i1; i2++) a[i1][i2] = a[i2][i1];
    for (i1=0; i1<nv1; i1++) {
        for (i2=i1+1; i2<(nv1+ns1); i2++) a[i1][i2] /= a[i1][i1];
        for (i3=0; i3<nv1; i3++) {
            if (i3 == i1) continue;
            for (i2=i1+1; i2<(nv1+ns1); i2++) a[i3][i2] -= a[i3][i1]*a[i1][i2];
        }
    }
    for (i2=0; i2<ns1; i2++) for (i1=0; i1<nv1; i1++) *coef++ = a[i1][nv1+i2];
    kf1 = 0; return(0);
}

int sm1rv(double vr[], double sr[])
{   int i1, i2;
    double tmp;
    if (kf1 != 0)
        { fprintf(stderr, "sm1rv: illegal call sequence\n"); return(-1); }
    for (i2=0; i2<ns1; i2++) {
        tmp = a[0][nv1+i2];
        for (i1=1; i1<nv1; i1++) tmp += a[i1][nv1+i2] * vr[i1-1];
        sr[i2] = tmp;
    }
    return(0);
}

int sm2opn(int ndeg, int nsrc)
{   static int nelm2[MxDG2] = { 3, 6, 10, 15, 21, 28, 36, 45 };
    int i1, i2;
    if (kf2 > 0)
        { fprintf(stderr, "sm2opn: illegal sequence\n"); return(-1); }
    if ((ndeg<1) || (nsrc<1))
        { fprintf(stderr, "sm2opn: erroneous argument\n"); return(-2); }
    if (ndeg > MxDG2)
        { fprintf(stderr, "sm2opn: too large degree\n"); return(-2); }
    if (nsrc > MxSR2)
        { fprintf(stderr, "sm2opn: too many variables\n"); return(-2); }
    nd2 = ndeg; ns2 = nsrc; ne2 = nelm2[nd2-1];
    for (i1=0; i1<NELM2; i1++) {
        for (i2=0; i2<(NELM2+MxSR2); i2++) b[i1][i2] = 0.;
    }
    uu[0] = vv[0] = 1.; kf2 = 1; return(0);
}

int sm2ex(double u, double v, double src[])
{   int i1, i2, i3, i4, j, k;
    if (kf2 < 1)
        { fprintf(stderr, "sm2ex: illegal call sequence\n"); return(-1); }
    for (j=0; j<(nd2*2); j++) { uu[j+1] = uu[j]*u; vv[j+1] = vv[j]*v; }
    for (j=0,i1=0; i1<=nd2; i1++) {
        for (i2=0; i2<=i1; i2++,j++) {
            for (i3=0; i3<ns2; i3++) b[j][ne2+i3] += src[i3]*uu[i1-i2]*vv[i2];
            for (k=0,i3=0; i3<=nd2; i3++) {
                for (i4=0; i4<=i3; i4++,k++) {
                    b[j][k] += uu[i1-i2+i3-i4]*vv[i2+i4];
                }
            }
        }
    }
    kf2 = 2; return(0);
}

int sm2cls(double *coef)
{   int i1, i2, i3;
    if (kf2 < 0)
        { fprintf(stderr, "sm2cls: illegal call sequence\n"); return(-1); }
    if (b[0][0] < (double)ne2)
        { fprintf(stderr, "sm2cls: not enough data\n"); return(-2); }
    for (i1=0; i1<ne2; i1++) {
        for (i2=i1+1; i2<(ne2+ns2); i2++) b[i1][i2] /= b[i1][i1];
        for (i3=0; i3<ne2; i3++) {
            if (i3 == i1) continue;
            for (i2=i1+1; i2<(ne2+ns2); i2++) b[i3][i2] -= b[i3][i1]*b[i1][i2];
        }
    }
    for (i2=0; i2<ns2; i2++) for (i1=0; i1<ne2; i1++) *coef++ = b[i1][ne2+i2];
    kf2 = 0; return(0);
}

int sm2rv(double ur, double vr, double sr[])
{   int i1, i2, i3, j;
    double tmp;
    if (kf2 != 0)
        { fprintf(stderr, "sm2rv: illegal call sequence\n"); return(-1); }
    for (j=0; j<(nd2*2); j++) { uu[j+1] = uu[j]*ur; vv[j+1] = vv[j]*vr; }
    for (i3=0; i3<ns2; i3++) {
        for (tmp=0.,j=0,i1=0; i1<=nd2; i1++) {
            for (i2=0; i2<=i1; i2++,j++) tmp += b[j][ne2+i3]*uu[i1-i2]*vv[i2];
        }
        sr[i3] = tmp;
    }
    return(0);
}

int sm3opn(int ndeg, int nsrc)
{   static int nelm3[MxDG3] = { 4, 10, 20, 35, 56 };
    int i1, i2;
    if (kf3 > 0)
        { fprintf(stderr, "sm3opn: illegal sequence\n"); return(-1); }
    if ((ndeg<1) || (nsrc<1))
        { fprintf(stderr, "sm3opn: erroneous argument\n"); return(-2); }
    if (ndeg > MxDG3)
        { fprintf(stderr, "sm3opn: too large degree\n"); return(-2); }
    if (nsrc > MxSR3)
        { fprintf(stderr, "sm3opn: too many variables\n"); return(-2); }
    nd3 = ndeg; ns3 = nsrc; ne3 = nelm3[nd3-1];
    for (i1=0; i1<NELM3; i1++) {
        for (i2=0; i2<(NELM3+MxSR3); i2++) c[i1][i2] = 0.;
    }
    xx[0] = yy[0] = zz[0] = 1.; kf3 = 1; return(0);
}

int sm3ex(double x, double y, double z, double src[])
{   int i1, i2, i3, i4, i5, i6, j, k;
    if (kf3 < 1)
        { fprintf(stderr, "sm3ex: illegal call sequence\n"); return(-1); }
    for (j=0; j<(nd3*2); j++) {
        xx[j+1] = xx[j]*x; yy[j+1] = yy[j]*y; zz[j+1] = zz[j]*z;
    }
    for (j=0,i1=0; i1<=nd3; i1++) {
        for (i2=0; i2<=i1; i2++) for (i3=0; i3<=i2; i3++,j++) {
            for (i4=0; i4<ns3; i4++) {
                c[j][ne3+i4] += src[i4]*xx[i1-i2]*yy[i2-i3]*zz[i3];
            }
            for (k=0,i4=0; i4<=nd3; i4++) {
                for (i5=0; i5<=i4; i5++) for (i6=0; i6<=i5; i6++,k++) {
                    c[j][k] += xx[i1-i2+i4-i5]*yy[i2-i3+i5-i6]*zz[i3+i6];
                }
            }
        }
    }
    kf3 = 2; return(0);
}

int sm3cls(double *coef)
{   int i1, i2, i3;
    if (kf3 < 0)
        { fprintf(stderr, "sm3cls: illegal call sequence\n"); return(-1); }
    if (c[0][0] < (double)ne3)
        { fprintf(stderr, "sm3cls: not enough data\n"); return(-2); }
    for (i1=0; i1<ne3; i1++) {
        for (i2=i1+1; i2<(ne3+ns3); i2++) c[i1][i2] /= c[i1][i1];
        for (i3=0; i3<ne3; i3++) {
            if (i3 == i1) continue;
            for (i2=i1+1; i2<(ne3+ns3); i2++) c[i3][i2] -= c[i3][i1]*c[i1][i2];
        }
    }
    for (i2=0; i2<ns3; i2++) for (i1=0; i1<ne3; i1++) *coef++ = c[i1][ne3+i2];
    kf3 = 0; return(0);
}

int sm3rv(double xr, double yr, double zr, double sr[])
{   int i1, i2, i3, i4, j;
    double tmp;
    if (kf3 != 0)
        { fprintf(stderr, "sm3rv: illegal call sequence\n"); return(-1); }
    for (j=0; j<(nd3*2); j++) {
        xx[j+1] = xx[j]*xr; yy[j+1] = yy[j]*yr; zz[j+1] = zz[j]*zr;
    }
    for (i4=0; i4<ns3; i4++) {
        for (tmp=0.,j=0,i1=0; i1<=nd3; i1++) {
            for (i2=0; i2<=i1; i2++) for (i3=0; i3<=i2; i3++,j++) {
                tmp += c[j][ne3+i4]*xx[i1-i2]*yy[i2-i3]*zz[i3];
            }
        }
        sr[i4] = tmp;
    }
    return(0);
}

/*******************************************/
/*   FORTRAN Language Interface Routines   */
/*******************************************/

void sm1opn_(int *nvr, int *nsr)
{   if (sm1opn(*nvr,*nsr) < 0) exit(1); }

void sm1ex_(float fv[], float fs[])
{   double dv[MxVAR], ds[MxSR1];
    int i;
    for (i=0; i<(nv1-1); i++) dv[i] = fv[i];
    for (i=0; i<ns1; i++) ds[i] = fs[i];
    if (sm1ex(dv,ds) < 0) exit(1);
}

void sm1cls_(float *fc, int *mc, int *nc)
{   double dc[(MxVAR+1)*MxSR1];
    int i, i1, i2;
    if (sm1cls(dc) < 0) exit(1);
    for (i=0,i2=0; ((i2<ns1) && (i2<*nc)); i2++) {
        for (i1=0; i1<nv1; i1++, i++) if (i1 < *mc) *fc++ = dc[i];
        while (i1 < *mc) { i1++; fc++; }
    }
}

void sm1rv_(float fvr[], float fsr[])
{   double dvr[MxVAR], dsr[MxSR1];
    int i;
    for (i=0; i<(nv1-1); i++) dvr[i] = fvr[i];
    if (sm1rv(dvr,dsr) < 0) exit(1);
    for (i=0; i<ns1; i++) fsr[i] = dsr[i];
}

void sm2opn_(int *ndg, int *nsr)
{   if (sm2opn(*ndg,*nsr) < 0) exit(1); }

void sm2ex_(float *fu, float *fv, float fs[])
{   double ds[MxSR2];
    int i;
    for (i=0; i<ns2; i++) ds[i] = fs[i];
    if (sm2ex((double)*fu,(double)*fv,ds) < 0) exit(1);
}

void sm2cls_(float *fc, int *mc, int *nc)
{   double dc[NELM2*MxSR2];
    int i, i1, i2;
    if (sm2cls(dc) < 0) exit(1);
    for (i=0,i2=0; ((i2<ns2) && (i2<*nc)); i2++) {
        for (i1=0; i1<ne2; i1++, i++) if (i1 < *mc) *fc++ = dc[i];
        while (i1 < *mc) { i1++; fc++; }
    }
}

void sm2rv_(float *fur, float *fvr, float fsr[])
{   double dsr[MxSR2];
    int i;
    if (sm2rv((double)*fur,(double)*fvr,dsr) < 0) exit(1);
    for (i=0; i<ns2; i++) fsr[i] = dsr[i];
}

void sm3opn_(int *ndg, int *nsr)
{   if (sm3opn(*ndg,*nsr) < 0) exit(1); }

void sm3ex_(float *fx, float *fy, float *fz, float fs[])
{   double ds[MxSR3];
    int i;
    for (i=0; i<ns3; i++) ds[i] = fs[i];
    if (sm3ex((double)*fx,(double)*fy,(double)*fz,ds) < 0) exit(1);
}

void sm3cls_(float *fc, int *mc, int *nc)
{   double dc[NELM3*MxSR3];
    int i, i1, i2;
    if (sm3cls(dc) < 0) exit(1);
    for (i=0,i2=0; ((i2<ns3) && (i2<*nc)); i2++) {
        for (i1=0; i1<ne3; i1++, i++) if (i1 < *mc) *fc++ = dc[i];
        while (i1 < *mc) { i1++; fc++; }
    }
}

void sm3rv_(float *fxr, float *fyr, float *fzr, float fsr[])
{   double dsr[MxSR3];
    int i;
    if (sm3rv((double)*fxr,(double)*fyr,(double)*fzr,dsr) < 0) exit(1);
    for (i=0; i<ns3; i++) fsr[i] = dsr[i];
}