1 /* GTS - Library for the manipulation of triangulated surfaces
2 * Copyright (C) 1999 St�phane Popinet
3 *
4 * This library is free software; you can redistribute it and/or
5 * modify it under the terms of the GNU Library General Public
6 * License as published by the Free Software Foundation; either
7 * version 2 of the License, or (at your option) any later version.
8 *
9 * This library 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 GNU
12 * Library General Public License for more details.
13 *
14 * You should have received a copy of the GNU Library General Public
15 * License along with this library; if not, write to the
16 * Free Software Foundation, Inc., 59 Temple Place - Suite 330,
17 * Boston, MA 02111-1307, USA.
18 */
19
20 /* This program creates a header file defining the various constants needed by
21 * predicates.c. These constant are machine dependent.
22 * adapted from predicates.c by Jonathan Richard Shewchuk
23 */
24
25 #include <stdio.h>
26
27 /* FPU control. We MUST have only double precision (not extended precision) */
28 #include "rounding.h"
29
main(int argc,char * argv[])30 int main (int argc, char * argv[])
31 {
32 double half = 0.5;
33 double check = 1.0, lastcheck;
34 int every_other = 1;
35 /* epsilon = 2^(-p). Used to estimate roundoff errors. */
36 double epsilon = 1.0;
37 /* splitter = 2^ceiling(p / 2) + 1. Used to split floats in half. */
38 double splitter = 1.0;
39 /* A set of coefficients used to calculate maximum roundoff errors. */
40 double resulterrbound;
41 double ccwerrboundA, ccwerrboundB, ccwerrboundC;
42 double o3derrboundA, o3derrboundB, o3derrboundC;
43 double iccerrboundA, iccerrboundB, iccerrboundC;
44 double isperrboundA, isperrboundB, isperrboundC;
45
46 FPU_ROUND_DOUBLE;
47
48 epsilon = 1.0;
49 splitter = 1.0;
50 /* Repeatedly divide `epsilon' by two until it is too small to add to */
51 /* one without causing roundoff. (Also check if the sum is equal to */
52 /* the previous sum, for machines that round up instead of using exact */
53 /* rounding. Not that this library will work on such machines anyway). */
54 do {
55 lastcheck = check;
56 epsilon *= half;
57 if (every_other) {
58 splitter *= 2.0;
59 }
60 every_other = !every_other;
61 check = 1.0 + epsilon;
62 } while ((check != 1.0) && (check != lastcheck));
63 splitter += 1.0;
64 /* Error bounds for orientation and incircle tests. */
65 resulterrbound = (3.0 + 8.0 * epsilon) * epsilon;
66 ccwerrboundA = (3.0 + 16.0 * epsilon) * epsilon;
67 ccwerrboundB = (2.0 + 12.0 * epsilon) * epsilon;
68 ccwerrboundC = (9.0 + 64.0 * epsilon) * epsilon * epsilon;
69 o3derrboundA = (7.0 + 56.0 * epsilon) * epsilon;
70 o3derrboundB = (3.0 + 28.0 * epsilon) * epsilon;
71 o3derrboundC = (26.0 + 288.0 * epsilon) * epsilon * epsilon;
72 iccerrboundA = (10.0 + 96.0 * epsilon) * epsilon;
73 iccerrboundB = (4.0 + 48.0 * epsilon) * epsilon;
74 iccerrboundC = (44.0 + 576.0 * epsilon) * epsilon * epsilon;
75 isperrboundA = (16.0 + 224.0 * epsilon) * epsilon;
76 isperrboundB = (5.0 + 72.0 * epsilon) * epsilon;
77 isperrboundC = (71.0 + 1408.0 * epsilon) * epsilon * epsilon;
78
79 puts ("/* This file was generated automatically by predicates_init\n"
80 " *\n"
81 " * This file is free software; you can redistribute it and/or\n"
82 " * modify it under the terms of the GNU Library General Public\n"
83 " * License as published by the Free Software Foundation; either\n"
84 " * version 2 of the License, or (at your option) any later version.\n"
85 " *\n"
86 " * This file is distributed in the hope that it will be useful,\n"
87 " * but WITHOUT ANY WARRANTY; without even the implied warranty of\n"
88 " * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.\n"
89 " */\n");
90 printf ("static double splitter = %f;\n", splitter);
91 printf ("static double resulterrbound = %.16g;\n", resulterrbound);
92 printf ("static double ccwerrboundA = %.16g;\n", ccwerrboundA);
93 printf ("static double ccwerrboundB = %.16g;\n", ccwerrboundB);
94 printf ("static double ccwerrboundC = %.16g;\n", ccwerrboundC);
95 printf ("static double o3derrboundA = %.16g;\n", o3derrboundA);
96 printf ("static double o3derrboundB = %.16g;\n", o3derrboundB);
97 printf ("static double o3derrboundC = %.16g;\n", o3derrboundC);
98 printf ("static double iccerrboundA = %.16g;\n", iccerrboundA);
99 printf ("static double iccerrboundB = %.16g;\n", iccerrboundB);
100 printf ("static double iccerrboundC = %.16g;\n", iccerrboundC);
101 printf ("static double isperrboundA = %.16g;\n", isperrboundA);
102 printf ("static double isperrboundB = %.16g;\n", isperrboundB);
103 printf ("static double isperrboundC = %.16g;\n", isperrboundC);
104
105 FPU_RESTORE;
106
107 return 0;
108 }
109