biveval.c
1.9 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
/* procedures for evaluating Tseries */
# include <projects.h>
# define NEAR_ONE 1.00001
static double ceval(struct PW_COEF *C, int n, projUV w, projUV w2) {
double d=0, dd=0, vd, vdd, tmp, *c;
int j;
for (C += n ; n-- ; --C ) {
if ((j = C->m)) {
vd = vdd = 0.;
for (c = C->c + --j; j ; --j ) {
vd = w2.v * (tmp = vd) - vdd + *c--;
vdd = tmp;
}
d = w2.u * (tmp = d) - dd + w.v * vd - vdd + 0.5 * *c;
} else
d = w2.u * (tmp = d) - dd;
dd = tmp;
}
if ((j = C->m)) {
vd = vdd = 0.;
for (c = C->c + --j; j ; --j ) {
vd = w2.v * (tmp = vd) - vdd + *c--;
vdd = tmp;
}
return (w.u * d - dd + 0.5 * ( w.v * vd - vdd + 0.5 * *c ));
} else
return (w.u * d - dd);
}
projUV /* bivariate Chebyshev polynomial entry point */
bcheval(projUV in, Tseries *T) {
projUV w2, w;
projUV out;
/* scale to +-1 */
w.u = ( in.u + in.u - T->a.u ) * T->b.u;
w.v = ( in.v + in.v - T->a.v ) * T->b.v;
if (fabs(w.u) > NEAR_ONE || fabs(w.v) > NEAR_ONE) {
out.u = out.v = HUGE_VAL;
pj_errno = -36;
} else { /* double evaluation */
w2.u = w.u + w.u;
w2.v = w.v + w.v;
out.u = ceval(T->cu, T->mu, w, w2);
out.v = ceval(T->cv, T->mv, w, w2);
}
return out;
}
projUV /* bivariate power polynomial entry point */
bpseval(projUV in, Tseries *T) {
projUV out;
double *c, row;
int i, m;
out.u = out.v = 0.;
for (i = T->mu; i >= 0; --i) {
row = 0.;
if ((m = T->cu[i].m)) {
c = T->cu[i].c + m;
while (m--)
row = *--c + in.v * row;
}
out.u = row + in.u * out.u;
}
for (i = T->mv; i >= 0; --i) {
row = 0.;
if ((m = T->cv[i].m)) {
c = T->cv[i].c + m;
while (m--)
row = *--c + in.v * row;
}
out.v = row + in.u * out.v;
}
return out;
}
projUV /* general entry point selecting evaluation mode */
biveval(projUV in, Tseries *T) {
if (T->power) {
return bpseval(in, T);
} else {
return bcheval(in, T);
}
}