1 /* Copyright (C) 2003 Free Software Foundation, Inc.
2
3 This file is free software; you can redistribute it and/or modify it
4 under the terms of the GNU General Public License as published by the
5 Free Software Foundation; either version 2, or (at your option) any
6 later version.
7
8 In addition to the permissions in the GNU General Public License, the
9 Free Software Foundation gives you unlimited permission to link the
10 compiled version of this file into combinations with other programs,
11 and to distribute those combinations without any restriction coming
12 from the use of this file. (The General Public License restrictions
13 do apply in other respects; for example, they cover modification of
14 the file, and distribution when not linked into a combine
15 executable.)
16
17 This file is distributed in the hope that it will be useful, but
18 WITHOUT ANY WARRANTY; without even the implied warranty of
19 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 General Public License for more details.
21
22 You should have received a copy of the GNU General Public License
23 along with this program; see the file COPYING. If not, write to
24 the Free Software Foundation, 51 Franklin Street, Fifth Floor,
25 Boston, MA 02110-1301, USA. */
26
27 /* Calculate division table for SH5Media integer division
28 Contributed by Joern Rennecke
29 joern.rennecke@superh.com */
30
31 #include <stdio.h>
32 #include <math.h>
33
34 #define BITS 5
35 #define N_ENTRIES (1 << BITS)
36 #define CUTOFF_BITS 20
37
38 #define BIAS (-330)
39
40 double max_defect = 0.;
41 double max_defect_x;
42
43 double min_defect = 1e9;
44 double min_defect_x;
45
46 double max_defect2 = 0.;
47 double max_defect2_x;
48
49 double min_defect2 = 0.;
50 double min_defect2_x;
51
52 double min_defect3 = 01e9;
53 double min_defect3_x;
54 int min_defect3_val;
55
56 double max_defect3 = 0.;
57 double max_defect3_x;
58 int max_defect3_val;
59
note_defect3(int val,double d2,double y2d,double x)60 static double note_defect3 (int val, double d2, double y2d, double x)
61 {
62 int cutoff_val = val >> CUTOFF_BITS;
63 double cutoff;
64 double defect;
65
66 if (val < 0)
67 cutoff_val++;
68 cutoff = (cutoff_val * (1<<CUTOFF_BITS) - val) * y2d;
69 defect = cutoff + val * d2;
70 if (val < 0)
71 defect = - defect;
72 if (defect > max_defect3)
73 {
74 max_defect3 = defect;
75 max_defect3_x = x;
76 max_defect3_val = val;
77 }
78 if (defect < min_defect3)
79 {
80 min_defect3 = defect;
81 min_defect3_x = x;
82 min_defect3_val = val;
83 }
84 }
85
86 /* This function assumes 32 bit integers. */
87 static double
calc_defect(double x,int constant,int factor)88 calc_defect (double x, int constant, int factor)
89 {
90 double y0 = (constant - (int) floor ((x * factor * 64.))) / 16384.;
91 double y1 = 2 * y0 -y0 * y0 * (x + BIAS / (1.*(1LL<<30)));
92 double y2d0, y2d;
93 int y2d1;
94 double d, d2;
95
96 y1 = floor (y1 * (1024 * 1024 * 1024)) / (1024 * 1024 * 1024);
97 d = y1 - 1 / x;
98 if (d > max_defect)
99 {
100 max_defect = d;
101 max_defect_x = x;
102 }
103 if (d < min_defect)
104 {
105 min_defect = d;
106 min_defect_x = x;
107 }
108 y2d0 = floor (y1 * x * (1LL << 60-16));
109 y2d1 = (int) (long long) y2d0;
110 y2d = - floor ((y1 - y0 / (1<<30-14)) * y2d1) / (1LL<<44);
111 d2 = y1 + y2d - 1/x;
112 if (d2 > max_defect2)
113 {
114 max_defect2 = d2;
115 max_defect2_x = x;
116 }
117 if (d2 < min_defect2)
118 {
119 min_defect2 = d2;
120 min_defect2_x = x;
121 }
122 /* zero times anything is trivially zero. */
123 note_defect3 ((1 << CUTOFF_BITS) - 1, d2, y2d, x);
124 note_defect3 (1 << CUTOFF_BITS, d2, y2d, x);
125 note_defect3 ((1U << 31) - (1 << CUTOFF_BITS), d2, y2d, x);
126 note_defect3 ((1U << 31) - 1, d2, y2d, x);
127 note_defect3 (-1, d2, y2d, x);
128 note_defect3 (-(1 << CUTOFF_BITS), d2, y2d, x);
129 note_defect3 ((1U << 31) - (1 << CUTOFF_BITS) + 1, d2, y2d, x);
130 note_defect3 (-(1U << 31), d2, y2d, x);
131 return d;
132 }
133
134 int
main()135 main ()
136 {
137 int i;
138 unsigned char factors[N_ENTRIES];
139 short constants[N_ENTRIES];
140 int steps = N_ENTRIES / 2;
141 double step = 1. / steps;
142 double eps30 = 1. / (1024 * 1024 * 1024);
143
144 for (i = 0; i < N_ENTRIES; i++)
145 {
146 double x_low = (i < steps ? 1. : -3.) + i * step;
147 double x_high = x_low + step - eps30;
148 double x_med;
149 int factor, constant;
150 double low_defect, med_defect, high_defect, max_defect;
151
152 factor = (1./x_low- 1./x_high) / step * 256. + 0.5;
153 if (factor == 256)
154 factor = 255;
155 factors[i] = factor;
156 /* Use minimum of error function for x_med. */
157 x_med = sqrt (256./factor);
158 if (x_low < 0)
159 x_med = - x_med;
160 low_defect = 1. / x_low + x_low * factor / 256.;
161 high_defect = 1. / x_high + x_high * factor / 256.;
162 med_defect = 1. / x_med + x_med * factor / 256.;
163 max_defect
164 = ((low_defect > high_defect) ^ (x_med < 0)) ? low_defect : high_defect;
165 constant = (med_defect + max_defect) * 0.5 * 16384. + 0.5;
166 if (constant < -32768 || constant > 32767)
167 abort ();
168 constants[i] = constant;
169 calc_defect (x_low, constant, factor);
170 calc_defect (x_med, constant, factor);
171 calc_defect (x_high, constant, factor);
172 }
173 printf ("/* This table has been generated by divtab.c .\n");
174 printf ("Defects for bias %d:\n", BIAS);
175 printf (" Max defect: %e at %e\n", max_defect, max_defect_x);
176 printf (" Min defect: %e at %e\n", min_defect, min_defect_x);
177 printf (" Max 2nd step defect: %e at %e\n", max_defect2, max_defect2_x);
178 printf (" Min 2nd step defect: %e at %e\n", min_defect2, min_defect2_x);
179 printf (" Max div defect: %e at %d:%e\n", max_defect3, max_defect3_val, max_defect3_x);
180 printf (" Min div defect: %e at %d:%e\n", min_defect3, min_defect3_val, min_defect3_x);
181 printf (" Defect at 1: %e\n",
182 calc_defect (1., constants[0], factors[0]));
183 printf (" Defect at -2: %e */\n",
184 calc_defect (-2., constants[steps], factors[steps]));
185 printf ("\t.section\t.rodata\n");
186 printf ("\t.balign 2\n");
187 printf ("/* negative division constants */\n");
188 for (i = steps; i < 2 * steps; i++)
189 printf ("\t.word\t%d\n", constants[i]);
190 printf ("/* negative division factors */\n");
191 for (i = steps; i < 2*steps; i++)
192 printf ("\t.byte\t%d\n", factors[i]);
193 printf ("\t.skip %d\n", steps);
194 printf ("\t.global GLOBAL(div_table):\n");
195 printf ("GLOBAL(div_table):\n");
196 printf ("\t.skip %d\n", steps);
197 printf ("/* positive division factors */\n");
198 for (i = 0; i < steps; i++)
199 printf ("\t.byte\t%d\n", factors[i]);
200 printf ("/* positive division constants */\n");
201 for (i = 0; i < steps; i++)
202 printf ("\t.word\t%d\n", constants[i]);
203 exit (0);
204 }
205