OpenCores
URL https://opencores.org/ocsvn/openrisc/openrisc/trunk

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-dev/] [or1k-gcc/] [libgcc/] [config/] [libbid/] [bid_round.c] - Blame information for rev 801

Go to most recent revision | Details | Compare with Previous | View Log

Line No. Rev Author Line
1 734 jeremybenn 
/* Copyright (C) 2007, 2009  Free Software Foundation, Inc.
2 
 
3 
This file is part of GCC.
4 
 
5 
GCC is free software; you can redistribute it and/or modify it under
6 
the terms of the GNU General Public License as published by the Free
7 
Software Foundation; either version 3, or (at your option) any later
8 
version.
9 
 
10 
GCC is distributed in the hope that it will be useful, but WITHOUT ANY
11 
WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
13 
for more details.
14 
 
15 
Under Section 7 of GPL version 3, you are granted additional
16 
permissions described in the GCC Runtime Library Exception, version
17 
3.1, as published by the Free Software Foundation.
18 
 
19 
You should have received a copy of the GNU General Public License and
20 
a copy of the GCC Runtime Library Exception along with this program;
21 
see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
22 
<http://www.gnu.org/licenses/>.  */
23 
 
24 
/*****************************************************************************
25 
 *
26 
 *  BID64 encoding:
27 
 * ****************************************
28 
 *  63  62              53 52           0
29 
 * |---|------------------|--------------|
30 
 * | S |  Biased Exp (E)  |  Coeff (c)   |
31 
 * |---|------------------|--------------|
32 
 *
33 
 * bias = 398
34 
 * number = (-1)^s * 10^(E-398) * c
35 
 * coefficient range - 0 to (2^53)-1
36 
 * COEFF_MAX = 2^53-1 = 9007199254740991
37 
 *
38 
 *****************************************************************************
39 
 *
40 
 * BID128 encoding:
41 
 *   1-bit sign
42 
 *   14-bit biased exponent in [0x21, 0x3020] = [33, 12320]
43 
 *         unbiased exponent in [-6176, 6111]; exponent bias = 6176
44 
 *   113-bit unsigned binary integer coefficient (49-bit high + 64-bit low)
45 
 *   Note: 10^34-1 ~ 2^112.945555... < 2^113 => coefficient fits in 113 bits
46 
 *
47 
 *   Note: assume invalid encodings are not passed to this function
48 
 *
49 
 * Round a number C with q decimal digits, represented as a binary integer
50 
 * to q - x digits. Six different routines are provided for different values
51 
 * of q. The maximum value of q used in the library is q = 3 * P - 1 where
52 
 * P = 16 or P = 34 (so q <= 111 decimal digits).
53 
 * The partitioning is based on the following, where Kx is the scaled
54 
 * integer representing the value of 10^(-x) rounded up to a number of bits
55 
 * sufficient to ensure correct rounding:
56 
 *
57 
 * --------------------------------------------------------------------------
58 
 * q    x           max. value of  a            max number      min. number
59 
 *                                              of bits in C    of bits in Kx
60 
 * --------------------------------------------------------------------------
61 
 *
62 
 *                          GROUP 1: 64 bits
63 
 *                          round64_2_18 ()
64 
 *
65 
 * 2    [1,1]       10^1 - 1 < 2^3.33            4              4
66 
 * ...  ...         ...                         ...             ...
67 
 * 18   [1,17]      10^18 - 1 < 2^59.80         60              61
68 
 *
69 
 *                        GROUP 2: 128 bits
70 
 *                        round128_19_38 ()
71 
 *
72 
 * 19   [1,18]      10^19 - 1 < 2^63.11         64              65
73 
 * 20   [1,19]      10^20 - 1 < 2^66.44         67              68
74 
 * ...  ...         ...                         ...             ...
75 
 * 38   [1,37]      10^38 - 1 < 2^126.24        127             128
76 
 *
77 
 *                        GROUP 3: 192 bits
78 
 *                        round192_39_57 ()
79 
 *
80 
 * 39   [1,38]      10^39 - 1 < 2^129.56        130             131
81 
 * ...  ...         ...                         ...             ...
82 
 * 57   [1,56]      10^57 - 1 < 2^189.35        190             191
83 
 *
84 
 *                        GROUP 4: 256 bits
85 
 *                        round256_58_76 ()
86 
 *
87 
 * 58   [1,57]      10^58 - 1 < 2^192.68        193             194
88 
 * ...  ...         ...                         ...             ...
89 
 * 76   [1,75]      10^76 - 1 < 2^252.47        253             254
90 
 *
91 
 *                        GROUP 5: 320 bits
92 
 *                        round320_77_96 ()
93 
 *
94 
 * 77   [1,76]      10^77 - 1 < 2^255.79        256             257
95 
 * 78   [1,77]      10^78 - 1 < 2^259.12        260             261
96 
 * ...  ...         ...                         ...             ...
97 
 * 96   [1,95]      10^96 - 1 < 2^318.91        319             320
98 
 *
99 
 *                        GROUP 6: 384 bits
100 
 *                        round384_97_115 ()
101 
 *
102 
 * 97   [1,96]      10^97 - 1 < 2^322.23        323             324
103 
 * ...  ...         ...                         ...             ...
104 
 * 115  [1,114]     10^115 - 1 < 2^382.03       383             384
105 
 *
106 
 ****************************************************************************/
107 
 
108 
#include "bid_internal.h"
109 
 
