1 /*
2 Copyright (C) 2015 Fredrik Johansson
3
4 This file is part of Arb.
5
6 Arb is free software: you can redistribute it and/or modify it under
7 the terms of the GNU Lesser General Public License (LGPL) as published
8 by the Free Software Foundation; either version 2.1 of the License, or
9 (at your option) any later version. See <http://www.gnu.org/licenses/>.
10 */
11
12 #include <string.h>
13 #include <ctype.h>
14 #include "arb.h"
15
16 #define RADIUS_DIGITS 3
17
18 char *
_arb_condense_digits(char * s,slong n)19 _arb_condense_digits(char * s, slong n)
20 {
21 slong i, j, run, out;
22 char * res;
23
24 res = flint_malloc(strlen(s) + 128); /* space for some growth */
25 out = 0;
26
27 for (i = 0; s[i] != '\0'; )
28 {
29 if (isdigit(s[i]))
30 {
31 run = 0;
32
33 for (j = 0; isdigit(s[i + j]); j++)
34 run++;
35
36 if (run > 3 * n)
37 {
38 for (j = 0; j < n; j++)
39 {
40 res[out] = s[i + j];
41 out++;
42 }
43
44 out += flint_sprintf(res + out, "{...%wd digits...}", run - 2 * n);
45
46 for (j = run - n; j < run; j++)
47 {
48 res[out] = s[i + j];
49 out++;
50 }
51 }
52 else
53 {
54 for (j = 0; j < run; j++)
55 {
56 res[out] = s[i + j];
57 out++;
58 }
59 }
60
61 i += run;
62 }
63 else
64 {
65 res[out] = s[i];
66 i++;
67 out++;
68 }
69 }
70
71 res[out] = '\0';
72 res = flint_realloc(res, strlen(res) + 1);
73
74 flint_free(s);
75 return res;
76 }
77
78 /* Format (digits=d, exponent=e) as floating-point or fixed-point.
79 Reallocates the input and mutates the exponent. */
80 void
_arb_digits_as_float_str(char ** d,fmpz_t e,slong minfix,slong maxfix)81 _arb_digits_as_float_str(char ** d, fmpz_t e, slong minfix, slong maxfix)
82 {
83 slong i, n, alloc, dotpos;
84
85 /* do nothing with 0 or something non-numerical */
86 if (!((*d)[0] >= '1' && (*d)[0] <= '9'))
87 return;
88
89 n = strlen(*d);
90
91 fmpz_add_ui(e, e, n - 1);
92
93 /* fixed-point or integer format */
94 /* we require e < n - 1; otherwise we would have to insert trailing zeros
95 [todo: could allow e < n, if printing integers without radix point] */
96 if (fmpz_cmp_si(e, minfix) >= 0 && fmpz_cmp_si(e, maxfix) <= 0 &&
97 fmpz_cmp_si(e, n - 1) < 0)
98 {
99 slong exp = *e;
100
101 /* 0.000xxx */
102 if (exp < 0)
103 {
104 /* 0. + (-1-exp) zeros + digits + null terminator */
105 alloc = 2 + (-1-exp) + n + 1;
106
107 *d = flint_realloc(*d, alloc);
108
109 /* copy in reverse order, including null terminator */
110 for (i = n; i >= 0; i--)
111 (*d)[2 + (-1-exp) + i] = (*d)[i];
112
113 for (i = 0; i < 2 + (-1-exp); i++)
114 (*d)[i] = (i == 1) ? '.' : '0';
115 }
116 else /* xxx.yyy --- must have dotpos < n - 1 */
117 {
118 dotpos = exp + 1;
119 alloc = n + 2; /* space for . and null terminator */
120
121 (*d) = flint_realloc(*d, alloc);
122
123 /* copy fractional part in reverse order, including null */
124 for (i = n; i >= dotpos; i--)
125 (*d)[i + 1] = (*d)[i];
126
127 (*d)[dotpos] = '.';
128 }
129 }
130 else
131 {
132 /* format as xe+zzz or x.yyye+zzz */
133 alloc = n + 1 + 2 + fmpz_sizeinbase(e, 10) + 1;
134 *d = flint_realloc(*d, alloc);
135
136 /* insert . */
137 if (n > 1)
138 {
139 /* copy fractional part in reverse order */
140 for (i = n; i >= 1; i--)
141 (*d)[i + 1] = (*d)[i];
142
143 (*d)[1] = '.';
144 }
145
146 (*d)[n + (n > 1)] = 'e';
147
148 if (fmpz_sgn(e) >= 0)
149 {
150 (*d)[n + (n > 1) + 1] = '+';
151 }
152 else
153 {
154 (*d)[n + (n > 1) + 1] = '-';
155 fmpz_neg(e, e);
156 }
157
158 fmpz_get_str((*d) + n + (n > 1) + 2, 10, e); /* writes null byte */
159 }
160 }
161
162 /* Rounds a string of decimal digits (null-terminated).
163 to length at most n. The rounding mode
164 can be ARF_RND_DOWN, ARF_RND_UP or ARF_RND_NEAR.
165 The string is overwritten in-place, truncating it as necessary.
166 The input should not have a leading sign or leading zero digits,
167 but can have trailing zero digits.
168
169 Computes shift and error such that
170
171 int(input) = int(output) * 10^shift + error
172
173 exactly.
174 */
175 void
_arb_digits_round_inplace(char * s,flint_bitcnt_t * shift,fmpz_t error,slong n,arf_rnd_t rnd)176 _arb_digits_round_inplace(char * s, flint_bitcnt_t * shift, fmpz_t error, slong n, arf_rnd_t rnd)
177 {
178 slong i, m;
179 int up;
180
181 if (n < 1)
182 {
183 flint_printf("_arb_digits_round_inplace: require n >= 1\n");
184 flint_abort();
185 }
186
187 m = strlen(s);
188
189 if (m <= n)
190 {
191 *shift = 0;
192 fmpz_zero(error);
193 return;
194 }
195
196 /* always round down */
197 if (rnd == ARF_RND_DOWN)
198 {
199 up = 0;
200 }
201 else if (rnd == ARF_RND_UP) /* round up if tail is nonzero */
202 {
203 up = 0;
204
205 for (i = n; i < m; i++)
206 {
207 if (s[i] != '0')
208 {
209 up = 1;
210 break;
211 }
212 }
213 }
214 else /* round to nearest (up on tie -- todo: round-to-even?) */
215 {
216 up = (s[n] >= '5' && s[n] <= '9');
217 }
218
219 if (!up)
220 {
221 /* simply truncate */
222 fmpz_set_str(error, s + n, 10);
223 s[n] = '\0';
224 *shift = m - n;
225 }
226 else
227 {
228 int digit, borrow, carry;
229
230 /* error = 10^(m-n) - s[n:], where s[n:] is nonzero */
231 /* i.e. 10s complement the truncated digits */
232 borrow = 0;
233
234 for (i = m - 1; i >= n; i--)
235 {
236 digit = 10 - (s[i] - '0') - borrow;
237
238 if (digit == 10)
239 {
240 digit = 0;
241 borrow = 0;
242 }
243 else
244 {
245 borrow = 1;
246 }
247
248 s[i] = digit + '0';
249 }
250
251 if (!borrow)
252 {
253 flint_printf("expected borrow!\n");
254 flint_abort();
255 }
256
257 fmpz_set_str(error, s + n, 10);
258 fmpz_neg(error, error);
259
260 /* add 1 ulp to the leading digits */
261 carry = 1;
262
263 for (i = n - 1; i >= 0; i--)
264 {
265 digit = (s[i] - '0') + carry;
266
267 if (digit > 9)
268 {
269 digit = 0;
270 carry = 1;
271 }
272 else
273 {
274 carry = 0;
275 }
276
277 s[i] = digit + '0';
278 }
279
280 /* carry-out -- only possible if we started with all 9s,
281 so now the rest will be 0s which we don't have to shift explicitly */
282 if (carry)
283 {
284 s[0] = '1';
285 *shift = m - n + 1;
286 }
287 else
288 {
289 *shift = m - n;
290 }
291
292 s[n] = '\0'; /* truncate */
293 }
294 }
295
296 void
arb_get_str_parts(int * negative,char ** mid_digits,fmpz_t mid_exp,char ** rad_digits,fmpz_t rad_exp,const arb_t x,slong n,int more)297 arb_get_str_parts(int * negative, char **mid_digits, fmpz_t mid_exp,
298 char **rad_digits, fmpz_t rad_exp,
299 const arb_t x, slong n, int more)
300 {
301 fmpz_t mid, rad, exp, err;
302 slong good;
303 flint_bitcnt_t shift;
304
305 if (!arb_is_finite(x))
306 {
307 *negative = 0;
308
309 fmpz_zero(mid_exp);
310 *mid_digits = flint_malloc(4);
311 if (arf_is_nan(arb_midref(x)))
312 strcpy(*mid_digits, "nan");
313 else
314 strcpy(*mid_digits, "0");
315
316 fmpz_zero(rad_exp);
317 *rad_digits = flint_malloc(4);
318 strcpy(*rad_digits, "inf");
319
320 return;
321 }
322
323 fmpz_init(mid);
324 fmpz_init(rad);
325 fmpz_init(exp);
326 fmpz_init(err);
327
328 /* heuristic part */
329 if (!more)
330 {
331 good = arb_rel_accuracy_bits(x) * 0.30102999566398119521 + 2;
332 n = FLINT_MIN(n, good);
333 }
334
335 arb_get_fmpz_mid_rad_10exp(mid, rad, exp, x, FLINT_MAX(n, 1));
336 *negative = arf_sgn(arb_midref(x)) < 0;
337 fmpz_abs(mid, mid);
338
339 *mid_digits = fmpz_get_str(NULL, 10, mid);
340 *rad_digits = NULL;
341
342 /* Truncate further so that 1 ulp error can be guaranteed (rigorous part)
343 Note: mid cannot be zero here if n >= 1 and rad != 0. */
344 if (n >= 1 && !(more || fmpz_is_zero(rad)))
345 {
346 slong lenmid, lenrad, rem;
347
348 *rad_digits = fmpz_get_str(NULL, 10, rad);
349
350 lenmid = strlen(*mid_digits);
351 lenrad = strlen(*rad_digits);
352
353 if (lenmid > lenrad)
354 {
355 /* we will truncate at n or n-1 */
356 good = lenmid - lenrad;
357
358 /* rounding to nearest can add at most 0.5 ulp */
359 /* look at first omitted digit */
360 rem = ((*mid_digits)[good]) - '0';
361 if (rem < 5)
362 rem = rem + 1;
363 else
364 rem = 10 - rem;
365
366 /* and include the leading digit of the radius */
367 rem = rem + ((*rad_digits)[0] - '0') + 1;
368
369 /* if error is <= 1.0 ulp, we get to keep the extra digit */
370 if (rem > 10)
371 good -= 1;
372
373 n = FLINT_MIN(n, good);
374 }
375 else
376 {
377 n = 0;
378 }
379
380 /* todo: avoid recomputing? */
381 flint_free(*rad_digits);
382 }
383
384 /* no accurate digits -- output 0 +/- rad */
385 if (n < 1)
386 {
387 fmpz_add(rad, rad, mid);
388 fmpz_zero(mid);
389 strcpy(*mid_digits, "0"); /* must have space already! */
390 }
391 else
392 {
393 _arb_digits_round_inplace(*mid_digits, &shift, err, n, ARF_RND_NEAR);
394 fmpz_add_ui(mid_exp, exp, shift);
395 fmpz_abs(err, err);
396 fmpz_add(rad, rad, err);
397 }
398
399 /* write radius */
400 if (fmpz_is_zero(rad))
401 {
402 *rad_digits = fmpz_get_str(NULL, 10, rad);
403 fmpz_zero(rad_exp);
404 }
405 else
406 {
407 *rad_digits = fmpz_get_str(NULL, 10, rad);
408 _arb_digits_round_inplace(*rad_digits, &shift, err, RADIUS_DIGITS, ARF_RND_UP);
409 fmpz_add_ui(rad_exp, exp, shift);
410 }
411
412 fmpz_clear(mid);
413 fmpz_clear(rad);
414 fmpz_clear(exp);
415 fmpz_clear(err);
416 }
417
arb_get_str(const arb_t x,slong n,ulong flags)418 char * arb_get_str(const arb_t x, slong n, ulong flags)
419 {
420 char * res;
421 char * mid_digits;
422 char * rad_digits;
423 int negative, more, skip_rad, skip_mid;
424 fmpz_t mid_exp;
425 fmpz_t rad_exp;
426 slong condense;
427
428 if (arb_is_zero(x))
429 {
430 res = flint_malloc(2);
431 strcpy(res, "0");
432 return res;
433 }
434
435 more = flags & ARB_STR_MORE;
436 condense = flags / ARB_STR_CONDENSE;
437
438 if (!arb_is_finite(x))
439 {
440 res = flint_malloc(10);
441
442 if (arf_is_nan(arb_midref(x)))
443 strcpy(res, "nan");
444 else
445 strcpy(res, "[+/- inf]");
446
447 return res;
448 }
449
450 fmpz_init(mid_exp);
451 fmpz_init(rad_exp);
452
453 arb_get_str_parts(&negative, &mid_digits, mid_exp, &rad_digits, rad_exp, x, n, more);
454
455 if ((flags & ARB_STR_NO_RADIUS) && mid_digits[0] == '0')
456 {
457 fmpz_add_ui(rad_exp, rad_exp, strlen(rad_digits));
458
459 res = flint_malloc(fmpz_sizeinbase(rad_exp, 10) + 4);
460 res[0] = '0';
461 res[1] = 'e';
462 if (fmpz_sgn(rad_exp) >= 0)
463 {
464 res[2] = '+';
465 fmpz_get_str(res + 3, 10, rad_exp);
466 }
467 else
468 {
469 fmpz_get_str(res + 2, 10, rad_exp);
470 }
471 }
472 else
473 {
474 skip_mid = mid_digits[0] == '0';
475 skip_rad = (rad_digits[0] == '0') || ((flags & ARB_STR_NO_RADIUS) && !skip_mid);
476
477 _arb_digits_as_float_str(&mid_digits, mid_exp, -4, FLINT_MAX(6, n - 1));
478 _arb_digits_as_float_str(&rad_digits, rad_exp, -2, 2);
479
480 if (skip_rad)
481 {
482 res = flint_malloc(strlen(mid_digits) + 2);
483
484 if (negative)
485 strcpy(res, "-");
486 else
487 strcpy(res, "");
488
489 strcat(res, mid_digits);
490 }
491 else if (skip_mid)
492 {
493 res = flint_malloc(strlen(rad_digits) + 7);
494
495 strcpy(res, "[+/- ");
496 strcat(res, rad_digits);
497 strcat(res, "]");
498 }
499 else
500 {
501 res = flint_malloc(strlen(mid_digits) + strlen(rad_digits) + 9);
502
503 strcpy(res, "[");
504
505 if (negative)
506 strcat(res, "-");
507
508 strcat(res, mid_digits);
509 strcat(res, " +/- ");
510 strcat(res, rad_digits);
511 strcat(res, "]");
512 }
513 }
514
515 if (condense)
516 res = _arb_condense_digits(res, condense);
517
518 flint_free(mid_digits);
519 flint_free(rad_digits);
520
521 fmpz_clear(mid_exp);
522 fmpz_clear(rad_exp);
523
524 return res;
525 }
526
527