add arb_get_str and arb_printn

This commit is contained in:
Fredrik Johansson 2015-01-27 18:09:38 +01:00
parent 5cc1e0d286
commit 8da850c378
8 changed files with 842 additions and 37 deletions

27
arb.h
View file

@ -189,8 +189,17 @@ void arb_neg_round(arb_t x, const arb_t y, long prec);
void arb_abs(arb_t x, const arb_t y);
void _arb_digits_round_inplace(char * s, mp_bitcnt_t * shift, fmpz_t error, long n, arf_rnd_t rnd);
int arb_set_str(arb_t res, const char * inp, long prec);
#define ARB_STR_MORE 1UL
#define ARB_STR_NO_RADIUS 2UL
#define ARB_STR_CONDENSE 16UL
char * arb_get_str(const arb_t x, long n, ulong flags);
ARB_INLINE void
arb_set_arf(arb_t x, const arf_t y)
{
@ -242,21 +251,11 @@ arb_one(arb_t f)
arb_set_ui(f, 1UL);
}
ARB_INLINE void
arb_print(const arb_t x)
{
arf_print(arb_midref(x));
printf(" +/- ");
mag_print(arb_radref(x));
}
void arb_print(const arb_t x);
ARB_INLINE void
arb_printd(const arb_t x, long digits)
{
arf_printd(arb_midref(x), FLINT_MAX(digits, 1));
printf(" +/- ");
mag_printd(arb_radref(x), 5);
}
void arb_printd(const arb_t x, long digits);
void arb_printn(const arb_t x, long digits, ulong flags);
ARB_INLINE void
arb_mul_2exp_si(arb_t y, const arb_t x, long e)

520
arb/get_str.c Normal file
View file

