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
|
/* procedures for evaluating Tseries */
#ifndef lint
static const char SCCSID[]="@(#)biveval.c 4.4 93/06/12 GIE REL";
#endif
# include <projects.h>
# define NEAR_ONE 1.00001
static UV
w2, w;
static double /* basic bivariate Chebyshev evaluation */
ceval(C, n) struct PW_COEF *C; {
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);
}
UV /* bivariate Chebyshev polynomial entry point */
bcheval(UV in, Tseries *T) {
UV 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);
out.v = ceval(T->cv, T->mv);
}
return out;
}
UV /* bivariate power polynomial entry point */
bpseval(UV in, Tseries *T) {
UV 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;
}
UV /* general entry point selecting evaluation mode */
biveval(UV in, Tseries *T) {
return (T->power ? bpseval(in, T) : bcheval(in, T));
}
|
