xref: /dflybsd-src/contrib/mpc/src/get_x.c (revision 04febcfb30580676d3e95f58a16c5137ee478b32)
1*d30dc8cbSJohn Marino /* mpc_get_dc, mpc_get_ldc -- Transform mpc number into C complex number
2*d30dc8cbSJohn Marino    mpc_get_str -- Convert a complex number into a string.
3*d30dc8cbSJohn Marino 
4*d30dc8cbSJohn Marino Copyright (C) 2009, 2010, 2011 INRIA
5*d30dc8cbSJohn Marino 
6*d30dc8cbSJohn Marino This file is part of GNU MPC.
7*d30dc8cbSJohn Marino 
8*d30dc8cbSJohn Marino GNU MPC is free software; you can redistribute it and/or modify it under
9*d30dc8cbSJohn Marino the terms of the GNU Lesser General Public License as published by the
10*d30dc8cbSJohn Marino Free Software Foundation; either version 3 of the License, or (at your
11*d30dc8cbSJohn Marino option) any later version.
12*d30dc8cbSJohn Marino 
13*d30dc8cbSJohn Marino GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY
14*d30dc8cbSJohn Marino WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
15*d30dc8cbSJohn Marino FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
16*d30dc8cbSJohn Marino more details.
17*d30dc8cbSJohn Marino 
18*d30dc8cbSJohn Marino You should have received a copy of the GNU Lesser General Public License
19*d30dc8cbSJohn Marino along with this program. If not, see http://www.gnu.org/licenses/ .
20*d30dc8cbSJohn Marino */
21*d30dc8cbSJohn Marino 
22*d30dc8cbSJohn Marino #include "config.h"
23*d30dc8cbSJohn Marino 
24*d30dc8cbSJohn Marino #ifdef HAVE_COMPLEX_H
25*d30dc8cbSJohn Marino #include <complex.h>
26*d30dc8cbSJohn Marino #endif
27*d30dc8cbSJohn Marino 
28*d30dc8cbSJohn Marino #ifdef HAVE_LOCALE_H
29*d30dc8cbSJohn Marino #include <locale.h>
30*d30dc8cbSJohn Marino #endif
31*d30dc8cbSJohn Marino 
32*d30dc8cbSJohn Marino #include <stdio.h> /* for sprintf, fprintf */
33*d30dc8cbSJohn Marino #include <ctype.h>
34*d30dc8cbSJohn Marino #include <string.h>
35*d30dc8cbSJohn Marino #include "mpc-impl.h"
36*d30dc8cbSJohn Marino 
37*d30dc8cbSJohn Marino #ifdef HAVE_COMPLEX_H
38*d30dc8cbSJohn Marino double _Complex
mpc_get_dc(mpc_srcptr op,mpc_rnd_t rnd)39*d30dc8cbSJohn Marino mpc_get_dc (mpc_srcptr op, mpc_rnd_t rnd) {
40*d30dc8cbSJohn Marino    return I * mpfr_get_d (mpc_imagref (op), MPC_RND_IM (rnd))
41*d30dc8cbSJohn Marino           + mpfr_get_d (mpc_realref (op), MPC_RND_RE (rnd));
42*d30dc8cbSJohn Marino }
43*d30dc8cbSJohn Marino 
44*d30dc8cbSJohn Marino long double _Complex
mpc_get_ldc(mpc_srcptr op,mpc_rnd_t rnd)45*d30dc8cbSJohn Marino mpc_get_ldc (mpc_srcptr op, mpc_rnd_t rnd) {
46*d30dc8cbSJohn Marino    return I * mpfr_get_ld (mpc_imagref (op), MPC_RND_IM (rnd))
47*d30dc8cbSJohn Marino           + mpfr_get_ld (mpc_realref (op), MPC_RND_RE (rnd));
48*d30dc8cbSJohn Marino }
49*d30dc8cbSJohn Marino #endif
50*d30dc8cbSJohn Marino 
51*d30dc8cbSJohn Marino 
52*d30dc8cbSJohn Marino /* Code for mpc_get_str. The output format is "(real imag)", the decimal point
53*d30dc8cbSJohn Marino    of the locale is used. */
54*d30dc8cbSJohn Marino 
55*d30dc8cbSJohn Marino /* mpfr_prec_t can be either int or long int */
56*d30dc8cbSJohn Marino #if (__GMP_MP_SIZE_T_INT == 1)
57*d30dc8cbSJohn Marino #define MPC_EXP_FORMAT_SPEC "i"
58*d30dc8cbSJohn Marino #elif (__GMP_MP_SIZE_T_INT == 0)
59*d30dc8cbSJohn Marino #define MPC_EXP_FORMAT_SPEC "li"
60*d30dc8cbSJohn Marino #else
61*d30dc8cbSJohn Marino #error "mpfr_exp_t size not supported"
62*d30dc8cbSJohn Marino #endif
63*d30dc8cbSJohn Marino 
64*d30dc8cbSJohn Marino static char *
pretty_zero(mpfr_srcptr zero)65*d30dc8cbSJohn Marino pretty_zero (mpfr_srcptr zero)
66*d30dc8cbSJohn Marino {
67*d30dc8cbSJohn Marino   char *pretty;
68*d30dc8cbSJohn Marino 
69*d30dc8cbSJohn Marino   pretty = mpc_alloc_str (3);
70*d30dc8cbSJohn Marino 
71*d30dc8cbSJohn Marino   pretty[0] = mpfr_signbit (zero) ? '-' : '+';
72*d30dc8cbSJohn Marino   pretty[1] = '0';
73*d30dc8cbSJohn Marino   pretty[2] = '\0';
74*d30dc8cbSJohn Marino 
75*d30dc8cbSJohn Marino   return pretty;
76*d30dc8cbSJohn Marino }
77*d30dc8cbSJohn Marino 
78*d30dc8cbSJohn Marino static char *
prettify(const char * str,const mp_exp_t expo,int base,int special)79*d30dc8cbSJohn Marino prettify (const char *str, const mp_exp_t expo, int base, int special)
80*d30dc8cbSJohn Marino {
81*d30dc8cbSJohn Marino   size_t sz;
82*d30dc8cbSJohn Marino   char *pretty;
83*d30dc8cbSJohn Marino   char *p;
84*d30dc8cbSJohn Marino   const char *s;
85*d30dc8cbSJohn Marino   mp_exp_t x;
86*d30dc8cbSJohn Marino   int sign;
87*d30dc8cbSJohn Marino 
88*d30dc8cbSJohn Marino   sz = strlen (str) + 1; /* + terminal '\0' */
89*d30dc8cbSJohn Marino 
90*d30dc8cbSJohn Marino   if (special)
91*d30dc8cbSJohn Marino     {
92*d30dc8cbSJohn Marino       /* special number: nan or inf */
93*d30dc8cbSJohn Marino       pretty = mpc_alloc_str (sz);
94*d30dc8cbSJohn Marino       strcpy (pretty, str);
95*d30dc8cbSJohn Marino 
96*d30dc8cbSJohn Marino       return pretty;
97*d30dc8cbSJohn Marino     }
98*d30dc8cbSJohn Marino 
99*d30dc8cbSJohn Marino   /* regular number */
100*d30dc8cbSJohn Marino 
101*d30dc8cbSJohn Marino   sign = (str[0] == '-' || str[0] == '+');
102*d30dc8cbSJohn Marino 
103*d30dc8cbSJohn Marino   x = expo - 1; /* expo is the exponent value with decimal point BEFORE
104*d30dc8cbSJohn Marino                    the first digit, we wants decimal point AFTER the first
105*d30dc8cbSJohn Marino                    digit */
106*d30dc8cbSJohn Marino   if (base == 16)
107*d30dc8cbSJohn Marino     x <<= 2; /* the output exponent is a binary exponent */
108*d30dc8cbSJohn Marino 
109*d30dc8cbSJohn Marino   ++sz; /* + decimal point */
110*d30dc8cbSJohn Marino 
111*d30dc8cbSJohn Marino   if (x != 0)
112*d30dc8cbSJohn Marino     {
113*d30dc8cbSJohn Marino       /* augment sz with the size needed for an exponent written in base
114*d30dc8cbSJohn Marino          ten */
115*d30dc8cbSJohn Marino       mp_exp_t xx;
116*d30dc8cbSJohn Marino 
117*d30dc8cbSJohn Marino       sz += 3; /* + exponent char + sign + 1 digit */
118*d30dc8cbSJohn Marino 
119*d30dc8cbSJohn Marino       if (x < 0)
120*d30dc8cbSJohn Marino         {
121*d30dc8cbSJohn Marino           /* avoid overflow when changing sign (assuming that, for the
122*d30dc8cbSJohn Marino              mp_exp_t type, (max value) is greater than (- min value / 10)) */
123*d30dc8cbSJohn Marino           if (x < -10)
124*d30dc8cbSJohn Marino             {
125*d30dc8cbSJohn Marino               xx = - (x / 10);
126*d30dc8cbSJohn Marino               sz++;
127*d30dc8cbSJohn Marino             }
128*d30dc8cbSJohn Marino           else
129*d30dc8cbSJohn Marino             xx = -x;
130*d30dc8cbSJohn Marino         }
131*d30dc8cbSJohn Marino       else
132*d30dc8cbSJohn Marino         xx = x;
133*d30dc8cbSJohn Marino 
134*d30dc8cbSJohn Marino       /* compute sz += floor(log(expo)/log(10)) without using libm
135*d30dc8cbSJohn Marino          functions */
136*d30dc8cbSJohn Marino       while (xx > 9)
137*d30dc8cbSJohn Marino         {
138*d30dc8cbSJohn Marino           sz++;
139*d30dc8cbSJohn Marino           xx /= 10;
140*d30dc8cbSJohn Marino         }
141*d30dc8cbSJohn Marino     }
142*d30dc8cbSJohn Marino 
143*d30dc8cbSJohn Marino   pretty = mpc_alloc_str (sz);
144*d30dc8cbSJohn Marino   p = pretty;
145*d30dc8cbSJohn Marino 
146*d30dc8cbSJohn Marino   /* 1. optional sign plus first digit */
147*d30dc8cbSJohn Marino   s = str;
148*d30dc8cbSJohn Marino   *p++ = *s++;
149*d30dc8cbSJohn Marino   if (sign)
150*d30dc8cbSJohn Marino     *p++ = *s++;
151*d30dc8cbSJohn Marino 
152*d30dc8cbSJohn Marino   /* 2. decimal point */
153*d30dc8cbSJohn Marino #ifdef HAVE_LOCALECONV
154*d30dc8cbSJohn Marino   *p++ = *localeconv ()->decimal_point;
155*d30dc8cbSJohn Marino #else
156*d30dc8cbSJohn Marino   *p++ = '.';
157*d30dc8cbSJohn Marino #endif
158*d30dc8cbSJohn Marino   *p = '\0';
159*d30dc8cbSJohn Marino 
160*d30dc8cbSJohn Marino   /* 3. other significant digits */
161*d30dc8cbSJohn Marino   strcat (pretty, s);
162*d30dc8cbSJohn Marino 
163*d30dc8cbSJohn Marino   /* 4. exponent (in base ten) */
164*d30dc8cbSJohn Marino   if (x == 0)
165*d30dc8cbSJohn Marino     return pretty;
166*d30dc8cbSJohn Marino 
167*d30dc8cbSJohn Marino   p = pretty + strlen (str) + 1;
168*d30dc8cbSJohn Marino 
169*d30dc8cbSJohn Marino   switch (base)
170*d30dc8cbSJohn Marino     {
171*d30dc8cbSJohn Marino     case 10:
172*d30dc8cbSJohn Marino       *p++ = 'e';
173*d30dc8cbSJohn Marino       break;
174*d30dc8cbSJohn Marino     case 2:
175*d30dc8cbSJohn Marino     case 16:
176*d30dc8cbSJohn Marino       *p++ = 'p';
177*d30dc8cbSJohn Marino       break;
178*d30dc8cbSJohn Marino     default:
179*d30dc8cbSJohn Marino       *p++ = '@';
180*d30dc8cbSJohn Marino     }
181*d30dc8cbSJohn Marino 
182*d30dc8cbSJohn Marino   *p = '\0';
183*d30dc8cbSJohn Marino 
184*d30dc8cbSJohn Marino   sprintf (p, "%+"MPC_EXP_FORMAT_SPEC, x);
185*d30dc8cbSJohn Marino 
186*d30dc8cbSJohn Marino   return pretty;
187*d30dc8cbSJohn Marino }
188*d30dc8cbSJohn Marino 
189*d30dc8cbSJohn Marino static char *
get_pretty_str(const int base,const size_t n,mpfr_srcptr x,mpfr_rnd_t rnd)190*d30dc8cbSJohn Marino get_pretty_str (const int base, const size_t n, mpfr_srcptr x, mpfr_rnd_t rnd)
191*d30dc8cbSJohn Marino {
192*d30dc8cbSJohn Marino   mp_exp_t expo;
193*d30dc8cbSJohn Marino   char *ugly;
194*d30dc8cbSJohn Marino   char *pretty;
195*d30dc8cbSJohn Marino 
196*d30dc8cbSJohn Marino   if (mpfr_zero_p (x))
197*d30dc8cbSJohn Marino     return pretty_zero (x);
198*d30dc8cbSJohn Marino 
199*d30dc8cbSJohn Marino   ugly = mpfr_get_str (NULL, &expo, base, n, x, rnd);
200*d30dc8cbSJohn Marino   MPC_ASSERT (ugly != NULL);
201*d30dc8cbSJohn Marino   pretty = prettify (ugly, expo, base, !mpfr_number_p (x));
202*d30dc8cbSJohn Marino   mpfr_free_str (ugly);
203*d30dc8cbSJohn Marino 
204*d30dc8cbSJohn Marino   return pretty;
205*d30dc8cbSJohn Marino }
206*d30dc8cbSJohn Marino 
207*d30dc8cbSJohn Marino char *
mpc_get_str(int base,size_t n,mpc_srcptr op,mpc_rnd_t rnd)208*d30dc8cbSJohn Marino mpc_get_str (int base, size_t n, mpc_srcptr op, mpc_rnd_t rnd)
209*d30dc8cbSJohn Marino {
210*d30dc8cbSJohn Marino   size_t needed_size;
211*d30dc8cbSJohn Marino   char *real_str;
212*d30dc8cbSJohn Marino   char *imag_str;
213*d30dc8cbSJohn Marino   char *complex_str = NULL;
214*d30dc8cbSJohn Marino 
215*d30dc8cbSJohn Marino   if (base < 2 || base > 36)
216*d30dc8cbSJohn Marino     return NULL;
217*d30dc8cbSJohn Marino 
218*d30dc8cbSJohn Marino   real_str = get_pretty_str (base, n, mpc_realref (op), MPC_RND_RE (rnd));
219*d30dc8cbSJohn Marino   imag_str = get_pretty_str (base, n, mpc_imagref (op), MPC_RND_IM (rnd));
220*d30dc8cbSJohn Marino 
221*d30dc8cbSJohn Marino   needed_size = strlen (real_str) + strlen (imag_str) + 4;
222*d30dc8cbSJohn Marino 
223*d30dc8cbSJohn Marino   complex_str = mpc_alloc_str (needed_size);
224*d30dc8cbSJohn Marino MPC_ASSERT (complex_str != NULL);
225*d30dc8cbSJohn Marino 
226*d30dc8cbSJohn Marino   strcpy (complex_str, "(");
227*d30dc8cbSJohn Marino   strcat (complex_str, real_str);
228*d30dc8cbSJohn Marino   strcat (complex_str, " ");
229*d30dc8cbSJohn Marino   strcat (complex_str, imag_str);
230*d30dc8cbSJohn Marino   strcat (complex_str, ")");
231*d30dc8cbSJohn Marino 
232*d30dc8cbSJohn Marino   mpc_free_str (real_str);
233*d30dc8cbSJohn Marino   mpc_free_str (imag_str);
234*d30dc8cbSJohn Marino 
235*d30dc8cbSJohn Marino   return complex_str;
236*d30dc8cbSJohn Marino }
237