1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
2 * gmpy2_math.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, 2014, *
11 * 2015, 2016, 2017, 2018, 2019, 2020 Case Van Horsen *
12 * *
13 * This file is part of GMPY2. *
14 * *
15 * GMPY2 is free software: you can redistribute it and/or modify it under *
16 * the terms of the GNU Lesser General Public License as published by the *
17 * Free Software Foundation, either version 3 of the License, or (at your *
18 * option) any later version. *
19 * *
20 * GMPY2 is distributed in the hope that it will be useful, but WITHOUT *
21 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or *
22 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public *
23 * License for more details. *
24 * *
25 * You should have received a copy of the GNU Lesser General Public *
26 * License along with GMPY2; if not, see <http://www.gnu.org/licenses/> *
27 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
28
29 PyDoc_STRVAR(GMPy_doc_context_sin,
30 "context.sin(x) -> number\n\n"
31 "Return sine of x; x in radians.");
32
33 PyDoc_STRVAR(GMPy_doc_function_sin,
34 "sin(x) -> number\n\n"
35 "Return sine of x; x in radians.");
36
37 GMPY_MPFR_MPC_UNIOP_EX(Sin, sin)
38
39 PyDoc_STRVAR(GMPy_doc_context_cos,
40 "context.cos(x) -> number\n\n"
41 "Return cosine of x; x in radians.");
42
43 PyDoc_STRVAR(GMPy_doc_function_cos,
44 "cos(x) -> number\n\n"
45 "Return cosine of x; x in radians.");
46
47 GMPY_MPFR_MPC_UNIOP_EX(Cos, cos)
48
49 PyDoc_STRVAR(GMPy_doc_context_tan,
50 "context.tan(x) -> number\n\n"
51 "Return tangent of x; x in radians.");
52
53 PyDoc_STRVAR(GMPy_doc_function_tan,
54 "tan(x) -> number\n\n"
55 "Return tangent of x; x in radians.");
56
57 GMPY_MPFR_MPC_UNIOP_EX(Tan, tan)
58
59 PyDoc_STRVAR(GMPy_doc_context_atan,
60 "context.atan(x) -> number\n\n"
61 "Return inverse tangent of x; result in radians.");
62
63 PyDoc_STRVAR(GMPy_doc_function_atan,
64 "atan(x) -> number\n\n"
65 "Return inverse tangent of x; result in radians.");
66
67 GMPY_MPFR_MPC_UNIOP_EX(Atan, atan)
68
69 PyDoc_STRVAR(GMPy_doc_context_sinh,
70 "context.sinh(x) -> number\n\n"
71 "Return hyperbolic sine of x.");
72
73 PyDoc_STRVAR(GMPy_doc_function_sinh,
74 "sinh(x) -> number\n\n"
75 "Return hyperbolic sine of x.");
76
77 GMPY_MPFR_MPC_UNIOP_EX(Sinh, sinh)
78
79 PyDoc_STRVAR(GMPy_doc_context_cosh,
80 "context.cosh(x) -> number\n\n"
81 "Return hyperbolic cosine of x.");
82
83 PyDoc_STRVAR(GMPy_doc_function_cosh,
84 "cosh(x) -> number\n\n"
85 "Return hyperbolic cosine of x.");
86
87 GMPY_MPFR_MPC_UNIOP_EX(Cosh, cosh)
88
89 PyDoc_STRVAR(GMPy_doc_context_tanh,
90 "context.tanh(x) -> number\n\n"
91 "Return hyperbolic tangent of x.");
92
93 PyDoc_STRVAR(GMPy_doc_function_tanh,
94 "tanh(x) -> number\n\n"
95 "Return hyperbolic tangent of x.");
96
97 GMPY_MPFR_MPC_UNIOP_EX(Tanh, tanh)
98
99 PyDoc_STRVAR(GMPy_doc_context_asinh,
100 "context.asinh(x) -> number\n\n"
101 "Return inverse hyperbolic sine of x.");
102
103 PyDoc_STRVAR(GMPy_doc_function_asinh,
104 "asinh(x) -> number\n\n"
105 "Return inverse hyperbolic sine of x.");
106
107 GMPY_MPFR_MPC_UNIOP_EX(Asinh, asinh)
108
109 PyDoc_STRVAR(GMPy_doc_context_acosh,
110 "context.acosh(x) -> number\n\n"
111 "Return inverse hyperbolic cosine of x.");
112
113 PyDoc_STRVAR(GMPy_doc_function_acosh,
114 "acosh(x) -> number\n\n"
115 "Return inverse hyperbolic cosine of x.");
116
117 GMPY_MPFR_MPC_UNIOP_EX(Acosh, acosh)
118
119 /* Section 2:
120 * These functions accept a single argument and return an mpfr result.
121 *
122 * GMPY_MPFR_UNIOP(NAME, FUNC) creates the following functions:
123 * GMPy_Real_NAME(x, context)
124 * GMPy_Number_NAME(x, context)
125 * GMPy_Context_NAME(self, other)
126 * - called with METH_O
127 */
128
129 PyDoc_STRVAR(GMPy_doc_context_sec,
130 "context.sec(x) -> number\n\n"
131 "Return secant of x; x in radians.");
132
133 PyDoc_STRVAR(GMPy_doc_function_sec,
134 "sec(x) -> number\n\n"
135 "Return secant of x; x in radians.");
136
137 GMPY_MPFR_UNIOP_EX(Sec, sec)
138
139 PyDoc_STRVAR(GMPy_doc_context_csc,
140 "context.csc(x) -> number\n\n"
141 "Return cosecant of x; x in radians.");
142
143 PyDoc_STRVAR(GMPy_doc_function_csc,
144 "csc(x) -> number\n\n"
145 "Return cosecant of x; x in radians.");
146
147 GMPY_MPFR_UNIOP_EX(Csc, csc)
148
149 PyDoc_STRVAR(GMPy_doc_context_cot,
150 "context.cot(x) -> number\n\n"
151 "Return cotangent of x; x in radians.");
152
153 PyDoc_STRVAR(GMPy_doc_function_cot,
154 "cot(x) -> number\n\n"
155 "Return cotangent of x; x in radians.");
156
157 GMPY_MPFR_UNIOP_EX(Cot, cot)
158
159 PyDoc_STRVAR(GMPy_doc_context_sech,
160 "context.sech(x) -> number\n\n"
161 "Return hyperbolic secant of x.");
162
163 PyDoc_STRVAR(GMPy_doc_function_sech,
164 "sech(x) -> number\n\n"
165 "Return hyperbolic secant of x.");
166
167 GMPY_MPFR_UNIOP_EX(Sech, sech)
168
169 PyDoc_STRVAR(GMPy_doc_context_csch,
170 "context.csch(x) -> number\n\n"
171 "Return hyperbolic cosecant of x.");
172
173 PyDoc_STRVAR(GMPy_doc_function_csch,
174 "csch(x) -> number\n\n"
175 "Return hyperbolic cosecant of x.");
176
177 GMPY_MPFR_UNIOP_EX(Csch, csch)
178
179 PyDoc_STRVAR(GMPy_doc_context_coth,
180 "context.coth(x) -> number\n\n"
181 "Return hyperbolic cotangent of x.");
182
183 PyDoc_STRVAR(GMPy_doc_function_coth,
184 "coth(x) -> number\n\n"
185 "Return hyperbolic cotangent of x.");
186
187 GMPY_MPFR_UNIOP_EX(Coth, coth)
188
189 PyDoc_STRVAR(GMPy_doc_context_rec_sqrt,
190 "context.rec_sqrt(x) -> number\n\n"
191 "Return the reciprocal of the square root of x.");
192
193 PyDoc_STRVAR(GMPy_doc_function_rec_sqrt,
194 "rec_sqrt(x) -> number\n\n"
195 "Return the reciprocal of the square root of x.");
196
197 GMPY_MPFR_UNIOP_EX(RecSqrt, rec_sqrt)
198
199 PyDoc_STRVAR(GMPy_doc_context_rint,
200 "context.rint(x) -> number\n\n"
201 "Return x rounded to the nearest integer using the context rounding\n"
202 "mode.");
203
204 PyDoc_STRVAR(GMPy_doc_function_rint,
205 "rint(x) -> number\n\n"
206 "Return x rounded to the nearest integer using the current rounding\n"
207 "mode.");
208
209 GMPY_MPFR_UNIOP_EX(Rint, rint)
210
211 PyDoc_STRVAR(GMPy_doc_context_rint_ceil,
212 "context.rint_ceil(x) -> number\n\n"
213 "Return x rounded to the nearest integer by first rounding to the\n"
214 "next higher or equal integer and then, if needed, using the context\n"
215 "rounding mode.");
216
217 PyDoc_STRVAR(GMPy_doc_function_rint_ceil,
218 "rint_ceil(x) -> number\n\n"
219 "Return x rounded to the nearest integer by first rounding to the\n"
220 "next higher or equal integer and then, if needed, using the current\n"
221 "rounding mode.");
222
223 GMPY_MPFR_UNIOP_EX(RintCeil, rint_ceil)
224
225 PyDoc_STRVAR(GMPy_doc_context_rint_floor,
226 "context.rint_floor(x) -> number\n\n"
227 "Return x rounded to the nearest integer by first rounding to the\n"
228 "next lower or equal integer and then, if needed, using the context\n"
229 "rounding mode.");
230
231 PyDoc_STRVAR(GMPy_doc_function_rint_floor,
232 "rint_floor(x) -> number\n\n"
233 "Return x rounded to the nearest integer by first rounding to the\n"
234 "next lower or equal integer and then, if needed, using the current\n"
235 "rounding mode.");
236
237 GMPY_MPFR_UNIOP_EX(RintFloor, rint_floor)
238
239 PyDoc_STRVAR(GMPy_doc_context_rint_round,
240 "context.rint_round(x) -> number\n\n"
241 "Return x rounded to the nearest integer by first rounding to the\n"
242 "nearest integer (ties away from 0) and then, if needed, using\n"
243 "the context rounding mode.");
244
245 PyDoc_STRVAR(GMPy_doc_function_rint_round,
246 "rint_round(x) -> number\n\n"
247 "Return x rounded to the nearest integer by first rounding to the\n"
248 "nearest integer (ties away from 0) and then, if needed, using\n"
249 "the current rounding mode.");
250
251 GMPY_MPFR_UNIOP_EX(RintRound, rint_round)
252
253 PyDoc_STRVAR(GMPy_doc_context_rint_trunc,
254 "context.rint_trunc(x) -> number\n\n"
255 "Return x rounded to the nearest integer by first rounding towards\n"
256 "zero and then, if needed, using the context rounding mode.");
257
258 PyDoc_STRVAR(GMPy_doc_function_rint_trunc,
259 "rint_trunc(x) -> number\n\n"
260 "Return x rounded to the nearest integer by first rounding towards\n"
261 "zero and then, if needed, using the current rounding mode.");
262
263 GMPY_MPFR_UNIOP_EX(RintTrunc, rint_trunc)
264
265 PyDoc_STRVAR(GMPy_doc_context_frac,
266 "context.frac(x) -> number\n\n"
267 "Return fractional part of x.");
268
269 PyDoc_STRVAR(GMPy_doc_function_frac,
270 "frac(x) -> number\n\n"
271 "Return fractional part of x.");
272
273 GMPY_MPFR_UNIOP_EX(Frac, frac)
274
275 PyDoc_STRVAR(GMPy_doc_context_cbrt,
276 "context.cbrt(x) -> number\n\n"
277 "Return the cube root of x.");
278
279 PyDoc_STRVAR(GMPy_doc_function_cbrt,
280 "cbrt(x) -> number\n\n"
281 "Return the cube root of x.");
282
283 GMPY_MPFR_UNIOP_EX(Cbrt, cbrt)
284
285 PyDoc_STRVAR(GMPy_doc_context_log2,
286 "context.log2(x) -> number\n\n"
287 "Return base-2 logarithm of x.");
288
289 PyDoc_STRVAR(GMPy_doc_function_log2,
290 "log2(x) -> number\n\n"
291 "Return base-2 logarithm of x.");
292
293 GMPY_MPFR_UNIOP_EX(Log2, log2)
294
295 PyDoc_STRVAR(GMPy_doc_context_exp2,
296 "context.exp2(x) -> number\n\n"
297 "Return 2**x.");
298
299 PyDoc_STRVAR(GMPy_doc_function_exp2,
300 "exp2(x) -> number\n\n"
301 "Return 2**x.");
302
303 GMPY_MPFR_UNIOP_EX(Exp2, exp2)
304
305 PyDoc_STRVAR(GMPy_doc_context_exp10,
306 "context.exp10(x) -> number\n\n"
307 "Return 10**x.");
308
309 PyDoc_STRVAR(GMPy_doc_function_exp10,
310 "exp10(x) -> number\n\n"
311 "Return 10**x.");
312
313 GMPY_MPFR_UNIOP_EX(Exp10, exp10)
314
315 PyDoc_STRVAR(GMPy_doc_context_log1p,
316 "context.log1p(x) -> number\n\n"
317 "Return natural logarithm of (1+x).");
318
319 PyDoc_STRVAR(GMPy_doc_function_log1p,
320 "log1p(x) -> number\n\n"
321 "Return natural logarithm of (1+x).");
322
323 GMPY_MPFR_UNIOP_EX(Log1p, log1p)
324
325 PyDoc_STRVAR(GMPy_doc_context_expm1,
326 "context.expm1(x) -> number\n\n"
327 "Return exp(x) - 1.");
328
329 PyDoc_STRVAR(GMPy_doc_function_expm1,
330 "expm1(x) -> number\n\n"
331 "Return exp(x) - 1.");
332
333 GMPY_MPFR_UNIOP_EX(Expm1, expm1)
334
335 PyDoc_STRVAR(GMPy_doc_context_eint,
336 "context.eint(x) -> number\n\n"
337 "Return exponential integral of x.");
338
339 PyDoc_STRVAR(GMPy_doc_function_eint,
340 "eint(x) -> number\n\n"
341 "Return exponential integral of x.");
342
343 GMPY_MPFR_UNIOP_EX(Eint, eint)
344
345 PyDoc_STRVAR(GMPy_doc_context_li2,
346 "context.li2(x) -> number\n\n"
347 "Return real part of dilogarithm of x.");
348
349 PyDoc_STRVAR(GMPy_doc_function_li2,
350 "li2(x) -> number\n\n"
351 "Return real part of dilogarithm of x.");
352
353 GMPY_MPFR_UNIOP_EX(Li2, li2)
354
355 PyDoc_STRVAR(GMPy_doc_context_gamma,
356 "context.gamma(x) -> number\n\n"
357 "Return gamma of x.");
358
359 PyDoc_STRVAR(GMPy_doc_function_gamma,
360 "gamma(x) -> number\n\n"
361 "Return gamma of x.");
362
363 GMPY_MPFR_UNIOP_EX(Gamma, gamma)
364
365 PyDoc_STRVAR(GMPy_doc_context_lngamma,
366 "context.lngamma(x) -> number\n\n"
367 "Return natural logarithm of gamma(x).");
368
369 PyDoc_STRVAR(GMPy_doc_function_lngamma,
370 "lngamma(x) -> number\n\n"
371 "Return natural logarithm of gamma(x).");
372
373 GMPY_MPFR_UNIOP_EX(Lngamma, lngamma)
374
375 PyDoc_STRVAR(GMPy_doc_context_digamma,
376 "context.digamma(x) -> number\n\n"
377 "Return digamma of x.");
378
379 PyDoc_STRVAR(GMPy_doc_function_digamma,
380 "digamma(x) -> number\n\n"
381 "Return digamma of x.");
382
383 GMPY_MPFR_UNIOP_EX(Digamma, digamma)
384
385 PyDoc_STRVAR(GMPy_doc_context_zeta,
386 "context.zeta(x) -> number\n\n"
387 "Return Riemann zeta of x.");
388
389 PyDoc_STRVAR(GMPy_doc_function_zeta,
390 "zeta(x) -> number\n\n"
391 "Return Riemann zeta of x.");
392
393 GMPY_MPFR_UNIOP_EX(Zeta, zeta)
394
395 PyDoc_STRVAR(GMPy_doc_context_erf,
396 "context.erf(x) -> number\n\n"
397 "Return error function of x.");
398
399 PyDoc_STRVAR(GMPy_doc_function_erf,
400 "erf(x) -> number\n\n"
401 "Return error function of x.");
402
403 GMPY_MPFR_UNIOP_EX(Erf, erf)
404
405 PyDoc_STRVAR(GMPy_doc_context_erfc,
406 "context.erfc(x) -> number\n\n"
407 "Return complementary error function of x.");
408
409 PyDoc_STRVAR(GMPy_doc_function_erfc,
410 "erfc(x) -> number\n\n"
411 "Return complementary error function of x.");
412
413 GMPY_MPFR_UNIOP_EX(Erfc, erfc)
414
415 PyDoc_STRVAR(GMPy_doc_context_j0,
416 "context.j0(x) -> number\n\n"
417 "Return first kind Bessel function of order 0 of x.");
418
419 PyDoc_STRVAR(GMPy_doc_function_j0,
420 "j0(x) -> number\n\n"
421 "Return first kind Bessel function of order 0 of x.");
422
423 GMPY_MPFR_UNIOP_EX(J0, j0)
424
425 PyDoc_STRVAR(GMPy_doc_context_j1,
426 "context.j1(x) -> number\n\n"
427 "Return first kind Bessel function of order 1 of x.");
428
429 PyDoc_STRVAR(GMPy_doc_function_j1,
430 "j1(x) -> number\n\n"
431 "Return first kind Bessel function of order 1 of x.");
432
433 GMPY_MPFR_UNIOP_EX(J1, j1)
434
435 PyDoc_STRVAR(GMPy_doc_context_y0,
436 "context.y0(x) -> number\n\n"
437 "Return second kind Bessel function of order 0 of x.");
438
439 PyDoc_STRVAR(GMPy_doc_function_y0,
440 "y0(x) -> number\n\n"
441 "Return second kind Bessel function of order 0 of x.");
442
443 GMPY_MPFR_UNIOP_EX(Y0, y0)
444
445 PyDoc_STRVAR(GMPy_doc_context_y1,
446 "context.y1(x) -> number\n\n"
447 "Return second kind Bessel function of order 1 of x.");
448
449 PyDoc_STRVAR(GMPy_doc_function_y1,
450 "y1(x) -> number\n\n"
451 "Return second kind Bessel function of order 1 of x.");
452
453 GMPY_MPFR_UNIOP_EX(Y1, y1)
454
455 PyDoc_STRVAR(GMPy_doc_context_ai,
456 "context.ai(x) -> number\n\n"
457 "Return Airy function of x.");
458
459 PyDoc_STRVAR(GMPy_doc_function_ai,
460 "ai(x) -> number\n\n"
461 "Return Airy function of x.");
462
463 GMPY_MPFR_UNIOP_EX(Ai, ai)
464
465 /* Section 3:
466 * The following functions may return an mpc result for certain mpfr arguments.
467 * Since the expectional values vary between functions, the 'Real' and 'Complex'
468 * functions do not use macros. However, they do use a macro to create the
469 * higher-level functions.
470 *
471 * GMPY_MPFR_MPC_UNIOP_TEMPLATE(NAME, FUNC) creates the following functions:
472 * GMPy_Number_NAME(x, context)
473 * - assumes GMPy_Real_NAME & GMPy_Complex_NAME exist
474 * GMPy_Context_NAME(self, other)
475 * - called with METH_O
476 */
477
478 PyDoc_STRVAR(GMPy_doc_context_acos,
479 "context.acos(x) -> number\n\n"
480 "Return inverse cosine of x; result in radians.");
481
482 PyDoc_STRVAR(GMPy_doc_function_acos,
483 "acos(x) -> number\n\n"
484 "Return inverse cosine of x; result in radians.");
485
486 static PyObject *
_GMPy_MPFR_Acos(PyObject * x,CTXT_Object * context)487 _GMPy_MPFR_Acos(PyObject *x, CTXT_Object *context)
488 {
489 MPFR_Object *result;
490
491 if (!mpfr_nan_p(MPFR(x)) &&
492 (mpfr_cmp_si(MPFR(x), 1) > 0 || mpfr_cmp_si(MPFR(x), -1) < 0) &&
493 context->ctx.allow_complex
494 ) {
495 return GMPy_Complex_Acos(x, context);
496 }
497
498 if (!(result = GMPy_MPFR_New(0, context))) {
499 return NULL;
500 }
501
502 mpfr_clear_flags();
503
504 result->rc = mpfr_acos(result->f, MPFR(x), GET_MPFR_ROUND(context));
505 _GMPy_MPFR_Cleanup(&result, context);
506 return (PyObject*)result;
507 }
508
509 static PyObject *
GMPy_Real_Acos(PyObject * x,CTXT_Object * context)510 GMPy_Real_Acos(PyObject *x, CTXT_Object *context)
511 {
512 PyObject *result, *tempx;
513
514 CHECK_CONTEXT(context);
515
516 if (!(tempx = (PyObject*)GMPy_MPFR_From_Real(x, 1, context))) {
517 return NULL;
518 }
519
520 result = _GMPy_MPFR_Acos(tempx, context);
521 Py_DECREF(tempx);
522 return result;
523 }
524
525 static PyObject *
_GMPy_MPC_Acos(PyObject * x,CTXT_Object * context)526 _GMPy_MPC_Acos(PyObject *x, CTXT_Object *context)
527 {
528 MPC_Object *result = NULL;
529
530 CHECK_CONTEXT(context);
531
532 if (!(result = GMPy_MPC_New(0, 0, context))) {
533 return NULL;
534 }
535
536 result->rc = mpc_acos(result->c, MPC(x), GET_MPC_ROUND(context));
537 _GMPy_MPC_Cleanup(&result, context);
538 return (PyObject*)result;
539 }
540
541 static PyObject *
GMPy_Complex_Acos(PyObject * x,CTXT_Object * context)542 GMPy_Complex_Acos(PyObject *x, CTXT_Object *context)
543 {
544 PyObject *result, *tempx;
545
546 CHECK_CONTEXT(context);
547
548 if (!(tempx = (PyObject*)GMPy_MPC_From_Complex(x, 1, 1, context))) {
549 return NULL;
550 }
551
552 result = _GMPy_MPC_Acos(tempx, context);
553 Py_DECREF(tempx);
554 return result;
555 }
556
557 GMPY_MPFR_MPC_UNIOP_TEMPLATE_EX(Acos, acos)
558
559 PyDoc_STRVAR(GMPy_doc_context_asin,
560 "context.asin(x) -> number\n\n"
561 "Return inverse sine of x; result in radians.");
562
563 PyDoc_STRVAR(GMPy_doc_function_asin,
564 "asin(x) -> number\n\n"
565 "Return inverse sine of x; result in radians.");
566
567 static PyObject *
_GMPy_MPFR_Asin(PyObject * x,CTXT_Object * context)568 _GMPy_MPFR_Asin(PyObject *x, CTXT_Object *context)
569 {
570 MPFR_Object *result;
571
572 CHECK_CONTEXT(context);
573
574 if (!mpfr_nan_p(MPFR(x)) &&
575 (mpfr_cmp_si(MPFR(x), 1) > 0 || mpfr_cmp_si(MPFR(x), -1) < 0) &&
576 context->ctx.allow_complex
577 ) {
578 return GMPy_Complex_Asin(x, context);
579 }
580
581 if (!(result = GMPy_MPFR_New(0, context))) {
582 return NULL;
583 }
584
585 mpfr_clear_flags();
586
587 result->rc = mpfr_asin(result->f, MPFR(x), GET_MPFR_ROUND(context));
588 _GMPy_MPFR_Cleanup(&result, context);
589 return (PyObject*)result;
590 }
591
592 static PyObject *
GMPy_Real_Asin(PyObject * x,CTXT_Object * context)593 GMPy_Real_Asin(PyObject *x, CTXT_Object *context)
594 {
595 PyObject *result, *tempx;
596
597 CHECK_CONTEXT(context);
598
599 if (!(tempx = (PyObject*)GMPy_MPFR_From_Real(x, 1, context))) {
600 return NULL;
601 }
602
603 result = _GMPy_MPFR_Asin(tempx, context);
604 Py_DECREF(tempx);
605 return result;
606 }
607
608 static PyObject *
_GMPy_MPC_Asin(PyObject * x,CTXT_Object * context)609 _GMPy_MPC_Asin(PyObject *x, CTXT_Object *context)
610 {
611 MPC_Object *result;
612
613 CHECK_CONTEXT(context);
614
615 if (!(result = GMPy_MPC_New(0, 0, context))) {
616 return NULL;
617 }
618
619 result->rc = mpc_asin(result->c, MPC(x), GET_MPC_ROUND(context));
620 _GMPy_MPC_Cleanup(&result, context);
621 return (PyObject*)result;
622 }
623
624 static PyObject *
GMPy_Complex_Asin(PyObject * x,CTXT_Object * context)625 GMPy_Complex_Asin(PyObject *x, CTXT_Object *context)
626 {
627 PyObject *result, *tempx;
628
629 CHECK_CONTEXT(context);
630
631 if (!(tempx = (PyObject*)GMPy_MPC_From_Complex(x, 1, 1, context))) {
632 return NULL;
633 }
634
635 result = _GMPy_MPC_Asin(tempx, context);
636 Py_DECREF(tempx);
637 return result;
638 }
639
640 GMPY_MPFR_MPC_UNIOP_TEMPLATE_EX(Asin, asin)
641
642 PyDoc_STRVAR(GMPy_doc_context_atanh,
643 "context.atanh(x) -> number\n\n"
644 "Return inverse hyperbolic tanget of x.");
645
646 PyDoc_STRVAR(GMPy_doc_function_atanh,
647 "atanh(x) -> number\n\n"
648 "Return inverse hyperbolic tangent of x.");
649
650 static PyObject *
_GMPy_MPFR_Atanh(PyObject * x,CTXT_Object * context)651 _GMPy_MPFR_Atanh(PyObject *x, CTXT_Object *context)
652 {
653 MPFR_Object *result;
654
655 CHECK_CONTEXT(context);
656
657 if (!mpfr_nan_p(MPFR(x)) &&
658 (mpfr_cmp_si(MPFR(x), 1) > 0 || mpfr_cmp_si(MPFR(x), -1) < 0) &&
659 context->ctx.allow_complex
660 ) {
661 return GMPy_Complex_Atanh(x, context);
662 }
663
664 if (!(result = GMPy_MPFR_New(0, context))) {
665 return NULL;
666 }
667
668 mpfr_clear_flags();
669
670 result->rc = mpfr_atanh(result->f, MPFR(x), GET_MPFR_ROUND(context));
671 _GMPy_MPFR_Cleanup(&result, context);
672 return (PyObject*)result;
673 }
674
675 static PyObject *
GMPy_Real_Atanh(PyObject * x,CTXT_Object * context)676 GMPy_Real_Atanh(PyObject *x, CTXT_Object *context)
677 {
678 PyObject *result, *tempx;
679
680 CHECK_CONTEXT(context);
681
682 if (!(tempx = (PyObject*)GMPy_MPFR_From_Real(x, 1, context))) {
683 return NULL;
684 }
685
686 result = _GMPy_MPFR_Atanh(tempx, context);
687 Py_DECREF(tempx);
688 return result;
689 }
690
691 static PyObject *
_GMPy_MPC_Atanh(PyObject * x,CTXT_Object * context)692 _GMPy_MPC_Atanh(PyObject *x, CTXT_Object *context)
693 {
694 MPC_Object *result;
695
696 CHECK_CONTEXT(context);
697
698 if (!(result = GMPy_MPC_New(0, 0, context))) {
699 return NULL;
700 }
701
702 result->rc = mpc_atanh(result->c, MPC(x), GET_MPC_ROUND(context));
703 _GMPy_MPC_Cleanup(&result, context);
704 return (PyObject*)result;
705 }
706
707 static PyObject *
GMPy_Complex_Atanh(PyObject * x,CTXT_Object * context)708 GMPy_Complex_Atanh(PyObject *x, CTXT_Object *context)
709 {
710 PyObject *result, *tempx;
711
712 CHECK_CONTEXT(context);
713
714 if (!(tempx = (PyObject*)GMPy_MPC_From_Complex(x, 1, 1, context))) {
715 return NULL;
716 }
717
718 result = _GMPy_MPC_Atanh(tempx, context);
719 Py_DECREF(tempx);
720 return result;
721 }
722
723 GMPY_MPFR_MPC_UNIOP_TEMPLATE_EX(Atanh, atanh)
724
725 PyDoc_STRVAR(GMPy_doc_function_atan2,
726 "atan2(y, x) -> number\n\n"
727 "Return arc-tangent of (y/x); result in radians.");
728
729 PyDoc_STRVAR(GMPy_doc_context_atan2,
730 "context.atan2(y, x) -> number\n\n"
731 "Return arc-tangent of (y/x); result in radians.");
732
733 GMPY_MPFR_BINOP_EX(Atan2, atan2)
734
735 PyDoc_STRVAR(GMPy_doc_function_hypot,
736 "hypot(x, y) -> number\n\n"
737 "Return square root of (x**2 + y**2).");
738
739 PyDoc_STRVAR(GMPy_doc_context_hypot,
740 "context.hypot(x, y) -> number\n\n"
741 "Return square root of (x**2 + y**2).");
742
GMPY_MPFR_BINOP_EX(Hypot,hypot)743 GMPY_MPFR_BINOP_EX(Hypot, hypot)
744
745 static PyObject *
746 _GMPy_MPFR_Sin_Cos(PyObject *x, CTXT_Object *context)
747 {
748 MPFR_Object *s, *c;
749 PyObject *result;
750 int code;
751
752 CHECK_CONTEXT(context);
753
754 s = GMPy_MPFR_New(0, context);
755 c = GMPy_MPFR_New(0, context);
756 result = PyTuple_New(2);
757 if (!s || !c || !result) {
758 Py_XDECREF((PyObject*)s);
759 Py_XDECREF((PyObject*)c);
760 Py_XDECREF(result);
761 return NULL;
762 }
763
764 mpfr_clear_flags();
765
766 code = mpfr_sin_cos(s->f, c->f, MPFR(x), GET_MPFR_ROUND(context));
767
768 s->rc = code & 0x03;
769 c->rc = code >> 2;
770 if (s->rc == 2) s->rc = -1;
771 if (c->rc == 2) c->rc = -1;
772
773 _GMPy_MPFR_Cleanup(&s, context);
774 _GMPy_MPFR_Cleanup(&c, context);
775
776 if (!s || !c) {
777 Py_XDECREF((PyObject*)s);
778 Py_XDECREF((PyObject*)c);
779 Py_DECREF(result);
780 return NULL;
781 }
782
783 PyTuple_SET_ITEM(result, 0, (PyObject*)s);
784 PyTuple_SET_ITEM(result, 1, (PyObject*)c);
785 return result;
786 }
787
788 static PyObject *
GMPy_Real_Sin_Cos(PyObject * x,CTXT_Object * context)789 GMPy_Real_Sin_Cos(PyObject *x, CTXT_Object *context)
790 {
791 PyObject *result, *tempx;
792
793 CHECK_CONTEXT(context);
794
795 if (!(tempx = (PyObject*)GMPy_MPFR_From_Real(x, 1, context))) {
796 return NULL;
797 }
798
799 result = _GMPy_MPFR_Sin_Cos(tempx, context);
800 Py_DECREF(tempx);
801 return result;
802 }
803
804 static PyObject *
_GMPy_MPC_Sin_Cos(PyObject * x,CTXT_Object * context)805 _GMPy_MPC_Sin_Cos(PyObject *x, CTXT_Object *context)
806 {
807 MPC_Object *s, *c;
808 PyObject *result;
809 int code;
810
811 CHECK_CONTEXT(context);
812
813 s = GMPy_MPC_New(0, 0, context);
814 c = GMPy_MPC_New(0, 0, context);
815 result = PyTuple_New(2);
816 if (!s || !c || !result) {
817 Py_XDECREF((PyObject*)s);
818 Py_XDECREF((PyObject*)c);
819 Py_XDECREF(result);
820 return NULL;
821 }
822
823 code = mpc_sin_cos(s->c, c->c, MPC(x), GET_MPC_ROUND(context), GET_MPC_ROUND(context));
824
825 s->rc = MPC_INEX1(code);
826 c->rc = MPC_INEX2(code);
827
828 _GMPy_MPC_Cleanup(&s, context);
829 _GMPy_MPC_Cleanup(&c, context);
830
831 if (!s || !c) {
832 Py_XDECREF((PyObject*)s);
833 Py_XDECREF((PyObject*)c);
834 Py_XDECREF(result);
835 return NULL;
836 }
837
838 PyTuple_SET_ITEM(result, 0, (PyObject*)s);
839 PyTuple_SET_ITEM(result, 1, (PyObject*)c);
840 return result;
841 }
842
843 static PyObject *
GMPy_Complex_Sin_Cos(PyObject * x,CTXT_Object * context)844 GMPy_Complex_Sin_Cos(PyObject *x, CTXT_Object *context)
845 {
846 PyObject *result, *tempx;
847
848 CHECK_CONTEXT(context);
849
850 if (!(tempx = (PyObject*)GMPy_MPC_From_Complex(x, 1, 1, context))) {
851 return NULL;
852 }
853
854 result = _GMPy_MPC_Sin_Cos(tempx, context);
855 Py_DECREF(tempx);
856 return result;
857 }
858
859 PyDoc_STRVAR(GMPy_doc_context_sin_cos,
860 "context.sin_cos(x) -> (number, number)\n\n"
861 "Return a tuple containing the sine and cosine of x; x in radians.");
862
863 PyDoc_STRVAR(GMPy_doc_function_sin_cos,
864 "sin_cos(x) -> (number, number)\n\n"
865 "Return a tuple containing the sine and cosine of x; x in radians.");
866
GMPY_MPFR_MPC_UNIOP_TEMPLATE_EX(Sin_Cos,sin_cos)867 GMPY_MPFR_MPC_UNIOP_TEMPLATE_EX(Sin_Cos, sin_cos)
868
869 static PyObject *
870 _GMPy_MPFR_Sinh_Cosh(PyObject *x, CTXT_Object *context)
871 {
872 MPFR_Object *s, *c;
873 PyObject *result;
874 int code;
875
876 CHECK_CONTEXT(context);
877
878 s = GMPy_MPFR_New(0, context);
879 c = GMPy_MPFR_New(0, context);
880 result = PyTuple_New(2);
881 if (!s || !c || !result) {
882 Py_XDECREF((PyObject*)s);
883 Py_XDECREF((PyObject*)c);
884 Py_XDECREF(result);
885 return NULL;
886 }
887
888 mpfr_clear_flags();
889
890 code = mpfr_sinh_cosh(s->f, c->f, MPFR(x), GET_MPFR_ROUND(context));
891
892 s->rc = code & 0x03;
893 c->rc = code >> 2;
894 if (s->rc == 2) s->rc = -1;
895 if (c->rc == 2) c->rc = -1;
896
897 _GMPy_MPFR_Cleanup(&s, context);
898 _GMPy_MPFR_Cleanup(&c, context);
899
900 if (!s || !c) {
901 Py_XDECREF((PyObject*)s);
902 Py_XDECREF((PyObject*)c);
903 Py_XDECREF(result);
904 return NULL;
905 }
906
907 PyTuple_SET_ITEM(result, 0, (PyObject*)s);
908 PyTuple_SET_ITEM(result, 1, (PyObject*)c);
909 return result;
910 }
911
912 static PyObject *
GMPy_Real_Sinh_Cosh(PyObject * x,CTXT_Object * context)913 GMPy_Real_Sinh_Cosh(PyObject *x, CTXT_Object *context)
914 {
915 PyObject *result, *tempx;
916
917 CHECK_CONTEXT(context);
918
919 if (!(tempx = (PyObject*)GMPy_MPFR_From_Real(x, 1, context))) {
920 return NULL;
921 }
922
923 result = _GMPy_MPFR_Sinh_Cosh(tempx, context);
924 Py_DECREF(tempx);
925 return result;
926 }
927
928 PyDoc_STRVAR(GMPy_doc_context_sinh_cosh,
929 "context.sinh_cosh(x) -> (number, number)\n\n"
930 "Return a tuple containing the hyperbolic sine and cosine of x.");
931
932 PyDoc_STRVAR(GMPy_doc_function_sinh_cosh,
933 "sinh_cosh(x) -> (number, number)\n\n"
934 "Return a tuple containing the hyperbolic sine and cosine of x.");
935
936 GMPY_MPFR_UNIOP_TEMPLATE_EX(Sinh_Cosh, sinh_cosh)
937
938 PyDoc_STRVAR(GMPy_doc_function_degrees,
939 "degrees(x) -> mpfr\n\n"
940 "Convert angle x from radians to degrees.\n"
941 "Note: In rare cases the result may not be correctly rounded.");
942
943 PyDoc_STRVAR(GMPy_doc_context_degrees,
944 "context.degrees(x) -> mpfr\n\n"
945 "Convert angle x from radians to degrees.\n"
946 "Note: In rare cases the result may not be correctly rounded.");
947
948 static PyObject *
GMPy_Context_Degrees(PyObject * self,PyObject * other)949 GMPy_Context_Degrees(PyObject *self, PyObject *other)
950 {
951 MPFR_Object *result, *tempx, *temp;
952 CTXT_Object *context = NULL;
953
954 if (self && CTXT_Check(self)) {
955 context = (CTXT_Object*)self;
956 }
957 else {
958 CHECK_CONTEXT(context);
959 }
960
961 result = GMPy_MPFR_New(0, context);
962 temp = GMPy_MPFR_New(context->ctx.mpfr_prec + 100, context);
963 tempx = GMPy_MPFR_From_Real(other, 1, context);
964 if (!result || !temp || !tempx) {
965 Py_XDECREF((PyObject*)temp);
966 Py_XDECREF((PyObject*)tempx);
967 Py_XDECREF((PyObject*)result);
968 return NULL;
969 }
970
971 mpfr_const_pi(temp->f, MPFR_RNDN);
972 mpfr_ui_div(temp->f, 180, temp->f, MPFR_RNDN);
973
974 mpfr_clear_flags();
975
976 mpfr_mul(result->f, temp->f, tempx->f, MPFR_RNDN);
977
978 Py_DECREF((PyObject*)temp);
979 Py_DECREF((PyObject*)tempx);
980 _GMPy_MPFR_Cleanup(&result, context);
981 return (PyObject*)result;
982 }
983
984 PyDoc_STRVAR(GMPy_doc_function_radians,
985 "radians(x) -> mpfr\n\n"
986 "Convert angle x from degrees to radians.\n"
987 "Note: In rare cases the result may not be correctly rounded.");
988
989 PyDoc_STRVAR(GMPy_doc_context_radians,
990 "context.radians(x) -> mpfr\n\n"
991 "Convert angle x from degrees to radians.\n"
992 "Note: In rare cases the result may not be correctly rounded.");
993
994 static PyObject *
GMPy_Context_Radians(PyObject * self,PyObject * other)995 GMPy_Context_Radians(PyObject *self, PyObject *other)
996 {
997 MPFR_Object *result, *tempx, *temp;
998 CTXT_Object *context = NULL;
999
1000 if (self && CTXT_Check(self)) {
1001 context = (CTXT_Object*)self;
1002 }
1003 else {
1004 CHECK_CONTEXT(context);
1005 }
1006
1007 result = GMPy_MPFR_New(0, context);
1008 temp = GMPy_MPFR_New(context->ctx.mpfr_prec + 100, context);
1009 tempx = GMPy_MPFR_From_Real(other, 1, context);
1010 if (!result || !temp || !tempx) {
1011 Py_XDECREF((PyObject*)temp);
1012 Py_XDECREF((PyObject*)tempx);
1013 Py_XDECREF((PyObject*)result);
1014 return NULL;
1015 }
1016
1017 mpfr_const_pi(temp->f, MPFR_RNDN);
1018 mpfr_div_ui(temp->f, temp->f, 180, MPFR_RNDN);
1019
1020 mpfr_clear_flags();
1021
1022 mpfr_mul(result->f, tempx->f, temp->f, MPFR_RNDN);
1023
1024 Py_DECREF((PyObject*)temp);
1025 Py_DECREF((PyObject*)tempx);
1026 _GMPy_MPFR_Cleanup(&result, context);
1027 return (PyObject*)result;
1028 }
1029
1030 PyDoc_STRVAR(GMPy_doc_context_log10,
1031 "context.log10(x) -> number\n\n"
1032 "Return the base-10 logarithm of x.");
1033
1034 PyDoc_STRVAR(GMPy_doc_function_log10,
1035 "log10(x) -> number\n\n"
1036 "Return the base-10 logarithm of x.");
1037
1038 GMPY_MPFR_MPC_UNIOP_EX(Log10, log10)
1039
1040 PyDoc_STRVAR(GMPy_doc_context_log,
1041 "context.log(x) -> number\n\n"
1042 "Return the natural logarithm of x.");
1043
1044 PyDoc_STRVAR(GMPy_doc_function_log,
1045 "log(x) -> number\n\n"
1046 "Return the natural logarithm of x.");
1047
1048 GMPY_MPFR_MPC_UNIOP_EX(Log, log)
1049
1050 PyDoc_STRVAR(GMPy_doc_context_exp,
1051 "context.exp(x) -> number\n\n"
1052 "Return the exponential of x.");
1053
1054 PyDoc_STRVAR(GMPy_doc_function_exp,
1055 "exp(x) -> number\n\n"
1056 "Return the exponential of x.");
1057
1058 GMPY_MPFR_MPC_UNIOP_EX(Exp, exp)
1059
1060 PyDoc_STRVAR(GMPy_doc_context_sqrt,
1061 "context.sqrt(x) -> number\n\n"
1062 "Return the square root of x.");
1063
1064 PyDoc_STRVAR(GMPy_doc_function_sqrt,
1065 "sqrt(x) -> number\n\n"
1066 "Return the square root of x.");
1067
1068 static PyObject *
_GMPy_MPFR_Sqrt(PyObject * x,CTXT_Object * context)1069 _GMPy_MPFR_Sqrt(PyObject *x, CTXT_Object *context)
1070 {
1071 MPFR_Object *result;
1072
1073 CHECK_CONTEXT(context);
1074
1075 if (mpfr_sgn(MPFR(x)) < 0 && context->ctx.allow_complex) {
1076 return GMPy_Complex_Sqrt(x, context);
1077 }
1078
1079 if (!(result = GMPy_MPFR_New(0, context))) {
1080 return NULL;
1081 }
1082
1083 mpfr_clear_flags();
1084
1085 result->rc = mpfr_sqrt(result->f, MPFR(x), GET_MPFR_ROUND(context));
1086 _GMPy_MPFR_Cleanup(&result, context);
1087 return (PyObject*)result;
1088 }
1089
1090 static PyObject *
GMPy_Real_Sqrt(PyObject * x,CTXT_Object * context)1091 GMPy_Real_Sqrt(PyObject *x, CTXT_Object *context)
1092 {
1093 PyObject *result, *tempx;
1094
1095 CHECK_CONTEXT(context);
1096
1097 if (!(tempx = (PyObject*)GMPy_MPFR_From_Real(x, 1, context))) {
1098 return NULL;
1099 }
1100
1101 result = _GMPy_MPFR_Sqrt(tempx, context);
1102 Py_DECREF(tempx);
1103 return result;
1104 }
1105
1106 static PyObject *
_GMPy_MPC_Sqrt(PyObject * x,CTXT_Object * context)1107 _GMPy_MPC_Sqrt(PyObject *x, CTXT_Object *context)
1108 {
1109 MPC_Object *result;
1110
1111 CHECK_CONTEXT(context);
1112
1113 if (!(result = GMPy_MPC_New(0, 0, context))) {
1114 return NULL;
1115 }
1116
1117 result->rc = mpc_sqrt(result->c, MPC(x), GET_MPFR_ROUND(context));
1118 _GMPy_MPC_Cleanup(&result, context);
1119 return (PyObject*)result;
1120 }
1121
1122 static PyObject *
GMPy_Complex_Sqrt(PyObject * x,CTXT_Object * context)1123 GMPy_Complex_Sqrt(PyObject *x, CTXT_Object *context)
1124 {
1125 PyObject *result, *tempx;
1126
1127 CHECK_CONTEXT(context);
1128
1129 if (!(tempx = (PyObject*)GMPy_MPC_From_Complex(x, 1, 1, context))) {
1130 return NULL;
1131 }
1132
1133 result = _GMPy_MPC_Sqrt(tempx, context);
1134 Py_DECREF(tempx);
1135 return result;
1136 }
1137
1138 GMPY_MPFR_MPC_UNIOP_TEMPLATE_EX(Sqrt, sqrt)
1139
1140 PyDoc_STRVAR(GMPy_doc_function_root,
1141 "root(x, n) -> mpfr\n\n"
1142 "Return n-th root of x. The result always an 'mpfr'.\n"
1143 "Note: not IEEE 754-2008 compliant; result differs when\n"
1144 "x = -0 and n is even. See rootn().");
1145
1146 PyDoc_STRVAR(GMPy_doc_context_root,
1147 "context.root(x, n) -> mpfr\n\n"
1148 "Return n-th root of x. The result always an 'mpfr'.\n"
1149 "Note: not IEEE 754-2008 compliant; result differs when\n"
1150 "x = -0 and n is even. See rootn().");
1151
1152 PyDoc_STRVAR(GMPy_doc_function_rootn,
1153 "rootn(x, n) -> mpfr\n\n"
1154 "Return n-th root of x. The result always an 'mpfr'.\n"
1155 "Note: this is IEEE 754-2008 compliant version of root().");
1156
1157 PyDoc_STRVAR(GMPy_doc_context_rootn,
1158 "context.rootn(x, n) -> mpfr\n\n"
1159 "Return n-th root of x. The result always an 'mpfr'.\n"
1160 "Note: this is IEEE 754-2008 compliant version of root().");
1161
1162 #if MPFR_VERSION_MAJOR > 3
1163
1164 /* Since mpfr_root is deprecated in MPFR 4, we use mpfr_rootn_ui to
1165 * mimic the behavior of mpfr_root. And the converse will be true when
1166 * using MPFR 3.
1167 */
1168
1169 static PyObject *
GMPy_Real_Rootn(PyObject * x,PyObject * y,CTXT_Object * context)1170 GMPy_Real_Rootn(PyObject *x, PyObject *y, CTXT_Object *context)
1171 {
1172 MPFR_Object *result = NULL, *tempx = NULL;
1173 unsigned long n;
1174
1175 CHECK_CONTEXT(context);
1176
1177 result = GMPy_MPFR_New(0, context);
1178 tempx = GMPy_MPFR_From_Real(x, 1, context);
1179 n = c_ulong_From_Integer(y);
1180
1181 if (!result || !tempx || (n == (unsigned long)(-1) && PyErr_Occurred())) {
1182 Py_XDECREF((PyObject*)tempx);
1183 Py_XDECREF((PyObject*)result);
1184 return NULL;
1185 }
1186
1187 mpfr_clear_flags();
1188 result->rc = mpfr_rootn_ui(result->f, tempx->f, n, GET_MPFR_ROUND(context));
1189 Py_DECREF((PyObject*)tempx);
1190 _GMPy_MPFR_Cleanup(&result, context);
1191 return (PyObject*)result;
1192 }
1193
1194 static PyObject *
GMPy_Real_Root(PyObject * x,PyObject * y,CTXT_Object * context)1195 GMPy_Real_Root(PyObject *x, PyObject *y, CTXT_Object *context)
1196 {
1197 MPFR_Object *result = NULL, *tempx = NULL;
1198 unsigned long n;
1199
1200 CHECK_CONTEXT(context);
1201
1202 result = GMPy_MPFR_New(0, context);
1203 tempx = GMPy_MPFR_From_Real(x, 1, context);
1204 n = c_ulong_From_Integer(y);
1205
1206 if (!result || !tempx || (n == (unsigned long)(-1) && PyErr_Occurred())) {
1207 Py_XDECREF((PyObject*)tempx);
1208 Py_XDECREF((PyObject*)result);
1209 return NULL;
1210 }
1211
1212 mpfr_clear_flags();
1213
1214 /* Mimic the non-compliant IEEE 752-2008 behavior. */
1215
1216 if (mpfr_zero_p(tempx->f)) {
1217 mpfr_set(result->f, tempx->f, GET_MPFR_ROUND(context));
1218 }
1219 else {
1220 result->rc = mpfr_rootn_ui(result->f, tempx->f, n, GET_MPFR_ROUND(context));
1221 }
1222
1223 Py_DECREF((PyObject*)tempx);
1224 _GMPy_MPFR_Cleanup(&result, context);
1225 return (PyObject*)result;
1226 }
1227 #else
1228 static PyObject *
GMPy_Real_Rootn(PyObject * x,PyObject * y,CTXT_Object * context)1229 GMPy_Real_Rootn(PyObject *x, PyObject *y, CTXT_Object *context)
1230 {
1231 MPFR_Object *result = NULL, *tempx = NULL;
1232 unsigned long n;
1233
1234 CHECK_CONTEXT(context);
1235
1236 result = GMPy_MPFR_New(0, context);
1237 tempx = GMPy_MPFR_From_Real(x, 1, context);
1238 n = c_ulong_From_Integer(y);
1239
1240 if (!result || !tempx || (n == (unsigned long)(-1) && PyErr_Occurred())) {
1241 Py_XDECREF((PyObject*)tempx);
1242 Py_XDECREF((PyObject*)result);
1243 return NULL;
1244 }
1245
1246 mpfr_clear_flags();
1247
1248 /* Mimic the compliant IEEE 752-2008 behavior. */
1249
1250 if (mpfr_zero_p(tempx->f) && mpfr_signbit(tempx->f)) {
1251 if ((n & 1)) {
1252 /* Odd, so result is -0. */
1253 mpfr_set_zero(result->f, -1);
1254 }
1255 else {
1256 /* Even, so result is 0. */
1257 mpfr_set_zero(result->f, 1);
1258 }
1259 }
1260 else {
1261 result->rc = mpfr_root(result->f, tempx->f, n, GET_MPFR_ROUND(context));
1262 }
1263 Py_DECREF((PyObject*)tempx);
1264 _GMPy_MPFR_Cleanup(&result, context);
1265 return (PyObject*)result;
1266 }
1267
1268 static PyObject *
GMPy_Real_Root(PyObject * x,PyObject * y,CTXT_Object * context)1269 GMPy_Real_Root(PyObject *x, PyObject *y, CTXT_Object *context)
1270 {
1271 MPFR_Object *result = NULL, *tempx = NULL;
1272 unsigned long n;
1273
1274 CHECK_CONTEXT(context);
1275
1276 result = GMPy_MPFR_New(0, context);
1277 tempx = GMPy_MPFR_From_Real(x, 1, context);
1278 n = c_ulong_From_Integer(y);
1279
1280 if (!result || !tempx || (n == (unsigned long)(-1) && PyErr_Occurred())) {
1281 Py_XDECREF((PyObject*)tempx);
1282 Py_XDECREF((PyObject*)result);
1283 return NULL;
1284 }
1285
1286 mpfr_clear_flags();
1287 result->rc = mpfr_root(result->f, tempx->f, n, GET_MPFR_ROUND(context));
1288 Py_DECREF((PyObject*)tempx);
1289 _GMPy_MPFR_Cleanup(&result, context);
1290 return (PyObject*)result;
1291 }
1292 #endif
1293
1294 static PyObject *
GMPy_Number_Rootn(PyObject * x,PyObject * y,CTXT_Object * context)1295 GMPy_Number_Rootn(PyObject *x, PyObject *y, CTXT_Object *context)
1296 {
1297 if (IS_REAL(x) && PyIntOrLong_Check(y))
1298 return GMPy_Real_Rootn(x, y, context);
1299 TYPE_ERROR("rootn() argument type not supported");
1300 return NULL;
1301 }
1302
1303 static PyObject *
GMPy_Context_Rootn(PyObject * self,PyObject * args)1304 GMPy_Context_Rootn(PyObject *self, PyObject *args)
1305 {
1306 CTXT_Object *context = NULL;
1307 if (PyTuple_GET_SIZE(args) != 2) {
1308 TYPE_ERROR("rootn() requires 2 arguments");
1309 return NULL;
1310 }
1311 if (self && CTXT_Check(self)) {
1312 context = (CTXT_Object*)self;
1313 }
1314 else {
1315 CHECK_CONTEXT(context);
1316 }
1317 return GMPy_Number_Rootn(PyTuple_GET_ITEM(args, 0), PyTuple_GET_ITEM(args, 1), context);
1318 }
1319
1320 static PyObject *
GMPy_Number_Root(PyObject * x,PyObject * y,CTXT_Object * context)1321 GMPy_Number_Root(PyObject *x, PyObject *y, CTXT_Object *context)
1322 {
1323 if (IS_REAL(x) && PyIntOrLong_Check(y))
1324 return GMPy_Real_Root(x, y, context);
1325 TYPE_ERROR("root() argument type not supported");
1326 return NULL;
1327 }
1328
1329 static PyObject *
GMPy_Context_Root(PyObject * self,PyObject * args)1330 GMPy_Context_Root(PyObject *self, PyObject *args)
1331 {
1332 CTXT_Object *context = NULL;
1333 if (PyTuple_GET_SIZE(args) != 2) {
1334 TYPE_ERROR("root() requires 2 arguments");
1335 return NULL;
1336 }
1337 if (self && CTXT_Check(self)) {
1338 context = (CTXT_Object*)self;
1339 }
1340 else {
1341 CHECK_CONTEXT(context);
1342 }
1343 return GMPy_Number_Root(PyTuple_GET_ITEM(args, 0), PyTuple_GET_ITEM(args, 1), context);
1344 }
1345
1346
1347
1348
1349 PyDoc_STRVAR(GMPy_doc_function_jn,
1350 "jn(x,n) -> mpfr\n\n"
1351 "Return the first kind Bessel function of order n of x.");
1352
1353 PyDoc_STRVAR(GMPy_doc_context_jn,
1354 "context.jn(x,n) -> mpfr\n\n"
1355 "Return the first kind Bessel function of order n of x.");
1356
1357 GMPY_MPFR_BINOP_REAL_LONG(Jn, jn)
1358
1359 PyDoc_STRVAR(GMPy_doc_function_yn,
1360 "yn(x,n) -> mpfr\n\n"
1361 "Return the second kind Bessel function of order n of x.");
1362
1363 PyDoc_STRVAR(GMPy_doc_context_yn,
1364 "context.yn(x,n) -> mpfr\n\n"
1365 "Return the second kind Bessel function of order n of x.");
1366
1367 GMPY_MPFR_BINOP_REAL_LONG(Yn, yn)
1368
1369 PyDoc_STRVAR(GMPy_doc_function_agm,
1370 "agm(x, y) -> mpfr\n\n"
1371 "Return arithmetic-geometric mean of x and y.");
1372
1373 PyDoc_STRVAR(GMPy_doc_context_agm,
1374 "context.agm(x, y) -> mpfr\n\n"
1375 "Return arithmetic-geometric mean of x and y.");
1376
1377 GMPY_MPFR_BINOP(AGM, agm)
1378
1379 PyDoc_STRVAR(GMPy_doc_function_maxnum,
1380 "maxnum(x, y) -> mpfr\n\n"
1381 "Return the maximum number of x and y. If x and y are not 'mpfr', they are\n"
1382 "converted to 'mpfr'. The result is rounded to match the current context.\n"
1383 "If only one of x or y is a number, then that number is returned.");
1384
1385 PyDoc_STRVAR(GMPy_doc_context_maxnum,
1386 "context.maxnum(x, y) -> mpfr\n\n"
1387 "Return the maximum number of x and y. If x and y are not 'mpfr', they are\n"
1388 "converted to 'mpfr'. The result is rounded to match the specified context.\n"
1389 "If only one of x or y is a number, then that number is returned.");
1390
1391 GMPY_MPFR_BINOP(Maxnum, max)
1392
1393 PyDoc_STRVAR(GMPy_doc_function_minnum,
1394 "minnum(x, y) -> mpfr\n\n"
1395 "Return the minimum number of x and y. If x and y are not 'mpfr', they are\n"
1396 "converted to 'mpfr'. The result is rounded to match the current context.\n"
1397 "If only one of x or y is a number, then that number is returned.");
1398
1399 PyDoc_STRVAR(GMPy_doc_context_minnum,
1400 "context.minnum(x, y) -> mpfr\n\n"
1401 "Return the minimum number of x and y. If x and y are not 'mpfr', they are\n"
1402 "converted to 'mpfr'. The result is rounded to match the specified context.\n"
1403 "If only one of x or y is a number, then that number is returned.");
1404
1405 GMPY_MPFR_BINOP(Minnum, min)
1406
1407 PyDoc_STRVAR(GMPy_doc_function_remainder,
1408 "remainder(x, y) -> mpfr\n\n"
1409 "Return x - n*y where n is the integer quotient of x/y, rounded to\n"
1410 "the nearest integer and ties rounded to even.");
1411
1412 PyDoc_STRVAR(GMPy_doc_context_remainder,
1413 "context.remainder(x, y) -> mpfr\n\n"
1414 "Return x - n*y where n is the integer quotient of x/y, rounded to\n"
1415 "the nearest integer and ties rounded to even.");
1416
1417 GMPY_MPFR_BINOP(Remainder, remainder)
1418
1419 PyDoc_STRVAR(GMPy_doc_function_fmod,
1420 "fmod(x, y) -> mpfr\n\n"
1421 "Return x - n*y where n is the integer quotient of x/y, rounded to 0.");
1422
1423 PyDoc_STRVAR(GMPy_doc_context_fmod,
1424 "context.fmod(x, y) -> mpfr\n\n"
1425 "Return x - n*y where n is the integer quotient of x/y, rounded to 0.");
1426
1427 GMPY_MPFR_BINOP(Fmod, fmod)
1428
1429 PyDoc_STRVAR(GMPy_doc_function_round2,
1430 "round2(x[, n]) -> mpfr\n\n"
1431 "Return x rounded to n bits. Uses default precision if n is not\n"
1432 "specified. See round_away() to access the mpfr_round() function.");
1433
1434 PyDoc_STRVAR(GMPy_doc_context_round2,
1435 "context.round2(x[, n]) -> mpfr\n\n"
1436 "Return x rounded to n bits. Uses default precision if n is not\n"
1437 "specified. See round_away() to access the mpfr_round() function.");
1438
1439 static PyObject *
GMPy_Real_Round2(PyObject * x,PyObject * y,CTXT_Object * context)1440 GMPy_Real_Round2(PyObject *x, PyObject *y, CTXT_Object *context)
1441 {
1442 MPFR_Object *result, *tempx;
1443 long n;
1444
1445 CHECK_CONTEXT(context);
1446 n = GET_MPFR_PREC(context);
1447
1448 if (y) {
1449 n = PyIntOrLong_AsLong(y);
1450 if ( (n == -1 && PyErr_Occurred()) || n < MPFR_PREC_MIN || n > MPFR_PREC_MAX) {
1451 VALUE_ERROR("invalid precision");
1452 return NULL;
1453 }
1454 }
1455
1456 if (!(tempx = GMPy_MPFR_From_Real(x, 1, context))) {
1457 return NULL;
1458 }
1459 if (!(result = GMPy_MPFR_New(mpfr_get_prec(tempx->f), context))) {
1460 Py_DECREF((PyObject*)tempx);
1461 return NULL;
1462 }
1463
1464 mpfr_set(result->f, tempx->f, GET_MPFR_ROUND(context));
1465 Py_DECREF((PyObject*)tempx);
1466
1467 mpfr_clear_flags();
1468
1469 result->rc = mpfr_prec_round(result->f, n, GET_MPFR_ROUND(context));
1470 _GMPy_MPFR_Cleanup(&result, context);
1471 return (PyObject*)result;
1472 }
1473
1474 static PyObject *
GMPy_Number_Round2(PyObject * x,PyObject * y,CTXT_Object * context)1475 GMPy_Number_Round2(PyObject *x, PyObject *y, CTXT_Object *context)
1476 {
1477 if (IS_REAL(x) && (!y || PyIntOrLong_Check(y)))
1478 return GMPy_Real_Round2(x, y, context);
1479
1480 TYPE_ERROR("round2() argument type not supported");
1481 return NULL;
1482 }
1483
1484 static PyObject *
GMPy_Context_Round2(PyObject * self,PyObject * args)1485 GMPy_Context_Round2(PyObject *self, PyObject *args)
1486 {
1487 CTXT_Object *context = NULL;
1488
1489 if (PyTuple_GET_SIZE(args) < 1 || PyTuple_GET_SIZE(args) > 2) {
1490 TYPE_ERROR("round2() requires 1 or 2 arguments");
1491 return NULL;
1492 }
1493
1494 if (self && CTXT_Check(self)) {
1495 context = (CTXT_Object*)self;
1496 }
1497 else {
1498 CHECK_CONTEXT(context);
1499 }
1500
1501 if (PyTuple_GET_SIZE(args) == 1) {
1502 return GMPy_Number_Round2(PyTuple_GET_ITEM(args, 0), NULL, context);
1503 }
1504 else {
1505 return GMPy_Number_Round2(PyTuple_GET_ITEM(args, 0), PyTuple_GET_ITEM(args, 1), context);
1506 }
1507 }
1508
1509 PyDoc_STRVAR(GMPy_doc_function_reldiff,
1510 "reldiff(x, y) -> mpfr\n\n"
1511 "Return the relative difference between x and y. Result is equal to\n"
1512 "abs(x-y)/x.");
1513
1514 PyDoc_STRVAR(GMPy_doc_context_reldiff,
1515 "context.reldiff(x, y) -> mpfr\n\n"
1516 "Return the relative difference between x and y. Result is equal to\n"
1517 "abs(x-y)/x.");
1518
1519 static PyObject *
GMPy_Real_RelDiff(PyObject * x,PyObject * y,CTXT_Object * context)1520 GMPy_Real_RelDiff(PyObject *x, PyObject *y, CTXT_Object *context)
1521 {
1522 MPFR_Object *tempx, *tempy, *result;
1523
1524 CHECK_CONTEXT(context);
1525
1526 result = GMPy_MPFR_New(0, context);
1527 tempx = GMPy_MPFR_From_Real(x, 1, context);
1528 tempy = GMPy_MPFR_From_Real(y, 1, context);
1529 if (!result || !tempx || !tempy) {
1530 Py_XDECREF((PyObject*)result);
1531 Py_XDECREF((PyObject*)tempx);
1532 Py_XDECREF((PyObject*)tempy);
1533 return NULL;
1534 }
1535
1536 mpfr_clear_flags();
1537
1538 mpfr_reldiff(result->f, tempx->f, tempy->f, GET_MPFR_ROUND(context));
1539 result->rc = 0;
1540 _GMPy_MPFR_Cleanup(&result, context);
1541 Py_DECREF((PyObject*)tempx);
1542 Py_DECREF((PyObject*)tempy);
1543 return (PyObject*)result;
1544 }
1545
1546 GMPY_MPFR_BINOP_TEMPLATE(RelDiff, reldiff)
1547
1548 PyDoc_STRVAR(GMPy_doc_mpfr_ceil_method,
1549 "x.__ceil__() -> mpfr\n\n"
1550 "Return an 'mpfr' that is the smallest integer >= x.");
1551
1552 PyDoc_STRVAR(GMPy_doc_function_ceil,
1553 "ceil(x) ->mpfr\n\n"
1554 "Return an 'mpfr' that is the smallest integer >= x.");
1555
1556 PyDoc_STRVAR(GMPy_doc_context_ceil,
1557 "context.ceil(x) ->mpfr\n\n"
1558 "Return an 'mpfr' that is the smallest integer >= x.");
1559
1560 GMPY_MPFR_UNIOP_NOROUND(Ceil, ceil)
1561
1562 PyDoc_STRVAR(GMPy_doc_mpfr_floor_method,
1563 "x.__floor__() -> mpfr\n\n"
1564 "Return an 'mpfr' that is the smallest integer <= x.");
1565
1566 PyDoc_STRVAR(GMPy_doc_function_floor,
1567 "floor(x) -> mpfr\n\n"
1568 "Return an 'mpfr' that is the smallest integer <= x.");
1569
1570 PyDoc_STRVAR(GMPy_doc_context_floor,
1571 "context.floor(x) -> mpfr\n\n"
1572 "Return an 'mpfr' that is the smallest integer <= x.");
1573
1574 GMPY_MPFR_UNIOP_NOROUND(Floor, floor);
1575
1576 PyDoc_STRVAR(GMPy_doc_mpfr_trunc_method,
1577 "x.__trunc__() -> mpfr\n\n"
1578 "Return an 'mpfr' that is truncated towards 0. Same as\n"
1579 "x.floor() if x>=0 or x.ceil() if x<0.");
1580
1581 PyDoc_STRVAR(GMPy_doc_function_trunc,
1582 "trunc(x) -> mpfr\n\n"
1583 "Return an 'mpfr' that is x truncated towards 0. Same as\n"
1584 "x.floor() if x>=0 or x.ceil() if x<0.");
1585
1586 PyDoc_STRVAR(GMPy_doc_context_trunc,
1587 "context.trunc(x) -> mpfr\n\n"
1588 "Return an 'mpfr' that is x truncated towards 0. Same as\n"
1589 "x.floor() if x>=0 or x.ceil() if x<0.");
1590
1591 GMPY_MPFR_UNIOP_NOROUND(Trunc, trunc)
1592
1593 PyDoc_STRVAR(GMPy_doc_function_round_away,
1594 "round_away(x) -> mpfr\n\n"
1595 "Return an 'mpfr' that is x rounded to the nearest integer,\n"
1596 "with ties rounded away from 0.");
1597
1598 PyDoc_STRVAR(GMPy_doc_context_round_away,
1599 "context.round_away(x) -> mpfr\n\n"
1600 "Return an 'mpfr' that is x rounded to the nearest integer,\n"
1601 "with ties rounded away from 0.");
1602
1603 GMPY_MPFR_UNIOP_NOROUND_NOMETHOD(RoundAway, round)
1604
1605 PyDoc_STRVAR(GMPy_doc_function_modf,
1606 "modf(x) -> (mpfr, mpfr)\n\n"
1607 "Return a tuple containing the integer and fractional portions\n"
1608 "of x.");
1609
1610 PyDoc_STRVAR(GMPy_doc_context_modf,
1611 "context.modf(x) -> (mpfr, mpfr)\n\n"
1612 "Return a tuple containing the integer and fractional portions\n"
1613 "of x.");
1614
1615 static PyObject *
GMPy_Real_Modf(PyObject * x,CTXT_Object * context)1616 GMPy_Real_Modf(PyObject *x, CTXT_Object *context)
1617 {
1618 MPFR_Object *s, *c, *tempx;
1619 PyObject *result;
1620 int code;
1621
1622 CHECK_CONTEXT(context);
1623
1624 tempx = GMPy_MPFR_From_Real(x, 1, context);
1625 s = GMPy_MPFR_New(0, context);
1626 c = GMPy_MPFR_New(0, context);
1627 result = PyTuple_New(2);
1628 if (! tempx || !s || !c || !result) {
1629 Py_XDECREF((PyObject*)tempx);
1630 Py_XDECREF((PyObject*)s);
1631 Py_XDECREF((PyObject*)c);
1632 Py_XDECREF(result);
1633 return NULL;
1634 }
1635
1636 mpfr_clear_flags();
1637
1638 code = mpfr_modf(s->f, c->f, tempx->f, GET_MPFR_ROUND(context));
1639 Py_DECREF((PyObject*)tempx);
1640
1641 s->rc = code & 0x03;
1642 c->rc = code >> 2;
1643 if (s->rc == 2) s->rc = -1;
1644 if (c->rc == 2) c->rc = -1;
1645
1646 _GMPy_MPFR_Cleanup(&s, context);
1647 _GMPy_MPFR_Cleanup(&c, context);
1648
1649 if (!s || !c) {
1650 Py_XDECREF((PyObject*)s);
1651 Py_XDECREF((PyObject*)c);
1652 Py_DECREF(result);
1653 return NULL;
1654 }
1655
1656 PyTuple_SET_ITEM(result, 0, (PyObject*)s);
1657 PyTuple_SET_ITEM(result, 1, (PyObject*)c);
1658 return result;
1659 }
1660
1661 GMPY_MPFR_UNIOP_TEMPLATE(Modf, modf)
1662
1663 PyDoc_STRVAR(GMPy_doc_function_lgamma,
1664 "lgamma(x) -> (mpfr, int)\n\n"
1665 "Return a tuple containing the logarithm of the absolute value of\n"
1666 "gamma(x) and the sign of gamma(x)");
1667
1668 PyDoc_STRVAR(GMPy_doc_context_lgamma,
1669 "context.lgamma(x) -> (mpfr, int)\n\n"
1670 "Return a tuple containing the logarithm of the absolute value of\n"
1671 "gamma(x) and the sign of gamma(x)");
1672
1673 static PyObject *
GMPy_Real_Lgamma(PyObject * x,CTXT_Object * context)1674 GMPy_Real_Lgamma(PyObject *x, CTXT_Object *context)
1675 {
1676 PyObject *result;
1677 MPFR_Object *value, *tempx;
1678 int signp = 0;
1679
1680 CHECK_CONTEXT(context)
1681
1682 tempx = GMPy_MPFR_From_Real(x, 1, context);
1683 value = GMPy_MPFR_New(0, context);
1684 result = PyTuple_New(2);
1685 if (!tempx || !value || !result) {
1686 Py_XDECREF((PyObject*)tempx);
1687 Py_XDECREF((PyObject*)value);
1688 Py_XDECREF(result);
1689 return NULL;
1690 }
1691
1692 mpfr_clear_flags();
1693
1694 value->rc = mpfr_lgamma(value->f, &signp, tempx->f, GET_MPFR_ROUND(context));
1695 Py_DECREF((PyObject*)tempx);
1696
1697 _GMPy_MPFR_Cleanup(&value, context);
1698
1699 if (!value) {
1700 Py_DECREF(result);
1701 return NULL;
1702 }
1703
1704 PyTuple_SET_ITEM(result, 0, (PyObject*)value);
1705 PyTuple_SET_ITEM(result, 1, PyIntOrLong_FromLong((long)signp));
1706 return result;
1707 }
1708
1709 GMPY_MPFR_UNIOP_TEMPLATE(Lgamma, lgamma)
1710
1711 PyDoc_STRVAR(GMPy_doc_function_remquo,
1712 "remquo(x, y) -> (mpfr, int)\n\n"
1713 "Return a tuple containing the remainder(x,y) and the low bits of the\n"
1714 "quotient.");
1715
1716 PyDoc_STRVAR(GMPy_doc_context_remquo,
1717 "context.remquo(x, y) -> (mpfr, int)\n\n"
1718 "Return a tuple containing the remainder(x,y) and the low bits of the\n"
1719 "quotient.");
1720
1721 static PyObject *
GMPy_Real_RemQuo(PyObject * x,PyObject * y,CTXT_Object * context)1722 GMPy_Real_RemQuo(PyObject *x, PyObject *y, CTXT_Object *context)
1723 {
1724 PyObject *result;
1725 MPFR_Object *value, *tempx, *tempy;
1726 long quobits = 0;
1727
1728 CHECK_CONTEXT(context);
1729
1730 value = GMPy_MPFR_New(0, context);
1731 tempx = GMPy_MPFR_From_Real(x, 1, context);
1732 tempy = GMPy_MPFR_From_Real(y, 1, context);
1733 result = PyTuple_New(2);
1734 if (!value || !tempx || !tempx || !result) {
1735 Py_XDECREF((PyObject*)tempx);
1736 Py_XDECREF((PyObject*)tempy);
1737 Py_XDECREF((PyObject*)value);
1738 Py_XDECREF(result);
1739 return NULL;
1740 }
1741
1742 mpfr_clear_flags();
1743
1744 value->rc = mpfr_remquo(value->f, &quobits, tempx->f, tempy->f, GET_MPFR_ROUND(context));
1745 Py_DECREF((PyObject*)tempx);
1746 Py_DECREF((PyObject*)tempy);
1747 _GMPy_MPFR_Cleanup(&value, context);
1748
1749 PyTuple_SET_ITEM(result, 0, (PyObject*)value);
1750 PyTuple_SET_ITEM(result, 1, PyIntOrLong_FromLong(quobits));
1751 return result;
1752 }
1753
1754 GMPY_MPFR_BINOP_TEMPLATE(RemQuo, remquo);
1755
1756 PyDoc_STRVAR(GMPy_doc_function_frexp,
1757 "frexp(x) -> (int, mpfr)\n\n"
1758 "Return a tuple containing the exponent and mantissa of x.");
1759
1760 PyDoc_STRVAR(GMPy_doc_context_frexp,
1761 "context.frexp(x) -> (int, mpfr)\n\n"
1762 "Return a tuple containing the exponent and mantissa of x.");
1763
1764 static PyObject *
GMPy_Real_Frexp(PyObject * x,CTXT_Object * context)1765 GMPy_Real_Frexp(PyObject *x, CTXT_Object *context)
1766 {
1767 PyObject *result;
1768 MPFR_Object *value, *tempx;
1769 mpfr_exp_t exp = 0;
1770
1771 CHECK_CONTEXT(context);
1772
1773 value = GMPy_MPFR_New(0, context);
1774 tempx = GMPy_MPFR_From_Real(x, 1, context);
1775 result = PyTuple_New(2);
1776 if (!value || !result || !tempx) {
1777 Py_XDECREF((PyObject*)tempx);
1778 Py_XDECREF((PyObject*)value);
1779 Py_XDECREF(result);
1780 return NULL;
1781 }
1782
1783 mpfr_clear_flags();
1784
1785 value->rc = mpfr_frexp(&exp, value->f, tempx->f, GET_MPFR_ROUND(context));
1786 Py_DECREF((PyObject*)tempx);
1787 _GMPy_MPFR_Cleanup(&value, context);
1788
1789 PyTuple_SET_ITEM(result, 0, PyIntOrLong_FromSsize_t((Py_ssize_t)exp));
1790 PyTuple_SET_ITEM(result, 1, (PyObject*)value);
1791 return result;
1792 }
1793
1794 GMPY_MPFR_UNIOP_TEMPLATE(Frexp, frexp)
1795
1796 PyDoc_STRVAR(GMPy_doc_function_next_toward,
1797 "next_toward(x, y) -> mpfr\n\n"
1798 "Return the next 'mpfr' from x in the direction of y. The result has\n"
1799 "the same precision as x.");
1800
1801 PyDoc_STRVAR(GMPy_doc_context_next_toward,
1802 "context.next_toward(x, y) -> mpfr\n\n"
1803 "Return the next 'mpfr' from x in the direction of y. The result has\n"
1804 "the same precision as x.");
1805
1806 static PyObject *
GMPy_Context_NextToward(PyObject * self,PyObject * args)1807 GMPy_Context_NextToward(PyObject *self, PyObject *args)
1808 {
1809 MPFR_Object *result, *tempx, *tempy;
1810 CTXT_Object *context = NULL;
1811 int direction;
1812 mpfr_rnd_t temp_round;
1813
1814 if (self && CTXT_Check(self)) {
1815 context = (CTXT_Object*)self;
1816 }
1817 else {
1818 CHECK_CONTEXT(context);
1819 }
1820
1821 if (PyTuple_GET_SIZE(args) != 2) {
1822 TYPE_ERROR("next_toward() requires 2 arguments");
1823 return NULL;
1824 }
1825
1826 tempx = GMPy_MPFR_From_Real(PyTuple_GET_ITEM(args, 0), 1, context);
1827 tempy = GMPy_MPFR_From_Real(PyTuple_GET_ITEM(args, 1), 1, context);
1828 if (!tempx || !tempy) {
1829 TYPE_ERROR("next_toward() argument type not supported");
1830 Py_XDECREF((PyObject*)tempx);
1831 Py_XDECREF((PyObject*)tempy);
1832 return NULL;
1833 }
1834
1835 if (!(result = GMPy_MPFR_New(mpfr_get_prec(tempx->f), context))) {
1836 Py_DECREF((PyObject*)tempx);
1837 Py_DECREF((PyObject*)tempy);
1838 return NULL;
1839 }
1840
1841 mpfr_clear_flags();
1842
1843 mpfr_set(result->f, tempx->f, GET_MPFR_ROUND(context));
1844 mpfr_nexttoward(result->f, tempy->f);
1845 result->rc = 0;
1846 direction = mpfr_signbit(tempy->f);
1847 Py_DECREF((PyObject*)tempx);
1848 Py_DECREF((PyObject*)tempy);
1849 temp_round = GET_MPFR_ROUND(context);
1850 if (direction)
1851 context->ctx.mpfr_round = MPFR_RNDD;
1852 else
1853 context->ctx.mpfr_round = MPFR_RNDU;
1854 _GMPy_MPFR_Cleanup(&result, context);
1855 context->ctx.mpfr_round = temp_round;
1856 return (PyObject*)result;
1857 }
1858
1859 PyDoc_STRVAR(GMPy_doc_function_next_above,
1860 "next_above(x) -> mpfr\n\n"
1861 "Return the next 'mpfr' from x toward +Infinity.");
1862
1863 PyDoc_STRVAR(GMPy_doc_context_next_above,
1864 "context.next_above(x) -> mpfr\n\n"
1865 "Return the next 'mpfr' from x toward +Infinity.");
1866
1867 static PyObject *
GMPy_Context_NextAbove(PyObject * self,PyObject * other)1868 GMPy_Context_NextAbove(PyObject *self, PyObject *other)
1869 {
1870 MPFR_Object *result, *tempx;
1871 CTXT_Object *context = NULL;
1872 mpfr_rnd_t temp_round;
1873
1874 if (self && CTXT_Check(self)) {
1875 context = (CTXT_Object*)self;
1876 }
1877 else {
1878 CHECK_CONTEXT(context);
1879 }
1880
1881 if (!(tempx = GMPy_MPFR_From_Real(other, 1, context))) {
1882 TYPE_ERROR("next_above() argument type not supported");
1883 return NULL;
1884 }
1885
1886 if (!(result = GMPy_MPFR_New(mpfr_get_prec(tempx->f), context))) {
1887 Py_DECREF((PyObject*)tempx);
1888 return NULL;
1889 }
1890
1891 mpfr_clear_flags();
1892
1893 mpfr_set(result->f, tempx->f, GET_MPFR_ROUND(context));
1894 Py_DECREF((PyObject*)tempx);
1895 mpfr_nextabove(result->f);
1896 result->rc = 0;
1897 temp_round = GET_MPFR_ROUND(context);
1898 context->ctx.mpfr_round = MPFR_RNDU;
1899 _GMPy_MPFR_Cleanup(&result, context);
1900 context->ctx.mpfr_round = temp_round;
1901 return (PyObject*)result;
1902 }
1903
1904 PyDoc_STRVAR(GMPy_doc_function_next_below,
1905 "next_below(x) -> mpfr\n\n"
1906 "Return the next 'mpfr' from x toward -Infinity.");
1907
1908 PyDoc_STRVAR(GMPy_doc_context_next_below,
1909 "context.next_below(x) -> mpfr\n\n"
1910 "Return the next 'mpfr' from x toward -Infinity.");
1911
1912 static PyObject *
GMPy_Context_NextBelow(PyObject * self,PyObject * other)1913 GMPy_Context_NextBelow(PyObject *self, PyObject *other)
1914 {
1915 MPFR_Object *result, *tempx;
1916 CTXT_Object *context = NULL;
1917 mpfr_rnd_t temp_round;
1918
1919 if (self && CTXT_Check(self)) {
1920 context = (CTXT_Object*)self;
1921 }
1922 else {
1923 CHECK_CONTEXT(context);
1924 }
1925
1926 if (!(tempx = GMPy_MPFR_From_Real(other, 1, context))) {
1927 TYPE_ERROR("next_below() argument type not supported");
1928 return NULL;
1929 }
1930
1931 if (!(result = GMPy_MPFR_New(mpfr_get_prec(tempx->f), context))) {
1932 Py_DECREF((PyObject*)tempx);
1933 return NULL;
1934 }
1935
1936 mpfr_clear_flags();
1937
1938 mpfr_set(result->f, tempx->f, GET_MPFR_ROUND(context));
1939 Py_DECREF((PyObject*)tempx);
1940 mpfr_nextbelow(result->f);
1941 result->rc = 0;
1942 temp_round = GET_MPFR_ROUND(context);
1943 context->ctx.mpfr_round = MPFR_RNDD;
1944 _GMPy_MPFR_Cleanup(&result, context);
1945 context->ctx.mpfr_round = temp_round;
1946 return (PyObject*)result;
1947 }
1948
1949 PyDoc_STRVAR(GMPy_doc_function_factorial,
1950 "factorial(n) -> mpfr\n\n"
1951 "Return the floating-point approximation to the factorial of n.\n\n"
1952 "See fac(n) to get the exact integer result.");
1953
1954 PyDoc_STRVAR(GMPy_doc_context_factorial,
1955 "context.factorial(n) -> mpfr\n\n"
1956 "Return the floating-point approximation to the factorial of n.\n\n"
1957 "See fac(n) to get the exact integer result.");
1958
1959 static PyObject *
GMPy_Context_Factorial(PyObject * self,PyObject * other)1960 GMPy_Context_Factorial(PyObject *self, PyObject *other)
1961 {
1962 MPFR_Object *result;
1963 long n;
1964 int error;
1965 CTXT_Object *context = NULL;
1966
1967 if (self && CTXT_Check(self)) {
1968 context = (CTXT_Object*)self;
1969 }
1970 else {
1971 CHECK_CONTEXT(context);
1972 }
1973
1974 if (!(IS_INTEGER(other))) {
1975 TYPE_ERROR("factorial() requires non-negative integer argument");
1976 return NULL;
1977 }
1978
1979 n = GMPy_Integer_AsLongAndError(other, &error);
1980 if ((error == -1) || (n < 0)) {
1981 VALUE_ERROR("factorial() of negative number");
1982 return NULL;
1983 }
1984
1985 if (error == 1) {
1986 VALUE_ERROR("factorial() argument too large");
1987 return NULL;
1988 }
1989
1990 if (error == 2) {
1991 TYPE_ERROR("factorial() requires non-negative integer argument");
1992 return NULL;
1993 }
1994
1995 if (!(result = GMPy_MPFR_New(0, context))) {
1996 return NULL;
1997 }
1998
1999 /* Force result to be 'inf' if n >= 44787928. The matches the behavior of
2000 * MPFR compiled on 64-bit Linux. MPFR compiled on 32-bit Linux and both 32
2001 * and 64-bit versions of Windows will enter an infinite loop. The constant
2002 * 44787929 occurs in the MPFR file gamma.c.
2003 */
2004
2005 mpfr_clear_flags();
2006
2007 if (n >= 44787928) {
2008 mpfr_set_inf(result->f, 1);
2009 mpfr_set_overflow();
2010 }
2011 else {
2012 mpfr_fac_ui(result->f, n, GET_MPFR_ROUND(context));
2013 }
2014
2015 _GMPy_MPFR_Cleanup(&result, context);
2016 return (PyObject*)result;
2017 }
2018
2019 PyDoc_STRVAR(GMPy_doc_function_fsum,
2020 "fsum(iterable) -> mpfr\n\n"
2021 "Return an accurate sum of the values in the iterable.");
2022
2023 PyDoc_STRVAR(GMPy_doc_context_fsum,
2024 "fsum(iterable) -> mpfr\n\n"
2025 "Return an accurate sum of the values in the iterable.");
2026
2027 static PyObject *
GMPy_Context_Fsum(PyObject * self,PyObject * other)2028 GMPy_Context_Fsum(PyObject *self, PyObject *other)
2029 {
2030 MPFR_Object *temp, *result;
2031 mpfr_ptr *tab;
2032 int errcode;
2033 Py_ssize_t i, seq_length = 0;
2034 CTXT_Object *context = NULL;
2035
2036 if (self && CTXT_Check(self)) {
2037 context = (CTXT_Object*)self;
2038 }
2039 else {
2040 CHECK_CONTEXT(context);
2041 }
2042
2043 if (!(result = GMPy_MPFR_New(0, context))) {
2044 return NULL;
2045 }
2046
2047 if (!(other = PySequence_List(other))) {
2048 Py_DECREF((PyObject*)result);
2049 TYPE_ERROR("argument must be an iterable");
2050 return NULL;
2051 }
2052
2053 /* other contains a new list containing all the values from the
2054 * iterable. Now make sure each item in the list is an mpfr.
2055 */
2056
2057 seq_length = PyList_GET_SIZE(other);
2058 if (seq_length > LONG_MAX) {
2059 OVERFLOW_ERROR("temporary array is too large");
2060 Py_DECREF(other);
2061 Py_DECREF((PyObject*)result);
2062 return NULL;
2063 }
2064 for (i=0; i < seq_length; i++) {
2065 if (!(temp = GMPy_MPFR_From_Real(PyList_GET_ITEM(other, i), 1, context))) {
2066 Py_DECREF(other);
2067 Py_DECREF((PyObject*)result);
2068 TYPE_ERROR("all items in iterable must be real numbers");
2069 return NULL;
2070 }
2071
2072 errcode = PyList_SetItem(other, i,(PyObject*)temp);
2073 if (errcode < 0) {
2074 Py_DECREF(other);
2075 Py_DECREF((PyObject*)result);
2076 TYPE_ERROR("all items in iterable must be real numbers");
2077 return NULL;
2078 }
2079 }
2080
2081 /* create an array of pointers to the mpfr_t field of a Pympfr object */
2082
2083 if (!(tab = (mpfr_ptr *)malloc((sizeof(mpfr_srcptr) * seq_length)))) {
2084 Py_DECREF(other);
2085 Py_DECREF((PyObject*)result);
2086 return PyErr_NoMemory();
2087 }
2088 for (i=0; i < seq_length; i++) {
2089 temp = (MPFR_Object*)PyList_GET_ITEM(other, i);
2090 tab[i] = temp->f;
2091 }
2092
2093 mpfr_clear_flags();
2094
2095 /* The cast is safe since we have compared seq_length to LONG_MAX. */
2096 result->rc = mpfr_sum(result->f, tab, (unsigned long)seq_length, GET_MPFR_ROUND(context));
2097 Py_DECREF(other);
2098 free(tab);
2099
2100 _GMPy_MPFR_Cleanup(&result, context);
2101 return (PyObject*)result;
2102 }
2103