1 /** @file numeric.cpp
2 *
3 * This file contains the interface to the underlying numerics packages.
4 * Its most important design principle is to completely hide the inner
5 * working of that other package from the user of GiNaC. It must either
6 * provide implementation of arithmetic operators and numerical evaluation
7 * of special functions or implement the interface to the numericsd package. */
8
9 /*
10 * This is a modified version of code included with Ginac.
11 * The modifications and modified version is:
12 *
13 * GiNaC-SAGE Copyright (C) 2008 William Stein
14 * Copyright (C) 2009 Burcin Erocal <burcin@erocal.org>
15 * (C) 2015-2017 Ralf Stephan <ralf@ark.in-berlin.de>
16 *
17 * This program is free software; you can redistribute it and/or modify
18 * it under the terms of the GNU General Public License as published by
19 * the Free Software Foundation; either version 2 of the License, or
20 * (at your option) any later version.
21 *
22 * This program is distributed in the hope that it will be useful,
23 * but WITHOUT ANY WARRANTY; without even the implied warranty of
24 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25 * GNU General Public License for more details.
26 *
27 * You should have received a copy of the GNU General Public License
28 * along with this program; if not, write to the Free Software
29 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
30 */
31
32
33 /* The original copyright:
34 *
35 * GiNaC Copyright (C) 1999-2008 Johannes Gutenberg University Mainz, Germany
36 *
37 * This program is free software; you can redistribute it and/or modify
38 * it under the terms of the GNU General Public License as published by
39 * the Free Software Foundation; either version 2 of the License, or
40 * (at your option) any later version.
41 *
42 * This program is distributed in the hope that it will be useful,
43 * but WITHOUT ANY WARRANTY; without even the implied warranty of
44 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
45 * GNU General Public License for more details.
46 *
47 * You should have received a copy of the GNU General Public License
48 * along with this program; if not, write to the Free Software
49 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
50 */
51
52 #define register
53 #define PY_SSIZE_T_CLEAN
54 #include <Python.h>
55 #include <longintrepr.h>
56 #include "flint/fmpz.h"
57 #include "flint/fmpz_factor.h"
58
59 #include <cmath>
60 #include <utility>
61
62 #include "numeric.h"
63 #include "operators.h"
64 #include "ex.h"
65 #include "mul.h"
66 #include "power.h"
67 #include "function.h"
68 #include "archive.h"
69 #include "tostring.h"
70 #include "utils.h"
71
72 #pragma clang diagnostic push
73 #pragma clang diagnostic ignored "-Wdeprecated-register"
74 #include "factory/factory.h"
75 #pragma clang diagnostic pop
76
77 #ifdef PYNAC_HAVE_LIBGIAC
78 #undef _POSIX_C_SOURCE
79 #undef _XOPEN_SOURCE
80
81 #include <giac/global.h>
82 #include <giac/gausspol.h>
83 #include <giac/fraction.h>
84 #endif
85
86 //#define Logging_refctr
87 #if defined(Logging_refctr)
88 #undef Py_INCREF
89 #undef Py_DECREF
90
91 #define Py_INCREF(op) ( \
92 _Py_INC_REFTOTAL _Py_REF_DEBUG_COMMA \
93 ((PyObject*)(op))->ob_refcnt++) ; \
94 std::cerr << "+ " << long(op) << ", " << Py_REFCNT(op) << ", " << Py_TYPE(op)->tp_name << std::endl; std::cerr.flush();
95
96 #define Py_DECREF(op) \
97 do { \
98 std::cerr << "- " << long(op) << ", " << Py_REFCNT(op) << ", " << Py_TYPE(op)->tp_name << std::endl; std::cerr.flush(); \
99 if (_Py_DEC_REFTOTAL _Py_REF_DEBUG_COMMA \
100 --((PyObject*)(op))->ob_refcnt != 0) \
101 _Py_CHECK_REFCNT(op) \
102 else \
103 _Py_Dealloc((PyObject *)(op)); \
104 } while (0); std::cerr << "Done" << std::endl; std::cerr.flush();
105 #endif
106
to_numeric(GiNaC::ex & e)107 const GiNaC::numeric& to_numeric(GiNaC::ex& e)
108 {
109 return GiNaC::ex_to<GiNaC::numeric>(e);
110 }
111
112 // Call the Python function f on *this as input and return the result
113 // as a PyObject*.
114 #define PY_RETURN(f) PyObject *a = to_pyobject(); \
115 PyObject *ans = f(a); \
116 Py_DECREF(a); \
117 if (!ans) py_error("error calling function"); \
118 return ans;
119
120 // Call the Python function f on *this and return the result
121 // as a PyObject*.
122 #define PY_RETURN2(f, b) PyObject *aa = to_pyobject(); \
123 PyObject* bb = b.to_pyobject(); \
124 PyObject *ans = f(aa, bb); \
125 if (!ans) py_error("error calling function"); \
126 Py_DECREF(aa); Py_DECREF(bb); return ans;
127
128 // Call the Python function f on *this and b, and get back
129 // a 2-tuple (z,w). Set c = w, where c should be a
130 // reference in the caller, and return z. This is used
131 // to return two inputs from a Python function call. See
132 // its usage in code below.
133 #define PY_RETURN3(f, b, c) \
134 PyObject *aa = to_pyobject(); \
135 PyObject* bb = b.to_pyobject(); \
136 PyObject *ans = f(aa, bb); \
137 if (!ans) py_error("error calling function"); \
138 if (!PyTuple_CheckExact(ans) || PyTuple_GET_SIZE(ans) != 2) \
139 py_error("error calling function -- return not a 2-tuple."); \
140 PyObject* z = PyTuple_GET_ITEM(ans, 0); Py_INCREF(z); \
141 PyObject* w = PyTuple_GET_ITEM(ans, 1); Py_INCREF(w); \
142 c = w; \
143 Py_DECREF(aa); Py_DECREF(bb); Py_DECREF(ans); \
144 return z;
145
146 //#define DEBUG
147 //#define VERBOSE
148
149 #ifdef DEBUG
150 #define todo(s) std::cerr << "TODO: " << s << std::endl;
151 #define stub(s) { std::cerr << "Hit STUB: " << s << std::endl; throw std::runtime_error("stub"); }
152 #define fake(s) std::cerr << "fake: " << s << std::endl;
153 #define ASSERT(s, msg) if (!s) { std::cerr << "Failed assertion: " << msg << std::endl; }
154 #else
155 #define todo(s)
156 #define stub(s) { std::cerr << "** Hit STUB**: " << s << std::endl; throw std::runtime_error("stub"); }
157 #define fake(s)
158 #endif
159
160 #ifdef VERBOSE
161 #define verbose(s) std::cerr << s << std::endl;
162 #define verbose2(s,t) std::cerr << s << " " << t << std::endl;
163 #define verbose3(s,t,u) std::cerr << s << " " << t << ", " << u << std::endl;
164 #else
165 #define verbose(s)
166 #define verbose2(s,t)
167 #define verbose3(s,t,u)
168 #endif
169
170
171 //////////////////////////////////////////////////////////////
172 // Python Interface
173 //////////////////////////////////////////////////////////////
174
py_error(const char * errmsg)175 inline void py_error(const char* errmsg) {
176 throw std::runtime_error((PyErr_Occurred() != nullptr) ? errmsg:
177 "pyerror() called but no error occurred!");
178 }
179
180 #if PY_MAJOR_VERSION < 3
181 #define PyNumber_TrueDivide PyNumber_Divide
182
183 #else
184 #define PyString_FromString PyUnicode_FromString
185 #endif
186
187 // The following variable gets changed to true once
188 // this library has been imported by the Python
189 // interpreter. This is because the interpreter calls
190 // ginac_pyinit_I below, which sets this to true.
191 // Once this is done we can call all the py_* functions
192 // defined in Cython modules, which is often much faster
193 // than what we can do without those.
194 static bool initialized = false;
195
196 static PyObject* pyfunc_Float = nullptr;
197
ginac_pyinit_Float(PyObject * f)198 void ginac_pyinit_Float(PyObject* f) {
199 Py_INCREF(f);
200 pyfunc_Float = f;
201 }
202
ginac_pyinit_I(PyObject * z)203 void ginac_pyinit_I(PyObject* z) {
204 initialized = true;
205 Py_INCREF(z);
206 GiNaC::I = z; // I is a global constant defined below.
207 }
208
RR_get()209 PyObject* RR_get()
210 {
211 static PyObject* ptr = nullptr;
212 if (ptr == nullptr) {
213 PyObject* m = PyImport_ImportModule("sage.rings.all");
214 if (m == nullptr)
215 py_error("Error importing sage.rings.all");
216 ptr = PyObject_GetAttrString(m, "RR");
217 if (ptr == nullptr)
218 py_error("Error getting RR attribute");
219 Py_INCREF(ptr);
220 }
221 return ptr;
222 }
223
CC_get()224 PyObject* CC_get()
225 {
226 static PyObject* ptr = nullptr;
227 if (ptr)
228 return ptr;
229 PyObject* m = PyImport_ImportModule("sage.rings.all");
230 if (m == nullptr)
231 py_error("Error importing sage.rings.all");
232 ptr = PyObject_GetAttrString(m, "ComplexField");
233 if (ptr == nullptr)
234 py_error("Error getting ComplexField attribute");
235 ptr = PyObject_CallObject(ptr, NULL);
236 if (ptr == nullptr)
237 py_error("Error getting CC attribute");
238 Py_INCREF(ptr);
239 return ptr;
240 }
241
242 static PyObject* pyfunc_Integer = nullptr;
243
ginac_pyinit_Integer(PyObject * f)244 void ginac_pyinit_Integer(PyObject* f) {
245 Py_INCREF(f);
246 pyfunc_Integer = f;
247 }
248
Integer_pyclass()249 PyObject* Integer_pyclass()
250 {
251 if (not initialized)
252 throw std::runtime_error("can't happen");
253 static PyObject *ptr = nullptr;
254 if (ptr)
255 return ptr;
256 PyObject *m = PyImport_ImportModule("sage.rings.integer");
257 if (m == nullptr)
258 py_error("Error importing sage.rings.integer");
259 ptr = PyObject_GetAttrString(m, "Integer");
260 if (ptr == nullptr)
261 py_error("Error getting Integer attribute");
262 return ptr;
263 }
264
Integer(const long int & x)265 PyObject* Integer(const long int& x) {
266 if (initialized)
267 return GiNaC::py_funcs.py_integer_from_long(x);
268
269 PyObject *Integer = Integer_pyclass();
270 PyObject* ans = PyObject_CallFunction(Integer, const_cast<char*> ("l"), x);
271 return ans;
272 }
273
precision(const GiNaC::numeric & num,PyObject * & a_parent)274 int precision(const GiNaC::numeric& num, PyObject*& a_parent) {
275 int prec = 0;
276 PyObject* the_parent = a_parent;
277 if (a_parent == nullptr) {
278 PyObject* m = PyImport_ImportModule("sage.structure.element");
279 if (m == nullptr)
280 py_error("Error importing element");
281 PyObject* f = PyObject_GetAttrString(m, "parent");
282 if (f == nullptr)
283 py_error("Error getting parent attribute");
284 PyObject* obj = num.to_pyobject();
285 the_parent = PyObject_CallFunctionObjArgs(f, obj, NULL);
286 Py_DECREF(obj);
287 Py_DECREF(f);
288 Py_DECREF(m);
289 }
290 else if (PyDict_Check(a_parent)) {
291 PyObject* pkey = PyString_FromString(const_cast<char*>("parent"));
292 the_parent = PyDict_GetItem(a_parent, pkey);
293 Py_DECREF(pkey);
294 }
295 PyObject *mprec = nullptr;
296 if (the_parent != nullptr)
297 mprec = PyObject_CallMethod(the_parent, const_cast<char*>("precision"), NULL);
298 if (mprec == nullptr) {
299 prec = 53;
300 PyErr_Clear();
301 }
302 else {
303 prec = PyLong_AsLong(mprec);
304 Py_DECREF(mprec);
305 }
306 if (a_parent == nullptr) {
307 a_parent = PyDict_New();
308 PyDict_SetItemString(a_parent, "parent", the_parent);
309 }
310 return prec;
311 }
312
CBF(int res)313 PyObject* CBF(int res) {
314 PyObject* m = PyImport_ImportModule("sage.rings.all");
315 if (m == nullptr)
316 py_error("Error importing arb");
317 PyObject* f = PyObject_GetAttrString(m, "ComplexBallField");
318 if (f == nullptr)
319 py_error("Error getting ComplexBallField attribute");
320 PyObject* list = PyTuple_New(1);
321 if (list == nullptr)
322 throw(std::runtime_error("GiNaC::CBF(): PyTuple_New returned NULL"));
323 PyObject *iobj = Integer(res);
324 int ret = PyTuple_SetItem(list, 0, iobj);
325 if (ret != 0)
326 throw(std::runtime_error("GiNaC::CBF(): PyTuple_SetItem unsuccessful"));
327 PyObject *obj = PyObject_Call(f, list, NULL);
328 if (obj == nullptr)
329 throw(std::runtime_error("GiNaC::CBF(): PyObject_Call unsuccessful"));
330 Py_DECREF(m);
331 Py_DECREF(f);
332 Py_DECREF(list);
333 return obj;
334 }
335
336 // Convert a to field elt, return a.meth()
CallBallMethod0Arg(PyObject * field,const char * meth,const GiNaC::numeric & a)337 PyObject* CallBallMethod0Arg(PyObject* field, const char* meth, const GiNaC::numeric& a) {
338 PyObject* list1 = PyTuple_New(1);
339 if (list1 == nullptr)
340 throw(std::runtime_error("GiNaC::CallBallMethod1Arg(): PyTuple_New returned NULL"));
341 PyObject *aobj = a.to_pyobject();
342 int ret = PyTuple_SetItem(list1, 0, aobj);
343 if (ret != 0)
344 throw(std::runtime_error("GiNaC::CallBallMethod1Arg(): PyTuple_SetItem unsuccessful"));
345 PyObject *aball = PyObject_Call(field, list1, NULL);
346 if (aball == nullptr)
347 throw(std::runtime_error("GiNaC::CallBallMethod1Arg(): PyObject_Call unsuccessful"));
348 PyObject* name = PyString_FromString(meth);
349 PyObject* ret_ball = PyObject_CallMethodObjArgs(aball, name, NULL);
350 if (ret_ball == nullptr)
351 throw(std::runtime_error("GiNaC::CallBallMethod1Arg(): PyObject_CallMethodObjArgs unsuccessful"));
352 Py_DECREF(list1);
353 Py_DECREF(aball);
354 Py_DECREF(name);
355
356 return ret_ball;
357 }
358
359 // Convert a to field elt, return a.meth(b)
CallBallMethod1Arg(PyObject * field,const char * meth,const GiNaC::numeric & a,const GiNaC::numeric & b)360 PyObject* CallBallMethod1Arg(PyObject* field, const char* meth, const GiNaC::numeric& a, const GiNaC::numeric& b) {
361 PyObject* list1 = PyTuple_New(1);
362 if (list1 == nullptr)
363 throw(std::runtime_error("GiNaC::CallBallMethod1Arg(): PyTuple_New returned NULL"));
364 PyObject *aobj = a.to_pyobject();
365 int ret = PyTuple_SetItem(list1, 0, aobj);
366 if (ret != 0)
367 throw(std::runtime_error("GiNaC::CallBallMethod1Arg(): PyTuple_SetItem unsuccessful"));
368 PyObject *aball = PyObject_Call(field, list1, NULL);
369 if (aball == nullptr)
370 throw(std::runtime_error("GiNaC::CallBallMethod1Arg(): PyObject_Call unsuccessful"));
371 PyObject* list2 = PyTuple_New(1);
372 if (list2 == nullptr)
373 throw(std::runtime_error("GiNaC::CallBallMethod1Arg(): PyTuple_New returned NULL"));
374 PyObject *bobj = b.to_pyobject();
375 ret = PyTuple_SetItem(list2, 0, bobj);
376 if (ret != 0)
377 throw(std::runtime_error("GiNaC::CallBallMethod1Arg(): PyTuple_SetItem unsuccessful"));
378 PyObject *bball = PyObject_Call(field, list2, NULL);
379 if (bball == nullptr)
380 throw(std::runtime_error("GiNaC::CallBallMethod1Arg(): PyObject_Call unsuccessful"));
381
382 PyObject* name = PyString_FromString(meth);
383 PyObject* ret_ball = PyObject_CallMethodObjArgs(aball, name, bball, NULL);
384 if (ret_ball == nullptr)
385 throw(std::runtime_error("GiNaC::CallBallMethod1Arg(): PyObject_CallMethodObjArgs unsuccessful"));
386 Py_DECREF(list1);
387 Py_DECREF(list2);
388 Py_DECREF(aball);
389 Py_DECREF(bball);
390 Py_DECREF(name);
391
392 return ret_ball;
393 }
394
CoerceBall(PyObject * ball,int prec)395 PyObject* CoerceBall(PyObject* ball, int prec) {
396 PyObject* m = PyImport_ImportModule("sage.rings.all");
397 if (m == nullptr)
398 py_error("Error importing sage.rings.all");
399 PyObject* f = PyObject_GetAttrString(m, "ComplexField");
400 if (f == nullptr)
401 py_error("Error getting ComplexField attribute");
402 PyObject* list1 = PyTuple_New(1);
403 if (list1 == nullptr)
404 throw(std::runtime_error("GiNaC: PyTuple_New returned NULL"));
405 PyObject *iobj = Integer(prec);
406 int ret = PyTuple_SetItem(list1, 0, iobj);
407 if (ret != 0)
408 throw(std::runtime_error("GiNaC: PyTuple_SetItem unsuccessful"));
409 PyObject *field = PyObject_CallObject(f, list1);
410 if (field == nullptr)
411 throw(std::runtime_error("GiNaC: PyObject_Call unsuccessful"));
412 PyObject* list2 = PyTuple_New(1);
413 if (list2 == nullptr)
414 throw(std::runtime_error("GiNaC: PyTuple_New returned NULL"));
415 ret = PyTuple_SetItem(list2, 0, ball);
416 if (ret != 0)
417 throw(std::runtime_error("GiNaC: PyTuple_SetItem unsuccessful"));
418 PyObject *celt = PyObject_CallObject(field, list2);
419 if (celt == nullptr)
420 throw(std::runtime_error("GiNaC: PyObject_Call unsuccessful"));
421
422 Py_INCREF(ball);
423 Py_DECREF(list1);
424 Py_DECREF(list2);
425 Py_DECREF(field);
426 Py_DECREF(f);
427 Py_DECREF(m);
428
429 PyObject* bret = PyObject_CallMethod(celt, const_cast<char*>("is_real"), NULL);
430 if (PyObject_IsTrue(bret)) {
431 PyObject* rret = PyObject_CallMethod(celt, const_cast<char*>("real"), NULL);
432 Py_DECREF(bret);
433 Py_DECREF(celt);
434 return rret;
435 }
436 Py_DECREF(bret);
437 return celt;
438 }
439
440 namespace GiNaC {
441
442 numeric I;
443
common_parent(const numeric & x,const numeric & y)444 PyObject* common_parent(const numeric& x, const numeric& y)
445 {
446 PyObject* m = PyImport_ImportModule("sage.structure.element");
447 if (m == nullptr)
448 py_error("Error importing coerce");
449 PyObject* f = PyObject_GetAttrString(m, "coercion_model");
450 if (f == nullptr)
451 py_error("Error getting coercion_model attribute");
452
453 PyObject* mname = PyString_FromString("common_parent");
454 PyObject* xobj = x.to_pyobject();
455 PyObject* yobj = y.to_pyobject();
456 PyObject* ret = PyObject_CallMethodObjArgs(f, mname,
457 xobj, yobj, NULL);
458 if (ret == nullptr)
459 throw(std::runtime_error("GiNaC::common_parent: PyObject_CallMethodObjArgs unsuccessful"));
460 Py_DECREF(xobj);
461 Py_DECREF(yobj);
462 Py_DECREF(m);
463 Py_DECREF(f);
464 Py_DECREF(mname);
465 return ret;
466 }
467
arbfunc_0arg(const char * name,PyObject * parent) const468 const numeric numeric::arbfunc_0arg(const char* name, PyObject* parent) const
469 {
470 int prec = precision(*this, parent);
471 PyObject* field = CBF(prec+15);
472 PyObject* ret = CallBallMethod0Arg(field, name, *this);
473 Py_DECREF(field);
474
475 numeric rnum(ret);
476 return ex_to<numeric>(rnum.evalf(0, parent));
477 }
478
py_tuple_from_numvector(const std::vector<numeric> & vec)479 static PyObject* py_tuple_from_numvector(const std::vector<numeric>& vec)
480 {
481 PyObject* list = PyTuple_New(vec.size());
482 if (list == nullptr)
483 throw(std::runtime_error("py_list_from_numvector(): PyList_New returned NULL"));
484 int pos = 0;
485 for (const numeric& num : vec) {
486 PyObject *numobj = num.to_pyobject();
487 int ret = PyTuple_SetItem(list, pos++, numobj);
488 if (ret != 0)
489 throw(std::runtime_error("py_list_from_numvector(): PyList_Append unsuccessful"));
490 }
491 return list;
492 }
493
494 ///////////////////////////////////////////////////////////////////////////////
495 // class numeric
496 ///////////////////////////////////////////////////////////////////////////////
497
498 #if PY_MAJOR_VERSION < 3
499 PyObject* ZERO = PyInt_FromLong(0); // todo: never freed
500 PyObject* ONE = PyInt_FromLong(1); // todo: never freed
501 PyObject* TWO = PyInt_FromLong(2); // todo: never freed
502 #else
503 PyObject* ZERO = PyLong_FromLong(0); // todo: never freed
504 PyObject* ONE = PyLong_FromLong(1); // todo: never freed
505 PyObject* TWO = PyLong_FromLong(2); // todo: never freed
506 #endif
507
operator <<(std::ostream & os,const numeric & s)508 std::ostream& operator<<(std::ostream& os, const numeric& s) {
509 switch (s.t) {
510 case LONG:
511 return os << s.v._long;
512 case MPZ: {
513 std::vector<char> cp(2 + mpz_sizeinbase(s.v._bigint, 10));
514 mpz_get_str(&cp[0], 10, s.v._bigint);
515 return os << &cp[0];
516 }
517 case MPQ: {
518 size_t size = mpz_sizeinbase(mpq_numref(s.v._bigrat), 10)
519 + mpz_sizeinbase(mpq_denref(s.v._bigrat), 10) + 5;
520 std::vector<char> cp(size);
521 mpq_get_str(&cp[0], 10, s.v._bigrat);
522 return os << &cp[0];
523 }
524 case PYOBJECT:
525 return os << *py_funcs.py_repr(s.v._pyobject, 0);
526 default:
527 stub("operator <<: type not yet handled");
528 }
529 }
530
operator =(const numeric & x)531 const numeric& numeric::operator=(const numeric& x) {
532 switch (t) {
533 case LONG:
534 break;
535 case MPZ:
536 mpz_clear(v._bigint);
537 break;
538 case MPQ:
539 mpq_clear(v._bigrat);
540 break;
541 case PYOBJECT:
542 Py_DECREF(v._pyobject);
543 break;
544 }
545 t = x.t;
546 hash = x.hash;
547 switch (x.t) {
548 case LONG:
549 v._long = x.v._long;
550 break;
551 case MPZ:
552 mpz_init(v._bigint);
553 mpz_set(v._bigint, x.v._bigint);
554 break;
555 case MPQ:
556 mpq_init(v._bigrat);
557 mpq_set(v._bigrat, x.v._bigrat);
558 break;
559 case PYOBJECT:
560 v = x.v;
561 Py_INCREF(v._pyobject);
562 break;
563 }
564 return *this;
565 }
566
swap(numeric & other)567 void numeric::swap(numeric& other)
568 {
569 std::swap(t, other.t);
570 std::swap(v, other.v);
571 std::swap(hash, other.hash);
572 std::swap(is_hashable, other.is_hashable);
573 }
574
compare_same_type(const numeric & right) const575 int numeric::compare_same_type(const numeric& right) const {
576 verbose("compare_same_type");
577 if (t != right.t) {
578 int ret;
579 if (t == MPZ and right.t == MPQ) {
580 ret = -mpq_cmp_z(right.v._bigrat, v._bigint);
581 }
582 else if (t == MPQ and right.t == MPZ) {
583 ret = mpq_cmp_z(v._bigrat, right.v._bigint);
584 }
585 else if (t == LONG and right.t == MPZ) {
586 ret = -mpz_cmp_si(right.v._bigint, v._long);
587 }
588 else if (t == LONG and right.t == MPQ) {
589 ret = -mpq_cmp_si(right.v._bigrat, v._long, 1);
590 }
591 else if (t == MPZ and right.t == LONG) {
592 ret = mpz_cmp_si(v._bigint, right.v._long);
593 }
594 else if (t == MPQ and right.t == LONG) {
595 ret = mpq_cmp_si(v._bigrat, right.v._long, 1);
596 }
597 else {
598 numeric a, b;
599 coerce(a, b, *this, right);
600 return a.compare_same_type(b);
601 }
602 if (ret > 0)
603 ret = 1;
604 else if (ret < 0)
605 ret = -1;
606 return ret;
607 }
608 int ret;
609 switch (t) {
610 case LONG:
611 if (v._long > right.v._long)
612 return 1;
613 else if (v._long < right.v._long)
614 return -1;
615 else
616 return 0;
617 case MPZ:
618 ret = mpz_cmp(v._bigint, right.v._bigint);
619 if (ret > 0)
620 ret = 1;
621 else if (ret < 0)
622 ret = -1;
623 return ret;
624 case MPQ:
625 ret = mpq_cmp(v._bigrat, right.v._bigrat);
626 if (ret > 0)
627 ret = 1;
628 else if (ret < 0)
629 ret = -1;
630 return ret;
631 case PYOBJECT: {
632 int result = PyObject_RichCompareBool(v._pyobject,
633 right.v._pyobject, Py_LT);
634 if (result == 1)
635 return -1;
636 if (result == -1)
637 py_error("richcmp failed");
638 else { // result == 0
639 result = PyObject_RichCompareBool(v._pyobject,
640 right.v._pyobject, Py_GT);
641 if (result == -1)
642 py_error("richcmp failed");
643 }
644 return result;
645 }
646 default:
647 stub("invalid type: compare_same_type type not handled");
648 }
649 }
650
651 #if PY_MAJOR_VERSION < 3
652 #define hash_bits (8 * sizeof(void*))
653 #else
654 #define hash_bits _PyHASH_BITS
655 #endif
656 #define limb_bits (8 * sizeof(mp_limb_t))
657
658 /* By convention hashes of PyObjects must be identical
659 with their Python hashes, this applies to our MPZ
660 and MPQ objects too. This implementation is copied
661 from sage.libs.gmp.pylong.mpz_pythonhash. */
_mpz_pythonhash_raw(mpz_t the_int)662 static long _mpz_pythonhash_raw(mpz_t the_int)
663 {
664 if (mpz_sgn(the_int) == 0)
665 return 0;
666
667 mp_limb_t modulus = ((((mp_limb_t)(1) << (hash_bits - 1)) - 1) * 2) + 1;
668 mp_limb_t h=0, x, y;
669 size_t n = mpz_size(the_int);
670 unsigned int r;
671 for (unsigned i=0; i<n; ++i) {
672 x = mpz_getlimbn(the_int, i);
673 if (limb_bits == hash_bits)
674 y = x;
675 else {
676 r = (limb_bits * i) % hash_bits;
677 y = (x << r) & modulus;
678 y += (x >> (hash_bits - r)) & modulus;
679 if (r > 2 * hash_bits - limb_bits)
680 y += (x >> (2 * hash_bits - r));
681 if (y > modulus)
682 y -= modulus;
683 }
684 if (h < modulus - y)
685 h = h + y;
686 else
687 h = h - (modulus - y);
688 }
689 if (limb_bits == hash_bits && h == 0)
690 h = -1;
691
692 if (mpz_sgn(the_int) < 0)
693 return -h;
694 return h;
695 }
696
_mpz_pythonhash(mpz_t the_int)697 static long _mpz_pythonhash(mpz_t the_int)
698 {
699 long h = _mpz_pythonhash_raw(the_int);
700 if (h == -1)
701 return -2;
702 return h;
703 }
704
_mpq_pythonhash(mpq_t the_rat)705 static long _mpq_pythonhash(mpq_t the_rat)
706 {
707 mpq_t rat;
708 mpq_init(rat);
709 mpq_set(rat, the_rat);
710 long n = _mpz_pythonhash_raw(mpq_numref(rat));
711 long d = _mpz_pythonhash_raw(mpq_denref(rat));
712 if (d != 1L)
713 n = n + (d-1) * 7461864723258187525;
714 mpq_clear(rat);
715 if (n == -1)
716 return -2;
717 return n;
718 }
719
720
721 // Initialize an mpz_t from a Python long integer
_mpz_set_pylong(mpz_t z,PyLongObject * l)722 static void _mpz_set_pylong(mpz_t z, PyLongObject* l)
723 {
724 Py_ssize_t pylong_size = Py_SIZE(l);
725 int sign = 1;
726
727 if (pylong_size < 0) {
728 pylong_size = -pylong_size;
729 sign = -1;
730 }
731
732 mpz_import(z, pylong_size, -1, sizeof(digit), 0,
733 8*sizeof(digit) - PyLong_SHIFT, l->ob_digit);
734
735 if (sign < 0)
736 mpz_neg(z, z);
737 }
738
739
740
741 ///////////////////////////////////////////////////////////////////////////////
742 // class numeric
743 ///////////////////////////////////////////////////////////////////////////////
744
745 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(numeric, basic,
746 print_func<print_context>(&numeric::do_print).
747 print_func<print_latex>(&numeric::do_print_latex).
748 print_func<print_tree>(&numeric::do_print_tree).
749 print_func<print_python_repr>(&numeric::do_print_python_repr))
750
751
numeric()752 numeric::numeric() : basic(&numeric::tinfo_static), t(LONG), v(0) {
753 setflag(status_flags::evaluated | status_flags::expanded);
754 }
755
756 //////////
757 // other constructors
758 //////////
759
760 // public
761
numeric(const numeric & other)762 numeric::numeric(const numeric& other) : basic(&numeric::tinfo_static) {
763 t = other.t;
764 hash = other.hash;
765 switch (t) {
766 case LONG:
767 v._long = other.v._long;
768 return;
769 case PYOBJECT:
770 v = other.v;
771 Py_INCREF(v._pyobject);
772 return;
773 case MPZ:
774 mpz_init(v._bigint);
775 mpz_set(v._bigint, other.v._bigint);
776 return;
777 case MPQ:
778 mpq_init(v._bigrat);
779 mpq_set(v._bigrat, other.v._bigrat);
780 return;
781 }
782 }
783
is_Sage_Integer(PyObject * o)784 bool is_Sage_Integer(PyObject* o)
785 {
786 PyObject *Integer = Integer_pyclass();
787 int res = PyObject_IsInstance(o, Integer);
788 if (res < 0)
789 py_error("Error testing Integer attribute");
790 return res != 0;
791 }
792
set_from(Type & t,Value & v,long & hash,mpz_t bigint)793 void set_from(Type& t, Value& v, long& hash, mpz_t bigint)
794 {
795 if (mpz_fits_sint_p(bigint)) {
796 t = LONG;
797 v._long = mpz_get_si(bigint);
798 hash = (v._long==-1) ? -2 : v._long;
799 }
800 else {
801 t = MPZ;
802 mpz_init_set(v._bigint, bigint);
803 hash = _mpz_pythonhash(v._bigint);
804 }
805 }
806
set_from(Type & t,Value & v,long & hash,mpq_t bigrat)807 void set_from(Type& t, Value& v, long& hash, mpq_t bigrat)
808 {
809 t = MPQ;
810 mpq_init(v._bigrat);
811 mpq_set(v._bigrat, bigrat);
812 hash = _mpq_pythonhash(v._bigrat);
813 }
814
numeric(PyObject * o,bool force_py)815 numeric::numeric(PyObject* o, bool force_py) : basic(&numeric::tinfo_static) {
816 if (o == nullptr) py_error("Error");
817 if (not force_py) {
818 #if PY_MAJOR_VERSION < 3
819 if (PyInt_Check(o)) {
820 t = LONG;
821 v._long = PyInt_AsLong(o);
822 hash = (v._long==-1) ? -2 : v._long;
823 setflag(status_flags::evaluated | status_flags::expanded);
824 Py_DECREF(o);
825 return;
826 }
827 #endif
828 if (PyLong_Check(o)) {
829 t = MPZ;
830 mpz_init(v._bigint);
831 _mpz_set_pylong(v._bigint, reinterpret_cast<PyLongObject*>( o));
832 hash = _mpz_pythonhash(v._bigint);
833 setflag(status_flags::evaluated | status_flags::expanded);
834 Py_DECREF(o);
835 return;
836 }
837 if (initialized) {
838 if (is_Sage_Integer(o)) {
839 mpz_ptr mpz = py_funcs.py_mpz_from_integer(o);
840 set_from(t, v, hash, mpz);
841 Py_DECREF(o);
842 setflag(status_flags::evaluated
843 | status_flags::expanded);
844 return;
845 }
846 if (py_funcs.py_is_Rational(o) != 0) {
847 mpq_ptr mpq = py_funcs.py_mpq_from_rational(o);
848 set_from(t, v, hash, mpq);
849 Py_DECREF(o);
850 setflag(status_flags::evaluated
851 | status_flags::expanded);
852 return;
853 }
854 }
855 }
856
857 t = PYOBJECT;
858 hash = PyObject_Hash(o);
859 if (hash == -1 && (PyErr_Occurred() != nullptr)) {
860 // error is thrown on first hash request
861 PyErr_Clear();
862 is_hashable = false;
863 }
864 v._pyobject = o; // STEAL a reference
865 setflag(status_flags::evaluated | status_flags::expanded);
866
867 }
868
numeric(mpz_t bigint)869 numeric::numeric(mpz_t bigint) : basic(&numeric::tinfo_static)
870 {
871 set_from(t, v, hash, bigint);
872 mpz_clear(bigint);
873 setflag(status_flags::evaluated | status_flags::expanded);
874 }
875
numeric(mpq_t bigrat)876 numeric::numeric(mpq_t bigrat) : basic(&numeric::tinfo_static)
877 {
878 set_from(t, v, hash, bigrat);
879 mpq_clear(bigrat);
880 setflag(status_flags::evaluated | status_flags::expanded);
881 }
882
883 /** Constructor for rational numerics a/b.
884 *
885 * @exception overflow_error (division by zero) */
numeric(long num,long den)886 numeric::numeric(long num, long den) : basic(&numeric::tinfo_static) {
887 if (den == 0)
888 throw std::overflow_error("numeric::div(): division by zero");
889 if ((num%den) == 0) {
890 t = LONG;
891 v._long = num/den;
892 hash = (v._long==-1) ? -2 : v._long;
893 }
894 else
895 {
896 t = MPQ;
897 mpq_init(v._bigrat);
898 mpq_set_si(v._bigrat, num, den);
899 mpq_canonicalize(v._bigrat);
900 hash = _mpq_pythonhash(v._bigrat);
901 }
902 setflag(status_flags::evaluated | status_flags::expanded);
903 }
904
numeric(double d)905 numeric::numeric(double d) : basic(&numeric::tinfo_static) {
906 t = PYOBJECT;
907 if ((v._pyobject = PyFloat_FromDouble(d)) == nullptr)
908 py_error("Error creating double");
909 setflag(status_flags::evaluated | status_flags::expanded);
910 }
911
~numeric()912 numeric::~numeric() {
913 switch (t) {
914 case LONG:
915 return;
916 case PYOBJECT:
917 Py_DECREF(v._pyobject);
918 return;
919 case MPZ:
920 mpz_clear(v._bigint);
921 return;
922 case MPQ:
923 mpq_clear(v._bigrat);
924 return;
925 }
926 }
927
928 //////////
929 // archiving
930 //////////
931
numeric(const archive_node & n,lst & sym_lst)932 numeric::numeric(const archive_node &n, lst &sym_lst) :
933 inherited(n, sym_lst) {
934 // get type information
935 unsigned int t_tmp;
936 if (!n.find_unsigned(std::string("T"), t_tmp))
937 throw std::runtime_error("archive error: cannot read type info");
938 t = Type(t_tmp);
939 std::string str;
940 // read object to a string
941 if (!n.find_string("S", str))
942 throw (std::runtime_error("archive error: cannot read object data"));
943 PyObject *arg;
944 switch (t) {
945 case LONG:
946 v._long = std::stol(str);
947 hash = (v._long == -1) ? -2 : v._long;
948 return;
949 case MPZ:
950 mpz_init(v._bigint);
951 mpz_set_str(v._bigint, str.c_str(), 10);
952 hash = _mpz_pythonhash(v._bigint);
953 return;
954 case MPQ:
955 mpq_init(v._bigrat);
956 mpq_set_str(v._bigrat, str.c_str(), 10);
957 hash = _mpq_pythonhash(v._bigrat);
958 return;
959 case PYOBJECT:
960 // read pickled python object to a string
961 if (!n.find_string("S", str))
962 throw(std::runtime_error("archive error: cannot read pyobject data"));
963 arg = Py_BuildValue("s#", str.c_str(), str.size());
964 // unpickle
965 v._pyobject = py_funcs.py_loads(arg);
966 Py_DECREF(arg);
967 if (PyErr_Occurred() != nullptr) {
968 throw(std::runtime_error("archive error: caught exception in py_loads"));
969 }
970 hash = PyObject_Hash(v._pyobject);
971 if (hash == -1 && (PyErr_Occurred() != nullptr)) {
972 PyErr_Clear();
973 is_hashable = false;
974 }
975 Py_INCREF(v._pyobject);
976 return;
977 default:
978 stub("unarchiving numeric");
979 return;
980 }
981 setflag(status_flags::evaluated | status_flags::expanded);
982 }
983
archive(archive_node & n) const984 void numeric::archive(archive_node &n) const {
985 // store type information
986 n.add_unsigned("T", t);
987
988 // create a string representation of this object
989 std::string *tstr;
990 switch (t) {
991 case LONG:
992 tstr = new std::string(std::to_string(v._long));
993 break;
994 case MPZ: {
995 std::vector<char> cp(2 + mpz_sizeinbase(v._bigint, 10));
996 mpz_get_str(&cp[0], 10, v._bigint);
997 tstr = new std::string(&cp[0]);
998 break;
999 }
1000 case MPQ: {
1001 size_t size = mpz_sizeinbase(mpq_numref(v._bigrat), 10)
1002 + mpz_sizeinbase(mpq_denref(v._bigrat), 10) + 5;
1003 std::vector<char> cp(size);
1004 mpq_get_str(&cp[0], 10, v._bigrat);
1005 tstr = new std::string(&cp[0]);
1006 break;
1007 }
1008 case PYOBJECT:
1009 tstr = py_funcs.py_dumps(v._pyobject);
1010 if (PyErr_Occurred() != nullptr) {
1011 throw(std::runtime_error("archive error: exception in py_dumps"));
1012 }
1013 break;
1014 default:
1015 stub("archive numeric");
1016 }
1017
1018 n.add_string("S", *tstr);
1019 delete tstr;
1020 inherited::archive(n);
1021 }
1022
1023 //////////
1024 // functions overriding virtual functions from base classes
1025 //////////
1026
1027 template<typename T1, typename T2>
1028 static inline bool coerce(T1& dst, const T2& arg);
1029
print_numeric(const print_context & c,const char *,const char *,const char *,const char *,unsigned level,bool latex=false) const1030 void numeric::print_numeric(const print_context & c, const char* /*unused*/,
1031 const char* /*unused*/, const char* /*unused*/, const char* /*unused*/,
1032 unsigned level, bool latex = false) const {
1033 std::string* out;
1034 PyObject* obj = to_pyobject();
1035 if (latex) {
1036 out = py_funcs.py_latex(obj, level);
1037 } else {
1038 out = py_funcs.py_repr(obj, level);
1039 }
1040 c.s << *out;
1041 Py_DECREF(obj);
1042 delete out;
1043 }
1044
do_print(const print_context & c,unsigned level) const1045 void numeric::do_print(const print_context & c, unsigned level) const {
1046 print_numeric(c, "(", ")", "I", "*", level, false);
1047 }
1048
do_print_latex(const print_latex & c,unsigned level) const1049 void numeric::do_print_latex(const print_latex & c, unsigned level) const {
1050 print_numeric(c, "{(", ")}", "i", " ", level, true);
1051 }
1052
do_print_python_repr(const print_python_repr & c,unsigned level) const1053 void numeric::do_print_python_repr(const print_python_repr & c, unsigned level) const {
1054 c.s << class_name() << "('";
1055 print_numeric(c, "(", ")", "I", "*", level);
1056 c.s << "')";
1057 }
1058
dbgstring(const numeric & num)1059 static std::string dbgstring(const numeric& num)
1060 {
1061 std::string ts;
1062 switch (num.t) {
1063 case LONG:
1064 ts = "LONG";
1065 break;
1066 case MPZ:
1067 ts = "MPZ";
1068 break;
1069 case MPQ:
1070 ts = "MPQ";
1071 break;
1072 case PYOBJECT:
1073 {
1074 ts = "PYOBJECT: ";
1075 PyObject *to = PyObject_Type(num.v._pyobject);
1076 if (to == nullptr)
1077 ts.append("NULL");
1078 else {
1079 PyObject *ro = PyObject_Repr(to);
1080 if (ro == nullptr)
1081 ts.append("NULL");
1082 else {
1083 ts.append(*py_funcs.py_repr(ro, 0));
1084 Py_DECREF(ro);
1085 }
1086 Py_DECREF(to);
1087 }
1088 break;
1089 }
1090 default: stub("typestr()");
1091 }
1092
1093 std::stringstream ss;
1094 ss << num << " (numeric)" << " @" << &num << std::hex << ", hash=0x"
1095 << num.hash << ", flags=0x" << num.flags << std::dec
1096 << ", type " << ts;
1097 return ss.str();
1098 }
1099
do_print_tree(const print_tree & c,unsigned level) const1100 void numeric::do_print_tree(const print_tree & c, unsigned level) const
1101 {
1102 c.s << std::string(level, ' ') << dbgstring(*this) << std::endl;
1103 }
1104
dbgprint() const1105 void numeric::dbgprint() const
1106 {
1107 std::cerr << dbgstring(*this) << std::endl;
1108 }
1109
info(unsigned inf) const1110 bool numeric::info(unsigned inf) const {
1111 switch (inf) {
1112 case info_flags::numeric:
1113 case info_flags::polynomial:
1114 case info_flags::rational_function:
1115 case info_flags::expanded:
1116 return true;
1117 case info_flags::real:
1118 return is_real();
1119 case info_flags::rational:
1120 case info_flags::rational_polynomial:
1121 return is_rational();
1122 case info_flags::crational:
1123 case info_flags::crational_polynomial:
1124 return is_crational();
1125 case info_flags::integer:
1126 case info_flags::integer_polynomial:
1127 return is_integer();
1128 case info_flags::cinteger:
1129 case info_flags::cinteger_polynomial:
1130 return is_cinteger();
1131 case info_flags::positive:
1132 return is_positive();
1133 case info_flags::negative:
1134 return is_negative();
1135 case info_flags::nonnegative:
1136 return is_zero() || is_positive();
1137 case info_flags::nonzero:
1138 return not is_zero();
1139 case info_flags::posint:
1140 return is_pos_integer();
1141 case info_flags::negint:
1142 return is_integer() && is_negative();
1143 case info_flags::nonnegint:
1144 return is_nonneg_integer();
1145 case info_flags::even:
1146 return is_even();
1147 case info_flags::odd:
1148 return is_odd();
1149 case info_flags::prime:
1150 return is_prime();
1151 case info_flags::algebraic:
1152 return !is_real();
1153 case info_flags::infinity:
1154 return false;
1155 case info_flags::inexact:
1156 return not is_exact();
1157 default:
1158 throw(std::runtime_error("numeric::info()"));
1159 }
1160 return false;
1161 }
1162
is_polynomial(const ex & var) const1163 bool numeric::is_polynomial(const ex & var) const {
1164 return true;
1165 }
1166
degree(const ex & s) const1167 numeric numeric::degree(const ex & s) const {
1168 // In sage deg (0) != 0 !!!
1169 return *_num0_p;
1170 }
1171
ldegree(const ex & s) const1172 numeric numeric::ldegree(const ex & s) const {
1173 return *_num0_p;
1174 }
1175
coeff(const ex & s,const ex & n) const1176 ex numeric::coeff(const ex & s, const ex & n) const {
1177 return n.is_zero() ? * this : _ex0;
1178 }
1179
1180 /** Disassemble real part and imaginary part to scan for the occurrence of a
1181 * single number. Also handles the imaginary unit. It ignores the sign on
1182 * both this and the argument, which may lead to what might appear as funny
1183 * results: (2+I).has(-2) -> true. But this is consistent, since we also
1184 * would like to have (-2+I).has(2) -> true and we want to think about the
1185 * sign as a multiplicative factor. */
has(const ex & other,unsigned) const1186 bool numeric::has(const ex &other, unsigned /*options*/) const {
1187 if (!is_exactly_a<numeric>(other))
1188 return false;
1189 const numeric &o = ex_to<numeric>(other);
1190 if (this->is_equal(o) || this->is_equal(-o))
1191 return true;
1192 if (o.imag().is_zero()) { // e.g. scan for 3 in -3*I
1193 if (!this->real().is_equal(*_num0_p))
1194 if (this->real().is_equal(o) || this->real().is_equal(-o))
1195 return true;
1196 if (!this->imag().is_equal(*_num0_p))
1197 if (this->imag().is_equal(o) || this->imag().is_equal(-o))
1198 return true;
1199 return false;
1200 }
1201 if (o.is_equal(I)) // e.g scan for I in 42*I
1202 return !this->is_real();
1203 if (o.real().is_zero() // e.g. scan for 2*I in 2*I+1
1204 and not this->imag().is_equal(*_num0_p)
1205 and (this->imag().is_equal(o * I)
1206 or this->imag().is_equal(-o * I)))
1207 return true;
1208
1209 return false;
1210 }
1211
conj() const1212 numeric numeric::conj() const {
1213 switch (t) {
1214 case LONG:
1215 case MPZ:
1216 case MPQ: return *this;
1217 case PYOBJECT: {
1218 PyObject *obj = PyObject_GetAttrString(v._pyobject,
1219 "conjugate");
1220 if (obj == nullptr)
1221 return *this;
1222 obj = PyObject_CallObject(obj, NULL);
1223 if (obj == nullptr)
1224 py_error("Error calling Python conjugate");
1225 return obj;
1226 }
1227 default:
1228 stub("invalid type: ::conjugate() type not handled");
1229 }
1230 }
1231
real_part() const1232 ex numeric::real_part() const {
1233 return real();
1234 }
1235
imag_part() const1236 ex numeric::imag_part() const {
1237 return imag();
1238 }
1239
1240 // protected
1241
1242 /** This method establishes a canonical order on all numbers. For complex
1243 * numbers this is not possible in a mathematically consistent way but we need
1244 * to establish some order and it ought to be fast. So we simply define it
1245 * to be compatible with our method csgn.
1246 *
1247 * @return csgn(*this-other) */
compare_same_type(const basic & other) const1248 int numeric::compare_same_type(const basic &other) const {
1249 GINAC_ASSERT(is_exactly_a<numeric>(other));
1250 const numeric &o = static_cast<const numeric &> (other);
1251
1252 return compare_same_type(o);
1253 }
1254
is_equal_same_type(const basic & other) const1255 bool numeric::is_equal_same_type(const basic &other) const {
1256 GINAC_ASSERT(is_exactly_a<numeric>(other));
1257 const numeric &o = static_cast<const numeric &> (other);
1258
1259 return this->is_equal(o);
1260 }
1261
calchash() const1262 long numeric::calchash() const {
1263 switch (t) {
1264 case LONG:
1265 return (v._long == -1) ? -2 : v._long;
1266 case MPZ:
1267 case MPQ:
1268 case PYOBJECT:
1269 if (is_hashable)
1270 return hash;
1271 throw(std::runtime_error("Python object not hashable"));
1272 default:
1273 stub("invalid type: ::hash() type not handled");
1274 }
1275 }
1276
1277
1278 //////////
1279 // new virtual functions which can be overridden by derived classes
1280 //////////
1281
1282 // none
1283
1284 //////////
1285 // non-virtual functions in this class
1286 //////////
1287
1288 // public
1289
1290 /** Numerical addition method. Adds argument to *this and returns result as
1291 * a numeric object. */
add(const numeric & other) const1292 const numeric numeric::add(const numeric &other) const {
1293 verbose("operator+");
1294 if (other.is_zero())
1295 return *this;
1296 if (is_zero())
1297 return other;
1298 if (t != other.t) {
1299 if (t == MPZ and other.t == MPQ) {
1300 mpq_t bigrat;
1301 mpq_init(bigrat);
1302 mpq_set_z(bigrat, v._bigint);
1303 mpq_add(bigrat, bigrat, other.v._bigrat);
1304 return bigrat;
1305 }
1306 if (t == MPQ and other.t == MPZ) {
1307 mpq_t bigrat;
1308 mpq_init(bigrat);
1309 mpq_set_z(bigrat, other.v._bigint);
1310 mpq_add(bigrat, bigrat, v._bigrat);
1311 return bigrat;
1312 }
1313
1314 numeric a, b;
1315 coerce(a, b, *this, other);
1316 return a + b;
1317 }
1318 switch (t) {
1319 case LONG: {
1320 if ((v._long > 0
1321 and v._long < std::numeric_limits<long>::max() / 2
1322 and other.v._long < std::numeric_limits<long>::max() / 2)
1323 or (v._long < 0
1324 and v._long > std::numeric_limits<long>::min() / 2
1325 and other.v._long > std::numeric_limits<long>::min() / 2))
1326 return v._long + other.v._long;
1327 mpz_t bigint;
1328 mpz_init_set_si(bigint, v._long);
1329 if (other.v._long < 0)
1330 mpz_sub_ui(bigint, bigint, -other.v._long);
1331 else
1332 mpz_add_ui(bigint, bigint, other.v._long);
1333 return bigint;
1334 }
1335 case MPZ:
1336 mpz_t bigint;
1337 mpz_init(bigint);
1338 mpz_add(bigint, v._bigint, other.v._bigint);
1339 return bigint;
1340 case MPQ:
1341 mpq_t bigrat;
1342 mpq_init(bigrat);
1343 mpq_add(bigrat, v._bigrat, other.v._bigrat);
1344 return bigrat;
1345 case PYOBJECT:
1346 return PyNumber_Add(v._pyobject, other.v._pyobject);
1347 default:
1348 stub("invalid type: operator+() type not handled");
1349 }
1350 }
1351
1352 /** Numerical subtraction method. Subtracts argument from *this and returns
1353 * result as a numeric object. */
sub(const numeric & other) const1354 const numeric numeric::sub(const numeric &other) const {
1355 verbose("operator-");
1356 if (other.is_zero())
1357 return *this;
1358 if (is_zero())
1359 return other.negative();
1360 if (t != other.t) {
1361 if (t == MPZ and other.t == MPQ) {
1362 mpq_t bigrat;
1363 mpq_init(bigrat);
1364 mpq_set_z(bigrat, v._bigint);
1365 mpq_sub(bigrat, bigrat, other.v._bigrat);
1366 return bigrat;
1367 }
1368 if (t == MPQ and other.t == MPZ) {
1369 mpq_t bigrat, tmp;
1370 mpq_init(bigrat);
1371 mpq_init(tmp);
1372 mpq_set(bigrat, v._bigrat);
1373 mpq_set_z(tmp, other.v._bigint);
1374 mpq_sub(bigrat, bigrat, tmp);
1375 mpq_clear(tmp);
1376 return bigrat;
1377 }
1378
1379 numeric a, b;
1380 coerce(a, b, *this, other);
1381 return a - b;
1382 }
1383 switch (t) {
1384 case LONG: {
1385 if ((v._long > 0
1386 and v._long < std::numeric_limits<long>::max() / 2
1387 and -other.v._long < std::numeric_limits<long>::max() / 2)
1388 or (v._long < 0
1389 and v._long > std::numeric_limits<long>::min() / 2
1390 and -other.v._long > std::numeric_limits<long>::min() / 2))
1391 return v._long - other.v._long;
1392 mpz_t bigint;
1393 mpz_init_set_si(bigint, v._long);
1394 if (other.v._long < 0)
1395 mpz_add_ui(bigint, bigint, -other.v._long);
1396 else
1397 mpz_sub_ui(bigint, bigint, other.v._long);
1398 return bigint;
1399 }
1400 case MPZ:
1401 mpz_t bigint;
1402 mpz_init(bigint);
1403 mpz_sub(bigint, v._bigint, other.v._bigint);
1404 return bigint;
1405 case MPQ:
1406 mpq_t bigrat;
1407 mpq_init(bigrat);
1408 mpq_sub(bigrat, v._bigrat, other.v._bigrat);
1409 return bigrat;
1410 case PYOBJECT:
1411 return PyNumber_Subtract(v._pyobject, other.v._pyobject);
1412 default:
1413 stub("invalid type: operator-() type not handled");
1414 }
1415 }
1416
1417 /** Numerical multiplication method. Multiplies *this and argument and returns
1418 * result as a numeric object. */
mul(const numeric & other) const1419 const numeric numeric::mul(const numeric &other) const {
1420 verbose("operator*");
1421 if ((is_zero() and t != PYOBJECT)
1422 or (other.is_zero() and other.t != PYOBJECT))
1423 return *_num0_p;
1424 if (other.is_one())
1425 return *this;
1426 if (is_one())
1427 return other;
1428 if (t != other.t) {
1429 numeric a, b;
1430 coerce(a, b, *this, other);
1431 return a * b;
1432 }
1433 switch (t) {
1434 case LONG: {
1435 static long lsqrt = std::lround(std::sqrt(
1436 std::numeric_limits<long>::max()));
1437 if (std::abs(v._long) < lsqrt
1438 and std::abs(other.v._long) < lsqrt)
1439 return v._long * other.v._long;
1440 mpz_t bigint;
1441 mpz_init_set_si(bigint, v._long);
1442 mpz_mul_si(bigint, bigint, other.v._long);
1443 return bigint;
1444 }
1445 case MPZ:
1446 mpz_t bigint;
1447 mpz_init(bigint);
1448 mpz_mul(bigint, v._bigint, other.v._bigint);
1449 return bigint;
1450 case MPQ:
1451 mpq_t bigrat;
1452 mpq_init(bigrat);
1453 mpq_mul(bigrat, v._bigrat, other.v._bigrat);
1454 return bigrat;
1455 case PYOBJECT:
1456 return PyNumber_Multiply(v._pyobject, other.v._pyobject);
1457 default:
1458 stub("invalid type: operator*() type not handled");
1459 }
1460 }
1461
1462 /** Numerical division method. Divides *this by argument and returns result as
1463 * a numeric object.
1464 *
1465 * @exception overflow_error (division by zero) */
div(const numeric & other) const1466 const numeric numeric::div(const numeric &other) const {
1467 verbose("operator/");
1468 if (other.is_zero())
1469 throw std::overflow_error("numeric::div(): division by zero");
1470 if (is_zero())
1471 return *_num0_p;
1472 if (other.is_one())
1473 return *this;
1474 if (t != other.t) {
1475 numeric a, b;
1476 coerce(a, b, *this, other);
1477 return a / b;
1478 }
1479 switch (t) {
1480 case LONG: {
1481 if (v._long == 1 and other.v._long == 2)
1482 return *_num1_2_p;
1483 auto ld = std::div(v._long, other.v._long);
1484 if (ld.rem == 0)
1485 return ld.quot;
1486
1487 mpq_t bigrat, obigrat;
1488 mpq_init(bigrat);
1489 mpq_init(obigrat);
1490 mpq_set_si(bigrat, v._long, 1);
1491 mpq_set_si(obigrat, other.v._long, 1);
1492 mpq_div(bigrat, bigrat, obigrat);
1493 mpq_clear(obigrat);
1494 return bigrat;
1495 }
1496 case MPZ: {
1497 if (mpz_divisible_p(v._bigint, other.v._bigint)) {
1498 mpz_t bigint;
1499 mpz_init(bigint);
1500 mpz_divexact(bigint, v._bigint,
1501 other.v._bigint);
1502 return bigint;
1503 }
1504 mpq_t bigrat, obigrat;
1505 mpq_init(bigrat);
1506 mpq_init(obigrat);
1507 mpq_set_z(bigrat, v._bigint);
1508 mpq_set_z(obigrat, other.v._bigint);
1509 mpq_div(bigrat, bigrat, obigrat);
1510 mpq_clear(obigrat);
1511 return bigrat;
1512 }
1513 case MPQ: {
1514 mpq_t bigrat;
1515 mpq_init(bigrat);
1516 mpq_div(bigrat, v._bigrat, other.v._bigrat);
1517 return bigrat;
1518 }
1519 case PYOBJECT:
1520 #if PY_MAJOR_VERSION < 3
1521 if (PyObject_Compare(other.v._pyobject, ONE) == 0
1522 and py_funcs.py_is_integer(other.v._pyobject) != 0) {
1523 return *this;
1524 }
1525 if (PyInt_Check(v._pyobject)) {
1526 if (PyInt_Check(other.v._pyobject)) {
1527 // This branch happens at startup.
1528 PyObject *o = PyNumber_TrueDivide(Integer(PyInt_AsLong(v._pyobject)),
1529 Integer(PyInt_AsLong(other.v._pyobject)));
1530 // I don't 100% understand why I have to incref this,
1531 // but if I don't, Sage crashes on exit.
1532 Py_INCREF(o);
1533 return o;
1534 }
1535 if (PyLong_Check(other.v._pyobject)) {
1536 PyObject *d = py_funcs.
1537 py_integer_from_python_obj(other.v._pyobject);
1538 PyObject *ans = PyNumber_TrueDivide(v._pyobject,
1539 d);
1540 Py_DECREF(d);
1541 return ans;
1542 }
1543 }
1544 #endif
1545 if (PyLong_Check(v._pyobject)) {
1546 PyObject *n = py_funcs.
1547 py_integer_from_python_obj(v._pyobject);
1548 PyObject *ans = PyNumber_TrueDivide(n,
1549 other.v._pyobject);
1550 Py_DECREF(n);
1551 return ans;
1552 }
1553 return PyNumber_TrueDivide(v._pyobject, other.v._pyobject);
1554
1555 default:
1556 stub("invalid type: operator/() type not handled");
1557 }
1558 }
1559
1560 // Compute `a^b` as an integer, if it is integral, or return ``false``.
1561 // The nonnegative real root is taken for even denominators.
integer_rational_power(numeric & res,const numeric & a,const numeric & b)1562 bool numeric::integer_rational_power(numeric& res,
1563 const numeric& a, const numeric& b)
1564 {
1565 if (b.t != MPQ)
1566 throw std::runtime_error("integer_rational_power: bad input");
1567 if (mpz_sgn(mpq_numref(b.v._bigrat)) < 0)
1568 throw std::runtime_error("integer_rational_power: bad input");
1569 if (a.t == LONG) {
1570 if (a.v._long == 1
1571 or mpz_cmp_ui(mpq_numref(b.v._bigrat), 0) == 0) {
1572 res = 1;
1573 return true;
1574 }
1575 if (a.v._long == 0) {
1576 res = 0;
1577 return true;
1578 }
1579 if (a.v._long < 0
1580 and mpz_cmp_ui(mpq_denref(b.v._bigrat), 1))
1581 return false;
1582 long z;
1583 if (not mpz_fits_ulong_p(mpq_numref(b.v._bigrat))
1584 or not mpz_fits_ulong_p(mpq_denref(b.v._bigrat)))
1585 // too big to take roots/powers
1586 return false;
1587 if (b.is_equal(*_num1_2_p)) {
1588 z = std::lround(std::sqrt(a.v._long));
1589 if (a.v._long == z*z) {
1590 res = numeric(z);
1591 return true;
1592 }
1593 return false;
1594 }
1595 return integer_rational_power(res, a.to_bigint(), b);
1596 }
1597 if (a.t != MPZ)
1598 throw std::runtime_error("integer_rational_power: bad input");
1599 int sgn = mpz_sgn(a.v._bigint);
1600 mpz_t z;
1601 mpz_init(z);
1602 mpz_set_ui(z, 0);
1603 if (mpz_cmp_ui(a.v._bigint, 1) == 0
1604 or mpz_cmp_ui(mpq_numref(b.v._bigrat), 0) == 0)
1605 mpz_set_ui(z, 1);
1606 else if (sgn == 0) {
1607 res = *_num0_p;
1608 return true;
1609 }
1610 else if (sgn < 0 and mpz_cmp_ui(mpq_denref(b.v._bigrat), 1))
1611 return false;
1612 else {
1613 if (not mpz_fits_ulong_p(mpq_numref(b.v._bigrat))
1614 or not mpz_fits_ulong_p(mpq_denref(b.v._bigrat)))
1615 // too big to take roots/powers
1616 return false;
1617 if (mpz_cmp_ui(mpq_denref(b.v._bigrat), 2) == 0) {
1618 if (mpz_perfect_square_p(a.v._bigint))
1619 mpz_sqrt(z, a.v._bigint);
1620 else
1621 return false;
1622 }
1623 else {
1624 bool exact = mpz_root(z, a.v._bigint,
1625 mpz_get_ui(mpq_denref(b.v._bigrat)));
1626 if (not exact)
1627 return false;
1628 }
1629 mpz_pow_ui(z, z, mpz_get_ui(mpq_numref(b.v._bigrat)));
1630 }
1631 res = numeric(z);
1632 return true;
1633 }
1634
1635 // for a^b return c,d such that a^b = c*d^b
1636 // only for MPZ/MPQ base and MPQ exponent
rational_power_parts(const numeric & a_orig,const numeric & b_orig,numeric & c,numeric & d,bool & c_unit)1637 void rational_power_parts(const numeric& a_orig, const numeric& b_orig,
1638 numeric& c, numeric& d, bool& c_unit)
1639 {
1640 if (b_orig.t == LONG or b_orig.t == MPZ) {
1641 c = a_orig.pow_intexp(b_orig);
1642 d = *_num1_p;
1643 c_unit = c.is_one();
1644 return;
1645 }
1646
1647 if ((a_orig.t != LONG
1648 and a_orig.t != MPZ
1649 and a_orig.t != MPQ)
1650 or b_orig.t != MPQ) {
1651 d = a_orig;
1652 c = *_num1_p;
1653 c_unit = true;
1654 return;
1655 }
1656 bool b_negative = b_orig.is_negative();
1657 const numeric& a = b_negative? a_orig.inverse() : a_orig;
1658 const numeric& b = b_negative? b_orig.negative() : b_orig;
1659 if (a.t == MPQ) {
1660 numeric c1, c2, d1, d2;
1661 rational_power_parts(a.numer(), b, c1, d1, c_unit);
1662 rational_power_parts(a.denom(), b, c2, d2, c_unit);
1663 c = c1 / c2;
1664 if (b_negative)
1665 d = d2 / d1;
1666 else
1667 d = d1 / d2;
1668 c_unit = c.is_one();
1669 return;
1670 }
1671
1672 // remaining a types are LONG and MPZ
1673 if (numeric::integer_rational_power(c, a, b)) {
1674 c_unit = c.is_one();
1675 d = *_num1_p;
1676 return;
1677 }
1678 numeric numer = b.numer();
1679 numeric denom = b.denom();
1680 if (denom.t == MPZ
1681 and not mpz_fits_slong_p(denom.v._bigint)) {
1682 c = *_num1_p;
1683 c_unit = true;
1684 d = a;
1685 return;
1686 }
1687 long denoml = denom.to_long();
1688 if (denoml > 1 and a.is_minus_one()) {
1689 c = *_num1_p;
1690 c_unit = true;
1691 d = *_num_1_p;
1692 return;
1693 }
1694 std::vector<std::pair<numeric, int>> factors;
1695 long max_prime_idx;
1696 static numeric maxnum = numeric(10).power(200);
1697 if (a.abs() < maxnum) {
1698 double ad;
1699 if (a.t == MPZ)
1700 ad = mpz_get_d(a.v._bigint);
1701 else
1702 ad = a.v._long;
1703 ad = std::abs(ad);
1704 double adrootlog = std::log(ad) / denoml;
1705 double adroot = std::exp(adrootlog);
1706 double mpid = adroot*1.25506/adrootlog;
1707 max_prime_idx = mpid>2000? 2000L : long(mpid) + 1;
1708 }
1709 else // can't convert a to double
1710 max_prime_idx = 2000L;
1711 a.factor(factors, max_prime_idx);
1712 c = *_num1_p;
1713 d = *_num1_p;
1714 for (auto p : factors) {
1715 c = c * p.first.pow_intexp(numer * numeric(p.second /
1716 denoml));
1717 d = d * p.first.pow_intexp(numeric(p.second % denoml));
1718 }
1719 if (a.is_negative() and numer.is_odd())
1720 d = d.negative();
1721 if (b_negative)
1722 d = d.inverse();
1723 c_unit = c.is_one();
1724 }
1725
power(signed long exp_si) const1726 const numeric numeric::power(signed long exp_si) const
1727 {
1728 PyObject *o, *r;
1729 if (exp_si == 0)
1730 return *_num1_p;
1731 if (exp_si == 1)
1732 return *this;
1733 switch (t) {
1734 case LONG:
1735 if (exp_si >= 0) {
1736 mpz_t bigint;
1737 mpz_init(bigint);
1738 mpz_set_si(bigint, v._long);
1739 mpz_pow_ui(bigint, bigint, exp_si);
1740 return numeric(bigint);
1741 } else {
1742 mpz_t bigint;
1743 mpz_init_set_si(bigint, v._long);
1744 if (exp_si != -1)
1745 mpz_pow_ui(bigint, bigint, -exp_si);
1746 mpq_t bigrat;
1747 mpq_init(bigrat);
1748 mpq_set_z(bigrat, bigint);
1749 mpq_inv(bigrat, bigrat);
1750 mpz_clear(bigint);
1751 return numeric(bigrat);
1752 }
1753 case MPZ:
1754 if (exp_si >= 0) {
1755 mpz_t bigint;
1756 mpz_init(bigint);
1757 mpz_pow_ui(bigint, v._bigint, exp_si);
1758 return numeric(bigint);
1759 } else {
1760 mpz_t bigint;
1761 mpz_init_set(bigint, v._bigint);
1762 mpz_pow_ui(bigint, bigint, -exp_si);
1763 mpq_t bigrat;
1764 mpq_init(bigrat);
1765 mpq_set_z(bigrat, bigint);
1766 mpq_inv(bigrat, bigrat);
1767 mpz_clear(bigint);
1768 return numeric(bigrat);
1769 }
1770 case MPQ:
1771 mpz_t bigint;
1772 mpq_t bigrat, obigrat;
1773 mpz_init(bigint);
1774 mpq_init(bigrat);
1775 mpq_init(obigrat);
1776 if (exp_si >= 0) {
1777 mpz_pow_ui(bigint, mpq_numref(v._bigrat), exp_si);
1778 mpq_set_z(bigrat, bigint);
1779 mpz_pow_ui(bigint, mpq_denref(v._bigrat), exp_si);
1780 mpq_set_z(obigrat, bigint);
1781 mpq_div(bigrat, bigrat, obigrat);
1782 } else {
1783 mpz_pow_ui(bigint, mpq_denref(v._bigrat), -exp_si);
1784 mpq_set_z(bigrat, bigint);
1785 mpz_pow_ui(bigint, mpq_numref(v._bigrat), -exp_si);
1786 mpq_set_z(obigrat, bigint);
1787 mpq_div(bigrat, bigrat, obigrat);
1788 }
1789 mpz_clear(bigint);
1790 mpq_clear(obigrat);
1791 return numeric(bigrat);
1792 case PYOBJECT:
1793 o = Integer(exp_si);
1794 r = PyNumber_Power(v._pyobject, o, Py_None);
1795 Py_DECREF(o);
1796 return numeric(r);
1797 default:
1798 stub("invalid type: pow_intexp numeric");
1799 }
1800 }
1801
pow_intexp(const numeric & exponent) const1802 const numeric numeric::pow_intexp(const numeric &exponent) const
1803 {
1804 if (not exponent.is_integer())
1805 throw std::runtime_error("numeric::pow_intexp: exponent not integer");
1806 // Because Sage may have simplified a fraction to an integer,
1807 // we should check whether the type is MPQ as well as MPZ
1808 if (exponent.t == MPZ) {
1809 if (not mpz_fits_sint_p(exponent.v._bigint))
1810 throw std::runtime_error("size of exponent exceeds signed long size");
1811 return power(mpz_get_si(exponent.v._bigint));
1812 }
1813 if (exponent.t == MPQ) {
1814 if (not mpz_fits_sint_p(mpq_numref(exponent.v._bigrat)))
1815 throw std::runtime_error("size of exponent exceeds signed long size");
1816 return power(mpz_get_si(mpq_numref(exponent.v._bigrat)));
1817 }
1818 return power(exponent.v._long);
1819 }
1820
1821 /** Numerical exponentiation.
1822 * Raises *this to the power given as argument and returns the result as ex,
1823 * because it is possible that the result is converted to symbolic by Sage,
1824 * and we don't want symbolic as PyObject inside numeric. */
power(const numeric & exponent) const1825 const ex numeric::power(const numeric &exponent) const {
1826 verbose("pow");
1827 numeric expo(exponent);
1828
1829 // any PyObjects castable to long are casted
1830 if (exponent.t == PYOBJECT) {
1831 #if PY_MAJOR_VERSION < 3
1832 if (PyInt_Check(exponent.v._pyobject)) {
1833 long si = PyInt_AsLong(exponent.v._pyobject);
1834 if (si == -1 and PyErr_Occurred())
1835 PyErr_Clear();
1836 else {
1837 expo.t = MPZ;
1838 mpz_set_si(expo.v._bigint, si);
1839 }
1840 } else
1841 #endif
1842 if (PyLong_Check(exponent.v._pyobject)) {
1843 expo.t = MPZ;
1844 _mpz_set_pylong(expo.v._bigint,
1845 reinterpret_cast<PyLongObject *>(exponent.v._pyobject));
1846 }
1847 }
1848
1849 // inexact PyObjects in base or exponent
1850 if (t == PYOBJECT and not is_exact()) {
1851 if (expo.t == PYOBJECT)
1852 return numeric(PyNumber_Power(v._pyobject,
1853 exponent.v._pyobject,
1854 Py_None));
1855 else {
1856 PyObject* obj = exponent.to_pyobject();
1857 const numeric& ret = numeric(PyNumber_Power(v._pyobject,
1858 obj,
1859 Py_None));
1860 Py_DECREF(obj);
1861 return ret;
1862 }
1863 }
1864 if (expo.t == PYOBJECT and not expo.is_exact()) {
1865 if (t == PYOBJECT)
1866 return numeric(PyNumber_Power(v._pyobject,
1867 exponent.v._pyobject,
1868 Py_None));
1869 else {
1870 PyObject* obj = to_pyobject();
1871 const numeric& ret = numeric(PyNumber_Power(obj,
1872 exponent.v._pyobject,
1873 Py_None));
1874 Py_DECREF(obj);
1875 return ret;
1876 }
1877 }
1878
1879 // handle all integer exponents
1880 if (expo.t == LONG or expo.t == MPZ)
1881 return pow_intexp(expo);
1882 if (expo.t == MPQ and expo.is_integer())
1883 return power(exponent.to_long());
1884
1885 using POW = class power;
1886 using MUL = class mul;
1887 // return any complex power unchanged
1888 if (not is_rational()
1889 or not exponent.is_rational())
1890 return (new POW(*this, expo))->
1891 setflag(status_flags::dynallocated
1892 | status_flags::evaluated);
1893 // rational exponent
1894 if ((t == LONG or t == MPZ or t == MPQ) and expo.t == MPQ) {
1895 numeric c, d;
1896 bool c_unit;
1897 rational_power_parts(*this, expo, c, d, c_unit);
1898 if (d.is_one())
1899 return c;
1900 if (d.is_minus_one()
1901 and expo.denom().is_equal(*_num2_p)) {
1902 switch (expo.numer().to_long() % 4) {
1903 case 0: return c;
1904 case 1: return I * c;
1905 case 2: return -c;
1906 case 3: return -I * c;
1907 };
1908 }
1909 else if (not c_unit)
1910 return d.power(expo) * c;
1911 if (exponent.is_negative()) {
1912 long int_exp = -(expo.to_long());
1913 numeric nexp = expo + numeric(int_exp);
1914 ex p = (new POW(*this, nexp))->
1915 setflag(status_flags::dynallocated
1916 | status_flags::evaluated);
1917 return p * pow_intexp(int_exp).inverse();
1918 }
1919 if (expo < *_num1_p)
1920 return (new POW(*this, expo))->
1921 setflag(status_flags::dynallocated
1922 | status_flags::evaluated);
1923 // rational expo > 1
1924 const numeric m = expo.denom();
1925 numeric rem;
1926 numeric quo = expo.numer().iquo(m, rem);
1927 if (rem.is_negative()) {
1928 rem += m;
1929 quo -= *_num1_p;
1930 }
1931 if (quo.is_zero()) {
1932 ex res = (new POW(d, expo))->
1933 setflag(status_flags::dynallocated
1934 | status_flags::evaluated);
1935 if (c_unit)
1936 return (new MUL(res, c))->
1937 setflag(status_flags::dynallocated
1938 | status_flags::evaluated);
1939 }
1940 if (rem.is_zero())
1941 return this->power(quo);
1942 return (new MUL(this->power(rem / m), this->power(quo)))->
1943 setflag(status_flags::dynallocated
1944 | status_flags::evaluated);
1945 }
1946
1947 // exact PyObjects in base or exponent
1948 if (t == PYOBJECT and exponent.t == PYOBJECT)
1949 return numeric(PyNumber_Power(v._pyobject,
1950 exponent.v._pyobject,
1951 Py_None));
1952 if (t == PYOBJECT) {
1953 PyObject* obj = exponent.to_pyobject();
1954 const numeric& ret = numeric(PyNumber_Power(v._pyobject,
1955 obj,
1956 Py_None));
1957 Py_DECREF(obj);
1958 return ret;
1959 }
1960 if (expo.t == PYOBJECT) {
1961 PyObject* obj = to_pyobject();
1962 const numeric& ret = numeric(PyNumber_Power(to_pyobject(),
1963 obj,
1964 Py_None));
1965 Py_DECREF(obj);
1966 return ret;
1967 }
1968
1969 throw std::runtime_error("numeric::power: can't happen");
1970 }
1971
1972
1973 /** Numerical addition method. Adds argument to *this and returns result as
1974 * a numeric object on the heap. Use internally only for direct wrapping into
1975 * an ex object, where the result would end up on the heap anyways. */
add_dyn(const numeric & other) const1976 const numeric &numeric::add_dyn(const numeric &other) const {
1977 // Efficiency shortcut: trap the neutral element by pointer. This hack
1978 // is supposed to keep the number of distinct numeric objects low.
1979 if (this == _num0_p)
1980 return other;
1981 if (&other == _num0_p)
1982 return *this;
1983
1984 return static_cast<const numeric &> ((new numeric(*this + other))->
1985 setflag(status_flags::dynallocated));
1986 }
1987
1988 /** Numerical subtraction method. Subtracts argument from *this and returns
1989 * result as a numeric object on the heap. Use internally only for direct
1990 * wrapping into an ex object, where the result would end up on the heap
1991 * anyways. */
sub_dyn(const numeric & other) const1992 const numeric &numeric::sub_dyn(const numeric &other) const {
1993 // Efficiency shortcut: trap the neutral exponent (first by pointer). This
1994 // hack is supposed to keep the number of distinct numeric objects low.
1995 if (&other == _num0_p || (other.is_zero()))
1996 return *this;
1997
1998 return static_cast<const numeric &> ((new numeric(*this - other))->
1999 setflag(status_flags::dynallocated));
2000 }
2001
2002 /** Numerical multiplication method. Multiplies *this and argument and returns
2003 * result as a numeric object on the heap. Use internally only for direct
2004 * wrapping into an ex object, where the result would end up on the heap
2005 * anyways. */
mul_dyn(const numeric & other) const2006 const numeric &numeric::mul_dyn(const numeric &other) const {
2007 // Efficiency shortcut: trap the neutral element by pointer. This hack
2008 // is supposed to keep the number of distinct numeric objects low.
2009 if (this == _num1_p)
2010 return other;
2011 if (&other == _num1_p)
2012 return *this;
2013
2014 return static_cast<const numeric &> ((new numeric(*this * other))->
2015 setflag(status_flags::dynallocated));
2016 }
2017
2018 /** Numerical division method. Divides *this by argument and returns result as
2019 * a numeric object on the heap. Use internally only for direct wrapping
2020 * into an ex object, where the result would end up on the heap
2021 * anyways.
2022 *
2023 * @exception overflow_error (division by zero) */
div_dyn(const numeric & other) const2024 const numeric &numeric::div_dyn(const numeric &other) const {
2025 // Efficiency shortcut: trap the neutral element by pointer. This hack
2026 // is supposed to keep the number of distinct numeric objects low.
2027 if (&other == _num1_p)
2028 return *this;
2029 if (other.is_zero())
2030 throw std::overflow_error("division by zero");
2031 return static_cast<const numeric &> ((new numeric(*this / other))->
2032 setflag(status_flags::dynallocated));
2033 }
2034
2035 /** Numerical exponentiation. Raises *this to the power given as argument and
2036 * returns result as a numeric object on the heap. Use internally only for
2037 * direct wrapping into an ex object, where the result would end up on the
2038 * heap anyways. */
2039 // NOTE: Use only if sure to return a non-expression (see power()).
power_dyn(const numeric & other) const2040 const numeric &numeric::power_dyn(const numeric &other) const {
2041 // Efficiency shortcut: trap the neutral exponent (first try by pointer, then
2042 // try harder, since calls to cln::expt() below may return amazing results for
2043 // floating point exponent 1.0).
2044 if (&other == _num1_p || (other == *_num1_p))
2045 return *this;
2046
2047 return static_cast<const numeric &> ((new numeric(ex_to<numeric>(pow(*this, other))))->setflag(status_flags::dynallocated));
2048 }
2049
operator =(int i)2050 const numeric &numeric::operator=(int i) {
2051 return operator=(numeric(i));
2052 }
2053
operator =(unsigned int i)2054 const numeric &numeric::operator=(unsigned int i) {
2055 return operator=(numeric(i));
2056 }
2057
operator =(long i)2058 const numeric &numeric::operator=(long i) {
2059 return operator=(numeric(i));
2060 }
2061
operator =(double d)2062 const numeric &numeric::operator=(double d) {
2063 return operator=(numeric(d));
2064 }
2065
negative() const2066 const numeric numeric::negative() const {
2067 verbose("operator-");
2068 switch (t) {
2069 case LONG:
2070 return -v._long;
2071 case MPZ:
2072 mpz_t bigint;
2073 mpz_init_set(bigint, v._bigint);
2074 mpz_neg(bigint, bigint);
2075 return bigint;
2076 case MPQ:
2077 mpq_t bigrat;
2078 mpq_init(bigrat);
2079 mpq_set(bigrat, v._bigrat);
2080 mpq_neg(bigrat, bigrat);
2081 return bigrat;
2082 case PYOBJECT:
2083 return PyNumber_Negative(v._pyobject);
2084 default:
2085 stub("invalid type: operator-() type not handled");
2086 }
2087 }
2088
2089 // binary arithmetic assignment operators with numeric
2090
operator +=(numeric & lh,const numeric & rh)2091 numeric & operator+=(numeric & lh, const numeric & rh)
2092 {
2093 if (rh.is_zero())
2094 return lh;
2095 if (lh.is_zero()) {
2096 lh = rh;
2097 return lh;
2098 }
2099 if (lh.t != rh.t) {
2100 if (lh.t == MPZ and rh.t == MPQ) {
2101 mpz_t bigint;
2102 mpz_init_set(bigint, lh.v._bigint);
2103 mpz_clear(lh.v._bigint);
2104 lh.t = MPQ;
2105 mpq_init(lh.v._bigrat);
2106 mpq_set_z(lh.v._bigrat, bigint);
2107 mpq_add(lh.v._bigrat, lh.v._bigrat, rh.v._bigrat);
2108 lh.hash = _mpq_pythonhash(lh.v._bigrat);
2109 mpz_clear(bigint);
2110 return lh;
2111 }
2112 if (lh.t == MPQ and rh.t == MPZ) {
2113 mpq_t tmp;
2114 mpq_init(tmp);
2115 mpq_set_z(tmp, rh.v._bigint);
2116 mpq_add(lh.v._bigrat, lh.v._bigrat, tmp);
2117 lh.hash = _mpq_pythonhash(lh.v._bigrat);
2118 mpq_clear(tmp);
2119 return lh;
2120 }
2121
2122 numeric a, b;
2123 coerce(a, b, lh, rh);
2124 lh = a + b;
2125 return lh;
2126 }
2127 switch (lh.t) {
2128 case LONG: {
2129 if ((lh.v._long > 0
2130 and lh.v._long < std::numeric_limits<long>::max() / 2
2131 and rh.v._long < std::numeric_limits<long>::max() / 2)
2132 or (lh.v._long < 0
2133 and lh.v._long > std::numeric_limits<long>::min() / 2
2134 and rh.v._long > std::numeric_limits<long>::min() / 2)) {
2135 lh.v._long += rh.v._long;
2136 lh.hash = (lh.v._long == -1) ? -2 : lh.v._long;
2137 return lh;
2138 }
2139 lh.t = MPZ;
2140 mpz_init_set_si(lh.v._bigint, lh.v._long);
2141 if (rh.v._long < 0)
2142 mpz_sub_ui(lh.v._bigint,
2143 lh.v._bigint, -rh.v._long);
2144 else
2145 mpz_add_ui(lh.v._bigint, lh.v._bigint,
2146 rh.v._long);
2147 lh.hash = _mpz_pythonhash(lh.v._bigint);
2148 return lh;
2149 }
2150 case MPZ:
2151 mpz_add(lh.v._bigint, lh.v._bigint, rh.v._bigint);
2152 lh.hash = _mpz_pythonhash(lh.v._bigint);
2153 return lh;
2154 case MPQ:
2155 mpq_add(lh.v._bigrat, lh.v._bigrat, rh.v._bigrat);
2156 lh.hash = _mpq_pythonhash(lh.v._bigrat);
2157 return lh;
2158 case PYOBJECT: {
2159 PyObject *p = lh.v._pyobject;
2160 lh.v._pyobject = PyNumber_Add(p, rh.v._pyobject);
2161 if (lh.v._pyobject == nullptr) {
2162 lh.v._pyobject = p;
2163 py_error("numeric operator+=");
2164 }
2165 lh.hash = PyObject_Hash(lh.v._pyobject);
2166 Py_DECREF(p);
2167 return lh;
2168 }
2169 default:
2170 stub("invalid type: operator+=() type not handled");
2171 }
2172 }
2173
operator -=(numeric & lh,const numeric & rh)2174 numeric & operator-=(numeric & lh, const numeric & rh)
2175 {
2176 if (rh.is_zero())
2177 return lh;
2178 if (lh.is_zero()) {
2179 lh = rh.negative();
2180 return lh;
2181 }
2182 if (lh.t != rh.t) {
2183 if (lh.t == MPZ and rh.t == MPQ) {
2184 mpz_t bigint;
2185 mpz_init_set(bigint, lh.v._bigint);
2186 mpz_clear(lh.v._bigint);
2187 lh.t = MPQ;
2188 mpq_init(lh.v._bigrat);
2189 mpq_set_z(lh.v._bigrat, bigint);
2190 mpq_sub(lh.v._bigrat, lh.v._bigrat, rh.v._bigrat);
2191 lh.hash = _mpq_pythonhash(lh.v._bigrat);
2192 mpz_clear(bigint);
2193 return lh;
2194 }
2195 if (lh.t == MPQ and rh.t == MPZ) {
2196 mpq_t tmp;
2197 mpq_init(tmp);
2198 mpq_set_z(tmp, rh.v._bigint);
2199 mpq_sub(lh.v._bigrat, lh.v._bigrat, tmp);
2200 lh.hash = _mpq_pythonhash(lh.v._bigrat);
2201 mpq_clear(tmp);
2202 return lh;
2203 }
2204
2205 numeric a, b;
2206 coerce(a, b, lh, rh);
2207 lh = a - b;
2208 return lh;
2209 }
2210 switch (lh.t) {
2211 case LONG: {
2212 if ((lh.v._long > 0
2213 and lh.v._long < std::numeric_limits<long>::max() / 2
2214 and -rh.v._long < std::numeric_limits<long>::max() / 2)
2215 or (lh.v._long < 0
2216 and lh.v._long > std::numeric_limits<long>::min() / 2
2217 and -rh.v._long > std::numeric_limits<long>::min() / 2)) {
2218 lh.v._long -= rh.v._long;
2219 lh.hash = (lh.v._long == -1) ? -2 : lh.v._long;
2220 return lh;
2221 }
2222 lh.t = MPZ;
2223 mpz_init_set_si(lh.v._bigint, lh.v._long);
2224 if (rh.v._long < 0)
2225 mpz_add_ui(lh.v._bigint,
2226 lh.v._bigint, -rh.v._long);
2227 else
2228 mpz_sub_ui(lh.v._bigint, lh.v._bigint,
2229 rh.v._long);
2230 lh.hash = _mpz_pythonhash(lh.v._bigint);
2231 return lh;
2232 }
2233 case MPZ:
2234 mpz_sub(lh.v._bigint, lh.v._bigint, rh.v._bigint);
2235 lh.hash = _mpz_pythonhash(lh.v._bigint);
2236 return lh;
2237 case MPQ:
2238 mpq_sub(lh.v._bigrat, lh.v._bigrat, rh.v._bigrat);
2239 lh.hash = _mpq_pythonhash(lh.v._bigrat);
2240 return lh;
2241 case PYOBJECT: {
2242 PyObject *p = lh.v._pyobject;
2243 lh.v._pyobject = PyNumber_Subtract(p, rh.v._pyobject);
2244 if (lh.v._pyobject == nullptr) {
2245 lh.v._pyobject = p;
2246 py_error("numeric operator-=");
2247 }
2248 lh.hash = PyObject_Hash(lh.v._pyobject);
2249 Py_DECREF(p);
2250 return lh;
2251 }
2252 default:
2253 stub("invalid type: operator-() type not handled");
2254 }
2255 }
2256
operator *=(numeric & lh,const numeric & rh)2257 numeric & operator*=(numeric & lh, const numeric & rh)
2258 {
2259 if (rh.is_one())
2260 return lh;
2261 if (lh.is_one()) {
2262 lh = rh;
2263 return lh;
2264 }
2265 if ((lh.is_zero() and lh.t != PYOBJECT)
2266 or (rh.is_zero() and rh.t != PYOBJECT)) {
2267 lh = *_num0_p;
2268 return lh;
2269 }
2270 if (lh.t != rh.t) {
2271 if (lh.t == MPZ and rh.t == MPQ) {
2272 mpq_t tmp;
2273 mpq_init(tmp);
2274 mpq_set_z(tmp, lh.v._bigint);
2275 mpq_mul(tmp, tmp, rh.v._bigrat);
2276 if (mpz_cmp_ui(mpq_denref(tmp),1) == 0) {
2277 mpz_set(lh.v._bigint, mpq_numref(tmp));
2278 lh.hash = _mpz_pythonhash(lh.v._bigint);
2279 mpq_clear(tmp);
2280 return lh;
2281 }
2282 mpz_clear(lh.v._bigint);
2283 lh.t = MPQ;
2284 mpq_init(lh.v._bigrat);
2285 mpq_set(lh.v._bigrat, tmp);
2286 lh.hash = _mpq_pythonhash(lh.v._bigrat);
2287 mpq_clear(tmp);
2288 return lh;
2289 }
2290 if (lh.t == MPQ and rh.t == MPZ) {
2291 mpq_t tmp;
2292 mpq_init(tmp);
2293 mpq_set_z(tmp, rh.v._bigint);
2294 mpq_mul(tmp, tmp, lh.v._bigrat);
2295 if (mpz_cmp_ui(mpq_denref(tmp),1) != 0) {
2296 mpq_set(lh.v._bigrat, tmp);
2297 lh.hash = _mpq_pythonhash(lh.v._bigrat);
2298 mpq_clear(tmp);
2299 return lh;
2300 }
2301 mpq_clear(lh.v._bigrat);
2302 lh.t = MPZ;
2303 mpz_init(lh.v._bigint);
2304 mpz_set(lh.v._bigint, mpq_numref(tmp));
2305 lh.hash = _mpz_pythonhash(lh.v._bigint);
2306 mpq_clear(tmp);
2307 return lh;
2308 }
2309
2310 numeric a, b;
2311 coerce(a, b, lh, rh);
2312 lh = a * b;
2313 return lh;
2314 }
2315 switch (lh.t) {
2316 case LONG: {
2317 static long lsqrt = std::lround(std::sqrt(
2318 std::numeric_limits<long>::max()));
2319 if (std::abs(lh.v._long) < lsqrt
2320 and std::abs(rh.v._long) < lsqrt) {
2321 lh.v._long *= rh.v._long;
2322 return lh;
2323 }
2324 lh.t = MPZ;
2325 mpz_init_set_si(lh.v._bigint, lh.v._long);
2326 mpz_mul_si(lh.v._bigint, lh.v._bigint, rh.v._long);
2327 lh.hash = _mpz_pythonhash(lh.v._bigint);
2328 return lh;
2329 }
2330 case MPZ:
2331 mpz_mul(lh.v._bigint, lh.v._bigint, rh.v._bigint);
2332 lh.hash = _mpz_pythonhash(lh.v._bigint);
2333 return lh;
2334 case MPQ:
2335 mpq_mul(lh.v._bigrat, lh.v._bigrat, rh.v._bigrat);
2336 lh.hash = _mpq_pythonhash(lh.v._bigrat);
2337 return lh;
2338 case PYOBJECT: {
2339 PyObject *p = lh.v._pyobject;
2340 lh.v._pyobject = PyNumber_Multiply(p, rh.v._pyobject);
2341 if (lh.v._pyobject == nullptr) {
2342 lh.v._pyobject = p;
2343 py_error("numeric operator*=");
2344 }
2345 lh.hash = PyObject_Hash(lh.v._pyobject);
2346 Py_DECREF(p);
2347 return lh;
2348 }
2349 default:
2350 stub("invalid type: operator*=() type not handled");
2351 }
2352 }
2353
sgn(T val)2354 template <typename T> int sgn(T val) {
2355 return (T(0) < val) - (val < T(0));
2356 }
operator /=(numeric & lh,const numeric & rh)2357 numeric & operator/=(numeric & lh, const numeric & rh)
2358 {
2359 if (rh.is_zero())
2360 throw std::overflow_error("numeric::/=(): division by zero");
2361 if (rh.is_one())
2362 return lh;
2363 if (lh.t != rh.t) {
2364 if (lh.t == MPZ and rh.t == MPQ) {
2365 mpq_t tmp;
2366 mpq_init(tmp);
2367 mpq_set_z(tmp, lh.v._bigint);
2368 mpq_div(tmp, tmp, rh.v._bigrat);
2369 if (mpz_cmp_ui(mpq_denref(tmp),1) == 0) {
2370 mpz_set(lh.v._bigint, mpq_numref(tmp));
2371 lh.hash = _mpz_pythonhash(lh.v._bigint);
2372 mpq_clear(tmp);
2373 return lh;
2374 }
2375 mpz_clear(lh.v._bigint);
2376 lh.t = MPQ;
2377 mpq_init(lh.v._bigrat);
2378 mpq_set(lh.v._bigrat, tmp);
2379 lh.hash = _mpq_pythonhash(lh.v._bigrat);
2380 mpq_clear(tmp);
2381 return lh;
2382 }
2383 if (lh.t == MPQ and rh.t == MPZ) {
2384 mpq_t tmp;
2385 mpq_init(tmp);
2386 mpq_set_z(tmp, rh.v._bigint);
2387 mpq_div(tmp, lh.v._bigrat, tmp);
2388 if (mpz_cmp_ui(mpq_denref(tmp),1) != 0) {
2389 mpq_set(lh.v._bigrat, tmp);
2390 lh.hash = _mpq_pythonhash(lh.v._bigrat);
2391 mpq_clear(tmp);
2392 return lh;
2393 }
2394 mpq_clear(lh.v._bigrat);
2395 lh.t = MPZ;
2396 mpz_init(lh.v._bigint);
2397 mpz_set(lh.v._bigint, mpq_numref(tmp));
2398 lh.hash = _mpz_pythonhash(lh.v._bigint);
2399 mpq_clear(tmp);
2400 return lh;
2401 }
2402
2403 numeric a, b;
2404 coerce(a, b, lh, rh);
2405 lh = a / b;
2406 return lh;
2407 }
2408 switch (lh.t) {
2409 case LONG: {
2410 auto ld = std::div(lh.v._long, rh.v._long);
2411 if (ld.rem == 0) {
2412 lh.v._long = ld.quot;
2413 lh.hash = (ld.quot == -1) ? -2 : ld.quot;
2414 return lh;
2415 }
2416
2417 mpq_t obigrat;
2418 lh.t = MPQ;
2419 unsigned long l = std::labs(lh.v._long);
2420 unsigned long r = std::labs(rh.v._long);
2421 int sign = sgn(lh.v._long) * sgn(rh.v._long);
2422 mpq_init(obigrat);
2423 mpq_init(lh.v._bigrat);
2424 mpq_set_ui(lh.v._bigrat, l, 1);
2425 mpq_set_ui(obigrat, r, 1);
2426 mpq_div(lh.v._bigrat, lh.v._bigrat, obigrat);
2427 if (sign == -1)
2428 mpq_neg(lh.v._bigrat, lh.v._bigrat);
2429 mpq_clear(obigrat);
2430 lh.hash = _mpq_pythonhash(lh.v._bigrat);
2431 return lh;
2432 }
2433 case MPZ: {
2434 if (mpz_divisible_p(lh.v._bigint, rh.v._bigint)) {
2435 mpz_divexact(lh.v._bigint,
2436 lh.v._bigint,
2437 rh.v._bigint);
2438 lh.hash = _mpz_pythonhash(lh.v._bigint);
2439 return lh;
2440 }
2441 mpq_t bigrat, obigrat;
2442 mpq_init(bigrat);
2443 mpq_init(obigrat);
2444 mpq_set_z(bigrat, lh.v._bigint);
2445 mpq_set_z(obigrat, rh.v._bigint);
2446 mpz_clear(lh.v._bigint);
2447 lh.t = MPQ;
2448 mpq_init(lh.v._bigrat);
2449 mpq_div(lh.v._bigrat, bigrat, obigrat);
2450 lh.hash = _mpq_pythonhash(lh.v._bigrat);
2451 mpq_clear(bigrat);
2452 mpq_clear(obigrat);
2453 return lh;
2454 }
2455 case MPQ:
2456 mpq_div(lh.v._bigrat, lh.v._bigrat, rh.v._bigrat);
2457 lh.hash = _mpq_pythonhash(lh.v._bigrat);
2458 return lh;
2459 case PYOBJECT: {
2460 PyObject *p = lh.v._pyobject;
2461 #if PY_MAJOR_VERSION < 3
2462 {
2463 if (PyInt_Check(p)) {
2464 if (PyInt_Check(rh.v._pyobject)) {
2465 // This branch happens at startup.
2466 lh.v._pyobject = PyNumber_TrueDivide(Integer(PyInt_AsLong(p)),
2467 Integer(PyInt_AsLong(rh.v._pyobject)));
2468 // I don't 100% understand why I have to incref this,
2469 // but if I don't, Sage crashes on exit.
2470 if (lh.v._pyobject == nullptr) {
2471 lh.v._pyobject = p;
2472 py_error("numeric operator/=");
2473 }
2474 lh.hash = PyObject_Hash(lh.v._pyobject);
2475 Py_DECREF(p);
2476 return lh;
2477 }
2478 if (PyLong_Check(rh.v._pyobject)) {
2479 PyObject *d = py_funcs.py_integer_from_python_obj(rh.v._pyobject);
2480 lh.v._pyobject = PyNumber_TrueDivide(p, d);
2481 if (lh.v._pyobject == nullptr) {
2482 lh.v._pyobject = p;
2483 py_error("numeric operator/=");
2484 }
2485 lh.hash = PyObject_Hash(lh.v._pyobject);
2486 Py_DECREF(d);
2487 Py_DECREF(p);
2488 return lh;
2489 }
2490 } else if (PyLong_Check(p)) {
2491 PyObject *n = py_funcs.py_integer_from_python_obj(p);
2492 lh.v._pyobject = PyNumber_TrueDivide(n, rh.v._pyobject);
2493 if (lh.v._pyobject == nullptr) {
2494 lh.v._pyobject = p;
2495 py_error("numeric operator/=");
2496 }
2497 lh.hash = PyObject_Hash(lh.v._pyobject);
2498 Py_DECREF(n);
2499 Py_DECREF(p);
2500 return lh;
2501 }
2502 }
2503 #else
2504 {
2505 if (PyLong_Check(p)) {
2506 PyObject *n = py_funcs.py_integer_from_python_obj(p);
2507 lh.v._pyobject = PyNumber_TrueDivide(n, rh.v._pyobject);
2508 if (lh.v._pyobject == nullptr) {
2509 lh.v._pyobject = p;
2510 py_error("numeric operator/=");
2511 }
2512 lh.hash = (long)PyObject_Hash(lh.v._pyobject);
2513 Py_DECREF(n);
2514 Py_DECREF(p);
2515 return lh;
2516 }
2517 }
2518 #endif
2519 lh.v._pyobject = PyNumber_TrueDivide(p, rh.v._pyobject);
2520 if (lh.v._pyobject == nullptr) {
2521 lh.v._pyobject = p;
2522 py_error("numeric operator/=");
2523 }
2524 lh.hash = PyObject_Hash(lh.v._pyobject);
2525 Py_DECREF(p);
2526 return lh;
2527 }
2528 default:
2529 stub("invalid type: operator/=() type not handled");
2530 }
2531 }
2532
2533
2534 /** Inverse of a number. */
inverse() const2535 const numeric numeric::inverse() const {
2536 if (is_zero())
2537 throw std::overflow_error("numeric::inverse(): division by zero");
2538 return numeric(1) / *this;
2539 }
2540
2541 /** Return the step function of a numeric. The imaginary part of it is
2542 * ignored because the step function is generally considered real but
2543 * a numeric may develop a small imaginary part due to rounding errors.
2544 */
step() const2545 const numeric numeric::step() const {
2546 switch (t) {
2547 case LONG: return v._long > 0;
2548 case MPZ:
2549 case MPQ:
2550 if (is_positive())
2551 return 1;
2552 else
2553 return 0;
2554 case PYOBJECT:
2555 return py_funcs.py_step(v._pyobject);
2556 default:
2557 stub("invalid type: step() type not handled");
2558 }
2559 }
2560
2561 /** Return the complex half-plane (left or right) in which the number lies.
2562 * csgn(x)==0 for x==0, csgn(x)==1 for Re(x)>0 or Re(x)=0 and Im(x)>0,
2563 * csgn(x)==-1 for Re(x)<0 or Re(x)=0 and Im(x)<0.
2564 *
2565 * @see numeric::compare(const numeric &other) */
csgn() const2566 int numeric::csgn() const {
2567 verbose("csgn");
2568 switch (t) {
2569 case LONG:
2570 if (v._long == 0)
2571 return 0;
2572 return v._long < 0 ? -1 : 1;
2573 case MPZ:
2574 return mpz_sgn(v._bigint);
2575 case MPQ:
2576 return mpq_sgn(v._bigrat);
2577 case PYOBJECT: {
2578 int result;
2579 if (is_real()) {
2580 numeric z(ZERO);
2581 Py_INCREF(ZERO);
2582 return compare_same_type(z);
2583 }
2584 numeric re = real();
2585 numeric z(ZERO);
2586 Py_INCREF(ZERO);
2587 result = re.compare_same_type(z);
2588 if (result == 0) {
2589 numeric im = imag();
2590 result = im.compare_same_type(z);
2591 }
2592 return result;
2593 }
2594 default:
2595 stub("invalid type: csgn() type not handled");
2596 }
2597 }
2598
is_equal(const numeric & other) const2599 bool numeric::is_equal(const numeric &other) const {
2600 return *this == other;
2601 }
2602
2603 /** True if object is zero. */
is_zero() const2604 bool numeric::is_zero() const {
2605 verbose("is_zero");
2606 int a;
2607 switch (t) {
2608 case LONG:
2609 return v._long == 0;
2610 case MPZ:
2611 return mpz_cmp_si(v._bigint, 0) == 0;
2612 case MPQ:
2613 return mpq_cmp_si(v._bigrat, 0, 1) == 0;
2614 case PYOBJECT:
2615 a = PyObject_Not(v._pyobject);
2616 if (a == -1)
2617 py_error("is_zero");
2618 return a == 1;
2619 default:
2620 std::cerr << "type = " << t << "\n";
2621 stub("invalid type: is_zero() type not handled");
2622 }
2623 }
2624
is_inexact_one() const2625 bool numeric::is_inexact_one() const {
2626 verbose("is_one");
2627 switch (t) {
2628 case LONG:
2629 return v._long == 1;
2630 case MPZ:
2631 return mpz_cmp_si(v._bigint, 1) == 0;
2632 case MPQ:
2633 return mpq_cmp_si(v._bigrat, 1, 1) == 0;
2634 case PYOBJECT:
2635 return is_equal(*_num1_p);
2636 default:
2637 std::cerr << "type = " << t << "\n";
2638 stub("invalid type: is_one() type not handled");
2639 }
2640 }
2641
is_one() const2642 bool numeric::is_one() const {
2643 verbose("is_one");
2644 switch (t) {
2645 case LONG:
2646 return v._long == 1;
2647 case MPZ:
2648 return mpz_cmp_si(v._bigint, 1) == 0;
2649 case MPQ:
2650 return mpq_cmp_si(v._bigrat, 1, 1) == 0;
2651 case PYOBJECT:
2652 return is_exact() and is_equal(*_num1_p);
2653 default:
2654 std::cerr << "type = " << t << "\n";
2655 stub("invalid type: is_one() type not handled");
2656 }
2657 }
2658
is_minus_one() const2659 bool numeric::is_minus_one() const {
2660 verbose("is_one");
2661 switch (t) {
2662 case LONG:
2663 return v._long == -1;
2664 case MPZ:
2665 return mpz_cmp_si(v._bigint, -1) == 0;
2666 case MPQ:
2667 return mpq_cmp_si(v._bigrat, -1, 1) == 0;
2668 case PYOBJECT:
2669 return is_exact() and is_equal(*_num_1_p);
2670 default:
2671 std::cerr << "type = " << t << "\n";
2672 stub("invalid type: is_minus_one() type not handled");
2673 }
2674 }
2675
2676 /** True if object is not complex and greater than zero. */
is_positive() const2677 bool numeric::is_positive() const {
2678 verbose("is_positive");
2679 switch (t) {
2680 case LONG:
2681 return v._long > 0;
2682 case MPZ:
2683 return mpz_cmp_si(v._bigint, 0) > 0;
2684 case MPQ:
2685 return mpq_cmp_si(v._bigrat, 0, 1) > 0;
2686 case PYOBJECT:
2687 if (is_real()) {
2688 int result;
2689 result = PyObject_RichCompareBool(v._pyobject,
2690 ZERO, Py_GT);
2691 if (result == 1)
2692 return true;
2693 if (result == -1)
2694 PyErr_Clear();
2695 }
2696 return false;
2697 default:
2698 stub("invalid type: is_positive() type not handled");
2699 }
2700 }
2701
2702 /** True if object is not complex and less than zero. */
is_negative() const2703 bool numeric::is_negative() const {
2704 verbose("is_negative");
2705 switch (t) {
2706 case LONG:
2707 return v._long < 0;
2708 case MPZ:
2709 return mpz_cmp_si(v._bigint, 0) < 0;
2710 case MPQ:
2711 return mpq_cmp_si(v._bigrat, 0, 1) < 0;
2712 case PYOBJECT:
2713 if (is_real()) {
2714 int result;
2715 result = PyObject_RichCompareBool(v._pyobject,
2716 ZERO, Py_LT);
2717 if (result == 1)
2718 return true;
2719 if (result == -1)
2720 PyErr_Clear();
2721 }
2722 return false;
2723 default:
2724 stub("invalid type: is_negative() type not handled");
2725 }
2726 }
2727
2728 /** True if object is a non-complex integer. */
is_integer() const2729 bool numeric::is_integer() const {
2730 verbose2("is_integer", *this);
2731
2732 bool ret;
2733 switch (t) {
2734 case LONG:
2735 case MPZ:
2736 return true;
2737 case MPQ: {
2738 mpq_t bigrat;
2739 mpq_init(bigrat);
2740 mpq_set(bigrat, v._bigrat);
2741 mpq_canonicalize(bigrat);
2742 ret = mpz_cmp_ui(mpq_denref(bigrat), 1) == 0;
2743 mpq_clear(bigrat);
2744 return ret;
2745 }
2746 case PYOBJECT:
2747 return py_funcs.py_is_integer(v._pyobject) != 0;
2748 default:
2749 stub("invalid type: is_integer() type not handled");
2750 }
2751 }
2752
2753 /** True if object is an exact integer greater than zero. */
is_pos_integer() const2754 bool numeric::is_pos_integer() const {
2755 verbose("is_pos_integer");
2756 switch (t) {
2757 case LONG:
2758 return v._long > 0;
2759 case MPZ:
2760 return is_positive();
2761 case MPQ:
2762 return (is_integer() && is_positive());
2763 case PYOBJECT:
2764 return (is_integer() && is_positive());
2765 default:
2766 stub("invalid type: is_pos_integer() type not handled");
2767 }
2768 }
2769
2770 /** True if object is an exact integer greater or equal zero. */
is_nonneg_integer() const2771 bool numeric::is_nonneg_integer() const {
2772 verbose("is_nonneg_integer");
2773 switch (t) {
2774 case LONG:
2775 return v._long >= 0;
2776 case MPZ:
2777 return is_positive() or is_zero();
2778 case MPQ:
2779 return (is_integer() and (is_positive() or is_zero()));
2780 case PYOBJECT:
2781 if (is_integer()) {
2782 int result;
2783 result = PyObject_RichCompareBool(v._pyobject,
2784 ZERO, Py_GE);
2785 if (result == 1)
2786 return true;
2787 if (result == -1)
2788 PyErr_Clear();
2789 }
2790 return false;
2791 default:
2792 stub("invalid type: is_nonneg_integer() type not handled");
2793 }
2794 }
2795
2796 /** True if object is an exact even integer. */
is_even() const2797 bool numeric::is_even() const {
2798 verbose("is_even");
2799
2800 if (!is_integer())
2801 return false;
2802
2803 switch (t) {
2804 case LONG:
2805 return v._long % 2 == 0;
2806 case MPZ:
2807 return mpz_tstbit(v._bigint, 0) == 0;
2808 case MPQ:
2809 return is_integer()
2810 and mpz_tstbit(mpq_numref(v._bigrat), 0) == 0;
2811 case PYOBJECT:
2812 return py_funcs.py_is_even(v._pyobject) != 0;
2813 default:
2814 stub("invalid type: is_even() type not handled");
2815 }
2816 }
2817
2818 /** True if object is an exact odd integer. */
is_odd() const2819 bool numeric::is_odd() const {
2820 switch (t) {
2821 case LONG:
2822 return v._long & 1;
2823 case MPZ:
2824 return mpz_tstbit(v._bigint, 0) == 1;
2825 case MPQ:
2826 return is_integer()
2827 and mpz_tstbit(mpq_numref(v._bigrat), 0) == 1;
2828 case PYOBJECT:
2829 return !is_even();
2830 default:
2831 stub("invalid type: is_odd() type not handled");
2832 }
2833 }
2834
2835 /** Probabilistic primality test.
2836 *
2837 * @return true if object is exact integer and prime. */
is_prime() const2838 bool numeric::is_prime() const {
2839 verbose("is_prime");
2840 switch (t) {
2841 case LONG: {
2842 mpz_t bigint;
2843 mpz_init_set_si(bigint, v._long);
2844 bool ret = mpz_probab_prime_p(bigint, 25) > 0;
2845 mpz_clear(bigint);
2846 return ret;
2847 }
2848 case MPZ:
2849 return mpz_probab_prime_p(v._bigint, 25) > 0;
2850 case MPQ:
2851 return is_integer()
2852 and mpz_probab_prime_p(mpq_numref(v._bigrat), 25) > 0;
2853 case PYOBJECT:
2854 return py_funcs.py_is_prime(v._pyobject) != 0;
2855 default:
2856 stub("invalid type: is_prime() type not handled");
2857 }
2858 }
2859
2860 /** True if object is an exact rational number, may even be complex
2861 * (denominator may be unity). */
is_rational() const2862 bool numeric::is_rational() const {
2863 verbose("is_rational");
2864 switch (t) {
2865 case LONG:
2866 case MPZ:
2867 case MPQ:
2868 return true;
2869 case PYOBJECT:
2870 return false;
2871 default:
2872 stub("invalid type -- is_rational() type not handled");
2873 }
2874 }
2875
2876 /** True if object is a real integer, rational or float (but not complex). */
is_real() const2877 bool numeric::is_real() const {
2878 verbose("is_real");
2879 switch (t) {
2880 case LONG:
2881 case MPZ:
2882 case MPQ:
2883 return true;
2884 case PYOBJECT:
2885 return py_funcs.py_is_real(v._pyobject) != 0;
2886 default:
2887 stub("invalid type -- is_real() type not handled");
2888 }
2889 }
2890
2891 /** True if object is an exact rational number, may even be complex
2892 * (denominator may be unity). */
is_exact() const2893 bool numeric::is_exact() const {
2894 verbose("is_exact");
2895 switch (t) {
2896 case LONG:
2897 case MPZ:
2898 case MPQ:
2899 return true;
2900 case PYOBJECT:
2901 return py_funcs.py_is_exact(v._pyobject) != 0;
2902 default:
2903 stub("invalid type -- is_exact() type not handled");
2904 }
2905 }
2906
2907 static std::map<long,std::pair<int,int>> small_powers;
2908
fill_small_powers()2909 static void fill_small_powers()
2910 {
2911 static int lim[] = {30, 18, 15, 12, 11, 10, 10, 9, 9};
2912 for (size_t b = sizeof(lim)/sizeof(int) + 1; b >= 2; --b) {
2913 long p = b*b;
2914 int c = 2;
2915 while (c <= lim[b-2]) {
2916 small_powers[p] = std::make_pair(int(b), c);
2917 p *= b;
2918 ++c;
2919 }
2920 }
2921 }
2922
is_small_power(std::pair<int,int> & p) const2923 bool numeric::is_small_power(std::pair<int,int>& p) const
2924 {
2925 int i;
2926 switch (t) {
2927 case LONG:
2928 if (v._long < 2)
2929 return false;
2930 i = v._long;
2931 break;
2932 case MPZ:
2933 if (not mpz_fits_sint_p(v._bigint))
2934 return false;
2935 i = mpz_get_si(v._bigint);
2936 if (i < 2)
2937 return false;
2938 break;
2939 case MPQ:
2940 case PYOBJECT:
2941 return false;
2942 default:
2943 stub("invalid type -- is_small_power() type not handled");
2944 }
2945 if (small_powers.empty())
2946 fill_small_powers();
2947 auto it = small_powers.find(i);
2948 if (it == small_powers.end())
2949 return false;
2950 p = it->second;
2951 return true;
2952 }
2953
operator ==(const numeric & right) const2954 bool numeric::operator==(const numeric &right) const {
2955 verbose3("operator==", *this, right);
2956
2957 if (this == &right)
2958 return true;
2959 if (t != right.t) {
2960 if (t == LONG and right.t == MPZ)
2961 return mpz_cmp_si(right.v._bigint, v._long) ==0;
2962 if (right.t == LONG and t == MPZ)
2963 return mpz_cmp_si(v._bigint, right.v._long) ==0;
2964 if (t == MPZ and right.t == MPQ) {
2965 if (mpz_cmp_ui(mpq_denref(right.v._bigrat), 1) != 0)
2966 return false;
2967 return mpz_cmp(v._bigint, mpq_numref(right.v._bigrat)) ==0;
2968 }
2969 if (t == MPQ and right.t == MPZ) {
2970 if (mpz_cmp_ui(mpq_denref(v._bigrat), 1) != 0)
2971 return false;
2972 return mpz_cmp(right.v._bigint, mpq_numref(v._bigrat)) ==0;
2973 }
2974 numeric a, b;
2975 coerce(a, b, *this, right);
2976 return a == b;
2977 }
2978 switch (t) {
2979 case LONG:
2980 return v._long == right.v._long;
2981 case MPZ:
2982 return mpz_cmp(v._bigint, right.v._bigint) == 0;
2983 case MPQ:
2984 return mpq_equal(v._bigrat, right.v._bigrat) != 0;
2985 case PYOBJECT:
2986 if (v._pyobject == right.v._pyobject)
2987 return true;
2988 return py_funcs.py_is_equal(v._pyobject,
2989 right.v._pyobject) != 0;
2990 default:
2991 stub("invalid type: operator== type not handled");
2992 }
2993 }
2994
operator !=(const numeric & right) const2995 bool numeric::operator!=(const numeric &right) const {
2996 verbose("operator!=");
2997 if (t != right.t) {
2998 if (t == LONG and right.t == MPZ)
2999 return mpz_cmp_si(right.v._bigint, v._long) !=0;
3000 if (right.t == LONG and t == MPZ)
3001 return mpz_cmp_si(v._bigint, right.v._long) !=0;
3002 if (t == MPZ and right.t == MPQ) {
3003 if (mpz_cmp_ui(mpq_denref(right.v._bigrat), 1) != 0)
3004 return true;
3005 return mpz_cmp(v._bigint, mpq_numref(right.v._bigrat)) !=0;
3006 }
3007 if (t == MPQ and right.t == MPZ) {
3008 if (mpz_cmp_ui(mpq_denref(v._bigrat), 1) != 0)
3009 return true;
3010 return mpz_cmp(right.v._bigint, mpq_numref(v._bigrat)) !=0;
3011 }
3012 numeric a, b;
3013 coerce(a, b, *this, right);
3014 return a != b;
3015 }
3016 switch (t) {
3017 case LONG:
3018 return v._long != right.v._long;
3019 case MPZ:
3020 return mpz_cmp(v._bigint, right.v._bigint) != 0;
3021 case MPQ:
3022 return mpq_equal(v._bigrat, right.v._bigrat) == 0;
3023 case PYOBJECT:
3024 return (py_funcs.py_is_equal(v._pyobject,
3025 right.v._pyobject) == 0);
3026 default:
3027 stub("invalid type: operator!= type not handled");
3028 }
3029 }
3030
3031 /** True if object is element of the domain of integers extended by I, i.e. is
3032 * of the form a+b*I, where a and b are integers. */
is_cinteger() const3033 bool numeric::is_cinteger() const {
3034 verbose("is_crational");
3035 switch (t) {
3036 case LONG:
3037 case MPZ:
3038 return true;
3039 case MPQ:
3040 return is_integer();
3041 case PYOBJECT:
3042 return real().is_integer()
3043 and imag().is_integer();
3044 default:
3045 stub("invalid type -- is_cinteger() type not handled");
3046 }
3047 }
3048
3049 /** True if object is an exact rational number, may even be complex
3050 * (denominator may be unity). */
is_crational() const3051 bool numeric::is_crational() const {
3052 verbose("is_crational");
3053 switch (t) {
3054 case LONG:
3055 case MPZ:
3056 case MPQ:
3057 return true;
3058 case PYOBJECT:
3059 return real().is_rational()
3060 and imag().is_rational();
3061 default:
3062 stub("invalid type -- is_crational() type not handled");
3063 }
3064 }
3065
3066 /** Numerical comparison: less.
3067 *
3068 * @exception invalid_argument (complex inequality) */
operator <(const numeric & right) const3069 bool numeric::operator<(const numeric &right) const {
3070 verbose("operator<");
3071 if (t == MPZ and right.t == LONG)
3072 return mpz_cmp_si(v._bigint, right.v._long) < 0;
3073 if (t == LONG and right.t == MPZ)
3074 return mpz_cmp_si(right.v._bigint, v._long) > 0;
3075 if (t != right.t) {
3076 numeric a, b;
3077 coerce(a, b, *this, right);
3078 return a < b;
3079 }
3080
3081 switch (t) {
3082 case LONG:
3083 return v._long < right.v._long;
3084 case MPZ:
3085 return mpz_cmp(v._bigint, right.v._bigint) < 0;
3086 case MPQ:
3087 return mpq_cmp(v._bigrat, right.v._bigrat) < 0;
3088 case PYOBJECT: {
3089 int result;
3090 result = PyObject_RichCompareBool(v._pyobject,
3091 right.v._pyobject, Py_LT);
3092 if (result == -1)
3093 py_error("richcmp failed");
3094
3095 return (result == 1);
3096 }
3097 default:
3098 stub("invalid type: operator< type not handled");
3099 }
3100 }
3101
3102 /** Numerical comparison: less or equal.
3103 *
3104 * @exception invalid_argument (complex inequality) */
operator <=(const numeric & right) const3105 bool numeric::operator<=(const numeric &right) const {
3106 verbose("operator<=");
3107 if (t == MPZ and right.t == LONG)
3108 return mpz_cmp_si(v._bigint, right.v._long) <= 0;
3109 if (t == LONG and right.t == MPZ)
3110 return mpz_cmp_si(right.v._bigint, v._long) >= 0;
3111 if (t != right.t) {
3112 numeric a, b;
3113 coerce(a, b, *this, right);
3114 return a <= b;
3115 }
3116 switch (t) {
3117 case LONG:
3118 return v._long <= right.v._long;
3119 case MPZ:
3120 return mpz_cmp(v._bigint, right.v._bigint) <= 0;
3121 case MPQ:
3122 return mpq_cmp(v._bigrat, right.v._bigrat) <= 0;
3123 case PYOBJECT:
3124 int result;
3125 result = PyObject_RichCompareBool(v._pyobject,
3126 right.v._pyobject, Py_LE);
3127 if (result == -1)
3128 py_error("richcmp failed");
3129
3130 return (result == 1);
3131 default:
3132 stub("invalid type: operator<= type not handled");
3133 }
3134 }
3135
3136 /** Numerical comparison: greater.
3137 *
3138 * @exception invalid_argument (complex inequality) */
operator >(const numeric & right) const3139 bool numeric::operator>(const numeric &right) const {
3140 verbose("operator>");
3141 if (t == MPZ and right.t == LONG)
3142 return mpz_cmp_si(v._bigint, right.v._long) > 0;
3143 if (t == LONG and right.t == MPZ)
3144 return mpz_cmp_si(right.v._bigint, v._long) < 0;
3145 if (t != right.t) {
3146 numeric a, b;
3147 coerce(a, b, *this, right);
3148 return a > b;
3149 }
3150 switch (t) {
3151 case LONG:
3152 return v._long > right.v._long;
3153 case MPZ:
3154 return mpz_cmp(v._bigint, right.v._bigint) > 0;
3155 case MPQ:
3156 return mpq_cmp(v._bigrat, right.v._bigrat) > 0;
3157 case PYOBJECT:
3158 int result;
3159 result = PyObject_RichCompareBool(v._pyobject,
3160 right.v._pyobject, Py_GT);
3161 if (result == -1)
3162 py_error("richcmp failed");
3163
3164 return (result == 1);
3165 default:
3166 stub("invalid type: operator> type not handled");
3167 }
3168 }
3169
3170 /** Numerical comparison: greater or equal.
3171 *
3172 * @exception invalid_argument (complex inequality) */
operator >=(const numeric & right) const3173 bool numeric::operator>=(const numeric &right) const {
3174 verbose("operator>=");
3175 if (t == MPZ and right.t == LONG)
3176 return mpz_cmp_si(v._bigint, right.v._long) >= 0;
3177 if (t == LONG and right.t == MPZ)
3178 return mpz_cmp_si(right.v._bigint, v._long) <= 0;
3179 if (t != right.t) {
3180 numeric a, b;
3181 coerce(a, b, *this, right);
3182 return a >= b;
3183 }
3184 switch (t) {
3185 case LONG:
3186 return v._long >= right.v._long;
3187 case MPZ:
3188 return mpz_cmp(v._bigint, right.v._bigint) >= 0;
3189 case MPQ:
3190 return mpq_cmp(v._bigrat, right.v._bigrat) >= 0;
3191 case PYOBJECT:
3192 int result;
3193 result = PyObject_RichCompareBool(v._pyobject,
3194 right.v._pyobject, Py_GE);
3195 if (result == -1)
3196 py_error("richcmp failed");
3197
3198 return (result == 1);
3199 default:
3200 stub("invalid type: operator!= type not handled");
3201 }
3202 }
3203
to_int() const3204 int numeric::to_int() const
3205 {
3206 switch (t) {
3207 case LONG:
3208 if (v._long < std::numeric_limits<int>::max()
3209 and v._long > std::numeric_limits<int>::min())
3210 return v._long;
3211 throw std::runtime_error("to_int");
3212 case MPZ:
3213 if (mpz_fits_sint_p(v._bigint))
3214 return mpz_get_si(v._bigint);
3215 throw conversion_error();
3216 case MPQ: {
3217 mpz_t bigint;
3218 mpz_init(bigint);
3219 mpz_fdiv_q(bigint, mpq_numref(v._bigrat),
3220 mpq_denref(v._bigrat));
3221 if (mpz_fits_sint_p(bigint)) {
3222 int n = mpz_get_si(bigint);
3223 mpz_clear(bigint);
3224 return n;
3225 }
3226 mpz_clear(bigint);
3227 throw conversion_error();
3228 }
3229 case PYOBJECT:
3230 return to_bigint().to_int();
3231 default:
3232 stub("invalid type: operator long int() type not handled");
3233 }
3234 }
3235 /** Converts numeric types to machine's long. You should check with
3236 * is_integer() if the number is really an integer before calling this method.
3237 * You may also consider checking the range first. */
to_long() const3238 long numeric::to_long() const {
3239 verbose("operator long int");
3240 switch (t) {
3241 case LONG:
3242 return v._long;
3243 case MPZ:
3244 if (mpz_fits_slong_p(v._bigint))
3245 return (long int)mpz_get_si(v._bigint);
3246 throw conversion_error();
3247 case MPQ: {
3248 mpz_t bigint;
3249 mpz_init(bigint);
3250 mpz_fdiv_q(bigint, mpq_numref(v._bigrat),
3251 mpq_denref(v._bigrat));
3252 if (mpz_fits_slong_p(v._bigint)) {
3253 long n = mpz_get_si(bigint);
3254 mpz_clear(bigint);
3255 return n;
3256 }
3257 mpz_clear(bigint);
3258 throw conversion_error();
3259 }
3260 case PYOBJECT:
3261 return to_bigint().to_long();
3262 default:
3263 stub("invalid type: operator long int() type not handled");
3264 }
3265 }
3266
3267 // Use this only if o was tested integer.
to_bigint() const3268 const numeric numeric::to_bigint() const {
3269 switch (t) {
3270 case LONG: {
3271 numeric nu;
3272 mpz_init_set_si(nu.v._bigint, v._long);
3273 nu.t = MPZ;
3274 nu.hash = _mpz_pythonhash(nu.v._bigint);
3275 return nu;
3276 }
3277 case MPZ: return *this;
3278 case MPQ:
3279 if (not denom().is_one())
3280 throw std::runtime_error("not integer in numeric::to_mpz_num()");
3281 return numer();
3282 case PYOBJECT: {
3283 PyObject *Integer = Integer_pyclass();
3284 PyObject *ans = PyObject_CallFunctionObjArgs(Integer,
3285 v._pyobject, NULL);
3286 return ans;
3287 }
3288 default:
3289 stub("invalid type: operator long int() type not handled");
3290 }
3291 }
3292
as_mpz() const3293 const mpz_t& numeric::as_mpz() const
3294 {
3295 if (t != MPZ)
3296 throw std::runtime_error("mpz_t requested from non-mpz numeric");
3297 return v._bigint;
3298 }
3299
as_mpq() const3300 const mpq_t& numeric::as_mpq() const
3301 {
3302 if (t != MPQ)
3303 throw std::runtime_error("mpq_t requested from non-mpq numeric");
3304 return v._bigrat;
3305 }
3306
canonicalize()3307 void numeric::canonicalize()
3308 {
3309 if (t == MPQ) {
3310 mpq_canonicalize(v._bigrat);
3311 if (mpz_cmp_ui(mpq_denref(v._bigrat), 1) == 0) {
3312 mpz_t tmp;
3313 mpz_init_set(tmp, mpq_numref(v._bigrat));
3314 mpq_clear(v._bigrat);
3315 set_from(t, v, hash, tmp);
3316 mpz_clear(tmp);
3317 }
3318 }
3319 }
3320
3321 #ifdef PYNAC_HAVE_LIBGIAC
to_giacgen(giac::context * cptr) const3322 giac::gen* numeric::to_giacgen(giac::context* cptr) const
3323 {
3324 if (t == LONG)
3325 return new giac::gen(v._long);
3326 if (t == MPZ) {
3327 mpz_t bigint;
3328 mpz_init_set(bigint, v._bigint);
3329 auto ret = new giac::gen(bigint);
3330 mpz_clear(bigint);
3331 return ret;
3332 }
3333 if (t == MPQ) {
3334 mpz_t bigint;
3335 mpz_init_set(bigint, mpq_numref(v._bigrat));
3336 giac::gen gn(bigint);
3337 mpz_set(bigint, mpq_denref(v._bigrat));
3338 giac::gen gd(bigint);
3339 giac::Tfraction<giac::gen> frac(gn, gd);
3340 mpz_clear(bigint);
3341 return new giac::gen(frac);
3342 }
3343 else
3344 return nullptr;
3345 }
3346 #endif
3347
to_canonical() const3348 CanonicalForm numeric::to_canonical() const
3349 {
3350 if (t == LONG)
3351 return CanonicalForm(v._long);
3352 if (t == MPZ) {
3353 if (mpz_fits_sint_p(v._bigint))
3354 return CanonicalForm(to_int());
3355 mpz_t bigint;
3356 mpz_init_set(bigint, v._bigint);
3357 return make_cf(mpz_ptr(const_cast<__mpz_struct*>(bigint)));
3358 }
3359 if (t == MPQ) {
3360 mpz_t num;
3361 mpz_init_set(num, mpz_ptr(mpq_numref(v._bigrat)));
3362 if (is_integer())
3363 return make_cf(mpz_ptr(const_cast<__mpz_struct*>(num)));
3364 mpz_t den;
3365 mpz_init_set(den, mpz_ptr(mpq_denref(v._bigrat)));
3366 return make_cf(num, den, false);
3367 }
3368
3369 throw (std::runtime_error("can't happen in numeric::to_canonical"));
3370 }
3371
3372 /* Return the underlying Python object corresponding to this
3373 numeric. If this numeric isn't implemented using a Python
3374 object, the corresponding Python object is constructed on
3375 the fly.
3376
3377 Returns a NEW REFERENCE.
3378 */
to_pyobject() const3379 PyObject* numeric::to_pyobject() const {
3380 // Returns a New Reference
3381 PyObject* o;
3382 switch (t) {
3383 case LONG: {
3384 mpz_t bigint;
3385 mpz_init_set_si(bigint, v._long);
3386 o = py_funcs.py_integer_from_mpz(bigint);
3387 mpz_clear(bigint);
3388 return o;
3389 }
3390 case MPZ: {
3391 mpz_t bigint;
3392 mpz_init_set(bigint, v._bigint);
3393 o = py_funcs.py_integer_from_mpz(bigint);
3394 mpz_clear(bigint);
3395 return o;
3396 }
3397 case MPQ: {
3398 mpq_t bigrat;
3399 mpq_init(bigrat);
3400 mpq_set(bigrat, v._bigrat);
3401 mpq_canonicalize(bigrat);
3402 o = py_funcs.py_rational_from_mpq(bigrat);
3403 mpq_clear(bigrat);
3404 return o;
3405 }
3406 case PYOBJECT:
3407 Py_INCREF(v._pyobject);
3408 return v._pyobject;
3409
3410 default:
3411 std::cout << t << std::endl;
3412 stub("numeric::to_pyobject -- not able to do conversion to pyobject; everything else will be nonsense");
3413 return nullptr;
3414 }
3415 }
3416
3417
3418 /** Converts numeric types to machine's double. You should check with is_real()
3419 * if the number is really not complex before calling this method. */
to_double() const3420 double numeric::to_double() const {
3421 GINAC_ASSERT(this->is_real());
3422 verbose("operator double");
3423 double d;
3424 switch (t) {
3425 case LONG:
3426 return v._long;
3427 case MPZ:
3428 return mpz_get_d(v._bigint);
3429 case MPQ:
3430 return mpq_get_d(v._bigrat);
3431 case PYOBJECT:
3432 d = PyFloat_AsDouble(v._pyobject);
3433 if (d == -1 && (PyErr_Occurred() != nullptr))
3434 py_error("Error converting to a double.");
3435 return d;
3436 default:
3437 std::cerr << "type = " << t << std::endl;
3438 stub("invalid type: operator double() type not handled");
3439 }
3440 }
3441
3442 /** Evaluation of numbers doesn't do anything at all. */
eval(int) const3443 ex numeric::eval(int /*level*/) const {
3444 // Warning: if this is ever gonna do something, the ex constructors from all kinds
3445 // of numbers should be checking for status_flags::evaluated.
3446 return this->hold();
3447 }
3448
3449 /** Cast numeric into a floating-point object. For example exact numeric(1) is
3450 * returned as a 1.0000000000000000000000 and so on according to how Digits is
3451 * currently set. In case the object already was a floating point number the
3452 * precision is trimmed to match the currently set default.
3453 *
3454 * @param level ignored, only needed for overriding basic::evalf.
3455 * @return an ex-handle to a numeric. */
evalf(int,PyObject * parent) const3456 ex numeric::evalf(int /*level*/, PyObject* parent) const {
3457 PyObject *ans, *a = to_pyobject();
3458 if (parent == nullptr)
3459 parent = RR_get();
3460 if (not PyDict_CheckExact(parent)) {
3461 PyObject* dict = PyDict_New();
3462 if (dict == nullptr)
3463 throw(std::runtime_error("PyDict_New returned NULL"));
3464 int r = PyDict_SetItemString(dict, "parent", parent);
3465 if (r<0)
3466 throw(std::runtime_error("PyDict_SetItemString failed"));
3467 ans = py_funcs.py_float(a, dict);
3468 Py_DECREF(dict);
3469 }
3470 else
3471 ans = py_funcs.py_float(a, parent);
3472 Py_DECREF(a);
3473 if (ans == nullptr)
3474 throw (std::runtime_error("numeric::evalf(): error calling py_float()"));
3475
3476 return ans;
3477 }
3478
try_py_method(const std::string & s) const3479 const numeric numeric::try_py_method(const std::string& s) const
3480 {
3481 PyObject *obj = to_pyobject();
3482 PyObject* ret = PyObject_CallMethod(obj,
3483 const_cast<char*>(s.c_str()), NULL);
3484 Py_DECREF(obj);
3485 if (ret == nullptr) {
3486 PyErr_Clear();
3487 throw std::logic_error("");
3488 }
3489
3490 return numeric(ret);
3491 }
3492
try_py_method(const std::string & s,const numeric & num2) const3493 const numeric numeric::try_py_method(const std::string& s,
3494 const numeric& num2) const
3495 {
3496 PyObject *obj1 = to_pyobject();
3497 PyObject *obj2 = num2.to_pyobject();
3498 PyObject *mstr = PyString_FromString(const_cast<char*>(s.c_str()));
3499 PyObject* ret = PyObject_CallMethodObjArgs(obj1, mstr, obj2, NULL);
3500 Py_DECREF(obj1);
3501 Py_DECREF(obj2);
3502 Py_DECREF(mstr);
3503 if (ret == nullptr) {
3504 PyErr_Clear();
3505 throw std::logic_error("");
3506 }
3507 return numeric(ret);
3508 }
3509
to_dict_parent(PyObject * obj) const3510 const numeric numeric::to_dict_parent(PyObject* obj) const
3511 {
3512 PyObject *ret = nullptr;
3513 PyObject *the_parent = nullptr;
3514 PyObject *the_arg = to_pyobject();
3515 if (obj != nullptr and PyDict_Check(obj)) {
3516 PyObject* pkey = PyString_FromString(const_cast<char*>("parent"));
3517 the_parent = PyDict_GetItem(obj, pkey);
3518 Py_DECREF(pkey);
3519 if (the_parent != nullptr
3520 and PyCallable_Check(the_parent)) {
3521 ret = PyObject_CallFunctionObjArgs(the_parent, the_arg, nullptr);
3522 Py_DECREF(the_arg);
3523 if (ret == nullptr) {
3524 PyErr_Clear();
3525 throw std::logic_error("");
3526 }
3527 return numeric(ret);
3528 }
3529 }
3530 ret = PyObject_CallFunctionObjArgs(RR_get(), the_arg, NULL);
3531 if (ret == nullptr) {
3532 PyErr_Clear();
3533 ret = PyObject_CallFunctionObjArgs(CC_get(), the_arg, NULL);
3534 Py_DECREF(the_arg);
3535 if (ret == nullptr) {
3536 PyErr_Clear();
3537 throw std::logic_error("");
3538 }
3539 }
3540 else
3541 Py_DECREF(the_arg);
3542 return numeric(ret);
3543 }
3544
subs(const exmap & m,unsigned) const3545 ex numeric::subs(const exmap & m, unsigned) const
3546 {
3547 const numeric& im = imag();
3548 if (im.is_zero()) {
3549 if (is_zero() or is_one() or is_minus_one())
3550 return *this;
3551 for (const auto & pair : m)
3552 if (is_exactly_a<numeric>(pair.first)
3553 and is_equal(ex_to<numeric>(pair.first)))
3554 return pair.second;
3555 return *this;
3556 }
3557
3558 const numeric& re = real();
3559 ex ret1 = re;
3560 ex ret2 = im;
3561 ex ret3 = I;
3562 bool changed1=false, changed2=false, changed3=false;
3563 for (const auto & pair : m) {
3564 if (not is_exactly_a<numeric>(pair.first))
3565 continue;
3566 const numeric& p = ex_to<numeric>(pair.first);
3567 const numeric& pim = p.imag();
3568 const numeric& pre = p.real();
3569 if (pim.is_zero()) {
3570 if (p.is_zero() or p.is_one() or p.is_minus_one())
3571 continue;
3572 if (re.is_equal(pre)) {
3573 ret1 = pair.second;
3574 changed1 = true;
3575 }
3576 if (im.is_equal(pre)) {
3577 ret2 = pair.second;
3578 changed2 = true;
3579 }
3580 continue;
3581 }
3582 if (pim.is_one() and pre.is_zero()) {
3583 ret3 = pair.second;
3584 changed3 = true;
3585 continue;
3586 }
3587 if (re.is_equal(pre) and im.is_equal(pim))
3588 return pair.second;
3589 }
3590 if (changed1 or changed2 or changed3)
3591 return ret1 + ret2*ret3;
3592 return *this;
3593 }
3594
3595 ///////////////////////////////////////////////////////////////////////
3596
3597 /** Real part of a number. */
real() const3598 const numeric numeric::real() const {
3599 switch (t) {
3600 case LONG:
3601 case MPZ:
3602 case MPQ:
3603 return *this;
3604 case PYOBJECT:
3605 {
3606 if (PyFloat_Check(v._pyobject))
3607 return *this;
3608 if (PyComplex_Check(v._pyobject))
3609 return PyComplex_RealAsDouble(v._pyobject);
3610 try {
3611 return try_py_method("real");
3612 }
3613 catch (std::logic_error) {}
3614 try {
3615 return try_py_method("real_part");
3616 }
3617 catch (std::logic_error) {}
3618 return *this;
3619 }
3620 default:
3621 stub("invalid type");
3622 }
3623 }
3624
3625 /** Imaginary part of a number. */
imag() const3626 const numeric numeric::imag() const {
3627 switch (t) {
3628 case LONG:
3629 case MPZ:
3630 case MPQ:
3631 return *_num0_p;
3632 case PYOBJECT:
3633 {
3634 if (PyFloat_Check(v._pyobject))
3635 return *_num0_p;
3636 if (PyComplex_Check(v._pyobject))
3637 return PyComplex_ImagAsDouble(v._pyobject);
3638 try {
3639 return try_py_method("imag");
3640 }
3641 catch (std::logic_error) {}
3642 try {
3643 return try_py_method("imag_part");
3644 }
3645 catch (std::logic_error) {}
3646 return *_num0_p;
3647 }
3648 default:
3649 stub("invalid type");
3650 }
3651 }
3652
3653 /** Numerator. Computes the numerator of rational numbers, rationalized
3654 * numerator of complex if real and imaginary part are both rational numbers
3655 * (i.e numer(4/3+5/6*I) == 8+5*I), the number carrying the sign in all other
3656 * cases. */
numer() const3657 const numeric numeric::numer() const {
3658
3659 switch (t) {
3660 case LONG:
3661 case MPZ:
3662 return *this;
3663 case MPQ: {
3664 mpz_t bigint;
3665 mpz_init_set(bigint, mpq_numref(v._bigrat));
3666 return bigint;
3667 }
3668 case PYOBJECT: {
3669 PyObject *a;
3670 a = py_funcs.py_numer(v._pyobject);
3671 if (a == nullptr) py_error("numer");
3672 return a;
3673 }
3674 default:
3675 stub("invalid type -- numer() type not handled");
3676 }
3677 }
3678
3679 /** Denominator. Computes the denominator of rational numbers, common integer
3680 * denominator of complex if real and imaginary part are both rational numbers
3681 * (i.e denom(4/3+5/6*I) == 6), one in all other cases. */
denom() const3682 const numeric numeric::denom() const {
3683
3684 switch (t) {
3685 case LONG:
3686 case MPZ:
3687 return 1;
3688 case MPQ: {
3689 mpz_t bigint;
3690 mpz_init_set(bigint, mpq_denref(v._bigrat));
3691 return bigint;
3692 }
3693 case PYOBJECT: {
3694 PyObject *a;
3695 a = py_funcs.py_denom(v._pyobject);
3696 if (a == nullptr) py_error("denom");
3697 return a;
3698 }
3699 default:
3700 stub("invalid type -- denom() type not handled");
3701 }
3702 }
3703
floor() const3704 const numeric numeric::floor() const {
3705 numeric d = denom();
3706 if (d.is_one())
3707 return *this;
3708 return numer().iquo(d);
3709 }
3710
frac() const3711 const numeric numeric::frac() const {
3712 numeric d = denom();
3713 if (d.is_one())
3714 return 0;
3715 return *this - numer().iquo(d);
3716 }
3717
fibonacci() const3718 const numeric numeric::fibonacci() const {
3719 PY_RETURN(py_funcs.py_fibonacci);
3720 }
3721
exp(PyObject * parent) const3722 const numeric numeric::exp(PyObject* parent) const {
3723 static numeric tentt20 = ex_to<numeric>(_num10_p->power(*_num20_p));
3724 // avoid Flint aborts
3725 if (real() > tentt20) {
3726 if (imag().is_zero())
3727 return py_funcs.py_eval_infinity();
3728 else
3729 return py_funcs.py_eval_unsigned_infinity();
3730 }
3731 if (real() < -tentt20)
3732 return ex_to<numeric>(_num0_p->evalf(0, parent));
3733
3734 return arbfunc_0arg("exp", parent);
3735 }
3736
log(PyObject * parent) const3737 const numeric numeric::log(PyObject* parent) const {
3738 return arbfunc_0arg("log", parent);
3739 }
3740
3741 // General log
log(const numeric & b,PyObject * parent) const3742 const numeric numeric::log(const numeric &b, PyObject* parent) const {
3743 if (b.is_one()) {
3744 if (is_one())
3745 throw (std::runtime_error("log(1,1) encountered"));
3746 else
3747 return py_funcs.py_eval_unsigned_infinity();
3748 }
3749 if (b.is_zero())
3750 return *_num0_p;
3751
3752 if ((t != LONG and t != MPZ and t != MPQ)
3753 or (b.t != LONG and b.t != MPZ and b.t != MPQ))
3754 return log(parent)/b.log(parent);
3755
3756 bool israt;
3757 numeric ret = ratlog(b, israt);
3758 if (not israt)
3759 return log(parent)/b.log(parent);
3760
3761 return ret;
3762 }
3763
3764 // General log
3765 // Handle special cases here that return MPZ/MPQ
ratlog(const numeric & b,bool & israt) const3766 const numeric numeric::ratlog(const numeric &b, bool& israt) const {
3767 israt = true;
3768 if (b.is_one()) {
3769 if (is_one())
3770 throw (std::runtime_error("log(1,1) encountered"));
3771 else
3772 return py_funcs.py_eval_unsigned_infinity();
3773 }
3774 if (b.is_zero())
3775 return *_num0_p;
3776
3777 if ((t != LONG and t != MPZ and t != MPQ)
3778 or (b.t != LONG and b.t != MPZ and b.t != MPQ)) {
3779 israt = false;
3780 return *_num0_p;
3781 }
3782
3783 if (t == LONG and b.t == LONG) {
3784 if (b.v._long <= 0) {
3785 israt = false;
3786 return *_num0_p;
3787 }
3788 int c = 0;
3789 std::ldiv_t ld;
3790 ld.quot = v._long;
3791 ld.rem = 0;
3792 do {
3793 ld = std::div(ld.quot, b.v._long);
3794 ++c;
3795 }
3796 while (ld.quot != 1 and ld.rem == 0);
3797 if (ld.quot != 1 or ld.rem != 0) {
3798 israt = false;
3799 return *_num0_p;
3800 }
3801 return c;
3802 }
3803 if (b.t == LONG)
3804 return ratlog(b.to_bigint(), israt);
3805 if (t == LONG)
3806 return to_bigint().ratlog(b, israt);
3807 if (t == MPZ and b.t == MPZ) {
3808 if (b > *_num0_p) {
3809 mpz_t ret;
3810 mpz_init(ret);
3811 mp_bitcnt_t r = mpz_remove(ret, v._bigint, b.v._bigint);
3812 if (mpz_cmp_ui(ret, 1) == 0) {
3813 mpz_clear(ret);
3814 return long(r);
3815 }
3816 mpz_clear(ret);
3817 }
3818 israt = false;
3819 return *_num0_p;
3820 }
3821
3822 if (t == MPZ) {
3823 if (b.numer().is_one())
3824 return -ratlog(b.denom(), israt);
3825
3826 israt = false;
3827 return *_num0_p;
3828 }
3829 if (b.t == MPZ) {
3830 if (numer().is_one())
3831 return -denom().ratlog(b, israt);
3832
3833 israt = false;
3834 return *_num0_p;
3835 }
3836
3837 // from here both MPQ
3838 numeric d = GiNaC::log(denom(), b.denom(), nullptr);
3839 if (numer().is_one() and b.numer().is_one())
3840 return d;
3841 numeric n = GiNaC::log(numer(), b.numer(), nullptr);
3842 if (n == d)
3843 return n;
3844 return log(nullptr)/b.log(nullptr);
3845 }
3846
sin() const3847 const numeric numeric::sin() const {
3848 PY_RETURN(py_funcs.py_sin);
3849 }
3850
cos() const3851 const numeric numeric::cos() const {
3852 PY_RETURN(py_funcs.py_cos);
3853 }
3854
tan() const3855 const numeric numeric::tan() const {
3856 PY_RETURN(py_funcs.py_tan);
3857 }
3858
asin(PyObject * parent) const3859 const numeric numeric::asin(PyObject* parent) const {
3860 return arbfunc_0arg("arcsin", parent);
3861 }
3862
acos(PyObject * parent) const3863 const numeric numeric::acos(PyObject* parent) const {
3864 return arbfunc_0arg("arccos", parent);
3865 }
3866
atan(PyObject * parent) const3867 const numeric numeric::atan(PyObject* parent) const {
3868 return arbfunc_0arg("arctan", parent);
3869 }
3870
atan(const numeric & y,PyObject * parent) const3871 const numeric numeric::atan(const numeric& y, PyObject* parent) const {
3872 PyObject *cparent = common_parent(*this, y);
3873 bool newdict = false;
3874 if (parent == nullptr) {
3875 parent = PyDict_New();
3876 PyDict_SetItemString(parent, "parent", cparent);
3877 newdict = true;
3878 }
3879 numeric rnum;
3880 if (not imag().is_zero()
3881 or not y.imag().is_zero()) {
3882 ex r = this->add(y.mul(I))
3883 / (mul(*this) + y*y).power(*_num1_2_p);
3884 rnum = -ex_to<numeric>(r.evalf(0, parent)).log(parent) * I;
3885 Py_DECREF(cparent);
3886 return rnum;
3887 }
3888 if (y.is_zero()) {
3889 if (is_zero())
3890 throw (std::runtime_error("atan2(): division by zero"));
3891 if (real().is_positive())
3892 rnum = *_num0_p;
3893 else
3894 rnum = ex_to<numeric>(Pi.evalf(0, parent));
3895 }
3896 else {
3897 if (is_zero())
3898 rnum = ex_to<numeric>(Pi.evalf(0, parent)) / *_num2_p;
3899 else if (real().is_positive())
3900 rnum = (y/(*this)).abs().atan(parent);
3901 else
3902 rnum = (ex_to<numeric>(Pi.evalf(0, parent))
3903 - (y/(*this)).abs().atan(parent));
3904 if (not y.real().is_positive())
3905 rnum = rnum.negative();
3906 }
3907 rnum = ex_to<numeric>(rnum.evalf(0, parent));
3908 Py_DECREF(cparent);
3909 if (newdict)
3910 Py_DECREF(parent);
3911 return rnum;
3912 }
3913
sinh(PyObject * parent) const3914 const numeric numeric::sinh(PyObject* parent) const {
3915 return (exp(parent) - negative().exp(parent)) / *_num2_p;
3916 }
3917
cosh(PyObject * parent) const3918 const numeric numeric::cosh(PyObject* parent) const {
3919 return (exp(parent) + negative().exp(parent)) / *_num2_p;
3920 }
3921
tanh(PyObject * parent) const3922 const numeric numeric::tanh(PyObject* parent) const {
3923 const numeric& e2x = exp(parent);
3924 const numeric& e2nx = negative().exp(parent);
3925 return (e2x - e2nx)/(e2x + e2nx);
3926 }
3927
asinh(PyObject * parent) const3928 const numeric numeric::asinh(PyObject* parent) const {
3929 return arbfunc_0arg("arcsinh", parent);
3930 }
3931
acosh(PyObject * parent) const3932 const numeric numeric::acosh(PyObject* parent) const {
3933 return arbfunc_0arg("arccosh", parent);
3934 }
3935
atanh(PyObject * parent) const3936 const numeric numeric::atanh(PyObject* parent) const {
3937 int prec = precision(*this, parent);
3938 PyObject* field = CBF(prec+15);
3939 PyObject* ret = CallBallMethod0Arg(field, const_cast<char*>("arctanh"), *this);
3940 Py_DECREF(field);
3941
3942 numeric rnum(ret);
3943 if ((is_real() or imag().is_zero())
3944 and abs()<(*_num1_p))
3945 return ex_to<numeric>(rnum.real().evalf(0, parent));
3946
3947 return ex_to<numeric>(rnum.evalf(0, parent));
3948 }
3949
Li2(const numeric & n,PyObject * parent) const3950 const numeric numeric::Li2(const numeric &n, PyObject* parent) const {
3951 PyObject *cparent = common_parent(*this, n);
3952 if (parent == nullptr)
3953 parent = cparent;
3954 int prec = precision(*this, parent);
3955 PyObject* field = CBF(prec+15);
3956 PyObject* ret = CallBallMethod1Arg(field, const_cast<char*>("polylog"), *this, n);
3957 Py_DECREF(field);
3958
3959 numeric rnum(ret), nret;
3960 if ((is_real() or imag().is_zero())
3961 and n.is_integer() and rnum.real()<(*_num1_p))
3962 nret = ex_to<numeric>(rnum.real().evalf(0, parent));
3963 else
3964 nret = ex_to<numeric>(rnum.evalf(0, parent));
3965 Py_DECREF(cparent);
3966 return nret;
3967 }
3968
lgamma(PyObject * parent) const3969 const numeric numeric::lgamma(PyObject* parent) const {
3970 int prec = precision(*this, parent);
3971 PyObject* field = CBF(prec+15);
3972 PyObject* ret = CallBallMethod0Arg(field, const_cast<char*>("log_gamma"), *this);
3973 Py_DECREF(field);
3974
3975 numeric rnum(ret);
3976 return ex_to<numeric>(rnum.evalf(0, parent));
3977 }
3978
gamma(PyObject * parent) const3979 const numeric numeric::gamma(PyObject* parent) const {
3980 return arbfunc_0arg("gamma", parent);
3981 }
3982
rgamma(PyObject * parent) const3983 const numeric numeric::rgamma(PyObject* parent) const {
3984 return arbfunc_0arg("rgamma", parent);
3985 }
3986
psi(PyObject * parent) const3987 const numeric numeric::psi(PyObject* parent) const {
3988 return arbfunc_0arg("psi", parent);
3989 }
3990
psi(const numeric & y) const3991 const numeric numeric::psi(const numeric& y) const {
3992 PY_RETURN2(py_funcs.py_psi2, y);
3993 }
3994
zeta() const3995 const numeric numeric::zeta() const {
3996 PY_RETURN(py_funcs.py_zeta);
3997 }
3998
stieltjes() const3999 const numeric numeric::stieltjes() const {
4000 PY_RETURN(py_funcs.py_stieltjes);
4001 }
4002
factorial() const4003 const numeric numeric::factorial() const {
4004 static long fac[] = {1, 1, 2, 6, 24, 120, 720,
4005 5040, 40320, 362880, 3628800, 39916800,
4006 479001600};
4007 if (is_integer()) {
4008 if (is_positive() and *this < 13)
4009 return fac[to_long()];
4010 mpz_t bigint;
4011 mpz_init(bigint);
4012 mpz_fac_ui(bigint, to_long());
4013 return bigint;
4014 }
4015 PY_RETURN(py_funcs.py_factorial);
4016 }
4017
doublefactorial() const4018 const numeric numeric::doublefactorial() const {
4019 PY_RETURN(py_funcs.py_doublefactorial);
4020 }
4021
binomial(const numeric & k) const4022 const numeric numeric::binomial(const numeric &k) const {
4023 if ((t == LONG or t == MPZ)
4024 and k.is_integer()) {
4025 if (is_positive() and k.is_positive()
4026 and *this < 13) {
4027 static long fac[] = {1, 1, 2, 6, 24, 120, 720,
4028 5040, 40320, 362880, 3628800, 39916800,
4029 479001600};
4030 long a = to_long();
4031 long b = k.to_long();
4032 if (b<=0 or b>12)
4033 return *_num0_p;
4034 return fac[a]/fac[b]/fac[a-b];
4035 }
4036 mpz_t bigint;
4037 mpz_init(bigint);
4038 if (t == MPZ) {
4039 mpz_bin_ui(bigint, v._bigint, k.to_long());
4040 }
4041 else {
4042 mpz_set_ui(bigint, v._long);
4043 mpz_bin_ui(bigint, bigint, k.to_long());
4044 }
4045 return bigint;
4046 }
4047 numeric res;
4048 PyObject *nobj, *kobj;
4049 nobj = this->to_pyobject();
4050 kobj = k.to_pyobject();
4051
4052 PyObject* m = PyImport_ImportModule("sage.arith.misc");
4053 if (m == nullptr)
4054 py_error("Error importing arith.misc");
4055 PyObject* binfunc = PyObject_GetAttrString(m, "binomial");
4056 if (binfunc == nullptr)
4057 py_error("Error getting binomial");
4058
4059 PyObject* pyresult = PyObject_CallFunctionObjArgs(binfunc, nobj, kobj, NULL);
4060 Py_DECREF(kobj);
4061 Py_DECREF(nobj);
4062 Py_DECREF(m);
4063 Py_DECREF(binfunc);
4064 if (pyresult == nullptr) {
4065 throw(std::runtime_error("numeric::binomial(): python function binomial raised exception"));
4066 }
4067 if ( pyresult == Py_None ) {
4068 throw(std::runtime_error("numeric::binomial: python function binomial returned None"));
4069 }
4070 // convert output Expression to an ex
4071 ex eval_result = py_funcs.pyExpression_to_ex(pyresult);
4072 Py_DECREF(pyresult);
4073 if (PyErr_Occurred() != nullptr) {
4074 throw(std::runtime_error("numeric::binomial(): python function (Expression_to_ex) raised exception"));
4075 }
4076 return ex_to<numeric>(eval_result);
4077 }
4078
binomial(unsigned long n,unsigned long k)4079 const numeric numeric::binomial(unsigned long n, unsigned long k)
4080 {
4081 if (n < 13) {
4082 static long fac[] = {1, 1, 2, 6, 24, 120, 720,
4083 5040, 40320, 362880, 3628800, 39916800,
4084 479001600};
4085 if (k == 0)
4086 return *_num1_p;
4087 if (k > n)
4088 return *_num0_p;
4089 return fac[n]/fac[k]/fac[n-k];
4090 }
4091 mpz_t bigint;
4092 mpz_init(bigint);
4093 mpz_bin_uiui(bigint, n, k);
4094 return bigint;
4095 }
4096
bernoulli() const4097 const numeric numeric::bernoulli() const {
4098 PY_RETURN(py_funcs.py_bernoulli);
4099 }
4100
hypergeometric_pFq(const std::vector<numeric> & a,const std::vector<numeric> & b,PyObject * parent) const4101 const numeric numeric::hypergeometric_pFq(const std::vector<numeric>& a, const std::vector<numeric>& b, PyObject *parent) const {
4102 PyObject *lista = py_tuple_from_numvector(a);
4103 PyObject *listb = py_tuple_from_numvector(b);
4104 PyObject *z = to_pyobject();
4105
4106 // call opt.evalf_f with this tuple
4107 PyObject* m = PyImport_ImportModule("sage.functions.hypergeometric");
4108 if (m == nullptr)
4109 py_error("Error importing hypergeometric");
4110 PyObject* hypfunc = PyObject_GetAttrString(m, "hypergeometric");
4111 if (hypfunc == nullptr)
4112 py_error("Error getting hypergeometric attribute");
4113
4114 if ((parent != nullptr) and PyDict_CheckExact(parent)) {
4115 Py_DECREF(z);
4116 z = ex_to<numeric>(this->evalf(0, parent)).to_pyobject();
4117 }
4118 PyObject* name = PyString_FromString(const_cast<char*>("_evalf_try_"));
4119 PyObject* pyresult = PyObject_CallMethodObjArgs(hypfunc, name, lista, listb, z, NULL);
4120 Py_DECREF(m);
4121 Py_DECREF(z);
4122 Py_DECREF(name);
4123 Py_DECREF(hypfunc);
4124 if (pyresult == nullptr) {
4125 throw(std::runtime_error("numeric::hypergeometric_pFq(): python function hypergeometric::_evalf_ raised exception"));
4126 }
4127 if ( pyresult == Py_None ) {
4128 throw(std::runtime_error("numeric::hypergeometric_pFq(): python function hypergeometric::_evalf_ returned None"));
4129 }
4130 // convert output Expression to an ex
4131 ex eval_result = py_funcs.pyExpression_to_ex(pyresult);
4132 Py_DECREF(pyresult);
4133 if (PyErr_Occurred() != nullptr) {
4134 throw(std::runtime_error("numeric::hypergeometric_pFq(): python function (Expression_to_ex) raised exception"));
4135 }
4136 return ex_to<numeric>(eval_result);
4137 }
4138
isqrt() const4139 const numeric numeric::isqrt() const {
4140 PY_RETURN(py_funcs.py_isqrt);
4141 }
4142
sqrt() const4143 const numeric numeric::sqrt() const {
4144 PY_RETURN(py_funcs.py_sqrt);
4145 }
4146
is_square() const4147 bool numeric::is_square() const
4148 {
4149 if (is_negative())
4150 return false;
4151 if (is_zero() or is_one())
4152 return true;
4153 switch (t) {
4154 case LONG: {
4155 long rt = std::lround(std::sqrt(v._long));
4156 return rt * rt == v._long;
4157 }
4158 case MPZ:
4159 return mpz_perfect_square_p(v._bigint);
4160 case MPQ:
4161 return (mpz_perfect_square_p(mpq_numref(v._bigrat))
4162 and mpz_perfect_square_p(mpq_denref(v._bigrat)));
4163 case PYOBJECT:
4164 default:
4165 stub("invalid type: type not handled");
4166 }
4167 return false;
4168 }
4169
4170 // Fast sqrt in case of square, symbolic sqrt else.
sqrt_as_ex() const4171 const ex numeric::sqrt_as_ex() const
4172 {
4173 if (is_negative())
4174 return I * (-*this).sqrt_as_ex();
4175 if (is_zero())
4176 return _ex0;
4177 if (is_one())
4178 return _ex1;
4179 switch (t) {
4180 case LONG: {
4181 long rt = std::lround(std::sqrt(v._long));
4182 if (rt * rt == v._long)
4183 return ex(numeric(rt));
4184 break;
4185 }
4186 case MPZ: {
4187 if (mpz_perfect_square_p(v._bigint)) {
4188 mpz_t bigint;
4189 mpz_init(bigint);
4190 mpz_sqrt(bigint, v._bigint);
4191 return ex(numeric(bigint));
4192 }
4193 break;
4194 }
4195 case MPQ: {
4196 if (mpz_perfect_square_p(mpq_numref(v._bigrat))
4197 and mpz_perfect_square_p(mpq_denref(v._bigrat))) {
4198 mpz_t bigint;
4199 mpq_t bigrat, obigrat;
4200 mpz_init(bigint);
4201 mpq_init(bigrat);
4202 mpq_init(obigrat);
4203 mpz_sqrt(bigint, mpq_numref(v._bigrat));
4204 mpq_set_z(bigrat, bigint);
4205 mpz_sqrt(bigint, mpq_denref(v._bigrat));
4206 mpq_set_z(obigrat, bigint);
4207 mpq_div(bigrat, bigrat, obigrat);
4208 mpz_clear(bigint);
4209 mpq_clear(obigrat);
4210 return ex(numeric(bigrat));
4211 }
4212 break;
4213 }
4214 case PYOBJECT:
4215 return sqrt();
4216 default:
4217 stub("invalid type: type not handled");
4218 }
4219 return (new GiNaC::power(ex(*this), _ex1_2))->setflag(status_flags::dynallocated | status_flags::evaluated);
4220 }
4221
abs() const4222 const numeric numeric::abs() const {
4223 switch (t) {
4224 case LONG:
4225 return std::labs(v._long);
4226 case MPZ: {
4227 mpz_t bigint;
4228 mpz_init(bigint);
4229 mpz_abs(bigint, v._bigint);
4230 return bigint;
4231 }
4232 case MPQ: {
4233 mpq_t bigrat;
4234 mpq_init(bigrat);
4235 mpq_abs(bigrat, v._bigrat);
4236 return bigrat;
4237 }
4238 case PYOBJECT: {
4239 PyObject *ret = PyNumber_Absolute(v._pyobject);
4240 if (ret == NULL) {
4241 PyErr_Clear();
4242 return *this;
4243 }
4244 return ret;
4245 }
4246 default:
4247 stub("invalid type: type not handled");
4248 }
4249 }
4250
mod(const numeric & b) const4251 const numeric numeric::mod(const numeric &b) const {
4252 switch (t) {
4253 case LONG:
4254 if (b.t == LONG)
4255 return v._long % b.v._long;
4256 if (b.t == MPZ)
4257 return to_bigint().mod(b);
4258 throw std::runtime_error("unsupported type in numeric::md");
4259 case MPZ:
4260 if (b.t == LONG)
4261 return mod(b.to_bigint());
4262 if (b.t == MPZ) {
4263 mpz_t bigint;
4264 mpz_init(bigint);
4265 mpz_mod(bigint, v._bigint, b.v._bigint);
4266 return bigint;
4267 }
4268 throw std::runtime_error("unsupported type in numeric::md");
4269 case MPQ:
4270 case PYOBJECT: {
4271 PY_RETURN2(py_funcs.py_mod, b);
4272 }
4273 default:
4274 stub("invalid type: type not handled");
4275 }
4276 }
4277
_smod(const numeric & b) const4278 const numeric numeric::_smod(const numeric &b) const {
4279 PY_RETURN2(py_funcs.py_smod, b);
4280 }
4281
irem(const numeric & b) const4282 const numeric numeric::irem(const numeric &b) const {
4283 switch (t) {
4284 case LONG:
4285 if (b.t == LONG)
4286 return std::remainder(v._long, b.v._long);
4287 if (b.t == MPZ)
4288 return to_bigint().irem(b);
4289 throw std::runtime_error("unsupported type in numeric::irem");
4290 case MPZ:
4291 if (b.t == LONG)
4292 return irem(b.to_bigint());
4293 if (b.t == MPZ) {
4294 mpz_t bigint;
4295 mpz_init(bigint);
4296 mpz_fdiv_r(bigint, v._bigint, b.v._bigint);
4297 return bigint;
4298 }
4299 throw std::runtime_error("unsupported type in numeric::irem");
4300 case MPQ:
4301 case PYOBJECT: {
4302 PY_RETURN2(py_funcs.py_irem, b);
4303 }
4304 default:
4305 stub("invalid type: type not handled");
4306 }
4307 }
4308
iquo(const numeric & b) const4309 const numeric numeric::iquo(const numeric &b) const {
4310 switch (t) {
4311 case LONG:
4312 if (b.t == LONG)
4313 return v._long / b.v._long;
4314 if (b.t == MPZ)
4315 return to_bigint().iquo(b);
4316 throw std::runtime_error("unsupported type in numeric::iquo");
4317 case MPZ:
4318 if (b.t == LONG) {
4319 mpz_t bigint;
4320 mpz_init(bigint);
4321 mpz_fdiv_q_ui(bigint, v._bigint,
4322 std::labs(b.v._long));
4323 if (b.v._long < 0)
4324 mpz_neg(bigint, bigint);
4325 return bigint;
4326 }
4327 if (b.t == MPZ) {
4328 mpz_t bigint;
4329 mpz_init(bigint);
4330 mpz_fdiv_q(bigint, v._bigint, b.v._bigint);
4331 return bigint;
4332 }
4333 throw std::runtime_error("unsupported type in numeric::iquo");
4334 case MPQ:
4335 case PYOBJECT:
4336 default:
4337 stub("invalid type: type not handled");
4338 }
4339 throw std::runtime_error("iquo: bad input");
4340 }
4341
iquo(const numeric & b,numeric & r) const4342 const numeric numeric::iquo(const numeric &b, numeric& r) const {
4343 switch (t) {
4344 case LONG:
4345 if (b.t == LONG) {
4346 auto ld = std::ldiv(v._long, b.v._long);
4347 r = ld.rem;
4348 return ld.quot;
4349 }
4350 if (b.t == MPZ)
4351 return to_bigint().iquo(b, r);
4352 throw std::runtime_error("unsupported type in numeric::iquo");
4353 case MPZ:
4354 if (b.t == LONG) {
4355 mpz_t bigint;
4356 mpz_init(bigint);
4357 r = long(mpz_fdiv_q_ui(bigint, v._bigint,
4358 std::labs(b.v._long)));
4359 return bigint;
4360 }
4361 if (b.t == MPZ) {
4362 mpz_t bigint, tmp;
4363 mpz_init(bigint);
4364 mpz_init(tmp);
4365 mpz_fdiv_q(bigint, v._bigint, b.v._bigint);
4366 mpz_mul(tmp, bigint, b.v._bigint);
4367 mpz_sub(tmp, v._bigint, tmp);
4368 r = numeric(tmp);
4369 return bigint;
4370 }
4371 throw std::runtime_error("unsupported type in numeric::iquo");
4372 case MPQ:
4373 case PYOBJECT:
4374 default:
4375 stub("invalid type: type not handled");
4376 }
4377 throw std::runtime_error("iquo2: bad input");
4378 }
4379
gcd(const numeric & B) const4380 const numeric numeric::gcd(const numeric &B) const
4381 {
4382 if (is_zero())
4383 return B.abs();
4384 if (B.is_zero())
4385 return abs();
4386 numeric a,b;
4387 if (is_negative())
4388 a = negative();
4389 else
4390 a = *this;
4391 if (B.is_negative())
4392 b = B.negative();
4393 else
4394 b = B;
4395 if (a.is_one() or b.is_one())
4396 return *_num1_p;
4397
4398 switch (a.t) {
4399 case LONG:
4400 switch (b.t) {
4401 case LONG: {
4402 // C++17 will have a function for this
4403 long al = a.v._long;
4404 long bl = b.v._long;
4405 long c;
4406 while (al != 0) {
4407 c = al;
4408 al = bl % al;
4409 bl = c;
4410 }
4411 return bl;
4412 }
4413 case MPZ:
4414 return a.to_bigint().gcd(b);
4415 case MPQ:
4416 return (a * b.denom()).gcd(b.numer()) / b.denom();
4417 case PYOBJECT:
4418 return *_num1_p;
4419 default: stub("invalid type: type not handled");
4420 }
4421 case MPZ:
4422 switch (b.t) {
4423 case LONG: {
4424 mpz_t bigint;
4425 mpz_init(bigint);
4426 long l = mpz_gcd_ui(bigint,
4427 a.v._bigint,
4428 std::labs(b.v._long));
4429 mpz_clear(bigint);
4430 return l;
4431 }
4432 case MPZ: {
4433 mpz_t bigint;
4434 mpz_init(bigint);
4435 mpz_gcd(bigint, a.v._bigint, b.v._bigint);
4436 return bigint;
4437 }
4438 case MPQ:
4439 return (a * b.denom()).gcd(b.numer()) / b.denom();
4440 case PYOBJECT:
4441 return *_num1_p;
4442 default:
4443 stub("invalid type: type not handled");
4444 }
4445 case MPQ:
4446 if (b.t == PYOBJECT)
4447 return *_num1_p;
4448 return (a.numer() * b.denom()).gcd(a.denom() * b.numer()) / (a.denom() * b.denom());
4449 stub("invalid type: type not handled");
4450 case PYOBJECT:
4451 return *_num1_p;
4452 default:
4453 stub("invalid type: type not handled");
4454 }
4455 }
4456
lcm(const numeric & b) const4457 const numeric numeric::lcm(const numeric &b) const
4458 {
4459 if (is_zero() or b.is_zero())
4460 return *_num0_p;
4461 if (is_one())
4462 return b;
4463 if (b.is_one())
4464 return *this;
4465 switch (t) {
4466 case LONG:
4467 if (b.t == LONG) {
4468 // C++17 will have a function for this
4469 long a = v._long;
4470 long bl = b.v._long;
4471 long c;
4472 while (a != 0) {
4473 c = a;
4474 a = bl % a;
4475 bl = c;
4476 }
4477 return (v._long / bl) * b.v._long;
4478 }
4479 if (b.t == MPZ)
4480 return to_bigint().lcm(b);
4481 throw std::runtime_error("unsupported type in numeric::lcm");
4482 case MPZ:
4483 if (b.t == LONG)
4484 return lcm(b.to_bigint());
4485 if (b.t == MPZ) {
4486 mpz_t bigint;
4487 mpz_init(bigint);
4488 mpz_lcm(bigint, v._bigint, b.v._bigint);
4489 return bigint;
4490 }
4491 throw std::runtime_error("unsupported type in numeric::lcm");
4492 case MPQ:
4493 case PYOBJECT: {
4494 PY_RETURN2(py_funcs.py_lcm, b);
4495 }
4496 default:
4497 stub("invalid type: type not handled");
4498 }
4499 }
4500
4501 // version of factor for factors fitting into long
factorsmall(std::vector<std::pair<long,int>> & factors,long range) const4502 void numeric::factorsmall(std::vector<std::pair<long, int>>& factors, long range) const
4503 {
4504 if (is_one() or is_zero() or is_minus_one())
4505 return;
4506 switch (t) {
4507 case LONG:
4508 return to_bigint().factorsmall(factors, range);
4509 case MPZ: {
4510 fmpz_t f;
4511 fmpz_init(f);
4512 mpz_t bigint;
4513 mpz_init(bigint);
4514 mpz_abs(bigint, v._bigint);
4515 fmpz_set_mpz(f, bigint);
4516 fmpz_factor_t fs;
4517 fmpz_factor_init(fs);
4518 if (range == 0)
4519 fmpz_factor(fs, f);
4520 else
4521 fmpz_factor_trial_range(fs, f, 0, range);
4522 for (slong i = 0; i < fs->num; ++i) {
4523 fmpz_get_mpz(bigint, fs->p + i);
4524 factors.push_back(std::make_pair(mpz_get_si(bigint),
4525 int(*(fs->exp + i))));
4526 }
4527 mpz_clear(bigint);
4528 fmpz_factor_clear(fs);
4529 fmpz_clear(f);
4530 return;
4531 }
4532 case MPQ:
4533 return to_bigint().factorsmall(factors);
4534 case PYOBJECT:
4535 default:
4536 stub("invalid type: type not handled");
4537 }
4538 }
4539
4540 // may take long time if no range given
factor(std::vector<std::pair<numeric,int>> & factors,long range) const4541 void numeric::factor(std::vector<std::pair<numeric, int>>& factors, long range) const
4542 {
4543 if (is_one() or is_minus_one())
4544 return;
4545 switch (t) {
4546 case LONG: {
4547 std::vector<std::pair<long, int>> f;
4548 factorsmall(f, range);
4549 for (const auto& p : f)
4550 factors.emplace_back(numeric(p.first), p.second);
4551 return;
4552 }
4553 case MPZ: {
4554 fmpz_t f;
4555 fmpz_init(f);
4556 mpz_t bigint, tmp;
4557 mpz_init(bigint);
4558 mpz_abs(bigint, v._bigint);
4559 fmpz_set_mpz(f, bigint);
4560 fmpz_factor_t fs;
4561 fmpz_factor_init(fs);
4562 if (range == 0)
4563 fmpz_factor(fs, f);
4564 else
4565 fmpz_factor_trial_range(fs, f, 0, range);
4566 for (slong i = 0; i < fs->num; ++i) {
4567 mpz_init(tmp);
4568 fmpz_get_mpz(tmp, fs->p + i);
4569 if (range != 0)
4570 for (int j = 0; j < int(*(fs->exp + i)); ++j) {
4571 mpz_divexact(bigint,
4572 bigint,
4573 tmp);
4574 }
4575 factors.emplace_back(numeric(tmp),
4576 int(*(fs->exp + i)));
4577 }
4578 fmpz_clear(f);
4579 fmpz_factor_clear(fs);
4580 if (range != 0 and mpz_cmp_ui(bigint, 1) != 0)
4581 factors.push_back(std::make_pair(numeric(bigint), 1));
4582 else
4583 mpz_clear(bigint);
4584 return;
4585 }
4586 case MPQ:
4587 to_bigint().factor(factors);
4588 return;
4589 case PYOBJECT:
4590 default:
4591 stub("invalid type: type not handled");
4592 }
4593 }
4594
setDivisors(long n,int i,std::set<int> & divisors,const std::vector<std::pair<long,int>> & factors)4595 static void setDivisors(long n, int i, std::set<int>& divisors,
4596 const std::vector<std::pair<long, int>>& factors)
4597 // author of this little gem:
4598 // http://stackoverflow.com/users/3864373/rocky-johnson
4599 {
4600 for (size_t j = i; j<factors.size(); j++) {
4601 long x = factors[j].first * n;
4602 for (int k = 0; k<factors[j].second; k++) {
4603 divisors.insert(x);
4604 setDivisors(x, j + 1, divisors, factors);
4605 x *= factors[j].first;
4606 }
4607 }
4608 }
4609
divisors(std::set<int> & divs) const4610 void numeric::divisors(std::set<int>& divs) const
4611 {
4612 divs.insert(1);
4613 if (is_one() or is_zero() or is_minus_one())
4614 return;
4615 switch (t) {
4616 case LONG:
4617 case MPZ: {
4618 std::vector<std::pair<long, int>> factors;
4619 factorsmall(factors);
4620 setDivisors(1, 0, divs, factors);
4621 return;
4622 }
4623 case MPQ:
4624 return to_bigint().divisors(divs);
4625 case PYOBJECT:
4626 default:
4627 stub("invalid type: type not handled");
4628 }
4629 }
4630
4631 //////////
4632 // global
4633 //////////
4634
coerce(numeric & new_left,numeric & new_right,const numeric & left,const numeric & right)4635 void coerce(numeric& new_left, numeric& new_right, const numeric& left, const numeric& right) {
4636 verbose("coerce");
4637 // Return a numeric
4638 if (left.t == right.t) {
4639 new_left = left;
4640 new_right = right;
4641 return;
4642 }
4643 PyObject *o;
4644 switch (left.t) {
4645 case LONG:
4646 switch (right.t) {
4647 case MPZ:
4648 if (mpz_fits_sint_p(right.v._bigint)) {
4649 new_right = numeric(mpz_get_si(right.v._bigint));
4650 new_left = left;
4651 } else {
4652 numeric n;
4653 mpz_init(n.v._bigint);
4654 n.t = MPZ;
4655 mpz_set_si(n.v._bigint, left.v._long);
4656 n.hash = (left.v._long == -1) ? -2 : left.v._long;
4657 new_left = n;
4658 new_right = right;
4659 }
4660 return;
4661 case MPQ: {
4662 mpq_t bigrat;
4663 mpq_init(bigrat);
4664 mpq_set_si(bigrat, left.v._long, 1);
4665 new_left = numeric(bigrat);
4666 new_right = right;
4667 return;
4668 }
4669 case PYOBJECT: {
4670 mpz_t bigint;
4671 mpz_init_set_si(bigint, left.v._long);
4672 o = py_funcs.py_integer_from_mpz(bigint);
4673 new_left = numeric(o, true);
4674 new_right = right;
4675 mpz_clear(bigint);
4676 return;
4677 }
4678 default:
4679 std::cerr << "type = " << right.t << "\n";
4680 stub("** invalid coercion -- left MPZ**");
4681 }
4682 case MPZ:
4683 switch (right.t) {
4684 case LONG:
4685 if (mpz_fits_sint_p(left.v._bigint)) {
4686 new_left = numeric(mpz_get_si(left.v._bigint));
4687 new_right = right;
4688 } else {
4689 numeric n;
4690 mpz_init(n.v._bigint);
4691 n.t = MPZ;
4692 mpz_set_si(n.v._bigint, right.v._long);
4693 n.hash = (right.v._long == -1) ? -2 : right.v._long;
4694 new_right = n;
4695 new_left = left;
4696 }
4697 return;
4698 case MPQ: {
4699 mpq_t bigrat;
4700 mpq_init(bigrat);
4701 mpq_set_z(bigrat, left.v._bigint);
4702 new_left = numeric(bigrat);
4703 new_right = right;
4704 return;
4705 }
4706 case PYOBJECT: {
4707 mpz_t bigint;
4708 mpz_init_set(bigint, left.v._bigint);
4709 o = py_funcs.py_integer_from_mpz(bigint);
4710 new_left = numeric(o, true);
4711 new_right = right;
4712 mpz_clear(bigint);
4713 return;
4714 }
4715 default:
4716 std::cerr << "type = " << right.t << "\n";
4717 stub("** invalid coercion -- left MPZ**");
4718 }
4719 case MPQ:
4720 switch (right.t) {
4721 case LONG: {
4722 mpq_t bigrat;
4723 mpq_init(bigrat);
4724 mpq_set_si(bigrat, right.v._long, 1);
4725 new_right = numeric(bigrat);
4726 new_left = left;
4727 return;
4728 }
4729 case MPZ: {
4730 mpq_t bigrat;
4731 mpq_init(bigrat);
4732 mpq_set_z(bigrat, right.v._bigint);
4733 new_left = left;
4734 new_right = numeric(bigrat);
4735 return;
4736 }
4737 case PYOBJECT: {
4738 mpq_t bigrat;
4739 mpq_init(bigrat);
4740 mpq_set(bigrat, left.v._bigrat);
4741 o = py_funcs.py_rational_from_mpq(bigrat);
4742 mpq_clear(bigrat);
4743 new_left = numeric(o, true);
4744 new_right = right;
4745 return;
4746 }
4747 default:
4748 std::cerr << "type = " << right.t << "\n";
4749 stub("** invalid coercion -- left MPQ**");
4750 }
4751 case PYOBJECT:
4752 new_left = left;
4753 switch (right.t) {
4754 case LONG: {
4755 mpz_t bigint;
4756 mpz_init_set_si(bigint, right.v._long);
4757 o = py_funcs.py_integer_from_mpz(bigint);
4758 mpz_clear(bigint);
4759 new_right = numeric(o, true);
4760 return;
4761 }
4762 case MPZ: {
4763 mpz_t bigint;
4764 mpz_init_set(bigint, right.v._bigint);
4765 o = py_funcs.py_integer_from_mpz(bigint);
4766 mpz_clear(bigint);
4767 new_right = numeric(o, true);
4768 return;
4769 }
4770 case MPQ: {
4771 mpq_t bigrat;
4772 mpq_init(bigrat);
4773 mpq_set(bigrat, right.v._bigrat);
4774 o = py_funcs.py_rational_from_mpq(bigrat);
4775 mpq_clear(bigrat);
4776 new_right = numeric(o, true);
4777 return;
4778 }
4779 default:
4780 std::cerr << "type = " << right.t << "\n";
4781 stub("** invalid coercion -- left PYOBJECT**");
4782 }
4783 default:
4784 std::cerr << "type = " << left.t << "\n";
4785 stub("** invalid coercion **");
4786 }
4787 }
4788
4789 /** Imaginary unit. This is not a constant but a numeric since we are
4790 * natively handing complex numbers anyways, so in each expression containing
4791 * an I it is automatically eval'ed away anyhow. */
4792
4793 /** Exponential function.
4794 *
4795 * @return arbitrary precision numerical exp(x). */
exp(const numeric & x,PyObject * parent)4796 const numeric exp(const numeric &x, PyObject* parent) {
4797 return x.exp(parent);
4798 }
4799
4800 /** Natural logarithm.
4801 *
4802 * @param x complex number
4803 * @return arbitrary precision numerical log(x).
4804 * @exception pole_error("log(): logarithmic pole",0) */
log(const numeric & x,PyObject * parent)4805 const numeric log(const numeric &x, PyObject* parent) {
4806 return x.log(parent);
4807 }
4808
4809 // General log
log(const numeric & x,const numeric & b,PyObject * parent)4810 const numeric log(const numeric &x, const numeric &b, PyObject* parent) {
4811 return x.log(b, parent);
4812 }
4813
4814 /** Numeric sine (trigonometric function).
4815 *
4816 * @return arbitrary precision numerical sin(x). */
sin(const numeric & x)4817 const numeric sin(const numeric &x) {
4818 return x.sin();
4819 }
4820
4821 /** Numeric cosine (trigonometric function).
4822 *
4823 * @return arbitrary precision numerical cos(x). */
cos(const numeric & x)4824 const numeric cos(const numeric &x) {
4825 return x.cos();
4826 }
4827
4828 /** Numeric tangent (trigonometric function).
4829 *
4830 * @return arbitrary precision numerical tan(x). */
tan(const numeric & x)4831 const numeric tan(const numeric &x) {
4832 return x.tan();
4833 }
4834
4835 /** Numeric inverse sine (trigonometric function).
4836 *
4837 * @return arbitrary precision numerical asin(x). */
asin(const numeric & x,PyObject * parent)4838 const numeric asin(const numeric &x, PyObject* parent) {
4839 return x.asin(parent);
4840 }
4841
4842 /** Numeric inverse cosine (trigonometric function).
4843 *
4844 * @return arbitrary precision numerical acos(x). */
acos(const numeric & x,PyObject * parent)4845 const numeric acos(const numeric &x, PyObject* parent) {
4846 return x.acos(parent);
4847 }
4848
4849 /** Numeric arcustangent.
4850 *
4851 * @param x complex number
4852 * @return atan(x)
4853 * @exception pole_error("atan(): logarithmic pole",0) if x==I or x==-I. */
atan(const numeric & x,PyObject * parent)4854 const numeric atan(const numeric &x, PyObject* parent) {
4855 if (!x.is_real() &&
4856 x.real().is_zero() &&
4857 abs(x.imag()).is_one())
4858 throw pole_error("atan(): logarithmic pole", 0);
4859 return x.atan(parent);
4860 }
4861
4862 /** Numeric arcustangent of two arguments, analytically continued in a suitable way.
4863 *
4864 * @param y complex number
4865 * @param x complex number
4866 * @return -I*log((x+I*y)/sqrt(x^2+y^2)), which is equal to atan(y/x) if y and
4867 * x are both real.
4868 * @exception pole_error("atan(): logarithmic pole",0) if y/x==+I or y/x==-I. */
atan(const numeric & y,const numeric & x,PyObject * parent)4869 const numeric atan(const numeric &y, const numeric &x, PyObject* parent) {
4870 return x.atan(y, parent);
4871 }
4872
4873 /** Numeric hyperbolic sine (trigonometric function).
4874 *
4875 * @return arbitrary precision numerical sinh(x). */
sinh(const numeric & x,PyObject * parent)4876 const numeric sinh(const numeric &x, PyObject* parent) {
4877 return x.sinh(parent);
4878 }
4879
4880 /** Numeric hyperbolic cosine (trigonometric function).
4881 *
4882 * @return arbitrary precision numerical cosh(x). */
cosh(const numeric & x,PyObject * parent)4883 const numeric cosh(const numeric &x, PyObject* parent) {
4884 return x.cosh(parent);
4885 }
4886
4887 /** Numeric hyperbolic tangent (trigonometric function).
4888 *
4889 * @return arbitrary precision numerical tanh(x). */
tanh(const numeric & x,PyObject * parent)4890 const numeric tanh(const numeric &x, PyObject* parent) {
4891 return x.tanh(parent);
4892 }
4893
4894 /** Numeric inverse hyperbolic sine (trigonometric function).
4895 *
4896 * @return arbitrary precision numerical asinh(x). */
asinh(const numeric & x,PyObject * parent)4897 const numeric asinh(const numeric &x, PyObject* parent) {
4898 return x.asinh(parent);
4899 }
4900
4901 /** Numeric inverse hyperbolic cosine (trigonometric function).
4902 *
4903 * @return arbitrary precision numerical acosh(x). */
acosh(const numeric & x,PyObject * parent)4904 const numeric acosh(const numeric &x, PyObject* parent) {
4905 return x.acosh(parent);
4906 }
4907
4908 /** Numeric inverse hyperbolic tangent (trigonometric function).
4909 *
4910 * @return arbitrary precision numerical atanh(x). */
atanh(const numeric & x,PyObject * parent)4911 const numeric atanh(const numeric &x, PyObject* parent) {
4912 return x.atanh(parent);
4913 }
4914
4915
Li2(const numeric & x,PyObject * parent)4916 const numeric Li2(const numeric &x, PyObject* parent) {
4917 try {
4918 return x.try_py_method("dilog");
4919 }
4920 catch (std::logic_error) {}
4921 try {
4922 return x.try_py_method("polylog", *_num2_p);
4923 }
4924 catch (std::logic_error) {}
4925
4926 return x.Li2(*_num2_p, parent);
4927 }
4928
Li2(const numeric & n,const numeric & x,PyObject * parent)4929 const numeric Li2(const numeric &n, const numeric &x, PyObject* parent) {
4930 return x.Li2(n, parent);
4931 }
4932
4933 /** Evaluation of Riemann's Zeta function. */
zeta(const numeric & x)4934 const numeric zeta(const numeric &x) {
4935 return x.zeta();
4936 }
4937
stieltjes(const numeric & x)4938 const numeric stieltjes(const numeric &x) {
4939 return x.stieltjes();
4940 }
4941
lgamma(const numeric & x,PyObject * parent)4942 const numeric lgamma(const numeric &x, PyObject* parent) {
4943 return x.lgamma(parent);
4944 }
4945
4946 /** The Gamma function. */
gamma(const numeric & x,PyObject * parent)4947 const numeric gamma(const numeric &x, PyObject* parent) {
4948 return x.gamma(parent);
4949 }
4950
4951 /** The psi function (aka polygamma function). */
psi(const numeric & x,PyObject * parent)4952 const numeric psi(const numeric &x, PyObject* parent) {
4953 return x.psi(parent);
4954 }
4955
4956 /** The psi functions (aka polygamma functions). */
psi(const numeric & n,const numeric & x)4957 const numeric psi(const numeric &n, const numeric &x) {
4958 return n.psi(x);
4959 }
beta(const numeric & x,const numeric & y,PyObject * parent)4960 const numeric beta(const numeric &x, const numeric &y, PyObject* parent)
4961 {
4962 PyObject *cparent = common_parent(x, y);
4963 if (parent == nullptr)
4964 parent = cparent;
4965 numeric ret = (x+y).rgamma(parent) * x.gamma(parent) * y.gamma(parent);
4966 Py_DECREF(cparent);
4967 return ret;
4968 }
4969
4970 /** Factorial combinatorial function.
4971 *
4972 * @param n integer argument >= 0
4973 * @exception range_error (argument must be integer >= 0) */
factorial(const numeric & n)4974 const numeric factorial(const numeric &n) {
4975 return n.factorial();
4976 }
4977
4978 /** The double factorial combinatorial function. (Scarcely used, but still
4979 * useful in cases, like for exact results of gamma(n+1/2) for instance.)
4980 *
4981 * @param n integer argument >= -1
4982 * @return n!! == n * (n-2) * (n-4) * ... * ({1|2}) with 0!! == (-1)!! == 1
4983 * @exception range_error (argument must be integer >= -1) */
doublefactorial(const numeric & n)4984 const numeric doublefactorial(const numeric &n) {
4985 return n.doublefactorial();
4986 }
4987
4988 /** The Binomial coefficients. It computes the binomial coefficients. For
4989 * integer n and k and positive n this is the number of ways of choosing k
4990 * objects from n distinct objects. If n is negative, the formula
4991 * binomial(n,k) == (-1)^k*binomial(k-n-1,k) is used to compute the result. */
binomial(const numeric & n,const numeric & k)4992 const numeric binomial(const numeric &n, const numeric &k) {
4993 return n.binomial(k);
4994 }
4995
4996 /** Bernoulli number. The nth Bernoulli number is the coefficient of x^n/n!
4997 * in the expansion of the function x/(e^x-1).
4998 *
4999 * @return the nth Bernoulli number (a rational number).
5000 * @exception range_error (argument must be integer >= 0) */
bernoulli(const numeric & n)5001 const numeric bernoulli(const numeric &n) {
5002 return n.bernoulli();
5003 }
5004
hypergeometric_pFq(const std::vector<numeric> & a,const std::vector<numeric> & b,const numeric & z,PyObject * parent)5005 const numeric hypergeometric_pFq(const std::vector<numeric>& a, const std::vector<numeric>& b, const numeric &z, PyObject* parent) {
5006 return z.hypergeometric_pFq(a, b, parent);
5007 }
5008
5009 /** Fibonacci number. The nth Fibonacci number F(n) is defined by the
5010 * recurrence formula F(n)==F(n-1)+F(n-2) with F(0)==0 and F(1)==1.
5011 *
5012 * @param n an integer
5013 * @return the nth Fibonacci number F(n) (an integer number)
5014 * @exception range_error (argument must be an integer) */
fibonacci(const numeric & n)5015 const numeric fibonacci(const numeric &n) {
5016 return n.fibonacci();
5017 }
5018
5019 /** Absolute value. */
abs(const numeric & x)5020 const numeric abs(const numeric& x) {
5021 return x.abs();
5022 }
5023
5024 /** Modulus (in positive representation).
5025 * In general, mod(a,b) has the sign of b or is zero, and rem(a,b) has the
5026 * sign of a or is zero. This is different from Maple's modp, where the sign
5027 * of b is ignored. It is in agreement with Mathematica's Mod.
5028 *
5029 * @return a mod b in the range [0,abs(b)-1] with sign of b if both are
5030 * integer, 0 otherwise. */
mod(const numeric & a,const numeric & b)5031 const numeric mod(const numeric &a, const numeric &b) {
5032 return a.mod(b);
5033 }
5034
5035 ///** Modulus (in symmetric representation).
5036 // * Equivalent to Maple's mods.
5037 // *
5038 // * @return a mod b in the range [-iquo(abs(b)-1,2), iquo(abs(b),2)]. */
smod(const numeric & a,const numeric & b)5039 const numeric smod(const numeric &a, const numeric &b) {
5040 return a._smod(b);
5041 }
5042
5043 /** Numeric integer remainder.
5044 * Equivalent to Maple's irem(a,b) as far as sign conventions are concerned.
5045 * In general, mod(a,b) has the sign of b or is zero, and irem(a,b) has the
5046 * sign of a or is zero.
5047 *
5048 * @return remainder of a/b if both are integer, 0 otherwise.
5049 * @exception overflow_error (division by zero) if b is zero. */
irem(const numeric & a,const numeric & b)5050 const numeric irem(const numeric &a, const numeric &b) {
5051 return a.irem(b);
5052 }
5053
5054 /** Numeric integer quotient.
5055 * Equivalent to Maple's iquo as far as sign conventions are concerned.
5056 *
5057 * @return truncated quotient of a/b if both are integer, 0 otherwise.
5058 * @exception overflow_error (division by zero) if b is zero. */
iquo(const numeric & a,const numeric & b)5059 const numeric iquo(const numeric &a, const numeric &b) {
5060 return a.iquo(b);
5061 }
5062
5063 /** Numeric integer quotient.
5064 * Equivalent to Maple's iquo(a,b,'r') it obeyes the relation
5065 * r == a - iquo(a,b,r)*b.
5066 *
5067 * @return truncated quotient of a/b and remainder stored in r if both are
5068 * integer, 0 otherwise.
5069 * @exception overflow_error (division by zero) if b is zero. */
iquo(const numeric & a,const numeric & b,numeric & r)5070 const numeric iquo(const numeric &a, const numeric &b, numeric &r) {
5071 return a.iquo(b, r);
5072 }
5073
5074 /** Greatest Common Divisor.
5075 *
5076 * @return The GCD of two numbers if both are integer, a numerical 1
5077 * if they are not. */
gcd(const numeric & a,const numeric & b)5078 const numeric gcd(const numeric &a, const numeric &b) {
5079 return a.gcd(b);
5080 }
5081
5082 /** Least Common Multiple.
5083 *
5084 * @return The LCM of two numbers if both are integer, the product of those
5085 * two numbers if they are not. */
lcm(const numeric & a,const numeric & b)5086 const numeric lcm(const numeric &a, const numeric &b) {
5087 return a.lcm(b);
5088 }
5089
5090 /** Numeric square root.
5091 * If possible, sqrt(x) should respect squares of exact numbers, i.e. sqrt(4)
5092 * should return integer 2.
5093 *
5094 * @param x numeric argument
5095 * @return square root of x. Branch cut along negative real axis, the negative
5096 * real axis itself where imag(x)==0 and real(x)<0 belongs to the upper part
5097 * where imag(x)>0. */
sqrt(const numeric & x)5098 const numeric sqrt(const numeric &x) {
5099 return x.sqrt();
5100 }
5101
5102 /** Integer numeric square root. */
isqrt(const numeric & x)5103 const numeric isqrt(const numeric &x) {
5104 return x.isqrt();
5105 }
5106
5107 /** Floating point evaluation of Sage's constants. */
ConstantEvalf(unsigned serial,PyObject * dict)5108 ex ConstantEvalf(unsigned serial, PyObject* dict) {
5109 if (dict == nullptr) {
5110 dict = PyDict_New();
5111 PyDict_SetItemString(dict, "parent", CC_get());
5112 }
5113 PyObject* x = py_funcs.py_eval_constant(serial, dict);
5114 if (x == nullptr) py_error("error getting digits of constant");
5115 return x;
5116 }
5117
UnsignedInfinityEvalf(unsigned serial,PyObject * parent)5118 ex UnsignedInfinityEvalf(unsigned serial, PyObject* parent) {
5119 PyObject* x = py_funcs.py_eval_unsigned_infinity();
5120 return x;
5121 }
5122
InfinityEvalf(unsigned serial,PyObject * parent)5123 ex InfinityEvalf(unsigned serial, PyObject* parent) {
5124 PyObject* x = py_funcs.py_eval_infinity();
5125 return x;
5126 }
5127
NegInfinityEvalf(unsigned serial,PyObject * parent)5128 ex NegInfinityEvalf(unsigned serial, PyObject* parent) {
5129 PyObject* x = py_funcs.py_eval_neg_infinity();
5130 return x;
5131 }
5132
5133
5134 } // namespace GiNaC
5135