@ -0,0 +1,520 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2015 Fredrik Johansson
******************************************************************************/
#include <string.h>
#include <ctype.h>
#include "arb.h"
#define RADIUS_DIGITS 3
char *
_arb_condense_digits(char * s, long n)
{
long i, j, run, out;
char * res;
res = flint_malloc(strlen(s) + 128); /* space for some growth */
out = 0;
for (i = 0; s[i] != '\0'; )
{
if (isdigit(s[i]))
{
run = 0;
for (j = 0; isdigit(s[i + j]); j++)
run++;
if (run > 3 * n)
{
for (j = 0; j < n; j++)
{
res[out] = s[i + j];
out++;
}
out += sprintf(res + out, "{...%ld digits...}", run - 2 * n);
for (j = run - n; j < run; j++)
{
res[out] = s[i + j];
out++;
}
}
else
{
for (j = 0; j < run; j++)
{
res[out] = s[i + j];
out++;
}
}
i += run;
}
else
{
res[out] = s[i];
i++;
out++;
}
}
res[out] = '\0';
res = flint_realloc(res, strlen(res) + 1);
flint_free(s);
return res;
}
/* Format (digits=d, exponent=e) as floating-point or fixed-point.
Reallocates the input and mutates the exponent. */
void
_arb_digits_as_float_str(char ** d, fmpz_t e, long minfix, long maxfix)
{
long i, n, alloc, dotpos;
/* do nothing with 0 or something non-numerical */
if (!((*d)[0] >= '1' && (*d)[0] <= '9'))
return;
n = strlen(*d);
fmpz_add_ui(e, e, n - 1);
/* fixed-point or integer format */
/* we require e < n - 1; otherwise we would have to insert trailing zeros
[todo: could allow e < n, if printing integers without radix point] */
if (fmpz_cmp_si(e, minfix) >= 0 && fmpz_cmp_si(e, maxfix) <= 0 &&
fmpz_cmp_si(e, n - 1) < 0)
{
long exp = *e;
/* 0.000xxx */
if (exp < 0)
{
/* 0. + (-1-exp) zeros + digits + null terminator */
alloc = 2 + (-1-exp) + n + 1;
*d = flint_realloc(*d, alloc);
/* copy in reverse order, including null terminator */
for (i = n; i >= 0; i--)
(*d)[2 + (-1-exp) + i] = (*d)[i];
for (i = 0; i < 2 + (-1-exp); i++)
(*d)[i] = (i == 1) ? '.' : '0';
}
else /* xxx.yyy --- must have dotpos < n - 1 */
{
dotpos = exp + 1;
alloc = n + 2; /* space for . and null terminator */
(*d) = flint_realloc(*d, alloc);
/* copy fractional part in reverse order, including null */
for (i = n; i >= dotpos; i--)
(*d)[i + 1] = (*d)[i];
(*d)[dotpos] = '.';
}
}
else
{
/* format as xe+zzz or x.yyye+zzz */
alloc = n + 1 + 2 + fmpz_sizeinbase(e, 10) + 1;
*d = flint_realloc(*d, alloc);
/* insert . */
if (n > 1)
{
/* copy fractional part in reverse order */
for (i = n; i >= 1; i--)
(*d)[i + 1] = (*d)[i];
(*d)[1] = '.';
}
(*d)[n + (n > 1)] = 'e';
if (fmpz_sgn(e) >= 0)
{
(*d)[n + (n > 1) + 1] = '+';
}
else
{
(*d)[n + (n > 1) + 1] = '-';
fmpz_neg(e, e);
}
fmpz_get_str((*d) + n + (n > 1) + 2, 10, e); /* writes null byte */
}
}
/* Rounds a string of decimal digits (null-terminated).
to length at most n. The rounding mode
can be ARF_RND_DOWN, ARF_RND_UP or ARF_RND_NEAR.
The string is overwritten in-place, truncating it as necessary.
The input should not have a leading sign or leading zero digits,
but can have trailing zero digits.
Computes shift and error such that
int(input) = int(output) * 10^shift + error
exactly.
*/
void
_arb_digits_round_inplace(char * s, mp_bitcnt_t * shift, fmpz_t error, long n, arf_rnd_t rnd)
{
long i, m;
int up;
if (n < 1)
{
printf("_arb_digits_round_inplace: require n >= 1\n");
abort();
}
m = strlen(s);
if (m <= n)
{
*shift = 0;
fmpz_zero(error);
return;
}
/* always round down */
if (rnd == ARF_RND_DOWN)
{
up = 0;
}
else if (rnd == ARF_RND_UP) /* round up if tail is nonzero */
{
up = 0;
for (i = n; i < m; i++)
{
if (s[i] != '0')
{
up = 1;
break;
}
}
}
else /* round to nearest (up on tie -- todo: round-to-even?) */
{
up = (s[n] >= '5' && s[n] <= '9');
}
if (!up)
{
/* simply truncate */
fmpz_set_str(error, s + n, 10);
s[n] = '\0';
*shift = m - n;
}
else
{
int digit, borrow, carry;
/* error = 10^(m-n) - s[n:], where s[n:] is nonzero */
/* i.e. 10s complement the truncated digits */
borrow = 0;
for (i = m - 1; i >= n; i--)
{
digit = 10 - (s[i] - '0') - borrow;
if (digit == 10)
{
digit = 0;
borrow = 0;
}
else
{
borrow = 1;
}
s[i] = digit + '0';
}
if (!borrow)
{
printf("expected borrow!\n");
abort();
}
fmpz_set_str(error, s + n, 10);
fmpz_neg(error, error);
/* add 1 ulp to the leading digits */
carry = 1;
for (i = n - 1; i >= 0; i--)
{
digit = (s[i] - '0') + carry;
if (digit > 9)
{
digit = 0;
carry = 1;
}
else
{
carry = 0;
}
s[i] = digit + '0';
}
/* carry-out -- only possible if we started with all 9s,
so now the rest will be 0s which we don't have to shift explicitly */
if (carry)
{
s[0] = '1';
*shift = m - n + 1;
}
else
{
*shift = m - n;
}
s[n] = '\0'; /* truncate */
}
}
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, long n, int more)
{
fmpz_t mid, rad, exp, err;
long good;
mp_bitcnt_t shift;
if (!arb_is_finite(x))
{
*negative = 0;
fmpz_zero(mid_exp);
*mid_digits = flint_malloc(4);
if (arf_is_nan(arb_midref(x)))
strcpy(*mid_digits, "nan");
else
strcpy(*mid_digits, "0");
fmpz_zero(rad_exp);
*rad_digits = flint_malloc(4);
strcpy(*rad_digits, "inf");
return;
}
fmpz_init(mid);
fmpz_init(rad);
fmpz_init(exp);
fmpz_init(err);
/* heuristic part */
if (!more)
{
good = arb_rel_accuracy_bits(x) * 0.30102999566398119521 + 2;
n = FLINT_MIN(n, good);
}
arb_get_fmpz_mid_rad_10exp(mid, rad, exp, x, FLINT_MAX(n, 1));
*negative = arf_sgn(arb_midref(x)) < 0;
fmpz_abs(mid, mid);
*mid_digits = fmpz_get_str(NULL, 10, mid);
*rad_digits = NULL;
/* Truncate further so that 1 ulp error can be guaranteed (rigorous part)
Note: mid cannot be zero here if n >= 1 and rad != 0. */
if (n >= 1 && !(more || fmpz_is_zero(rad)))
{
long lenmid, lenrad, rem;
*rad_digits = fmpz_get_str(NULL, 10, rad);
lenmid = strlen(*mid_digits);
lenrad = strlen(*rad_digits);
if (lenmid > lenrad)
{
/* we will truncate at n or n-1 */
good = lenmid - lenrad;
/* rounding to nearest can add at most 0.5 ulp */
/* look at first omitted digit */
rem = ((*mid_digits)[good]) - '0';
if (rem < 5)
rem = rem + 1;
else
rem = 10 - rem;
/* and include the leading digit of the radius */
rem = rem + ((*rad_digits)[0] - '0') + 1;
/* if error is <= 1.0 ulp, we get to keep the extra digit */
if (rem > 10)
good -= 1;
n = FLINT_MIN(n, good);
}
else
{
n = 0;
}
/* todo: avoid recomputing? */
flint_free(*rad_digits);
}
/* no accurate digits -- output 0 +/- rad */
if (n < 1)
{
fmpz_add(rad, rad, mid);
fmpz_zero(mid);
strcpy(*mid_digits, "0"); /* must have space already! */
}
else
{
_arb_digits_round_inplace(*mid_digits, &shift, err, n, ARF_RND_NEAR);
fmpz_add_ui(mid_exp, exp, shift);
fmpz_abs(err, err);
fmpz_add(rad, rad, err);
}
/* write radius */
if (fmpz_is_zero(rad))
{
*rad_digits = fmpz_get_str(NULL, 10, rad);
fmpz_zero(rad_exp);
}
else
{
*rad_digits = fmpz_get_str(NULL, 10, rad);
_arb_digits_round_inplace(*rad_digits, &shift, err, RADIUS_DIGITS, ARF_RND_UP);
fmpz_add_ui(rad_exp, exp, shift);
}
fmpz_clear(mid);
fmpz_clear(rad);
fmpz_clear(exp);
fmpz_clear(err);
}
char * arb_get_str(const arb_t x, long n, ulong flags)
{
char * res;
char * mid_digits;
char * rad_digits;
int negative, more, skip_rad, skip_mid;
fmpz_t mid_exp;
fmpz_t rad_exp;
long condense;
if (arb_is_zero(x))
{
res = flint_malloc(2);
strcpy(res, "0");
return res;
}
more = flags & ARB_STR_MORE;
condense = flags / ARB_STR_CONDENSE;
if (!arb_is_finite(x))
{
res = flint_malloc(10);
if (arf_is_nan(arb_midref(x)))
strcpy(res, "nan");
else
strcpy(res, "[+/- inf]");
return res;
}
fmpz_init(mid_exp);
fmpz_init(rad_exp);
arb_get_str_parts(&negative, &mid_digits, mid_exp, &rad_digits, rad_exp, x, n, more);
skip_rad = (rad_digits[0] == '0') || (flags & ARB_STR_NO_RADIUS);
skip_mid = mid_digits[0] == '0';
_arb_digits_as_float_str(&mid_digits, mid_exp, -4, FLINT_MAX(6, n - 1));
_arb_digits_as_float_str(&rad_digits, rad_exp, -2, 2);
if (skip_rad)
{
res = flint_malloc(strlen(mid_digits) + 2);
if (negative)
strcpy(res, "-");
else
strcpy(res, "");
strcat(res, mid_digits);
}
else if (skip_mid)
{
res = flint_malloc(strlen(rad_digits) + 7);
strcpy(res, "[+/- ");
strcat(res, rad_digits);
strcat(res, "]");
}
else
{
res = flint_malloc(strlen(mid_digits) + strlen(rad_digits) + 9);
strcpy(res, "[");
if (negative)
strcat(res, "-");
strcat(res, mid_digits);
strcat(res, " +/- ");
strcat(res, rad_digits);
strcat(res, "]");
}
if (condense)
res = _arb_condense_digits(res, condense);
flint_free(mid_digits);
flint_free(rad_digits);
fmpz_clear(mid_exp);
fmpz_clear(rad_exp);
return res;
}

