/src/gmpy_mpmath.c

  1/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
  2 * gmpy_mpmath.c                                                           *
  3 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
  4 * Python interface to the GMP or MPIR, MPFR, and MPC multiple precision   *
  5 * libraries.                                                              *
  6 *                                                                         *
  7 * Copyright 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007,               *
  8 *           2008, 2009 Alex Martelli                                      *
  9 *                                                                         *
 10 * Copyright 2008, 2009, 2010, 2011, 2012, 2013 Case Van Horsen            *
 11 *                                                                         *
 12 * This file is part of GMPY2.                                             *
 13 *                                                                         *
 14 * GMPY2 is free software: you can redistribute it and/or modify it under  *
 15 * the terms of the GNU Lesser General Public License as published by the  *
 16 * Free Software Foundation, either version 3 of the License, or (at your  *
 17 * option) any later version.                                              *
 18 *                                                                         *
 19 * GMPY2 is distributed in the hope that it will be useful, but WITHOUT    *
 20 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or   *
 21 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public    *
 22 * License for more details.                                               *
 23 *                                                                         *
 24 * You should have received a copy of the GNU Lesser General Public        *
 25 * License along with GMPY2; if not, see <http://www.gnu.org/licenses/>    *
 26 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 27
 28
 29/* Internal helper function for mpmath. */
 30
 31static PyObject *
 32mpmath_build_mpf(long sign, PympzObject *man, PyObject *exp, mpir_si bc)
 33{
 34    PyObject *tup, *tsign, *tbc;
 35
 36    if (!(tup = PyTuple_New(4))) {
 37        Py_DECREF((PyObject*)man);
 38        Py_DECREF(exp);
 39        return NULL;
 40    }
 41
 42    if (!(tsign = PyIntOrLong_FromLong(sign))) {
 43        Py_DECREF((PyObject*)man);
 44        Py_DECREF(exp);
 45        Py_DECREF(tup);
 46        return NULL;
 47    }
 48
 49    if (!(tbc = PyIntOrLong_FromSI(bc))) {
 50        Py_DECREF((PyObject*)man);
 51        Py_DECREF(exp);
 52        Py_DECREF(tup);
 53        Py_DECREF(tsign);
 54        return NULL;
 55    }
 56
 57    PyTuple_SET_ITEM(tup, 0, tsign);
 58    PyTuple_SET_ITEM(tup, 1, (PyObject*)man);
 59    PyTuple_SET_ITEM(tup, 2, (exp)?exp:PyIntOrLong_FromLong(0));
 60    PyTuple_SET_ITEM(tup, 3, tbc);
 61    return tup;
 62}
 63
 64PyDoc_STRVAR(doc_mpmath_normalizeg,
 65"_mpmath_normalize(...): helper function for mpmath.");
 66
 67static PyObject *
 68Pympz_mpmath_normalize(PyObject *self, PyObject *args)
 69{
 70    long sign = 0;
 71    mpir_si bc = 0, prec = 0, shift, zbits, carry = 0;
 72    PyObject *exp = 0, *newexp = 0, *newexp2 = 0, *tmp = 0, *rndstr = 0;
 73    PympzObject *man = 0, *upper = 0, *lower = 0;
 74    char rnd = 0;
 75
 76    if (PyTuple_GET_SIZE(args) == 6) {
 77        /* Need better error-checking here. Under Python 3.0, overflow into
 78           C-long is possible. */
 79        sign = clong_From_Integer(PyTuple_GET_ITEM(args, 0));
 80        man = (PympzObject *)PyTuple_GET_ITEM(args, 1);
 81        exp = PyTuple_GET_ITEM(args, 2);
 82        bc = SI_From_Integer(PyTuple_GET_ITEM(args, 3));
 83        prec = SI_From_Integer(PyTuple_GET_ITEM(args, 4));
 84        rndstr = PyTuple_GET_ITEM(args, 5);
 85        if (PyErr_Occurred()) {
 86            TYPE_ERROR("arguments long, PympzObject*, PyObject*, long, long, char needed");
 87            return NULL;
 88        }
 89    }
 90    else {
 91        TYPE_ERROR("6 arguments required");
 92        return NULL;
 93    }
 94
 95    if (!Pympz_Check(man)) {
 96        TYPE_ERROR("argument is not an mpz");
 97        return NULL;
 98    }
 99
100    /* If rndstr really is a string, extract the first character. */
101    if (Py2or3String_Check(rndstr)) {
102        rnd = Py2or3String_AsString(rndstr)[0];
103    }
104    else {
105        VALUE_ERROR("invalid rounding mode specified");
106        return NULL;
107    }
108
109    /* If the mantissa is 0, return the normalized representation. */
110    if (!mpz_sgn(man->z)) {
111        Py_INCREF((PyObject*)man);
112        return mpmath_build_mpf(0, man, 0, 0);
113    }
114
115    /* if bc <= prec and the number is odd return it */
116    if ((bc <= prec) && mpz_odd_p(man->z)) {
117        Py_INCREF((PyObject*)man);
118        Py_INCREF((PyObject*)exp);
119        return mpmath_build_mpf(sign, man, exp, bc);
120    }
121
122    if (!(upper = (PympzObject*)Pympz_new()) ||
123        !(lower = (PympzObject*)Pympz_new())) {
124        Py_XDECREF((PyObject*)upper);
125        Py_XDECREF((PyObject*)lower);
126    }
127
128    shift = bc - prec;
129    if (shift>0) {
130        switch (rnd) {
131            case 'f':
132                if(sign) {
133                    mpz_cdiv_q_2exp(upper->z, man->z, shift);
134                }
135                else {
136                    mpz_fdiv_q_2exp(upper->z, man->z, shift);
137                }
138                break;
139            case 'c':
140                if(sign) {
141                    mpz_fdiv_q_2exp(upper->z, man->z, shift);
142                }
143                else {
144                    mpz_cdiv_q_2exp(upper->z, man->z, shift);
145                }
146                break;
147            case 'd':
148                mpz_fdiv_q_2exp(upper->z, man->z, shift);
149                break;
150            case 'u':
151                mpz_cdiv_q_2exp(upper->z, man->z, shift);
152                break;
153            case 'n':
154            default:
155                mpz_tdiv_r_2exp(lower->z, man->z, shift);
156                mpz_tdiv_q_2exp(upper->z, man->z, shift);
157                if (mpz_sgn(lower->z)) {
158                    /* lower is not 0 so it must have at least 1 bit set */
159                    if (mpz_sizeinbase(lower->z, 2) == shift) {
160                        /* lower is >= 1/2 */
161                        if (mpz_scan1(lower->z, 0) == shift-1) {
162                            /* lower is exactly 1/2 */
163                            if (mpz_odd_p(upper->z))
164                                carry = 1;
165                        }
166                        else {
167                            carry = 1;
168                        }
169                    }
170                }
171                if (carry)
172                    mpz_add_ui(upper->z, upper->z, 1);
173        }
174
175        if (!(tmp = PyIntOrLong_FromSI(shift))) {
176            Py_DECREF((PyObject*)upper);
177            Py_DECREF((PyObject*)lower);
178            return NULL;
179        }
180
181        if (!(newexp = PyNumber_Add(exp, tmp))) {
182            Py_DECREF((PyObject*)upper);
183            Py_DECREF((PyObject*)lower);
184            Py_DECREF(tmp);
185            return NULL;
186        }
187        Py_DECREF(tmp);
188        bc = prec;
189    }
190    else {
191        mpz_set(upper->z, man->z);
192        newexp = exp;
193        Py_INCREF(newexp);
194    }
195
196    /* Strip trailing 0 bits. */
197    if ((zbits = mpz_scan1(upper->z, 0)))
198        mpz_tdiv_q_2exp(upper->z, upper->z, zbits);
199
200    if (!(tmp = PyIntOrLong_FromSI(zbits))) {
201        Py_DECREF((PyObject*)upper);
202        Py_DECREF((PyObject*)lower);
203        Py_DECREF(newexp);
204        return NULL;
205    }
206    if (!(newexp2 = PyNumber_Add(newexp, tmp))) {
207        Py_DECREF((PyObject*)upper);
208        Py_DECREF((PyObject*)lower);
209        Py_DECREF(tmp);
210        Py_DECREF(newexp);
211        return NULL;
212    }
213    Py_DECREF(newexp);
214    Py_DECREF(tmp);
215
216    bc -= zbits;
217    /* Check if one less than a power of 2 was rounded up. */
218    if (!mpz_cmp_ui(upper->z, 1))
219        bc = 1;
220
221    Py_DECREF((PyObject*)lower);
222    return mpmath_build_mpf(sign, upper, newexp2, bc);
223}
224
225PyDoc_STRVAR(doc_mpmath_createg,
226"_mpmath_create(...): helper function for mpmath.");
227
228static PyObject *
229Pympz_mpmath_create(PyObject *self, PyObject *args)
230{
231    long sign;
232    mpir_si bc, shift, zbits, carry = 0, prec = 0;
233    PyObject *exp = 0, *newexp = 0, *newexp2 = 0, *tmp = 0;
234    PympzObject *man = 0, *upper = 0, *lower = 0;
235
236    const char *rnd = "f";
237
238    if (PyTuple_GET_SIZE(args) < 2) {
239        TYPE_ERROR("mpmath_create() expects 'mpz','int'[,'int','str'] arguments");
240        return NULL;
241    }
242
243    switch (PyTuple_GET_SIZE(args)) {
244        case 4:
245            rnd = Py2or3String_AsString(PyTuple_GET_ITEM(args, 3));
246        case 3:
247            prec = SI_From_Integer(PyTuple_GET_ITEM(args, 2));
248            if (prec == -1 && PyErr_Occurred())
249                return NULL;
250            prec = ABS(prec);
251        case 2:
252            exp = PyTuple_GET_ITEM(args, 1);
253        case 1:
254            man = Pympz_From_Integer(PyTuple_GET_ITEM(args, 0));
255            if (!man) {
256                TYPE_ERROR("mpmath_create() expects 'mpz','int'[,'int','str'] arguments");
257                return NULL;
258            }
259    }
260
261    /* If the mantissa is 0, return the normalized representation. */
262    if (!mpz_sgn(man->z)) {
263        return mpmath_build_mpf(0, man, 0, 0);
264    }
265
266    upper = (PympzObject*)Pympz_new();
267    lower = (PympzObject*)Pympz_new();
268    if (!upper || !lower) {
269        Py_DECREF((PyObject*)man);
270        Py_XDECREF((PyObject*)upper);
271        Py_XDECREF((PyObject*)lower);
272        return NULL;
273    }
274
275    /* Extract sign, make man positive, and set bit count */
276    sign = (mpz_sgn(man->z) == -1);
277    mpz_abs(upper->z, man->z);
278    bc = mpz_sizeinbase(upper->z, 2);
279
280    if (!prec) prec = bc;
281
282    shift = bc - prec;
283    if (shift > 0) {
284        switch (rnd[0]) {
285            case 'f':
286                if(sign) {
287                    mpz_cdiv_q_2exp(upper->z, upper->z, shift);
288                }
289                else {
290                    mpz_fdiv_q_2exp(upper->z, upper->z, shift);
291                }
292                break;
293            case 'c':
294                if(sign) {
295                    mpz_fdiv_q_2exp(upper->z, upper->z, shift);
296                }
297                else {
298                    mpz_cdiv_q_2exp(upper->z, upper->z, shift);
299                }
300                break;
301            case 'd':
302                mpz_fdiv_q_2exp(upper->z, upper->z, shift);
303                break;
304            case 'u':
305                mpz_cdiv_q_2exp(upper->z, upper->z, shift);
306                break;
307            case 'n':
308            default:
309                mpz_tdiv_r_2exp(lower->z, upper->z, shift);
310                mpz_tdiv_q_2exp(upper->z, upper->z, shift);
311                if (mpz_sgn(lower->z)) {
312                    /* lower is not 0 so it must have at least 1 bit set */
313                    if (mpz_sizeinbase(lower->z, 2)==shift) {
314                        /* lower is >= 1/2 */
315                        if (mpz_scan1(lower->z, 0)==shift-1) {
316                            /* lower is exactly 1/2 */
317                            if (mpz_odd_p(upper->z))
318                                carry = 1;
319                        }
320                        else {
321                            carry = 1;
322                        }
323                    }
324                }
325                if (carry)
326                    mpz_add_ui(upper->z, upper->z, 1);
327        }
328        if (!(tmp = PyIntOrLong_FromSI(shift))) {
329            Py_DECREF((PyObject*)upper);
330            Py_DECREF((PyObject*)lower);
331            return NULL;
332        }
333        if (!(newexp = PyNumber_Add(exp, tmp))) {
334            Py_DECREF((PyObject*)man);
335            Py_DECREF((PyObject*)upper);
336            Py_DECREF((PyObject*)lower);
337            Py_DECREF(tmp);
338            return NULL;
339        }
340        Py_DECREF(tmp);
341        bc = prec;
342    }
343    else {
344        newexp = exp;
345        Py_INCREF(newexp);
346    }
347
348    /* Strip trailing 0 bits. */
349    if ((zbits = mpz_scan1(upper->z, 0)))
350        mpz_tdiv_q_2exp(upper->z, upper->z, zbits);
351
352    if (!(tmp = PyIntOrLong_FromSI(zbits))) {
353        Py_DECREF((PyObject*)man);
354        Py_DECREF((PyObject*)upper);
355        Py_DECREF((PyObject*)lower);
356        Py_DECREF(newexp);
357        return NULL;
358    }
359    if (!(newexp2 = PyNumber_Add(newexp, tmp))) {
360        Py_DECREF((PyObject*)man);
361        Py_DECREF((PyObject*)upper);
362        Py_DECREF((PyObject*)lower);
363        Py_DECREF(tmp);
364        Py_DECREF(newexp);
365        return NULL;
366    }
367    Py_DECREF(newexp);
368    Py_DECREF(tmp);
369
370    bc -= zbits;
371    /* Check if one less than a power of 2 was rounded up. */
372    if (!mpz_cmp_ui(upper->z, 1))
373        bc = 1;
374
375    Py_DECREF((PyObject*)lower);
376    Py_DECREF((PyObject*)man);
377    return mpmath_build_mpf(sign, upper, newexp2, bc);
378}