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
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
|
/* catdvi - get text from DVI files
Copyright (C) 2000-01 Bjoern Brill <brill@fs.math.uni-frankfurt.de>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
*/
#ifndef DENSITY_H
#define DENSITY_H
/* Implements a staircase (i.e. piecewise constant) density function
* on an interval [xmin, xmax].
*
* The domain can be integral or "real" (float, double), the range should
* be "real".
*
* The implementation is NOT numerically sophisticated, so don't expect
* miracles.
*/
/* There's no need to use these two typdefs in an application. They
* are here for logical clarity and easy customization and can be changed
* as needed.
*/
#include "bytesex.h" /* for sint32 */
typedef sint32 scdf_domain_t;
typedef double scdf_range_t;
/* convenience defs - c++ does this automatically */
#ifndef __cplusplus
typedef struct scdf_t scdf_t;
typedef struct scdf_step_t scdf_step_t;
#endif
/* The function is stored as a (singly) linked list of steps. Its value
* f(x) is step.f for x in the half-open interval [step.x, step.next->x) .
* For technical reasons, we keep a last step with last.x = xmax and
* last.next = NULL. last.f is not important since any value f(xmax) will
* give the same integral of f.
*
* Typical applications will traverse [xmin, xmax) as a union of subintervals
* [x0, x1) from left to right. We try to keep this direction efficient.
*/
struct scdf_step_t {
scdf_domain_t x;
scdf_range_t f;
scdf_step_t * next;
};
struct scdf_t {
scdf_domain_t xmin;
scdf_domain_t xmax;
scdf_step_t * head;
scdf_step_t * curr;
};
void scdf_init(
scdf_t * this,
scdf_domain_t xmin,
scdf_domain_t xmax,
scdf_range_t f /* The initial (constant) value of f - usually 0 */
);
void scdf_done(scdf_t * this);
/* Join neighbouring steps with the same f. This should be done at the
* end of a sequence of operations traversing [xmin, xmax] .
*/
void scdf_normalize(scdf_t * this);
/* Force the density function to have at least value fmin in the interval
* [x0, x1). Mathematically: let g have value fmin on [x0, x1) and value
* (-infinity) elsewhere, then f is replaced by the pointwise maximum of
* f and g.
*/
void scdf_force_min_value(
scdf_t * this,
scdf_domain_t x0,
scdf_domain_t x1,
scdf_range_t fmin
);
/* Force f to have at least integral Jmin on [x0, x1]. This is currently
* done by first chacking if the integral is large enough anyway, and
* forcing f to have value at least Jmin / (x1 - x0) if not. More
* sophisticated (keeping f smaller in some cases) methods are possible
* and may be implemented later.
*/
void scdf_force_min_integral(
scdf_t * this,
scdf_domain_t x0,
scdf_domain_t x1,
scdf_range_t Jmin
);
/* Find the value of f at x */
scdf_range_t scdf_eval(scdf_t * this, scdf_domain_t x);
/* Compute the integral of f on [x0, x1] */
scdf_range_t scdf_integral(scdf_t * this, scdf_domain_t x0, scdf_domain_t x1);
/* Solve the equation "integral of f on [x0, x1] = J" for x1;
* set errno = EDOM if this is not possible.
*
* The algorithm used has to do a conversion from scdf_range_t to
* scdf_domain_t, which is done by casting a _positive_ value of
* type scdf_range_t to scdf_domain_t. This should result in rounding
* the return value towards (-infinity) in cases where loss of precision
* occurs.
*/
scdf_domain_t scdf_solve_integral_for_x1(
scdf_t * this,
scdf_domain_t x0,
scdf_range_t J
);
/* Dump a textual representation of f to stderr */
void scdf_dump(scdf_t * this);
#endif
|
