1 /* catdvi - get text from DVI files
2 Copyright (C) 2000-01 Bjoern Brill <brill@fs.math.uni-frankfurt.de>
3
4 This program is free software; you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation; either version 2 of the License, or
7 (at your option) any later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program; if not, write to the Free Software
16 Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
17 */
18
19 #include "density.h"
20 #include <assert.h>
21 #include <stdlib.h>
22 #include "util.h" /* for enomem() and xmalloc() */
23 #include <stdio.h> /* used by scdf_dump only */
24 #include <errno.h>
25
26 /* Insert a step at x, with value f to the right. where is the nearest
27 * existing step to the left.
28 * Returns address of the newly created step.
29 */
30 scdf_step_t * scdf_insert_after(
31 scdf_t * this,
32 scdf_step_t * where,
33 scdf_domain_t x,
34 scdf_range_t f
35 );
36
37
38 /* Delete the step _after_ (i.e. _right from_ ) where */
39 void scdf_delete_after(scdf_t * this, scdf_step_t * where);
40
41
42 /* Let scdf.curr point to the step in which x lies. Returns new curr. */
43 scdf_step_t * scdf_set_curr_to_x(scdf_t * this, scdf_domain_t x);
44
45 /**************************************************************************
46 **************************************************************************/
47
48
scdf_init(scdf_t * this,scdf_domain_t xmin,scdf_domain_t xmax,scdf_range_t f)49 void scdf_init(
50 scdf_t * this,
51 scdf_domain_t xmin,
52 scdf_domain_t xmax,
53 scdf_range_t f
54 )
55 {
56 scdf_step_t * pa, * pb;
57
58 assert(xmin < xmax);
59
60 this->xmin = xmin;
61 this->xmax = xmax;
62
63 pa = malloc(sizeof(scdf_step_t));
64 if(pa == NULL) enomem();
65 pb = malloc(sizeof(scdf_step_t));
66 if(pb == NULL) enomem();
67
68 pa->x = xmin;
69 pa->f = f;
70 pa->next = pb;
71 pb->x = xmax;
72 pb->f = f;
73 pb->next = NULL;
74
75 this->curr = this->head = pa;
76 }
77
scdf_done(scdf_t * this)78 void scdf_done(scdf_t * this)
79 {
80 scdf_step_t * p0, * p1;
81
82 for(p0 = this->head; p0 != NULL; p0 = p1) {
83 p1 = p0->next;
84 free(p0);
85 }
86 this->curr = this->head = NULL;
87 this->xmin = this->xmax = 0;
88 }
89
90
scdf_insert_after(scdf_t * this,scdf_step_t * where,scdf_domain_t x,scdf_range_t f)91 scdf_step_t * scdf_insert_after(
92 scdf_t * this,
93 scdf_step_t * where,
94 scdf_domain_t x,
95 scdf_range_t f
96 )
97 {
98 scdf_step_t * p;
99
100 assert(where->next != NULL);
101 /* Can't insert after the end of the interval */
102 assert(where->x < x ); assert(x < where->next->x);
103
104 p = malloc(sizeof(scdf_step_t));
105 if(p == NULL) enomem();
106
107 p->x = x;
108 p->f = f;
109 p->next = where->next;
110 where->next = p;
111
112 return p;
113 }
114
115
scdf_delete_after(scdf_t * this,scdf_step_t * where)116 void scdf_delete_after(scdf_t * this, scdf_step_t * where)
117 {
118 scdf_step_t * p;
119
120 p = where->next;
121 assert(p != NULL);
122 /* Can't delete after the end of the interval */
123 assert(p->next != NULL);
124 /* Don't delete the "step" at xmax */
125
126 where->next = p->next;
127 if(this->curr == p) this->curr = where;
128 /* don't leave dangling curr */
129 free(p);
130 }
131
132
scdf_set_curr_to_x(scdf_t * this,scdf_domain_t x)133 scdf_step_t * scdf_set_curr_to_x(scdf_t * this, scdf_domain_t x)
134 {
135 scdf_step_t * p0, * p1;
136
137 assert(this->xmin <= x); assert(x <= this->xmax);
138
139 p0 = this->curr;
140 if(p0->x > x) p0 = this->head;
141 /* too far to the right, start at the left */
142
143 while(p1 = p0->next, p1 != NULL && p1->x <= x) p0 = p1;
144
145 this->curr = p0;
146 return p0;
147 }
148
149
scdf_force_min_value(scdf_t * this,scdf_domain_t x0,scdf_domain_t x1,scdf_range_t fmin)150 void scdf_force_min_value(
151 scdf_t * this,
152 scdf_domain_t x0,
153 scdf_domain_t x1,
154 scdf_range_t fmin
155 )
156 {
157 scdf_step_t * p0, * p1;
158
159 assert(this->xmin <= x0); assert(x0 < x1); assert(x1 <= this->xmax);
160
161 p0 = scdf_set_curr_to_x(this, x0);
162 p1 = p0->next;
163
164 for( ; x0 < x1; x0 = p1->x, p0 = p1, p1 = p1->next) {
165
166 if(p1->x > x1) p1 = scdf_insert_after(this, p0, x1, p0->f);
167 /* Avoid complications that would arise in the last iteration
168 * if x1 wasn't a step boundary by forcing it to be one.
169 * A final scdf_normalize() will clean up this crap.
170 */
171
172 if(fmin <= p0->f) continue;
173 /* the nothing-to-do case: f is already large enough */
174
175 if(p0->x < x0) {
176 /* If we change the function value in the middle of a step
177 * we have to subdivide. The newly created step is automatically
178 * skipped by the for loop iteration code.
179 */
180 scdf_insert_after(this, p0, x0, fmin);
181 continue;
182 }
183
184 /* Now we know p0->x == x0 && x1 <= p1->x && p0->f < fmin ,
185 * i.e. we have a step which is completely contained in [x0, x1]
186 * and its value is too small.
187 */
188 p0->f = fmin;
189 }
190
191 this->curr = p0;
192 /* p0 points to the step containing x1 now - this is where the
193 * next call will most probably start.
194 */
195 }
196
197
scdf_force_min_integral(scdf_t * this,scdf_domain_t x0,scdf_domain_t x1,scdf_range_t Jmin)198 void scdf_force_min_integral(
199 scdf_t * this,
200 scdf_domain_t x0,
201 scdf_domain_t x1,
202 scdf_range_t Jmin
203 )
204 {
205 if(Jmin > scdf_integral(this, x0, x1)) {
206 scdf_force_min_value(this, x0, x1, Jmin / (x1 - x0));
207 }
208 }
209
scdf_normalize(scdf_t * this)210 void scdf_normalize(scdf_t * this)
211 {
212 scdf_step_t * p0, * p1, * p2;
213
214 p0 = this->head;
215 p1 = p0->next; /* always != NULL */
216 p2 = p1->next; /* could be NULL */
217
218 this->curr = this->head;
219 /* May point to free()d memory otherwise. Besides, we expect to
220 * be called after processing the whole interval of definition
221 * anyway.
222 */
223
224 while(p2 != NULL) {
225 if(p0->f == p1->f) scdf_delete_after(this, p0);
226 else p0 = p1;
227
228 p1 = p2;
229 p2 = p2->next;
230 }
231 }
232
233
scdf_solve_integral_for_x1(scdf_t * this,scdf_domain_t x0,scdf_range_t J)234 scdf_domain_t scdf_solve_integral_for_x1(
235 scdf_t * this,
236 scdf_domain_t x0,
237 scdf_range_t J
238 )
239 {
240 scdf_step_t * p0, * p1;
241 scdf_range_t K;
242
243 assert(this->xmin <= x0); assert(x0 < this->xmax);
244 assert(J > 0);
245
246 p0 = scdf_set_curr_to_x(this, x0);
247
248 while(p1 = p0->next, p1 != NULL) {
249 K = p0->f * (p1->x - x0);
250 if(K > J) {
251 return x0 + (scdf_domain_t) (J / p0->f);
252 }
253
254 J -= K;
255 x0 = p1->x;
256 p0 = p1;
257 }
258 errno = EDOM;
259 return this->xmax;
260 }
261
scdf_eval(scdf_t * this,scdf_domain_t x)262 scdf_range_t scdf_eval(scdf_t * this, scdf_domain_t x)
263 {
264 return scdf_set_curr_to_x(this, x)->f;
265 }
266
267
scdf_integral(scdf_t * this,scdf_domain_t x0,scdf_domain_t x1)268 scdf_range_t scdf_integral(scdf_t * this, scdf_domain_t x0, scdf_domain_t x1)
269 {
270 scdf_step_t * p0, * p1;
271 scdf_domain_t d;
272 scdf_range_t J = 0;
273
274 assert(this->xmin <= x0); assert(x0 < x1); assert(x1 <= this->xmax);
275
276 p0 = scdf_set_curr_to_x(this, x0);
277 p1 = p0->next;
278
279 for( ; x0 < x1; x0 = p1->x, p0 = p1, p1 = p1->next) {
280
281 /* Sum it up step by step - and be careful with the last one. */
282 if(x1 < p1->x) d = x1 - x0;
283 else d = p1->x - x0;
284 J += d * p0->f;
285 }
286
287 return J;
288 }
289
scdf_dump(scdf_t * this)290 void scdf_dump(scdf_t * this)
291 {
292 scdf_step_t * p;
293
294 p = this->head;
295 fprintf(stderr, "scdf at %p {\n", (void *) this);
296 for(p = this->head; p != NULL; p = p->next) {
297 fprintf(stderr, "\tf(%10.8g) = %10.8g\n", (double) p->x, (double) p->f);
298 /* should work with any type of domain and range */
299 }
300 fputs("}\n", stderr);
301 }
302
303
scdf_floor_of_integral(scdf_t * this)304 scdf_t * scdf_floor_of_integral(scdf_t * this)
305 {
306 scdf_domain_t x0, x1;
307 scdf_range_t f = 0;
308 scdf_t * result;
309 scdf_step_t * p;
310
311 result = xmalloc(sizeof(scdf_t));
312 scdf_init(result, this->xmin, this->xmax, 0);
313 /* we start with result identically 0 */
314 p = result->head;
315 x0 = this->xmin;
316
317 for( ; ; ) {
318 /* Where should the next jump by 1 of result be? */
319 errno = 0;
320 x1 = scdf_solve_integral_for_x1(this, x0, 1);
321
322 if(errno == EDOM) {
323 /* EDOM means not enough density right from x0 to give integral 1,
324 * i.e. no more jumps and we're done.
325 */
326 break;
327 }
328
329 /* OK, result should jump by 1 at x1 */
330 f += 1;
331 if(x1 == this->xmax) {
332 /* Jump is at the right margin, result already has a pseudo step
333 * there and we're done.
334 */
335 break;
336 }
337
338 /* Insert a step at x1 with value f */
339 p = scdf_insert_after(result, p, x1, f);
340
341 x0 = x1;
342 }
343
344 /* The pseudo step at xmax still has to get its f value */
345 p = scdf_set_curr_to_x(result, this->xmax);
346 p->f = f;
347
348 return result;
349 }
350