110 
void
111 
round64_2_18 (int q,
112 
              int x,
113 
              UINT64 C,
114 
              UINT64 * ptr_Cstar,
115 
              int *incr_exp,
116 
              int *ptr_is_midpoint_lt_even,
117 
              int *ptr_is_midpoint_gt_even,
118 
              int *ptr_is_inexact_lt_midpoint,
119 
              int *ptr_is_inexact_gt_midpoint) {
120 
 
121 
  UINT128 P128;
122 
  UINT128 fstar;
123 
  UINT64 Cstar;
124 
  UINT64 tmp64;
125 
  int shift;
126 
  int ind;
127 
 
128 
  // Note:
129 
  //    In round128_2_18() positive numbers with 2 <= q <= 18 will be
130 
  //    rounded to nearest only for 1 <= x <= 3:
131 
  //     x = 1 or x = 2 when q = 17
132 
  //     x = 2 or x = 3 when q = 18
133 
  // However, for generality and possible uses outside the frame of IEEE 754R
134 
  // this implementation works for 1 <= x <= q - 1
135 
 
136 
  // assume *ptr_is_midpoint_lt_even, *ptr_is_midpoint_gt_even,
137 
  // *ptr_is_inexact_lt_midpoint, and *ptr_is_inexact_gt_midpoint are
138 
  // initialized to 0 by the caller
139 
 
140 
  // round a number C with q decimal digits, 2 <= q <= 18
141 
  // to q - x digits, 1 <= x <= 17
142 
  // C = C + 1/2 * 10^x where the result C fits in 64 bits
143 
  // (because the largest value is 999999999999999999 + 50000000000000000 =
144 
  // 0x0e92596fd628ffff, which fits in 60 bits)
145 
  ind = x - 1;  // 0 <= ind <= 16
146 
  C = C + midpoint64[ind];
147 
  // kx ~= 10^(-x), kx = Kx64[ind] * 2^(-Ex), 0 <= ind <= 16
148 
  // P128 = (C + 1/2 * 10^x) * kx * 2^Ex = (C + 1/2 * 10^x) * Kx
149 
  // the approximation kx of 10^(-x) was rounded up to 64 bits
150 
  __mul_64x64_to_128MACH (P128, C, Kx64[ind]);
151 
  // calculate C* = floor (P128) and f*
152 
  // Cstar = P128 >> Ex
153 
  // fstar = low Ex bits of P128
154 
  shift = Ex64m64[ind]; // in [3, 56]
155 
  Cstar = P128.w[1] >> shift;
156 
  fstar.w[1] = P128.w[1] & mask64[ind];
157 
  fstar.w[0] = P128.w[0];
158 
  // the top Ex bits of 10^(-x) are T* = ten2mxtrunc64[ind], e.g.
159 
  // if x=1, T*=ten2mxtrunc64[0]=0xcccccccccccccccc
160 
  // if (0 < f* < 10^(-x)) then the result is a midpoint
161 
  //   if floor(C*) is even then C* = floor(C*) - logical right
162 
  //       shift; C* has q - x decimal digits, correct by Prop. 1)
163 
  //   else if floor(C*) is odd C* = floor(C*)-1 (logical right
164 
  //       shift; C* has q - x decimal digits, correct by Pr. 1)
165 
  // else
166 
  //   C* = floor(C*) (logical right shift; C has q - x decimal digits,
167 
  //       correct by Property 1)
168 
  // in the caling function n = C* * 10^(e+x)
169 
 
170 
  // determine inexactness of the rounding of C*
171 
  // if (0 < f* - 1/2 < 10^(-x)) then
172 
  //   the result is exact
173 
  // else // if (f* - 1/2 > T*) then
174 
  //   the result is inexact
175 
  if (fstar.w[1] > half64[ind] ||
176 
      (fstar.w[1] == half64[ind] && fstar.w[0])) {
177 
    // f* > 1/2 and the result may be exact
178 
    // Calculate f* - 1/2
179 
    tmp64 = fstar.w[1] - half64[ind];
180 
    if (tmp64 || fstar.w[0] > ten2mxtrunc64[ind]) {      // f* - 1/2 > 10^(-x)
181 
      *ptr_is_inexact_lt_midpoint = 1;
182 
    }   // else the result is exact
183 
  } else {      // the result is inexact; f2* <= 1/2
184 
    *ptr_is_inexact_gt_midpoint = 1;
185 
  }
186 
  // check for midpoints (could do this before determining inexactness)
187 
  if (fstar.w[1] == 0 && fstar.w[0] <= ten2mxtrunc64[ind]) {
188 
    // the result is a midpoint
189 
    if (Cstar & 0x01) { // Cstar is odd; MP in [EVEN, ODD]
190 
      // if floor(C*) is odd C = floor(C*) - 1; the result may be 0
191 
      Cstar--;  // Cstar is now even
192 
      *ptr_is_midpoint_gt_even = 1;
193 
      *ptr_is_inexact_lt_midpoint = 0;
194 
      *ptr_is_inexact_gt_midpoint = 0;
195 
    } else {    // else MP in [ODD, EVEN]
196 
      *ptr_is_midpoint_lt_even = 1;
197 
      *ptr_is_inexact_lt_midpoint = 0;
198 
      *ptr_is_inexact_gt_midpoint = 0;
199 
    }
200 
  }
201 
  // check for rounding overflow, which occurs if Cstar = 10^(q-x)
202 
  ind = q - x;  // 1 <= ind <= q - 1
203 
  if (Cstar == ten2k64[ind]) {  // if  Cstar = 10^(q-x)
204 
    Cstar = ten2k64[ind - 1];   // Cstar = 10^(q-x-1)
205 
    *incr_exp = 1;
206 
  } else {      // 10^33 <= Cstar <= 10^34 - 1
207 
    *incr_exp = 0;
208 
  }
209 
  *ptr_Cstar = Cstar;
210 
}
211 
 
212 
 
