blob: 2755489d4fa199cce3924d33cbc44be13b84feef (
plain)
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
|
/*
Fluct.c
compute the fluctuation in the left and right half
this file is part of Suave
last modified 23 Aug 11 th
*/
#if defined(HAVE_LONG_DOUBLE) && defined(HAVE_POWL)
typedef long double realx;
#define XDBL_MAX_EXP LDBL_MAX_EXP
#define XDBL_MAX LDBL_MAX
#define powx powl
#define ldexpx ldexpl
#else
typedef double realx;
#define XDBL_MAX_EXP DBL_MAX_EXP
#define XDBL_MAX DBL_MAX
#define powx pow
#define ldexpx ldexp
#endif
typedef const realx crealx;
typedef struct {
realx fluct;
number n;
} Var;
/*********************************************************************/
static void Fluct(cThis *t, Var *var,
cBounds *b, creal *w, number n, ccount comp, creal avg, creal err)
{
creal *x = w + n;
creal *f = x + n*t->ndim + comp;
count nvar = 2*t->ndim;
creal norm = 1/(err*Max(fabs(avg), err));
creal flat = 2/3./t->flatness;
crealx max = ldexpx(1., (int)((XDBL_MAX_EXP - 2)/t->flatness));
Clear(var, nvar);
while( n-- ) {
count dim;
crealx arg = 1 + fabs(*w++)*Sq(*f - avg)*norm;
crealx ft = powx(arg < max ? arg : max, t->flatness);
f += t->ncomp;
for( dim = 0; dim < t->ndim; ++dim ) {
Var *v = &var[2*dim + (*x++ >= b[dim].mid)];
crealx f = v->fluct + ft;
v->fluct = (f > XDBL_MAX/2) ? XDBL_MAX/2 : f;
++v->n;
}
}
while( nvar-- ) {
var->fluct = powx(var->fluct, flat);
++var;
}
}
|