51
arb/print.c Normal file
View file

@ -0,0 +1,51 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2015 Fredrik Johansson
******************************************************************************/
#include "arb.h"
void
arb_print(const arb_t x)
{
arf_print(arb_midref(x));
printf(" +/- ");
mag_print(arb_radref(x));
}
void
arb_printd(const arb_t x, long digits)
{
arf_printd(arb_midref(x), FLINT_MAX(digits, 1));
printf(" +/- ");
mag_printd(arb_radref(x), 5);
}
void
arb_printn(const arb_t x, long digits, ulong flags)
{
char * s = arb_get_str(x, digits, flags);
printf("%s", s);
flint_free(s);
}

View file

@ -31,7 +31,6 @@ static int
arb_set_float_str(arb_t res, const char * inp, long prec)
{
char * emarker;
char * fracmarker;
char * buf;
int error;
long i;

View file

@ -0,0 +1,107 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2015 Fredrik Johansson
******************************************************************************/
#include <string.h>
#include "arb.h"
int main()
{
flint_rand_t state;
printf("digits_round_inplace....");
fflush(stdout);
flint_randinit(state);
{
char s[30];
long i, j, len, n;
mp_bitcnt_t shift;
fmpz_t inp, out, err, t;
arf_rnd_t rnd;
fmpz_init(inp);
fmpz_init(out);
fmpz_init(err);
fmpz_init(t);
for (i = 0; i < 100000; i++)
{
len = 1 + n_randint(state, 20);
n = 1 + n_randint(state, 20);
s[0] = (n_randint(state, 9) + '1');
for (j = 1; j < len; j++)
s[j] = (n_randint(state, 10) + '0');
s[len] = '\0';
fmpz_set_str(inp, s, 10);
switch (n_randint(state, 3))
{
case 0:
rnd = ARF_RND_DOWN;
break;
case 1:
rnd = ARF_RND_UP;
break;
default:
rnd = ARF_RND_NEAR;
break;
}
_arb_digits_round_inplace(s, &shift, err, n, rnd);
fmpz_set_str(out, s, 10);
fmpz_set_ui(t, 10);
fmpz_pow_ui(t, t, shift);
fmpz_mul(t, t, out);
fmpz_add(t, t, err);
if (!fmpz_equal(t, inp) || (rnd == ARF_RND_UP && fmpz_sgn(err) > 0))
{
printf("FAIL!\n");
printf("inp = "); fmpz_print(inp); printf("\n\n");
printf("shift = %ld\n\n", shift);
printf("err = "); fmpz_print(err); printf("\n\n");
printf("out = "); fmpz_print(out); printf("\n\n");
printf(" t = "); fmpz_print(t); printf("\n\n");
abort();
}
}
fmpz_clear(inp);
fmpz_clear(out);
fmpz_clear(err);
fmpz_clear(t);
}
flint_randclear(state);
flint_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}

97
arb/test/t-get_str.c Normal file
View file

@ -0,0 +1,97 @@
/*=============================================================================
This file is part of ARB.
ARB is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
ARB is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with ARB; if not, write to the Free Software
Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
=============================================================================*/
/******************************************************************************
Copyright (C) 2015 Fredrik Johansson
******************************************************************************/
#include <string.h>
#include "arb.h"
int main()
{
flint_rand_t state;
long iter;
printf("get_str....");
fflush(stdout);
flint_randinit(state);
/* just test no crashing... */
for (iter = 0; iter < 10000; iter++)
{
arb_t x;
char * s;
long n;
arb_init(x);
arb_randtest_special(x, state, 1 + n_randint(state, 1000), 1 + n_randint(state, 100));
n = 1 + n_randint(state, 300);
s = arb_get_str(x, n, (n_randint(state, 2) * ARB_STR_MORE)
| (n_randint(state, 2) * ARB_STR_NO_RADIUS)
| (ARB_STR_CONDENSE * n_randint(state, 50)));
flint_free(s);
arb_clear(x);
}
for (iter = 0; iter < 100000; iter++)
{
arb_t x, y;
char * s;
long n, prec;
int conversion_error;
arb_init(x);
arb_init(y);
arb_randtest_special(x, state, 1 + n_randint(state, 1000), 1 + n_randint(state, 100));
arb_randtest_special(y, state, 1 + n_randint(state, 1000), 1 + n_randint(state, 100));
n = 1 + n_randint(state, 300);
prec = 2 + n_randint(state, 1000);
s = arb_get_str(x, n, n_randint(state, 2) * ARB_STR_MORE);
conversion_error = arb_set_str(y, s, prec);
if (conversion_error || !arb_contains(y, x))
{
printf("FAIL (roundtrip) iter = %ld\n", iter);
printf("x = "); arb_printd(x, 50); printf("\n\n");
printf("s = %s", s); printf("\n\n");
printf("y = "); arb_printd(y, 50); printf("\n\n");
abort();
}
flint_free(s);
arb_clear(x);
arb_clear(y);
}
flint_randclear(state);
flint_cleanup();
printf("PASS\n");
return EXIT_SUCCESS;
}

View file

@ -89,23 +89,6 @@ Assignment and rounding
Sets *y* to the value of *x* without rounding.
.. function:: int arb_set_str(arb_t res, const char * inp, long prec)
Sets *res* to the value specified by the human-readable string *inp*.
The input may be a decimal floating-point literal,
such as "25", "0.001", "7e+141" or "-31.4159e-1", and may also consist
of two such literals separated by the symbol "+/-" and optionally
enclosed in brackets, e.g. "[3.25 +/- 0.0001]", or simply
"[+/- 10]" with an implicit zero midpoint.
The output is rounded to *prec* bits, and if the binary-to-decimal
conversion is inexact, the resulting error is added to the radius.
The symbols "inf" and "nan" are recognized (a nan midpoint results in an
indeterminate interval, with infinite radius).
Returns 0 if successful and nonzero if unsuccessful. If unsuccessful,
the result is set to an indeterminate interval
.. function:: void arb_set_fmpz_2exp(arb_t y, const fmpz_t x, const fmpz_t e)
Sets *y* to `x \cdot 2^e`.
@ -124,6 +107,48 @@ Assignment and rounding
Sets *y* to the rational number *x*, rounded to *prec* bits.
.. function:: int arb_set_str(arb_t res, const char * inp, long prec)
Sets *res* to the value specified by the human-readable string *inp*.
The input may be a decimal floating-point literal,
such as "25", "0.001", "7e+141" or "-31.4159e-1", and may also consist
of two such literals separated by the symbol "+/-" and optionally
enclosed in brackets, e.g. "[3.25 +/- 0.0001]", or simply
"[+/- 10]" with an implicit zero midpoint.
The output is rounded to *prec* bits, and if the binary-to-decimal
conversion is inexact, the resulting error is added to the radius.
The symbols "inf" and "nan" are recognized (a nan midpoint results in an
indeterminate interval, with infinite radius).
Returns 0 if successful and nonzero if unsuccessful. If unsuccessful,
the result is set to an indeterminate interval.
.. function:: char * arb_get_str(const arb_t x, long n, ulong flags)
Returns a nice human-readable representation of *x*, with at most *n*
digits of the midpoint printed.
With default flags, the output can be parsed back with :func:`arb_set_str`,
and this is guaranteed to produce an interval containing the original
interval *x*.
By default, the output is rounded so that the value given for the
midpoint is correct up to 1 ulp (unit in the last decimal place).
If *ARB_STR_MORE* is added to *flags*, more (possibly incorrect)
digits may be printed.
If *ARB_STR_NO_RADIUS* is added to *flags*, the radius is not
included in the output if at least 1 digit of the midpoint
can be printed.
By adding a multiple *m* of *ARB_STR_CONDENSE* to *flags*, strings
of more than three times *m* consecutive digits are condensed, only
printing the leading and trailing *m* digits along with
brackets indicating the number of digits omitted
(useful when computing values to extremely high precision).
Assignment of special values
-------------------------------------------------------------------------------
@ -165,6 +190,13 @@ Input and output
to compensate for the fact that the binary-to-decimal conversion
of both the midpoint and the radius introduces additional error.
.. function:: void arb_printn(const arb_t x, long digits, ulong flags)
Prints a nice decimal representation of *x*.
By default, the output is guaranteed to be correct to within one unit
in the last digit, and includes an error bound on top of that.
See :func:`arb_get_str` for details.
Random number generation
-------------------------------------------------------------------------------

View file

@ -6,20 +6,20 @@
int main(int argc, char *argv[])
{
arb_t x;
long prec, digits, digits_to_print;
long prec, digits, condense;
if (argc < 2)
{
printf("usage: build/examples/pi digits [digits_to_print = 10]\n");
printf("usage: build/examples/pi digits [condense = 20]\n");
return 1;
}
digits = atol(argv[1]);
if (argc > 2)
digits_to_print = atol(argv[2]);
condense = atol(argv[2]);
else
digits_to_print = 10;
condense = 20;
arb_init(x);
@ -33,7 +33,7 @@ int main(int argc, char *argv[])
SHOW_MEMORY_USAGE
arb_printd(x, digits_to_print);
arb_printn(x, digits, ARB_STR_CONDENSE * condense);
printf("\n");