213 
void
214 
round128_19_38 (int q,
215 
                int x,
216 
                UINT128 C,
217 
                UINT128 * ptr_Cstar,
218 
                int *incr_exp,
219 
                int *ptr_is_midpoint_lt_even,
220 
                int *ptr_is_midpoint_gt_even,
221 
                int *ptr_is_inexact_lt_midpoint,
222 
                int *ptr_is_inexact_gt_midpoint) {
223 
 
224 
  UINT256 P256;
225 
  UINT256 fstar;
226 
  UINT128 Cstar;
227 
  UINT64 tmp64;
228 
  int shift;
229 
  int ind;
230 
 
231 
  // Note:
232 
  //    In round128_19_38() positive numbers with 19 <= q <= 38 will be
233 
  //    rounded to nearest only for 1 <= x <= 23:
234 
  //     x = 3 or x = 4 when q = 19
235 
  //     x = 4 or x = 5 when q = 20
236 
  //     ...
237 
  //     x = 18 or x = 19 when q = 34
238 
  //     x = 1 or x = 2 or x = 19 or x = 20 when q = 35
239 
  //     x = 2 or x = 3 or x = 20 or x = 21 when q = 36
240 
  //     x = 3 or x = 4 or x = 21 or x = 22 when q = 37
241 
  //     x = 4 or x = 5 or x = 22 or x = 23 when q = 38
242 
  // However, for generality and possible uses outside the frame of IEEE 754R
243 
  // this implementation works for 1 <= x <= q - 1
244 
 
245 
  // assume *ptr_is_midpoint_lt_even, *ptr_is_midpoint_gt_even,
246 
  // *ptr_is_inexact_lt_midpoint, and *ptr_is_inexact_gt_midpoint are
247 
  // initialized to 0 by the caller
248 
 
249 
  // round a number C with q decimal digits, 19 <= q <= 38
250 
  // to q - x digits, 1 <= x <= 37
251 
  // C = C + 1/2 * 10^x where the result C fits in 128 bits
252 
  // (because the largest value is 99999999999999999999999999999999999999 +
253 
  // 5000000000000000000000000000000000000 =
254 
  // 0x4efe43b0c573e7e68a043d8fffffffff, which fits is 127 bits)
255 
 
256 
  ind = x - 1;  // 0 <= ind <= 36
257 
  if (ind <= 18) {      // if 0 <= ind <= 18
258 
    tmp64 = C.w[0];
259 
    C.w[0] = C.w[0] + midpoint64[ind];
260 
    if (C.w[0] < tmp64)
261 
      C.w[1]++;
262 
  } else {      // if 19 <= ind <= 37
263 
    tmp64 = C.w[0];
264 
    C.w[0] = C.w[0] + midpoint128[ind - 19].w[0];
265 
    if (C.w[0] < tmp64) {
266 
      C.w[1]++;
267 
    }
268 
    C.w[1] = C.w[1] + midpoint128[ind - 19].w[1];
269 
  }
270 
  // kx ~= 10^(-x), kx = Kx128[ind] * 2^(-Ex), 0 <= ind <= 36
271 
  // P256 = (C + 1/2 * 10^x) * kx * 2^Ex = (C + 1/2 * 10^x) * Kx
272 
  // the approximation kx of 10^(-x) was rounded up to 128 bits
273 
  __mul_128x128_to_256 (P256, C, Kx128[ind]);
274 
  // calculate C* = floor (P256) and f*
275 
  // Cstar = P256 >> Ex
276 
  // fstar = low Ex bits of P256
277 
  shift = Ex128m128[ind];       // in [2, 63] but have to consider two cases
278 
  if (ind <= 18) {      // if 0 <= ind <= 18
279 
    Cstar.w[0] = (P256.w[2] >> shift) | (P256.w[3] << (64 - shift));
280 
    Cstar.w[1] = (P256.w[3] >> shift);
281 
    fstar.w[0] = P256.w[0];
282 
    fstar.w[1] = P256.w[1];
283 
    fstar.w[2] = P256.w[2] & mask128[ind];
284 
    fstar.w[3] = 0x0ULL;
285 
  } else {      // if 19 <= ind <= 37
286 
    Cstar.w[0] = P256.w[3] >> shift;
287 
    Cstar.w[1] = 0x0ULL;
288 
    fstar.w[0] = P256.w[0];
289 
    fstar.w[1] = P256.w[1];
290 
    fstar.w[2] = P256.w[2];
291 
    fstar.w[3] = P256.w[3] & mask128[ind];
292 
  }
293 
  // the top Ex bits of 10^(-x) are T* = ten2mxtrunc64[ind], e.g.
294 
  // if x=1, T*=ten2mxtrunc128[0]=0xcccccccccccccccccccccccccccccccc
295 
  // if (0 < f* < 10^(-x)) then the result is a midpoint
296 
  //   if floor(C*) is even then C* = floor(C*) - logical right
297 
  //       shift; C* has q - x decimal digits, correct by Prop. 1)
298 
  //   else if floor(C*) is odd C* = floor(C*)-1 (logical right
299 
  //       shift; C* has q - x decimal digits, correct by Pr. 1)
300 
  // else
301 
  //   C* = floor(C*) (logical right shift; C has q - x decimal digits,
302 
  //       correct by Property 1)
303 
  // in the caling function n = C* * 10^(e+x)
304 
 
305 
  // determine inexactness of the rounding of C*
306 
  // if (0 < f* - 1/2 < 10^(-x)) then
307 
  //   the result is exact
308 
  // else // if (f* - 1/2 > T*) then
309 
  //   the result is inexact
310 
  if (ind <= 18) {      // if 0 <= ind <= 18
311 
    if (fstar.w[2] > half128[ind] ||
312 
        (fstar.w[2] == half128[ind] && (fstar.w[1] || fstar.w[0]))) {
313 
      // f* > 1/2 and the result may be exact
314 
      // Calculate f* - 1/2
315 
      tmp64 = fstar.w[2] - half128[ind];
316 
      if (tmp64 || fstar.w[1] > ten2mxtrunc128[ind].w[1] || (fstar.w[1] == ten2mxtrunc128[ind].w[1] && fstar.w[0] > ten2mxtrunc128[ind].w[0])) {  // f* - 1/2 > 10^(-x)
317 
        *ptr_is_inexact_lt_midpoint = 1;
318 
      } // else the result is exact
319 
    } else {    // the result is inexact; f2* <= 1/2
320 
      *ptr_is_inexact_gt_midpoint = 1;
321 
    }
322 
  } else {      // if 19 <= ind <= 37
323 
    if (fstar.w[3] > half128[ind] || (fstar.w[3] == half128[ind] &&
324 
                                      (fstar.w[2] || fstar.w[1]
325 
                                       || fstar.w[0]))) {
326 
      // f* > 1/2 and the result may be exact
327 
      // Calculate f* - 1/2
328 
      tmp64 = fstar.w[3] - half128[ind];
329 
      if (tmp64 || fstar.w[2] || fstar.w[1] > ten2mxtrunc128[ind].w[1] || (fstar.w[1] == ten2mxtrunc128[ind].w[1] && fstar.w[0] > ten2mxtrunc128[ind].w[0])) {    // f* - 1/2 > 10^(-x)
330 
        *ptr_is_inexact_lt_midpoint = 1;
331 
      } // else the result is exact
332 
    } else {    // the result is inexact; f2* <= 1/2
333 
      *ptr_is_inexact_gt_midpoint = 1;
334 
    }
335 
  }
336 
  // check for midpoints (could do this before determining inexactness)
337 
  if (fstar.w[3] == 0 && fstar.w[2] == 0 &&
338 
      (fstar.w[1] < ten2mxtrunc128[ind].w[1] ||
339 
       (fstar.w[1] == ten2mxtrunc128[ind].w[1] &&
340 
        fstar.w[0] <= ten2mxtrunc128[ind].w[0]))) {
341 
    // the result is a midpoint
342 
    if (Cstar.w[0] & 0x01) {     // Cstar is odd; MP in [EVEN, ODD]
343 
      // if floor(C*) is odd C = floor(C*) - 1; the result may be 0
344 
      Cstar.w[0]--;      // Cstar is now even
345 
      if (Cstar.w[0] == 0xffffffffffffffffULL) {
346 
        Cstar.w[1]--;
347 
      }
348 
      *ptr_is_midpoint_gt_even = 1;
349 
      *ptr_is_inexact_lt_midpoint = 0;
350 
      *ptr_is_inexact_gt_midpoint = 0;
351 
    } else {    // else MP in [ODD, EVEN]
352 
      *ptr_is_midpoint_lt_even = 1;
353 
      *ptr_is_inexact_lt_midpoint = 0;
354 
      *ptr_is_inexact_gt_midpoint = 0;
355 
    }
356 
  }
357 
  // check for rounding overflow, which occurs if Cstar = 10^(q-x)
358 
  ind = q - x;  // 1 <= ind <= q - 1
359 
  if (ind <= 19) {
360 
    if (Cstar.w[1] == 0x0ULL && Cstar.w[0] == ten2k64[ind]) {
361 
      // if  Cstar = 10^(q-x)
362 
      Cstar.w[0] = ten2k64[ind - 1];     // Cstar = 10^(q-x-1)
363 
      *incr_exp = 1;
364 
    } else {
365 
      *incr_exp = 0;
366 
    }
367 
  } else if (ind == 20) {
368 
    // if ind = 20
369 
    if (Cstar.w[1] == ten2k128[0].w[1]
370 
        && Cstar.w[0] == ten2k128[0].w[0]) {
371 
      // if  Cstar = 10^(q-x)
372 
      Cstar.w[0] = ten2k64[19];  // Cstar = 10^(q-x-1)
373 
      Cstar.w[1] = 0x0ULL;
374 
      *incr_exp = 1;
375 
    } else {
376 
      *incr_exp = 0;
377 
    }
378 
  } else {      // if 21 <= ind <= 37
379 
    if (Cstar.w[1] == ten2k128[ind - 20].w[1] &&
380 
        Cstar.w[0] == ten2k128[ind - 20].w[0]) {
381 
      // if  Cstar = 10^(q-x)
382 
      Cstar.w[0] = ten2k128[ind - 21].w[0];       // Cstar = 10^(q-x-1)
383 
      Cstar.w[1] = ten2k128[ind - 21].w[1];
384 
      *incr_exp = 1;
385 
    } else {
386 
      *incr_exp = 0;
387 
    }
388 
  }
389 
  ptr_Cstar->w[1] = Cstar.w[1];
390 
  ptr_Cstar->w[0] = Cstar.w[0];
391 
}
392 
 
393 
 
