Annotation of OpenXM_contrib/gmp/mpfr/get_str.c, Revision 1.1.1.1
1.1 maekawa 1: /* mpfr_get_str -- output a floating-point number to a string
2:
3: Copyright (C) 1999 PolKA project, Inria Lorraine and Loria
4:
5: This file is part of the MPFR Library.
6:
7: The MPFR Library is free software; you can redistribute it and/or modify
8: it under the terms of the GNU Library General Public License as published by
9: the Free Software Foundation; either version 2 of the License, or (at your
10: option) any later version.
11:
12: The MPFR Library is distributed in the hope that it will be useful, but
13: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
15: License for more details.
16:
17: You should have received a copy of the GNU Library General Public License
18: along with the MPFR Library; see the file COPYING.LIB. If not, write to
19: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
20: MA 02111-1307, USA. */
21:
22: #include <math.h>
23: #include <stdio.h>
24: #include <stdlib.h>
25: #include <string.h>
26: #include "gmp.h"
27: #include "gmp-impl.h"
28: #include "longlong.h"
29: #include "mpfr.h"
30:
31: /*
32: Convert op to a string in base 'base' with 'n' digits and writes the
33: mantissa in 'str', the exponent in 'expptr'.
34: The result is rounded wrt 'rnd_mode'.
35:
36: For op = 3.1416 we get str = "31416" and expptr=1.
37: */
38: #if __STDC__
39: char *mpfr_get_str(char *str, mp_exp_t *expptr, int base, size_t n,
40: mpfr_srcptr op, unsigned char rnd_mode)
41: #else
42: char *mpfr_get_str(str, expptr, base, n, op, rnd_mode)
43: char *str;
44: mp_exp_t *expptr;
45: int base;
46: size_t n;
47: mpfr_srcptr op;
48: unsigned char rnd_mode;
49: #endif
50: {
51: double d; long e, q, div, p, err, prec, sh; mpfr_t a, b; mpz_t bz;
52: char *str0; unsigned char rnd1; int f, pow2, ok=0, neg;
53:
54: if (base<2 || 36<base) {
55: fprintf(stderr, "Error: too small or too large base in mpfr_get_str: %d\n",
56: base);
57: exit(1);
58: }
59:
60: neg = (SIGN(op)<0) ? 1 : 0;
61:
62: if (!NOTZERO(op)) {
63: if (str==NULL) str0=str=(*_mp_allocate_func)(neg + n + 2);
64: if (SIGN(op)<0) *str++ = '-';
65: for (f=0;f<n;f++) *str++ = '0';
66: *expptr = 1;
67: return str0;
68: }
69:
70: count_leading_zeros(pow2, (mp_limb_t)base);
71: pow2 = BITS_PER_MP_LIMB - pow2 - 1;
72: if (base != (1<<pow2)) pow2=0;
73: /* if pow2 <> 0, then base = 2^pow2 */
74:
75: /* first determines the exponent */
76: e = EXP(op);
77: d = fabs(mpfr_get_d2(op, 0));
78: /* the absolute value of op is between 1/2*2^e and 2^e */
79: /* the output exponent f is such that base^(f-1) <= |op| < base^f
80: i.e. f = 1 + floor(log(|op|)/log(base))
81: = 1 + floor((log(|m|)+e*log(2))/log(base)) */
82: f = 1 + (int) floor((log(d)+(double)e*log(2.0))/log((double)base));
83: if (n==0) {
84: /* performs exact rounding, i.e. returns y such that for rnd_mode=RNDN
85: for example, we have:
86: y*base^(f-n) <= x*2^(e-p) < (x+1)*2^(e-p) <= (y+1)*base^(f-n)
87: which implies 2^(EXP(op)-PREC(op)) <= base^(f-n)
88: */
89: n = f + (int) ceil(((double)PREC(op)-e)*log(2.0)/log((double)base));
90: }
91: /* now the first n digits of the mantissa are obtained from
92: rnd(op*base^(n-f)) */
93: prec = (long) ceil((double)n*log((double)base)/log(2.0));
94: err = 5;
95: q = prec+err;
96: /* one has to use at least q bits */
97: q = (((q-1)/BITS_PER_MP_LIMB)+1)*BITS_PER_MP_LIMB;
98: mpfr_init2(a,q); mpfr_init2(b,q);
99:
100: do {
101: p = n-f; if ((div=(p<0))) p=-p;
102: rnd1 = rnd_mode;
103: if (div) {
104: /* if div we divide by base^p so we have to invert the rounding mode */
105: switch (rnd1) {
106: case GMP_RNDN: rnd1=GMP_RNDN; break;
107: case GMP_RNDZ: rnd1=GMP_RNDU; break;
108: case GMP_RNDU: rnd1=GMP_RNDZ; break;
109: case GMP_RNDD: rnd1=GMP_RNDZ; break;
110: }
111: }
112:
113: if (pow2) {
114: if (div) mpfr_div_2exp(b, op, pow2*p, rnd_mode);
115: else mpfr_mul_2exp(b, op, pow2*p, rnd_mode);
116: }
117: else {
118: /* compute base^p with q bits and rounding towards zero */
119: mpfr_set_prec(b, q);
120: if (p==0) { mpfr_set(b, op, rnd_mode); mpfr_set_ui(a, 1, rnd_mode); }
121: else {
122: mpfr_set_prec(a, q);
123: mpfr_ui_pow_ui(a, base, p, rnd1);
124: if (div) {
125: mpfr_set_ui(b, 1, rnd_mode);
126: mpfr_div(a, b, a, rnd_mode);
127: }
128: /* now a is an approximation by default of 1/base^(f-n) */
129: mpfr_mul(b, op, a, rnd_mode);
130: }
131: }
132: if (neg) CHANGE_SIGN(b); /* put b positive */
133:
134: if (q>2*prec+BITS_PER_MP_LIMB) {
135: /* happens when just in the middle between two digits */
136: n--; q-=BITS_PER_MP_LIMB;
137: if (n==0) {
138: fprintf(stderr, "cannot determine leading digit\n"); exit(1);
139: }
140: }
141: ok = pow2 || mpfr_can_round(b, q-err, rnd_mode, rnd_mode, prec);
142:
143: if (ok) {
144: if (pow2) {
145: sh = e-PREC(op) + pow2*(n-f); /* error at most 2^e */
146: ok = mpfr_can_round(b, EXP(b)-sh-1, rnd_mode, rnd_mode, n*pow2);
147: }
148: else {
149: /* check that value is the same at distance 2^(e-PREC(op))/base^(f-n)
150: in opposite from rounding direction */
151: if (e>=PREC(op)) mpfr_mul_2exp(a, a, e-PREC(op), rnd_mode);
152: else mpfr_div_2exp(a, a, PREC(op)-e, rnd_mode);
153: if (rnd_mode==GMP_RNDN) {
154: mpfr_div_2exp(a, a, 2, rnd_mode);
155: mpfr_sub(b, b, a, rnd_mode); /* b - a/2 */
156: mpfr_mul_2exp(a, a, 2, rnd_mode);
157: mpfr_add(a, b, a, rnd_mode); /* b + a/2 */
158: }
159: else if ((rnd_mode==GMP_RNDU && neg==0) || (rnd_mode==GMP_RNDD && neg))
160: mpfr_sub(a, b, a, rnd_mode);
161: else mpfr_add(a, b, a, rnd_mode);
162: /* check that a and b are rounded similarly */
163: prec=EXP(b);
164: if (EXP(a) != prec) ok=0;
165: else {
166: mpfr_round(b, rnd_mode, prec);
167: mpfr_round(a, rnd_mode, prec);
168: if (mpfr_cmp(a, b)) ok=0;
169: }
170: }
171: if (ok==0) { /* n is too large */
172: n--;
173: if (n==0) {
174: fprintf(stderr, "cannot determine leading digit\n"); exit(1);
175: }
176: q -= BITS_PER_MP_LIMB;
177: }
178: }
179: } while (ok==0 && (q+=BITS_PER_MP_LIMB) );
180: if (neg)
181: switch (rnd_mode) {
182: case GMP_RNDU: rnd_mode=GMP_RNDZ; break;
183: case GMP_RNDD: rnd_mode=GMP_RNDU; break;
184: }
185:
186: prec=EXP(b); /* may have changed due to rounding */
187:
188: /* now the mantissa is the integer part of b */
189: mpz_init(bz); q=1+(prec-1)/BITS_PER_MP_LIMB;
190: _mpz_realloc(bz, q);
191: sh = prec%BITS_PER_MP_LIMB;
192: e = 1 + (PREC(b)-1)/BITS_PER_MP_LIMB-q;
193: if (sh) mpn_rshift(PTR(bz), MANT(b)+e, q, BITS_PER_MP_LIMB-sh);
194: else MPN_COPY(PTR(bz), MANT(b)+e, q);
195: bz->_mp_size=q;
196:
197: /* computes the number of characters needed */
198: q = neg + n + 2; /* n+1 may not be enough for 100000... */
199: if (str==NULL) str0=str=(*_mp_allocate_func)(q);
200: if (neg) *str++='-';
201: mpz_get_str(str, base, bz); /* n digits of mantissa */
202: if (strlen(str)==n+1) f++; /* possible due to rounding */
203: *expptr = f;
204: mpfr_clear(a); mpfr_clear(b); mpz_clear(bz);
205: return str0;
206: }
207:
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to get_str.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpfr