proj/src/mk_cheby.c

   1 #include <projects.h>
   2 static void /* sum coefficients less than res */
   3 eval(projUV **w, int nu, int nv, double res, projUV *resid) {
   4     int i, j;
   5     double ab;
   6     projUV *s;
   7
   8     resid->u = resid->v = 0.;
   9     for (i = 0; i < nu; ++i)
  10         for (s = w[i], j = 0; j < nv; ++j, ++s) {
  11             if ((ab = fabs(s->u)) < res)
  12                 resid->u += ab;
  13             if ((ab = fabs(s->v)) < res)
  14                 resid->v += ab;
  15         }
  16 }
  17 static Tseries * /* create power series structure */
  18 makeT(int nru, int nrv) {
  19     Tseries *T;
  20     int i;
  21
  22     if ((T = (Tseries *)pj_malloc(sizeof(Tseries))) &&
  23         (T->cu = (struct PW_COEF *)pj_malloc(
  24             sizeof(struct PW_COEF) * nru)) &&
  25         (T->cv = (struct PW_COEF *)pj_malloc(
  26             sizeof(struct PW_COEF) * nrv))) {
  27         for (i = 0; i < nru; ++i)
  28             T->cu[i].c = 0;
  29         for (i = 0; i < nrv; ++i)
  30             T->cv[i].c = 0;
  31         return T;
  32     } else
  33         return 0;
  34 }
  35 Tseries *
  36 mk_cheby(projUV a, projUV b, double res, projUV *resid, projUV (*func)(projUV),
  37          int nu, int nv, int power) {
  38     int j, i, nru, nrv, *ncu, *ncv;
  39     Tseries *T;
  40     projUV **w;
  41     double cutres;
  42
  43     if (!(w = (projUV **)vector2(nu, nv, sizeof(projUV))) ||
  44         !(ncu = (int *)vector1(nu + nv, sizeof(int))))
  45         return 0;
  46     ncv = ncu + nu;
  47     if (!bchgen(a, b, nu, nv, w, func)) {
  48         projUV *s;
  49         double ab, *p;
  50
  51         /* analyse coefficients and adjust until residual OK */
  52         cutres = res;
  53         for (i = 4; i ; --i) {
  54             eval(w, nu, nv, cutres, resid);
  55             if (resid->u < res && resid->v < res)
  56                 break;
  57             cutres *= 0.5;
  58         }
  59         if (i <= 0) /* warn of too many tries */
  60             resid->u = - resid->u;
  61         /* apply cut resolution and set pointers */
  62         nru = nrv = 0;
  63         for (j = 0; j < nu; ++j) {
  64             ncu[j] = ncv[j] = 0; /* clear column maxes */
  65             for (s = w[j], i = 0; i < nv; ++i, ++s) {
  66                 if ((ab = fabs(s->u)) < cutres) /* < resolution ? */
  67                     s->u = 0.;          /* clear coefficient */
  68                 else
  69                     ncu[j] = i + 1;     /* update column max */
  70                 if ((ab = fabs(s->v)) < cutres) /* same for v coef's */
  71                     s->v = 0.;
  72                 else
  73                     ncv[j] = i + 1;
  74             }
  75             if (ncu[j]) nru = j + 1;    /* update row max */
  76             if (ncv[j]) nrv = j + 1;
  77         }
  78         if (power) { /* convert to bivariate power series */
  79             if (!bch2bps(a, b, w, nu, nv))
  80                 goto error;
  81             /* possible change in some row counts, so readjust */
  82             nru = nrv = 0;
  83             for (j = 0; j < nu; ++j) {
  84                 ncu[j] = ncv[j] = 0; /* clear column maxes */
  85                 for (s = w[j], i = 0; i < nv; ++i, ++s) {
  86                     if (s->u)
  87                         ncu[j] = i + 1; /* update column max */
  88                     if (s->v)
  89                         ncv[j] = i + 1;
  90                 }
  91                 if (ncu[j]) nru = j + 1;        /* update row max */
  92                 if (ncv[j]) nrv = j + 1;
  93             }
  94             if ((T = makeT(nru, nrv)) != NULL ) {
  95                 T->a = a;
  96                 T->b = b;
  97                 T->mu = nru - 1;
  98                 T->mv = nrv - 1;
  99                 T->power = 1;
 100                 for (i = 0; i < nru; ++i) /* store coefficient rows for u */
 101                 {
 102                     if ((T->cu[i].m = ncu[i]) != 0)
 103                     {
 104                         if ((p = T->cu[i].c =
 105                              (double *)pj_malloc(sizeof(double) * ncu[i])))
 106                             for (j = 0; j < ncu[i]; ++j)
 107                                 *p++ = (w[i] + j)->u;
 108                         else
 109                             goto error;
 110                     }
 111                 }
 112                 for (i = 0; i < nrv; ++i) /* same for v */
 113                 {
 114                     if ((T->cv[i].m = ncv[i]) != 0)
 115                     {
 116                         if ((p = T->cv[i].c =
 117                              (double *)pj_malloc(sizeof(double) * ncv[i])))
 118                             for (j = 0; j < ncv[i]; ++j)
 119                                 *p++ = (w[i] + j)->v;
 120                         else
 121                             goto error;
 122                     }
 123                 }
 124             }
 125         } else if ((T = makeT(nru, nrv)) != NULL) {
 126             /* else make returned Chebyshev coefficient structure */
 127             T->mu = nru - 1; /* save row degree */
 128             T->mv = nrv - 1;
 129             T->a.u = a.u + b.u; /* set argument scaling */
 130             T->a.v = a.v + b.v;
 131             T->b.u = 1. / (b.u - a.u);
 132             T->b.v = 1. / (b.v - a.v);
 133             T->power = 0;
 134             for (i = 0; i < nru; ++i) /* store coefficient rows for u */
 135             {
 136                 if ((T->cu[i].m = ncu[i]) != 0)
 137                 {
 138                     if ((p = T->cu[i].c =
 139                          (double *)pj_malloc(sizeof(double) * ncu[i])))
 140                         for (j = 0; j < ncu[i]; ++j)
 141                             *p++ = (w[i] + j)->u;
 142                     else
 143                         goto error;
 144                 }
 145             }
 146             for (i = 0; i < nrv; ++i) /* same for v */
 147             {
 148                 if ((T->cv[i].m = ncv[i]) != 0)
 149                 {
 150                     if ((p = T->cv[i].c =
 151                          (double *)pj_malloc(sizeof(double) * ncv[i])))
 152                         for (j = 0; j < ncv[i]; ++j)
 153                             *p++ = (w[i] + j)->v;
 154                     else
 155                         goto error;
 156                 }
 157             }
 158         } else
 159             goto error;
 160     }
 161     goto gohome;
 162   error:
 163     if (T) { /* pj_dalloc up possible allocations */
 164         for (i = 0; i <= T->mu; ++i)
 165             if (T->cu[i].c)
 166                 pj_dalloc(T->cu[i].c);
 167         for (i = 0; i <= T->mv; ++i)
 168             if (T->cv[i].c)
 169                 pj_dalloc(T->cv[i].c);
 170         pj_dalloc(T);
 171     }
 172     T = 0;
 173   gohome:
 174     freev2((void **) w, nu);
 175     pj_dalloc(ncu);
 176     return T;
 177 }