394 
void
395 
round192_39_57 (int q,
396 
                int x,
397 
                UINT192 C,
398 
                UINT192 * ptr_Cstar,
399 
                int *incr_exp,
400 
                int *ptr_is_midpoint_lt_even,
401 
                int *ptr_is_midpoint_gt_even,
402 
                int *ptr_is_inexact_lt_midpoint,
403 
                int *ptr_is_inexact_gt_midpoint) {
404 
 
405 
  UINT384 P384;
406 
  UINT384 fstar;
407 
  UINT192 Cstar;
408 
  UINT64 tmp64;
409 
  int shift;
410 
  int ind;
411 
 
412 
  // Note:
413 
  //    In round192_39_57() positive numbers with 39 <= q <= 57 will be
414 
  //    rounded to nearest only for 5 <= x <= 42:
415 
  //     x = 23 or x = 24 or x = 5 or x = 6 when q = 39
416 
  //     x = 24 or x = 25 or x = 6 or x = 7 when q = 40
417 
  //     ...
418 
  //     x = 41 or x = 42 or x = 23 or x = 24 when q = 57
419 
  // However, for generality and possible uses outside the frame of IEEE 754R
420 
  // this implementation works for 1 <= x <= q - 1
421 
 
422 
  // assume *ptr_is_midpoint_lt_even, *ptr_is_midpoint_gt_even,
423 
  // *ptr_is_inexact_lt_midpoint, and *ptr_is_inexact_gt_midpoint are
424 
  // initialized to 0 by the caller
425 
 
426 
  // round a number C with q decimal digits, 39 <= q <= 57
427 
  // to q - x digits, 1 <= x <= 56
428 
  // C = C + 1/2 * 10^x where the result C fits in 192 bits
429 
  // (because the largest value is
430 
  // 999999999999999999999999999999999999999999999999999999999 +
431 
  //  50000000000000000000000000000000000000000000000000000000 =
432 
  // 0x2ad282f212a1da846afdaf18c034ff09da7fffffffffffff, which fits in 190 bits)
433 
  ind = x - 1;  // 0 <= ind <= 55
434 
  if (ind <= 18) {      // if 0 <= ind <= 18
435 
    tmp64 = C.w[0];
436 
    C.w[0] = C.w[0] + midpoint64[ind];
437 
    if (C.w[0] < tmp64) {
438 
      C.w[1]++;
439 
      if (C.w[1] == 0x0) {
440 
        C.w[2]++;
441 
      }
442 
    }
443 
  } else if (ind <= 37) {       // if 19 <= ind <= 37
444 
    tmp64 = C.w[0];
445 
    C.w[0] = C.w[0] + midpoint128[ind - 19].w[0];
446 
    if (C.w[0] < tmp64) {
447 
      C.w[1]++;
448 
      if (C.w[1] == 0x0) {
449 
        C.w[2]++;
450 
      }
451 
    }
452 
    tmp64 = C.w[1];
453 
    C.w[1] = C.w[1] + midpoint128[ind - 19].w[1];
454 
    if (C.w[1] < tmp64) {
455 
      C.w[2]++;
456 
    }
457 
  } else {      // if 38 <= ind <= 57 (actually ind <= 55)
458 
    tmp64 = C.w[0];
459 
    C.w[0] = C.w[0] + midpoint192[ind - 38].w[0];
460 
    if (C.w[0] < tmp64) {
461 
      C.w[1]++;
462 
      if (C.w[1] == 0x0ull) {
463 
        C.w[2]++;
464 
      }
465 
    }
466 
    tmp64 = C.w[1];
467 
    C.w[1] = C.w[1] + midpoint192[ind - 38].w[1];
468 
    if (C.w[1] < tmp64) {
469 
      C.w[2]++;
470 
    }
471 
    C.w[2] = C.w[2] + midpoint192[ind - 38].w[2];
472 
  }
473 
  // kx ~= 10^(-x), kx = Kx192[ind] * 2^(-Ex), 0 <= ind <= 55
474 
  // P384 = (C + 1/2 * 10^x) * kx * 2^Ex = (C + 1/2 * 10^x) * Kx
475 
  // the approximation kx of 10^(-x) was rounded up to 192 bits
476 
  __mul_192x192_to_384 (P384, C, Kx192[ind]);
477 
  // calculate C* = floor (P384) and f*
478 
  // Cstar = P384 >> Ex
479 
  // fstar = low Ex bits of P384
480 
  shift = Ex192m192[ind];       // in [1, 63] but have to consider three cases
481 
  if (ind <= 18) {      // if 0 <= ind <= 18
482 
    Cstar.w[2] = (P384.w[5] >> shift);
483 
    Cstar.w[1] = (P384.w[5] << (64 - shift)) | (P384.w[4] >> shift);
484 
    Cstar.w[0] = (P384.w[4] << (64 - shift)) | (P384.w[3] >> shift);
485 
    fstar.w[5] = 0x0ULL;
486 
    fstar.w[4] = 0x0ULL;
487 
    fstar.w[3] = P384.w[3] & mask192[ind];
488 
    fstar.w[2] = P384.w[2];
489 
    fstar.w[1] = P384.w[1];
490 
    fstar.w[0] = P384.w[0];
491 
  } else if (ind <= 37) {       // if 19 <= ind <= 37
492 
    Cstar.w[2] = 0x0ULL;
493 
    Cstar.w[1] = P384.w[5] >> shift;
494 
    Cstar.w[0] = (P384.w[5] << (64 - shift)) | (P384.w[4] >> shift);
495 
    fstar.w[5] = 0x0ULL;
496 
    fstar.w[4] = P384.w[4] & mask192[ind];
497 
    fstar.w[3] = P384.w[3];
498 
    fstar.w[2] = P384.w[2];
499 
    fstar.w[1] = P384.w[1];
500 
    fstar.w[0] = P384.w[0];
501 
  } else {      // if 38 <= ind <= 57
502 
    Cstar.w[2] = 0x0ULL;
503 
    Cstar.w[1] = 0x0ULL;
504 
    Cstar.w[0] = P384.w[5] >> shift;
505 
    fstar.w[5] = P384.w[5] & mask192[ind];
506 
    fstar.w[4] = P384.w[4];
507 
    fstar.w[3] = P384.w[3];
508 
    fstar.w[2] = P384.w[2];
509 
    fstar.w[1] = P384.w[1];
510 
    fstar.w[0] = P384.w[0];
511 
  }
512 
 
513 
  // the top Ex bits of 10^(-x) are T* = ten2mxtrunc192[ind], e.g. if x=1,
514 
  // T*=ten2mxtrunc192[0]=0xcccccccccccccccccccccccccccccccccccccccccccccccc
515 
  // if (0 < f* < 10^(-x)) then the result is a midpoint
516 
  //   if floor(C*) is even then C* = floor(C*) - logical right
517 
  //       shift; C* has q - x decimal digits, correct by Prop. 1)
518 
  //   else if floor(C*) is odd C* = floor(C*)-1 (logical right
519 
  //       shift; C* has q - x decimal digits, correct by Pr. 1)
520 
  // else
521 
  //   C* = floor(C*) (logical right shift; C has q - x decimal digits,
522 
  //       correct by Property 1)
523 
  // in the caling function n = C* * 10^(e+x)
524 
 
525 
  // determine inexactness of the rounding of C*
526 
  // if (0 < f* - 1/2 < 10^(-x)) then
527 
  //   the result is exact
528 
  // else // if (f* - 1/2 > T*) then
529 
  //   the result is inexact
530 
  if (ind <= 18) {      // if 0 <= ind <= 18
531 
    if (fstar.w[3] > half192[ind] || (fstar.w[3] == half192[ind] &&
532 
                                      (fstar.w[2] || fstar.w[1]
533 
                                       || fstar.w[0]))) {
534 
      // f* > 1/2 and the result may be exact
535 
      // Calculate f* - 1/2
536 
      tmp64 = fstar.w[3] - half192[ind];
537 
      if (tmp64 || fstar.w[2] > ten2mxtrunc192[ind].w[2] || (fstar.w[2] == ten2mxtrunc192[ind].w[2] && fstar.w[1] > ten2mxtrunc192[ind].w[1]) || (fstar.w[2] == ten2mxtrunc192[ind].w[2] && fstar.w[1] == ten2mxtrunc192[ind].w[1] && fstar.w[0] > ten2mxtrunc192[ind].w[0])) {   // f* - 1/2 > 10^(-x)
538 
        *ptr_is_inexact_lt_midpoint = 1;
539 
      } // else the result is exact
540 
    } else {    // the result is inexact; f2* <= 1/2
541 
      *ptr_is_inexact_gt_midpoint = 1;
542 
    }
543 
  } else if (ind <= 37) {       // if 19 <= ind <= 37
544 
    if (fstar.w[4] > half192[ind] || (fstar.w[4] == half192[ind] &&
545 
                                      (fstar.w[3] || fstar.w[2]
546 
                                       || fstar.w[1] || fstar.w[0]))) {
547 
      // f* > 1/2 and the result may be exact
548 
      // Calculate f* - 1/2
549 
      tmp64 = fstar.w[4] - half192[ind];
550 
      if (tmp64 || fstar.w[3] || fstar.w[2] > ten2mxtrunc192[ind].w[2] || (fstar.w[2] == ten2mxtrunc192[ind].w[2] && fstar.w[1] > ten2mxtrunc192[ind].w[1]) || (fstar.w[2] == ten2mxtrunc192[ind].w[2] && fstar.w[1] == ten2mxtrunc192[ind].w[1] && fstar.w[0] > ten2mxtrunc192[ind].w[0])) {     // f* - 1/2 > 10^(-x)
551 
        *ptr_is_inexact_lt_midpoint = 1;
552 
      } // else the result is exact
553 
    } else {    // the result is inexact; f2* <= 1/2
554 
      *ptr_is_inexact_gt_midpoint = 1;
555 
    }
556 
  } else {      // if 38 <= ind <= 55
557 
    if (fstar.w[5] > half192[ind] || (fstar.w[5] == half192[ind] &&
558 
                                      (fstar.w[4] || fstar.w[3]
559 
                                       || fstar.w[2] || fstar.w[1]
560 
                                       || fstar.w[0]))) {
561 
      // f* > 1/2 and the result may be exact
562 
      // Calculate f* - 1/2
563 
      tmp64 = fstar.w[5] - half192[ind];
564 
      if (tmp64 || fstar.w[4] || fstar.w[3] || fstar.w[2] > ten2mxtrunc192[ind].w[2] || (fstar.w[2] == ten2mxtrunc192[ind].w[2] && fstar.w[1] > ten2mxtrunc192[ind].w[1]) || (fstar.w[2] == ten2mxtrunc192[ind].w[2] && fstar.w[1] == ten2mxtrunc192[ind].w[1] && fstar.w[0] > ten2mxtrunc192[ind].w[0])) {       // f* - 1/2 > 10^(-x)
565 
        *ptr_is_inexact_lt_midpoint = 1;
566 
      } // else the result is exact
567 
    } else {    // the result is inexact; f2* <= 1/2
568 
      *ptr_is_inexact_gt_midpoint = 1;
569 
    }
570 
  }
571 
  // check for midpoints (could do this before determining inexactness)
572 
  if (fstar.w[5] == 0 && fstar.w[4] == 0 && fstar.w[3] == 0 &&
573 
      (fstar.w[2] < ten2mxtrunc192[ind].w[2] ||
574 
       (fstar.w[2] == ten2mxtrunc192[ind].w[2] &&
575 
        fstar.w[1] < ten2mxtrunc192[ind].w[1]) ||
576 
       (fstar.w[2] == ten2mxtrunc192[ind].w[2] &&
577 
        fstar.w[1] == ten2mxtrunc192[ind].w[1] &&
578 
        fstar.w[0] <= ten2mxtrunc192[ind].w[0]))) {
579 
    // the result is a midpoint
580 
    if (Cstar.w[0] & 0x01) {     // Cstar is odd; MP in [EVEN, ODD]
581 
      // if floor(C*) is odd C = floor(C*) - 1; the result may be 0
582 
      Cstar.w[0]--;      // Cstar is now even
583 
      if (Cstar.w[0] == 0xffffffffffffffffULL) {
584 
        Cstar.w[1]--;
585 
        if (Cstar.w[1] == 0xffffffffffffffffULL) {
586 
          Cstar.w[2]--;
587 
        }
588 
      }
589 
      *ptr_is_midpoint_gt_even = 1;
590 
      *ptr_is_inexact_lt_midpoint = 0;
591 
      *ptr_is_inexact_gt_midpoint = 0;
592 
    } else {    // else MP in [ODD, EVEN]
593 
      *ptr_is_midpoint_lt_even = 1;
594 
      *ptr_is_inexact_lt_midpoint = 0;
595 
      *ptr_is_inexact_gt_midpoint = 0;
596 
    }
597 
  }
598 
  // check for rounding overflow, which occurs if Cstar = 10^(q-x)
599 
  ind = q - x;  // 1 <= ind <= q - 1
600 
  if (ind <= 19) {
601 
    if (Cstar.w[2] == 0x0ULL && Cstar.w[1] == 0x0ULL &&
602 
        Cstar.w[0] == ten2k64[ind]) {
603 
      // if  Cstar = 10^(q-x)
604 
      Cstar.w[0] = ten2k64[ind - 1];     // Cstar = 10^(q-x-1)
605 
      *incr_exp = 1;
606 
    } else {
607 
      *incr_exp = 0;
608 
    }
609 
  } else if (ind == 20) {
610 
    // if ind = 20
611 
    if (Cstar.w[2] == 0x0ULL && Cstar.w[1] == ten2k128[0].w[1] &&
612 
        Cstar.w[0] == ten2k128[0].w[0]) {
613 
      // if  Cstar = 10^(q-x)
614 
      Cstar.w[0] = ten2k64[19];  // Cstar = 10^(q-x-1)
615 
      Cstar.w[1] = 0x0ULL;
616 
      *incr_exp = 1;
617 
    } else {
618 
      *incr_exp = 0;
619 
    }
620 
  } else if (ind <= 38) {       // if 21 <= ind <= 38
621 
    if (Cstar.w[2] == 0x0ULL && Cstar.w[1] == ten2k128[ind - 20].w[1] &&
622 
        Cstar.w[0] == ten2k128[ind - 20].w[0]) {
623 
      // if  Cstar = 10^(q-x)
624 
      Cstar.w[0] = ten2k128[ind - 21].w[0];       // Cstar = 10^(q-x-1)
625 
      Cstar.w[1] = ten2k128[ind - 21].w[1];
626 
      *incr_exp = 1;
627 
    } else {
628 
      *incr_exp = 0;
629 
    }
630 
  } else if (ind == 39) {
631 
    if (Cstar.w[2] == ten2k256[0].w[2] && Cstar.w[1] == ten2k256[0].w[1]
632 
        && Cstar.w[0] == ten2k256[0].w[0]) {
633 
      // if  Cstar = 10^(q-x)
634 
      Cstar.w[0] = ten2k128[18].w[0];     // Cstar = 10^(q-x-1)
635 
      Cstar.w[1] = ten2k128[18].w[1];
636 
      Cstar.w[2] = 0x0ULL;
637 
      *incr_exp = 1;
638 
    } else {
639 
      *incr_exp = 0;
640 
    }
641 
  } else {      // if 40 <= ind <= 56
642 
    if (Cstar.w[2] == ten2k256[ind - 39].w[2] &&
643 
        Cstar.w[1] == ten2k256[ind - 39].w[1] &&
644 
        Cstar.w[0] == ten2k256[ind - 39].w[0]) {
645 
      // if  Cstar = 10^(q-x)
646 
      Cstar.w[0] = ten2k256[ind - 40].w[0];       // Cstar = 10^(q-x-1)
647 
      Cstar.w[1] = ten2k256[ind - 40].w[1];
648 
      Cstar.w[2] = ten2k256[ind - 40].w[2];
649 
      *incr_exp = 1;
650 
    } else {
651 
      *incr_exp = 0;
652 
    }
653 
  }
654 
  ptr_Cstar->w[2] = Cstar.w[2];
655 
  ptr_Cstar->w[1] = Cstar.w[1];
656 
  ptr_Cstar->w[0] = Cstar.w[0];
657 
}
658 
 
