1 #include <projects.h>
2 static void /* sum coefficients less than res */
eval(projUV ** w,int nu,int nv,double res,projUV * resid)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 */
makeT(int nru,int nrv)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 *
mk_cheby(projUV a,projUV b,double res,projUV * resid,projUV (* func)(projUV),int nu,int nv,int power)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 = NULL;
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 }
178