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