659 
 
660 
void
661 
round256_58_76 (int q,
662 
                int x,
663 
                UINT256 C,
664 
                UINT256 * ptr_Cstar,
665 
                int *incr_exp,
666 
                int *ptr_is_midpoint_lt_even,
667 
                int *ptr_is_midpoint_gt_even,
668 
                int *ptr_is_inexact_lt_midpoint,
669 
                int *ptr_is_inexact_gt_midpoint) {
670 
 
671 
  UINT512 P512;
672 
  UINT512 fstar;
673 
  UINT256 Cstar;
674 
  UINT64 tmp64;
675 
  int shift;
676 
  int ind;
677 
 
678 
  // Note:
679 
  //    In round256_58_76() positive numbers with 58 <= q <= 76 will be
680 
  //    rounded to nearest only for 24 <= x <= 61:
681 
  //     x = 42 or x = 43 or x = 24 or x = 25 when q = 58
682 
  //     x = 43 or x = 44 or x = 25 or x = 26 when q = 59
683 
  //     ...
684 
  //     x = 60 or x = 61 or x = 42 or x = 43 when q = 76
685 
  // However, for generality and possible uses outside the frame of IEEE 754R
686 
  // this implementation works for 1 <= x <= q - 1
687 
 
688 
  // assume *ptr_is_midpoint_lt_even, *ptr_is_midpoint_gt_even,
689 
  // *ptr_is_inexact_lt_midpoint, and *ptr_is_inexact_gt_midpoint are
690 
  // initialized to 0 by the caller
691 
 
692 
  // round a number C with q decimal digits, 58 <= q <= 76
693 
  // to q - x digits, 1 <= x <= 75
694 
  // C = C + 1/2 * 10^x where the result C fits in 256 bits
695 
  // (because the largest value is 9999999999999999999999999999999999999999
696 
  //     999999999999999999999999999999999999 + 500000000000000000000000000
697 
  //     000000000000000000000000000000000000000000000000 =
698 
  //     0x1736ca15d27a56cae15cf0e7b403d1f2bd6ebb0a50dc83ffffffffffffffffff,
699 
  // which fits in 253 bits)
700 
  ind = x - 1;  // 0 <= ind <= 74
701 
  if (ind <= 18) {      // if 0 <= ind <= 18
702 
    tmp64 = C.w[0];
703 
    C.w[0] = C.w[0] + midpoint64[ind];
704 
    if (C.w[0] < tmp64) {
705 
      C.w[1]++;
706 
      if (C.w[1] == 0x0) {
707 
        C.w[2]++;
708 
        if (C.w[2] == 0x0) {
709 
          C.w[3]++;
710 
        }
711 
      }
712 
    }
713 
  } else if (ind <= 37) {       // if 19 <= ind <= 37
714 
    tmp64 = C.w[0];
715 
    C.w[0] = C.w[0] + midpoint128[ind - 19].w[0];
716 
    if (C.w[0] < tmp64) {
717 
      C.w[1]++;
718 
      if (C.w[1] == 0x0) {
719 
        C.w[2]++;
720 
        if (C.w[2] == 0x0) {
721 
          C.w[3]++;
722 
        }
723 
      }
724 
    }
725 
    tmp64 = C.w[1];
726 
    C.w[1] = C.w[1] + midpoint128[ind - 19].w[1];
727 
    if (C.w[1] < tmp64) {
728 
      C.w[2]++;
729 
      if (C.w[2] == 0x0) {
730 
        C.w[3]++;
731 
      }
732 
    }
733 
  } else if (ind <= 57) {       // if 38 <= ind <= 57
734 
    tmp64 = C.w[0];
735 
    C.w[0] = C.w[0] + midpoint192[ind - 38].w[0];
736 
    if (C.w[0] < tmp64) {
737 
      C.w[1]++;
738 
      if (C.w[1] == 0x0ull) {
739 
        C.w[2]++;
740 
        if (C.w[2] == 0x0) {
741 
          C.w[3]++;
742 
        }
743 
      }
744 
    }
745 
    tmp64 = C.w[1];
746 
    C.w[1] = C.w[1] + midpoint192[ind - 38].w[1];
747 
    if (C.w[1] < tmp64) {
748 
      C.w[2]++;
749 
      if (C.w[2] == 0x0) {
750 
        C.w[3]++;
751 
      }
752 
    }
753 
    tmp64 = C.w[2];
754 
    C.w[2] = C.w[2] + midpoint192[ind - 38].w[2];
755 
    if (C.w[2] < tmp64) {
756 
      C.w[3]++;
757 
    }
758 
  } else {      // if 58 <= ind <= 76 (actually 58 <= ind <= 74)
759 
    tmp64 = C.w[0];
760 
    C.w[0] = C.w[0] + midpoint256[ind - 58].w[0];
761 
    if (C.w[0] < tmp64) {
762 
      C.w[1]++;
763 
      if (C.w[1] == 0x0ull) {
764 
        C.w[2]++;
765 
        if (C.w[2] == 0x0) {
766 
          C.w[3]++;
767 
        }
768 
      }
769 
    }
770 
    tmp64 = C.w[1];
771 
    C.w[1] = C.w[1] + midpoint256[ind - 58].w[1];
772 
    if (C.w[1] < tmp64) {
773 
      C.w[2]++;
774 
      if (C.w[2] == 0x0) {
775 
        C.w[3]++;
776 
      }
777 
    }
778 
    tmp64 = C.w[2];
779 
    C.w[2] = C.w[2] + midpoint256[ind - 58].w[2];
780 
    if (C.w[2] < tmp64) {
781 
      C.w[3]++;
782 
    }
783 
    C.w[3] = C.w[3] + midpoint256[ind - 58].w[3];
784 
  }
785 
  // kx ~= 10^(-x), kx = Kx256[ind] * 2^(-Ex), 0 <= ind <= 74
786 
  // P512 = (C + 1/2 * 10^x) * kx * 2^Ex = (C + 1/2 * 10^x) * Kx
787 
  // the approximation kx of 10^(-x) was rounded up to 192 bits
788 
  __mul_256x256_to_512 (P512, C, Kx256[ind]);
789 
  // calculate C* = floor (P512) and f*
790 
  // Cstar = P512 >> Ex
791 
  // fstar = low Ex bits of P512
792 
  shift = Ex256m256[ind];       // in [0, 63] but have to consider four cases
793 
  if (ind <= 18) {      // if 0 <= ind <= 18
794 
    Cstar.w[3] = (P512.w[7] >> shift);
795 
    Cstar.w[2] = (P512.w[7] << (64 - shift)) | (P512.w[6] >> shift);
796 
    Cstar.w[1] = (P512.w[6] << (64 - shift)) | (P512.w[5] >> shift);
797 
    Cstar.w[0] = (P512.w[5] << (64 - shift)) | (P512.w[4] >> shift);
798 
    fstar.w[7] = 0x0ULL;
799 
    fstar.w[6] = 0x0ULL;
800 
    fstar.w[5] = 0x0ULL;
801 
    fstar.w[4] = P512.w[4] & mask256[ind];
802 
    fstar.w[3] = P512.w[3];
803 
    fstar.w[2] = P512.w[2];
804 
    fstar.w[1] = P512.w[1];
805 
    fstar.w[0] = P512.w[0];
806 
  } else if (ind <= 37) {       // if 19 <= ind <= 37
807 
    Cstar.w[3] = 0x0ULL;
808 
    Cstar.w[2] = P512.w[7] >> shift;
809 
    Cstar.w[1] = (P512.w[7] << (64 - shift)) | (P512.w[6] >> shift);
810 
    Cstar.w[0] = (P512.w[6] << (64 - shift)) | (P512.w[5] >> shift);
811 
    fstar.w[7] = 0x0ULL;
812 
    fstar.w[6] = 0x0ULL;
813 
    fstar.w[5] = P512.w[5] & mask256[ind];
814 
    fstar.w[4] = P512.w[4];
815 
    fstar.w[3] = P512.w[3];
816 
    fstar.w[2] = P512.w[2];
817 
    fstar.w[1] = P512.w[1];
818 
    fstar.w[0] = P512.w[0];
819 
  } else if (ind <= 56) {       // if 38 <= ind <= 56
820 
    Cstar.w[3] = 0x0ULL;
821 
    Cstar.w[2] = 0x0ULL;
822 
    Cstar.w[1] = P512.w[7] >> shift;
823 
    Cstar.w[0] = (P512.w[7] << (64 - shift)) | (P512.w[6] >> shift);
824 
    fstar.w[7] = 0x0ULL;
825 
    fstar.w[6] = P512.w[6] & mask256[ind];
826 
    fstar.w[5] = P512.w[5];
827 
    fstar.w[4] = P512.w[4];
828 
    fstar.w[3] = P512.w[3];
829 
    fstar.w[2] = P512.w[2];
830 
    fstar.w[1] = P512.w[1];
831 
    fstar.w[0] = P512.w[0];
832 
  } else if (ind == 57) {
833 
    Cstar.w[3] = 0x0ULL;
834 
    Cstar.w[2] = 0x0ULL;
835 
    Cstar.w[1] = 0x0ULL;
836 
    Cstar.w[0] = P512.w[7];
837 
    fstar.w[7] = 0x0ULL;
838 
    fstar.w[6] = P512.w[6];
839 
    fstar.w[5] = P512.w[5];
840 
    fstar.w[4] = P512.w[4];
841 
    fstar.w[3] = P512.w[3];
842 
    fstar.w[2] = P512.w[2];
843 
    fstar.w[1] = P512.w[1];
844 
    fstar.w[0] = P512.w[0];
845 
  } else {      // if 58 <= ind <= 74
846 
    Cstar.w[3] = 0x0ULL;
847 
    Cstar.w[2] = 0x0ULL;
848 
    Cstar.w[1] = 0x0ULL;
849 
    Cstar.w[0] = P512.w[7] >> shift;
850 
    fstar.w[7] = P512.w[7] & mask256[ind];
851 
    fstar.w[6] = P512.w[6];
852 
    fstar.w[5] = P512.w[5];
853 
    fstar.w[4] = P512.w[4];
854 
    fstar.w[3] = P512.w[3];
855 
    fstar.w[2] = P512.w[2];
856 
    fstar.w[1] = P512.w[1];
857 
    fstar.w[0] = P512.w[0];
858 
  }
859 
 
860 
  // the top Ex bits of 10^(-x) are T* = ten2mxtrunc256[ind], e.g. if x=1,
861 
  // T*=ten2mxtrunc256[0]=
862 
  //     0xcccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
863 
  // if (0 < f* < 10^(-x)) then the result is a midpoint
864 
  //   if floor(C*) is even then C* = floor(C*) - logical right
865 
  //       shift; C* has q - x decimal digits, correct by Prop. 1)
866 
  //   else if floor(C*) is odd C* = floor(C*)-1 (logical right
867 
  //       shift; C* has q - x decimal digits, correct by Pr. 1)
868 
  // else
869 
  //   C* = floor(C*) (logical right shift; C has q - x decimal digits,
870 
  //       correct by Property 1)
871 
  // in the caling function n = C* * 10^(e+x)
872 
 
873 
  // determine inexactness of the rounding of C*
874 
  // if (0 < f* - 1/2 < 10^(-x)) then
875 
  //   the result is exact
876 
  // else // if (f* - 1/2 > T*) then
877 
  //   the result is inexact
878 
  if (ind <= 18) {      // if 0 <= ind <= 18
879 
    if (fstar.w[4] > half256[ind] || (fstar.w[4] == half256[ind] &&
880 
                                      (fstar.w[3] || fstar.w[2]
881 
                                       || fstar.w[1] || fstar.w[0]))) {
882 
      // f* > 1/2 and the result may be exact
883 
      // Calculate f* - 1/2
884 
      tmp64 = fstar.w[4] - half256[ind];
885 
      if (tmp64 || fstar.w[3] > ten2mxtrunc256[ind].w[2] || (fstar.w[3] == ten2mxtrunc256[ind].w[3] && fstar.w[2] > ten2mxtrunc256[ind].w[2]) || (fstar.w[3] == ten2mxtrunc256[ind].w[3] && fstar.w[2] == ten2mxtrunc256[ind].w[2] && fstar.w[1] > ten2mxtrunc256[ind].w[1]) || (fstar.w[3] == ten2mxtrunc256[ind].w[3] && fstar.w[2] == ten2mxtrunc256[ind].w[2] && fstar.w[1] == ten2mxtrunc256[ind].w[1] && fstar.w[0] > ten2mxtrunc256[ind].w[0])) {  // f* - 1/2 > 10^(-x)
886 
        *ptr_is_inexact_lt_midpoint = 1;
887 
      } // else the result is exact
888 
    } else {    // the result is inexact; f2* <= 1/2
889 
      *ptr_is_inexact_gt_midpoint = 1;
890 
    }
891 
  } else if (ind <= 37) {       // if 19 <= ind <= 37
892 
    if (fstar.w[5] > half256[ind] || (fstar.w[5] == half256[ind] &&
893 
                                      (fstar.w[4] || fstar.w[3]
894 
                                       || fstar.w[2] || fstar.w[1]
895 
                                       || fstar.w[0]))) {
896 
      // f* > 1/2 and the result may be exact
897 
      // Calculate f* - 1/2
898 
      tmp64 = fstar.w[5] - half256[ind];
899 
      if (tmp64 || fstar.w[4] || fstar.w[3] > ten2mxtrunc256[ind].w[3] || (fstar.w[3] == ten2mxtrunc256[ind].w[3] && fstar.w[2] > ten2mxtrunc256[ind].w[2]) || (fstar.w[3] == ten2mxtrunc256[ind].w[3] && fstar.w[2] == ten2mxtrunc256[ind].w[2] && fstar.w[1] > ten2mxtrunc256[ind].w[1]) || (fstar.w[3] == ten2mxtrunc256[ind].w[3] && fstar.w[2] == ten2mxtrunc256[ind].w[2] && fstar.w[1] == ten2mxtrunc256[ind].w[1] && fstar.w[0] > ten2mxtrunc256[ind].w[0])) {    // f* - 1/2 > 10^(-x)
900 
        *ptr_is_inexact_lt_midpoint = 1;
901 
      } // else the result is exact
902 
    } else {    // the result is inexact; f2* <= 1/2
903 
      *ptr_is_inexact_gt_midpoint = 1;
904 
    }
905 
  } else if (ind <= 57) {       // if 38 <= ind <= 57
906 
    if (fstar.w[6] > half256[ind] || (fstar.w[6] == half256[ind] &&
907 
                                      (fstar.w[5] || fstar.w[4]
908 
                                       || fstar.w[3] || fstar.w[2]
909 
                                       || fstar.w[1] || fstar.w[0]))) {
910 
      // f* > 1/2 and the result may be exact
911 
      // Calculate f* - 1/2
912 
      tmp64 = fstar.w[6] - half256[ind];
913 
      if (tmp64 || fstar.w[5] || fstar.w[4] || fstar.w[3] > ten2mxtrunc256[ind].w[3] || (fstar.w[3] == ten2mxtrunc256[ind].w[3] && fstar.w[2] > ten2mxtrunc256[ind].w[2]) || (fstar.w[3] == ten2mxtrunc256[ind].w[3] && fstar.w[2] == ten2mxtrunc256[ind].w[2] && fstar.w[1] > ten2mxtrunc256[ind].w[1]) || (fstar.w[3] == ten2mxtrunc256[ind].w[3] && fstar.w[2] == ten2mxtrunc256[ind].w[2] && fstar.w[1] == ten2mxtrunc256[ind].w[1] && fstar.w[0] > ten2mxtrunc256[ind].w[0])) {      // f* - 1/2 > 10^(-x)
914 
        *ptr_is_inexact_lt_midpoint = 1;
915 
      } // else the result is exact
916 
    } else {    // the result is inexact; f2* <= 1/2
917 
      *ptr_is_inexact_gt_midpoint = 1;
918 
    }
919 
  } else {      // if 58 <= ind <= 74
920 
    if (fstar.w[7] > half256[ind] || (fstar.w[7] == half256[ind] &&
921 
                                      (fstar.w[6] || fstar.w[5]
922 
                                       || fstar.w[4] || fstar.w[3]
923 
                                       || fstar.w[2] || fstar.w[1]
924 
                                       || fstar.w[0]))) {
925 
      // f* > 1/2 and the result may be exact
926 
      // Calculate f* - 1/2
927 
      tmp64 = fstar.w[7] - half256[ind];
928 
      if (tmp64 || fstar.w[6] || fstar.w[5] || fstar.w[4] || fstar.w[3] > ten2mxtrunc256[ind].w[3] || (fstar.w[3] == ten2mxtrunc256[ind].w[3] && fstar.w[2] > ten2mxtrunc256[ind].w[2]) || (fstar.w[3] == ten2mxtrunc256[ind].w[3] && fstar.w[2] == ten2mxtrunc256[ind].w[2] && fstar.w[1] > ten2mxtrunc256[ind].w[1]) || (fstar.w[3] == ten2mxtrunc256[ind].w[3] && fstar.w[2] == ten2mxtrunc256[ind].w[2] && fstar.w[1] == ten2mxtrunc256[ind].w[1] && fstar.w[0] > ten2mxtrunc256[ind].w[0])) {        // f* - 1/2 > 10^(-x)
929 
        *ptr_is_inexact_lt_midpoint = 1;
930 
      } // else the result is exact
931 
    } else {    // the result is inexact; f2* <= 1/2
932 
      *ptr_is_inexact_gt_midpoint = 1;
933 
    }
934 
  }
935 
  // check for midpoints (could do this before determining inexactness)
936 
  if (fstar.w[7] == 0 && fstar.w[6] == 0 &&
937 
      fstar.w[5] == 0 && fstar.w[4] == 0 &&
938 
      (fstar.w[3] < ten2mxtrunc256[ind].w[3] ||
939 
       (fstar.w[3] == ten2mxtrunc256[ind].w[3] &&
940 
        fstar.w[2] < ten2mxtrunc256[ind].w[2]) ||
941 
       (fstar.w[3] == ten2mxtrunc256[ind].w[3] &&
942 
        fstar.w[2] == ten2mxtrunc256[ind].w[2] &&
943 
        fstar.w[1] < ten2mxtrunc256[ind].w[1]) ||
944 
       (fstar.w[3] == ten2mxtrunc256[ind].w[3] &&
945 
        fstar.w[2] == ten2mxtrunc256[ind].w[2] &&
946 
        fstar.w[1] == ten2mxtrunc256[ind].w[1] &&
947 
        fstar.w[0] <= ten2mxtrunc256[ind].w[0]))) {
948 
    // the result is a midpoint
949 
    if (Cstar.w[0] & 0x01) {     // Cstar is odd; MP in [EVEN, ODD]
950 
      // if floor(C*) is odd C = floor(C*) - 1; the result may be 0
951 
      Cstar.w[0]--;      // Cstar is now even
952 
      if (Cstar.w[0] == 0xffffffffffffffffULL) {
953 
        Cstar.w[1]--;
954 
        if (Cstar.w[1] == 0xffffffffffffffffULL) {
955 
          Cstar.w[2]--;
956 
          if (Cstar.w[2] == 0xffffffffffffffffULL) {
957 
            Cstar.w[3]--;
958 
          }
959 
        }
960 
      }
961 
      *ptr_is_midpoint_gt_even = 1;
962 
      *ptr_is_inexact_lt_midpoint = 0;
963 
      *ptr_is_inexact_gt_midpoint = 0;
964 
    } else {    // else MP in [ODD, EVEN]
965 
      *ptr_is_midpoint_lt_even = 1;
966 
      *ptr_is_inexact_lt_midpoint = 0;
967 
      *ptr_is_inexact_gt_midpoint = 0;
968 
    }
969 
  }
970 
  // check for rounding overflow, which occurs if Cstar = 10^(q-x)
971 
  ind = q - x;  // 1 <= ind <= q - 1
972 
  if (ind <= 19) {
973 
    if (Cstar.w[3] == 0x0ULL && Cstar.w[2] == 0x0ULL &&
974 
        Cstar.w[1] == 0x0ULL && Cstar.w[0] == ten2k64[ind]) {
975 
      // if  Cstar = 10^(q-x)
976 
      Cstar.w[0] = ten2k64[ind - 1];     // Cstar = 10^(q-x-1)
977 
      *incr_exp = 1;
978 
    } else {
979 
      *incr_exp = 0;
980 
    }
981 
  } else if (ind == 20) {
982 
    // if ind = 20
983 
    if (Cstar.w[3] == 0x0ULL && Cstar.w[2] == 0x0ULL &&
984 
        Cstar.w[1] == ten2k128[0].w[1]
985 
        && Cstar.w[0] == ten2k128[0].w[0]) {
986 
      // if  Cstar = 10^(q-x)
987 
      Cstar.w[0] = ten2k64[19];  // Cstar = 10^(q-x-1)
988 
      Cstar.w[1] = 0x0ULL;
989 
      *incr_exp = 1;
990 
    } else {
991 
      *incr_exp = 0;
992 
    }
993 
  } else if (ind <= 38) {       // if 21 <= ind <= 38
994 
    if (Cstar.w[3] == 0x0ULL && Cstar.w[2] == 0x0ULL &&
995 
        Cstar.w[1] == ten2k128[ind - 20].w[1] &&
996 
        Cstar.w[0] == ten2k128[ind - 20].w[0]) {
997 
      // if  Cstar = 10^(q-x)
998 
      Cstar.w[0] = ten2k128[ind - 21].w[0];       // Cstar = 10^(q-x-1)
999 
      Cstar.w[1] = ten2k128[ind - 21].w[1];
1000 
      *incr_exp = 1;
1001 
    } else {
1002 
      *incr_exp = 0;
1003 
    }
1004 
  } else if (ind == 39) {
1005 
    if (Cstar.w[3] == 0x0ULL && Cstar.w[2] == ten2k256[0].w[2] &&
1006 
        Cstar.w[1] == ten2k256[0].w[1]
1007 
        && Cstar.w[0] == ten2k256[0].w[0]) {
1008 
      // if  Cstar = 10^(q-x)
1009 
      Cstar.w[0] = ten2k128[18].w[0];     // Cstar = 10^(q-x-1)
1010 
      Cstar.w[1] = ten2k128[18].w[1];
1011 
      Cstar.w[2] = 0x0ULL;
1012 
      *incr_exp = 1;
1013 
    } else {
1014 
      *incr_exp = 0;
1015 
    }
1016 
  } else if (ind <= 57) {       // if 40 <= ind <= 57
1017 
    if (Cstar.w[3] == 0x0ULL && Cstar.w[2] == ten2k256[ind - 39].w[2] &&
1018 
        Cstar.w[1] == ten2k256[ind - 39].w[1] &&
1019 
        Cstar.w[0] == ten2k256[ind - 39].w[0]) {
1020 
      // if  Cstar = 10^(q-x)
1021 
      Cstar.w[0] = ten2k256[ind - 40].w[0];       // Cstar = 10^(q-x-1)
1022 
      Cstar.w[1] = ten2k256[ind - 40].w[1];
1023 
      Cstar.w[2] = ten2k256[ind - 40].w[2];
1024 
      *incr_exp = 1;
1025 
    } else {
1026 
      *incr_exp = 0;
1027 
    }
1028 
    // else if (ind == 58) is not needed becauae we do not have ten2k192[] yet
1029 
  } else {      // if 58 <= ind <= 77 (actually 58 <= ind <= 74)
1030 
    if (Cstar.w[3] == ten2k256[ind - 39].w[3] &&
1031 
        Cstar.w[2] == ten2k256[ind - 39].w[2] &&
1032 
        Cstar.w[1] == ten2k256[ind - 39].w[1] &&
1033 
        Cstar.w[0] == ten2k256[ind - 39].w[0]) {
1034 
      // if  Cstar = 10^(q-x)
1035 
      Cstar.w[0] = ten2k256[ind - 40].w[0];       // Cstar = 10^(q-x-1)
1036 
      Cstar.w[1] = ten2k256[ind - 40].w[1];
1037 
      Cstar.w[2] = ten2k256[ind - 40].w[2];
1038 
      Cstar.w[3] = ten2k256[ind - 40].w[3];
1039 
      *incr_exp = 1;
1040 
    } else {
1041 
      *incr_exp = 0;
1042 
    }
1043 
  }
1044 
  ptr_Cstar->w[3] = Cstar.w[3];
1045 
  ptr_Cstar->w[2] = Cstar.w[2];
1046 
  ptr_Cstar->w[1] = Cstar.w[1];
1047 
  ptr_Cstar->w[0] = Cstar.w[0];
1048 
 
1049 
}

powered by: WebSVN 2.1.0

© copyright 1999-2024 OpenCores.org, equivalent to Oliscience, all rights reserved. OpenCores®, registered trademark.