GiNaC 1.8.7
numeric.cpp
Go to the documentation of this file.
1
9/*
10 * GiNaC Copyright (C) 1999-2023 Johannes Gutenberg University Mainz, Germany
11 *
12 * This program is free software; you can redistribute it and/or modify
13 * it under the terms of the GNU General Public License as published by
14 * the Free Software Foundation; either version 2 of the License, or
15 * (at your option) any later version.
16 *
17 * This program is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 * GNU General Public License for more details.
21 *
22 * You should have received a copy of the GNU General Public License
23 * along with this program; if not, write to the Free Software
24 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
25 */
26
27#include "numeric.h"
28#include "ex.h"
29#include "operators.h"
30#include "archive.h"
31#include "utils.h"
32
33#include <limits>
34#include <sstream>
35#include <stdexcept>
36#include <string>
37#include <vector>
38
39// CLN should pollute the global namespace as little as possible. Hence, we
40// include most of it here and include only the part needed for properly
41// declaring cln::cl_number in numeric.h. This can only be safely done in
42// namespaced versions of CLN, i.e. version > 1.1.0. Also, we only need a
43// subset of CLN, so we don't include the complete <cln/cln.h> but only the
44// essential stuff:
45#include <cln/output.h>
46#include <cln/integer_io.h>
47#include <cln/integer_ring.h>
48#include <cln/rational_io.h>
49#include <cln/rational_ring.h>
50#include <cln/lfloat_class.h>
51#include <cln/lfloat_io.h>
52#include <cln/real_io.h>
53#include <cln/real_ring.h>
54#include <cln/complex_io.h>
55#include <cln/complex_ring.h>
56#include <cln/numtheory.h>
57
58namespace GiNaC {
59
62 print_func<print_latex>(&numeric::do_print_latex).
63 print_func<print_csrc>(&numeric::do_print_csrc).
64 print_func<print_csrc_cl_N>(&numeric::do_print_csrc_cl_N).
65 print_func<print_tree>(&numeric::do_print_tree).
66 print_func<print_python_repr>(&numeric::do_print_python_repr))
67
68
69// default constructor
71
74{
75 value = cln::cl_I(0);
77}
78
80// other constructors
82
83// public
84
86{
87 // Not the whole int-range is available if we don't cast to long
88 // first. This is due to the behavior of the cl_I-ctor, which
89 // emphasizes efficiency. However, if the integer is small enough
90 // we save space and dereferences by using an immediate type.
91 // (C.f. <cln/object.h>)
92 // The #if clause prevents compiler warnings on 64bit machines where the
93 // comparision is always true.
94#if cl_value_len >= 32
95 value = cln::cl_I(i);
96#else
97 if (i < (1L << (cl_value_len-1)) && i >= -(1L << (cl_value_len-1)))
98 value = cln::cl_I(i);
99 else
100 value = cln::cl_I(static_cast<long>(i));
101#endif
103}
104
105
106numeric::numeric(unsigned int i)
107{
108 // Not the whole uint-range is available if we don't cast to ulong
109 // first. This is due to the behavior of the cl_I-ctor, which
110 // emphasizes efficiency. However, if the integer is small enough
111 // we save space and dereferences by using an immediate type.
112 // (C.f. <cln/object.h>)
113 // The #if clause prevents compiler warnings on 64bit machines where the
114 // comparision is always true.
115#if cl_value_len >= 32
116 value = cln::cl_I(i);
117#else
118 if (i < (1UL << (cl_value_len-1)))
119 value = cln::cl_I(i);
120 else
121 value = cln::cl_I(static_cast<unsigned long>(i));
122#endif
124}
125
126
128{
129 value = cln::cl_I(i);
131}
132
133
134numeric::numeric(unsigned long i)
135{
136 value = cln::cl_I(i);
138}
139
141{
142 value = cln::cl_I(i);
144}
145
146numeric::numeric(unsigned long long i)
147{
148 value = cln::cl_I(i);
150}
151
156{
157 if (!denom)
158 throw std::overflow_error("division by zero");
159 value = cln::cl_I(numer) / cln::cl_I(denom);
161}
162
163
165{
166 // We really want to explicitly use the type cl_LF instead of the
167 // more general cl_F, since that would give us a cl_DF only which
168 // will not be promoted to cl_LF if overflow occurs:
169 value = cln::cl_float(d, cln::default_float_format);
171}
172
173
176numeric::numeric(const char *s)
177{
178 cln::cl_N ctorval = 0;
179 // parse complex numbers (functional but not completely safe, unfortunately
180 // std::string does not understand regexpese):
181 // ss should represent a simple sum like 2+5*I
182 std::string ss = s;
183 std::string::size_type delim;
184
185 // make this implementation safe by adding explicit sign
186 if (ss.at(0) != '+' && ss.at(0) != '-' && ss.at(0) != '#')
187 ss = '+' + ss;
188
189 // We use 'E' as exponent marker in the output, but some people insist on
190 // writing 'e' at input, so let's substitute them right at the beginning:
191 while ((delim = ss.find("e"))!=std::string::npos)
192 ss.replace(delim,1,"E");
193
194 // main parser loop:
195 do {
196 // chop ss into terms from left to right
197 std::string term;
198 bool imaginary = false;
199 delim = ss.find_first_of(std::string("+-"),1);
200 // Do we have an exponent marker like "31.415E-1"? If so, hop on!
201 if (delim!=std::string::npos && ss.at(delim-1)=='E')
202 delim = ss.find_first_of(std::string("+-"),delim+1);
203 term = ss.substr(0,delim);
204 if (delim!=std::string::npos)
205 ss = ss.substr(delim);
206 // is the term imaginary?
207 if (term.find("I")!=std::string::npos) {
208 // erase 'I':
209 term.erase(term.find("I"),1);
210 // erase '*':
211 if (term.find("*")!=std::string::npos)
212 term.erase(term.find("*"),1);
213 // correct for trivial +/-I without explicit factor on I:
214 if (term.size()==1)
215 term += '1';
216 imaginary = true;
217 }
218 if (term.find('.')!=std::string::npos || term.find('E')!=std::string::npos) {
219 // CLN's short type cl_SF is not very useful within the GiNaC
220 // framework where we are mainly interested in the arbitrary
221 // precision type cl_LF. Hence we go straight to the construction
222 // of generic floats. In order to create them we have to convert
223 // our own floating point notation used for output and construction
224 // from char * to CLN's generic notation:
225 // 3.14 --> 3.14e0_<Digits>
226 // 31.4E-1 --> 31.4e-1_<Digits>
227 // and s on.
228 // No exponent marker? Let's add a trivial one.
229 if (term.find("E")==std::string::npos)
230 term += "E0";
231 // E to lower case
232 term = term.replace(term.find("E"),1,"e");
233 // append _<Digits> to term
234 term += "_" + std::to_string((unsigned)Digits);
235 // construct float using cln::cl_F(const char *) ctor.
236 if (imaginary)
237 ctorval = ctorval + cln::complex(cln::cl_I(0),cln::cl_F(term.c_str()));
238 else
239 ctorval = ctorval + cln::cl_F(term.c_str());
240 } else {
241 // this is not a floating point number...
242 if (imaginary)
243 ctorval = ctorval + cln::complex(cln::cl_I(0),cln::cl_R(term.c_str()));
244 else
245 ctorval = ctorval + cln::cl_R(term.c_str());
246 }
247 } while (delim != std::string::npos);
248 value = ctorval;
250}
251
252
255numeric::numeric(const cln::cl_N &z)
256{
257 value = z;
259}
260
261
263// archiving
265
269static const cln::cl_F make_real_float(const cln::cl_idecoded_float& dec)
270{
271 cln::cl_F x = cln::cl_float(dec.mantissa, cln::default_float_format);
272 x = cln::scale_float(x, dec.exponent);
273 cln::cl_F sign = cln::cl_float(dec.sign, cln::default_float_format);
274 x = cln::float_sign(sign, x);
275 return x;
276}
277
281static const cln::cl_F read_real_float(std::istream& s)
282{
283 cln::cl_idecoded_float dec;
284 s >> dec.sign >> dec.mantissa >> dec.exponent;
285 const cln::cl_F x = make_real_float(dec);
286 return x;
287}
288
290{
291 inherited::read_archive(n, sym_lst);
292 value = 0;
293
294 // Read number as string
295 std::string str;
296 if (n.find_string("number", str)) {
297 std::istringstream s(str);
298 cln::cl_R re, im;
299 char c;
300 s.get(c);
301 switch (c) {
302 case 'R':
303 // real FP (floating point) number
304 re = read_real_float(s);
305 value = re;
306 break;
307 case 'C':
308 // both real and imaginary part are FP numbers
309 re = read_real_float(s);
310 im = read_real_float(s);
311 value = cln::complex(re, im);
312 break;
313 case 'H':
314 // real part is a rational number,
315 // imaginary part is a FP number
316 s >> re;
317 im = read_real_float(s);
318 value = cln::complex(re, im);
319 break;
320 case 'J':
321 // real part is a FP number,
322 // imaginary part is a rational number
323 re = read_real_float(s);
324 s >> im;
325 value = cln::complex(re, im);
326 break;
327 default:
328 // both real and imaginary parts are rational
329 s.putback(c);
330 s >> value;
331 break;
332 }
333 }
335}
337
338static void write_real_float(std::ostream& s, const cln::cl_R& n)
339{
340 const cln::cl_idecoded_float dec = cln::integer_decode_float(cln::the<cln::cl_F>(n));
341 s << dec.sign << ' ' << dec.mantissa << ' ' << dec.exponent;
342}
343
345{
346 inherited::archive(n);
347
348 // Write number as string
349
350 const cln::cl_R re = cln::realpart(value);
351 const cln::cl_R im = cln::imagpart(value);
352 const bool re_rationalp = cln::instanceof(re, cln::cl_RA_ring);
353 const bool im_rationalp = cln::instanceof(im, cln::cl_RA_ring);
354
355 // Non-rational numbers are written in an integer-decoded format
356 // to preserve the precision
357 std::ostringstream s;
358 if (re_rationalp && im_rationalp)
359 s << value;
360 else if (zerop(im)) {
361 // real FP (floating point) number
362 s << 'R';
363 write_real_float(s, re);
364 } else if (re_rationalp) {
365 s << 'H'; // just any unique character
366 // real part is a rational number,
367 // imaginary part is a FP number
368 s << re << ' ';
369 write_real_float(s, im);
370 } else if (im_rationalp) {
371 s << 'J';
372 // real part is a FP number,
373 // imaginary part is a rational number
374 write_real_float(s, re);
375 s << ' ' << im;
376 } else {
377 // both real and imaginary parts are floating point
378 s << 'C';
379 write_real_float(s, re);
380 s << ' ';
381 write_real_float(s, im);
382 }
383 n.add_string("number", s.str());
384}
385
387// functions overriding virtual functions from base classes
389
397static void print_real_number(const print_context & c, const cln::cl_R & x)
398{
399 cln::cl_print_flags ourflags;
400 if (cln::instanceof(x, cln::cl_RA_ring)) {
401 // case 1: integer or rational
402 if (cln::instanceof(x, cln::cl_I_ring) ||
403 !is_a<print_latex>(c)) {
404 cln::print_real(c.s, ourflags, x);
405 } else { // rational output in LaTeX context
406 if (x < 0)
407 c.s << "-";
408 c.s << "\\frac{";
409 cln::print_real(c.s, ourflags, cln::abs(cln::numerator(cln::the<cln::cl_RA>(x))));
410 c.s << "}{";
411 cln::print_real(c.s, ourflags, cln::denominator(cln::the<cln::cl_RA>(x)));
412 c.s << '}';
413 }
414 } else {
415 // case 2: float
416 // make CLN believe this number has default_float_format, so it prints
417 // 'E' as exponent marker instead of 'L':
418 ourflags.default_float_format = cln::float_format(cln::the<cln::cl_F>(x));
419 cln::print_real(c.s, ourflags, x);
420 }
421}
422
426static void print_integer_csrc(const print_context & c, const cln::cl_I & x)
427{
428 // Print small numbers in compact float format, but larger numbers in
429 // scientific format
430 const int max_cln_int = 536870911; // 2^29-1
431 if (x >= cln::cl_I(-max_cln_int) && x <= cln::cl_I(max_cln_int))
432 c.s << cln::cl_I_to_int(x) << ".0";
433 else
434 c.s << cln::double_approx(x);
435}
436
440static void print_real_csrc(const print_context & c, const cln::cl_R & x)
441{
442 if (cln::instanceof(x, cln::cl_I_ring)) {
443
444 // Integer number
445 print_integer_csrc(c, cln::the<cln::cl_I>(x));
446
447 } else if (cln::instanceof(x, cln::cl_RA_ring)) {
448
449 // Rational number
450 const cln::cl_I numer = cln::numerator(cln::the<cln::cl_RA>(x));
451 const cln::cl_I denom = cln::denominator(cln::the<cln::cl_RA>(x));
452 if (cln::plusp(x)) {
453 c.s << "(";
455 } else {
456 c.s << "-(";
458 }
459 c.s << "/";
461 c.s << ")";
462
463 } else {
464
465 // Anything else
466 c.s << cln::double_approx(x);
467 }
468}
469
470template<typename T1, typename T2>
471static inline bool coerce(T1& dst, const T2& arg);
472
478template<>
479inline bool coerce<int, cln::cl_I>(int& dst, const cln::cl_I& arg)
480{
481 static const cln::cl_I cl_max_int =
482 (cln::cl_I)(long)(std::numeric_limits<int>::max());
483 static const cln::cl_I cl_min_int =
484 (cln::cl_I)(long)(std::numeric_limits<int>::min());
485 if ((arg >= cl_min_int) && (arg <= cl_max_int)) {
486 dst = cl_I_to_int(arg);
487 return true;
488 }
489 return false;
490}
491
492template<>
493inline bool coerce<unsigned int, cln::cl_I>(unsigned int& dst, const cln::cl_I& arg)
494{
495 static const cln::cl_I cl_max_uint =
496 (cln::cl_I)(unsigned long)(std::numeric_limits<unsigned int>::max());
497 if ((! minusp(arg)) && (arg <= cl_max_uint)) {
498 dst = cl_I_to_uint(arg);
499 return true;
500 }
501 return false;
502}
503
507static void print_real_cl_N(const print_context & c, const cln::cl_R & x)
508{
509 if (cln::instanceof(x, cln::cl_I_ring)) {
510
511 int dst;
512 // fixnum
513 if (coerce(dst, cln::the<cln::cl_I>(x))) {
514 // can be converted to native int
515 if (dst < 0)
516 c.s << '(' << dst << ')';
517 else
518 c.s << dst;
519 } else {
520 // bignum
521 c.s << "cln::cl_I(\"";
523 c.s << "\")";
524 }
525 } else if (cln::instanceof(x, cln::cl_RA_ring)) {
526
527 // Rational number
528 cln::cl_print_flags ourflags;
529 c.s << "cln::cl_RA(\"";
530 cln::print_rational(c.s, ourflags, cln::the<cln::cl_RA>(x));
531 c.s << "\")";
532
533 } else {
534
535 // Anything else
536 c.s << "cln::cl_F(\"";
537 print_real_number(c, cln::cl_float(1.0, cln::default_float_format) * x);
538 c.s << "_" << Digits << "\")";
539 }
540}
541
542void numeric::print_numeric(const print_context & c, const char *par_open, const char *par_close, const char *imag_sym, const char *mul_sym, unsigned level) const
543{
544 const cln::cl_R r = cln::realpart(value);
545 const cln::cl_R i = cln::imagpart(value);
546
547 if (cln::zerop(i)) {
548
549 // case 1, real: x or -x
550 if ((precedence() <= level) && (!this->is_nonneg_integer())) {
551 c.s << par_open;
553 c.s << par_close;
554 } else {
556 }
557
558 } else {
559 if (cln::zerop(r)) {
560
561 // case 2, imaginary: y*I or -y*I
562 if (i == 1)
563 c.s << imag_sym;
564 else {
565 if (precedence()<=level)
566 c.s << par_open;
567 if (i == -1)
568 c.s << "-" << imag_sym;
569 else {
571 c.s << mul_sym << imag_sym;
572 }
573 if (precedence()<=level)
574 c.s << par_close;
575 }
576
577 } else {
578
579 // case 3, complex: x+y*I or x-y*I or -x+y*I or -x-y*I
580 if (precedence() <= level)
581 c.s << par_open;
583 if (i < 0) {
584 if (i == -1) {
585 c.s << "-" << imag_sym;
586 } else {
588 c.s << mul_sym << imag_sym;
589 }
590 } else {
591 if (i == 1) {
592 c.s << "+" << imag_sym;
593 } else {
594 c.s << "+";
596 c.s << mul_sym << imag_sym;
597 }
598 }
599 if (precedence() <= level)
600 c.s << par_close;
601 }
602 }
603}
604
605void numeric::do_print(const print_context & c, unsigned level) const
606{
607 print_numeric(c, "(", ")", "I", "*", level);
608}
609
610void numeric::do_print_latex(const print_latex & c, unsigned level) const
611{
612 print_numeric(c, "{(", ")}", "i", " ", level);
613}
614
615void numeric::do_print_csrc(const print_csrc & c, unsigned level) const
616{
617 std::ios::fmtflags oldflags = c.s.flags();
618 c.s.setf(std::ios::scientific);
619 int oldprec = c.s.precision();
620
621 // Set precision
622 if (is_a<print_csrc_double>(c))
623 c.s.precision(std::numeric_limits<double>::digits10 + 1);
624 else
625 c.s.precision(std::numeric_limits<float>::digits10 + 1);
626
627 if (this->is_real()) {
628
629 // Real number
630 print_real_csrc(c, cln::the<cln::cl_R>(value));
631
632 } else {
633
634 // Complex number
635 c.s << "std::complex<";
636 if (is_a<print_csrc_double>(c))
637 c.s << "double>(";
638 else
639 c.s << "float>(";
640
641 print_real_csrc(c, cln::realpart(value));
642 c.s << ",";
643 print_real_csrc(c, cln::imagpart(value));
644 c.s << ")";
645 }
646
647 c.s.flags(oldflags);
648 c.s.precision(oldprec);
649}
650
651void numeric::do_print_csrc_cl_N(const print_csrc_cl_N & c, unsigned level) const
652{
653 if (this->is_real()) {
654
655 // Real number
656 print_real_cl_N(c, cln::the<cln::cl_R>(value));
657
658 } else {
659
660 // Complex number
661 c.s << "cln::complex(";
662 print_real_cl_N(c, cln::realpart(value));
663 c.s << ",";
664 print_real_cl_N(c, cln::imagpart(value));
665 c.s << ")";
666 }
667}
668
669void numeric::do_print_tree(const print_tree & c, unsigned level) const
670{
671 c.s << std::string(level, ' ') << value
672 << " (" << class_name() << ")" << " @" << this
673 << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
674 << std::endl;
675}
676
677void numeric::do_print_python_repr(const print_python_repr & c, unsigned level) const
678{
679 c.s << class_name() << "('";
680 print_numeric(c, "(", ")", "I", "*", level);
681 c.s << "')";
682}
683
684bool numeric::info(unsigned inf) const
685{
686 switch (inf) {
691 return true;
692 case info_flags::real:
693 return is_real();
696 return is_rational();
699 return is_crational();
702 return is_integer();
705 return is_cinteger();
707 return is_positive();
709 return is_negative();
711 return is_zero() || is_positive();
713 return is_pos_integer();
715 return is_integer() && is_negative();
717 return is_nonneg_integer();
718 case info_flags::even:
719 return is_even();
720 case info_flags::odd:
721 return is_odd();
723 return is_prime();
724 }
725 return false;
726}
727
728bool numeric::is_polynomial(const ex & var) const
729{
730 return true;
731}
732
733int numeric::degree(const ex & s) const
734{
735 return 0;
736}
737
738int numeric::ldegree(const ex & s) const
739{
740 return 0;
741}
742
743ex numeric::coeff(const ex & s, int n) const
744{
745 return n==0 ? *this : _ex0;
746}
747
754bool numeric::has(const ex &other, unsigned options) const
755{
756 if (!is_exactly_a<numeric>(other))
757 return false;
758 const numeric &o = ex_to<numeric>(other);
759 if (this->is_equal(o) || this->is_equal(-o))
760 return true;
761 if (o.imag().is_zero()) { // e.g. scan for 3 in -3*I
762 if (!this->real().is_equal(*_num0_p))
763 if (this->real().is_equal(o) || this->real().is_equal(-o))
764 return true;
765 if (!this->imag().is_equal(*_num0_p))
766 if (this->imag().is_equal(o) || this->imag().is_equal(-o))
767 return true;
768 return false;
769 }
770 else {
771 if (o.is_equal(I)) // e.g scan for I in 42*I
772 return !this->is_real();
773 if (o.real().is_zero()) // e.g. scan for 2*I in 2*I+1
774 if (!this->imag().is_equal(*_num0_p))
775 if (this->imag().is_equal(o*I) || this->imag().is_equal(-o*I))
776 return true;
777 }
778 return false;
779}
780
781
784{
785 return this->hold();
786}
787
788
796{
797 return numeric(cln::cl_float(1.0, cln::default_float_format) * value);
798}
799
801{
802 if (is_real()) {
803 return *this;
804 }
805 return numeric(cln::conjugate(this->value));
806}
807
809{
810 return numeric(cln::realpart(value));
811}
812
814{
815 return numeric(cln::imagpart(value));
816}
817
818// protected
819
820int numeric::compare_same_type(const basic &other) const
821{
822 GINAC_ASSERT(is_exactly_a<numeric>(other));
823 const numeric &o = static_cast<const numeric &>(other);
824
825 return this->compare(o);
826}
827
828
829bool numeric::is_equal_same_type(const basic &other) const
830{
831 GINAC_ASSERT(is_exactly_a<numeric>(other));
832 const numeric &o = static_cast<const numeric &>(other);
833
834 return this->is_equal(o);
835}
836
837
838unsigned numeric::calchash() const
839{
840 // Base computation of hashvalue on CLN's hashcode. Note: That depends
841 // only on the number's value, not its type or precision (i.e. a true
842 // equivalence relation on numbers). As a consequence, 3 and 3.0 share
843 // the same hashvalue. That shouldn't really matter, though.
845 hashvalue = golden_ratio_hash(cln::equal_hashcode(value));
846 return hashvalue;
847}
848
849
851// new virtual functions which can be overridden by derived classes
853
854// none
855
857// non-virtual functions in this class
859
860// public
861
864const numeric numeric::add(const numeric &other) const
865{
866 return numeric(value + other.value);
867}
868
869
872const numeric numeric::sub(const numeric &other) const
873{
874 return numeric(value - other.value);
875}
876
877
880const numeric numeric::mul(const numeric &other) const
881{
882 return numeric(value * other.value);
883}
884
885
890const numeric numeric::div(const numeric &other) const
891{
892 if (cln::zerop(other.value))
893 throw std::overflow_error("numeric::div(): division by zero");
894 return numeric(value / other.value);
895}
896
897
900const numeric numeric::power(const numeric &other) const
901{
902 // Shortcut for efficiency and numeric stability (as in 1.0 exponent):
903 // trap the neutral exponent.
904 if (&other==_num1_p || cln::equal(other.value,_num1_p->value))
905 return *this;
906
907 if (cln::zerop(value)) {
908 if (cln::zerop(other.value))
909 throw std::domain_error("numeric::eval(): pow(0,0) is undefined");
910 else if (cln::zerop(cln::realpart(other.value)))
911 throw std::domain_error("numeric::eval(): pow(0,I) is undefined");
912 else if (cln::minusp(cln::realpart(other.value)))
913 throw std::overflow_error("numeric::eval(): division by zero");
914 else
915 return *_num0_p;
916 }
917 return numeric(cln::expt(value, other.value));
918}
919
920
921
925const numeric &numeric::add_dyn(const numeric &other) const
926{
927 // Efficiency shortcut: trap the neutral element by pointer. This hack
928 // is supposed to keep the number of distinct numeric objects low.
929 if (this==_num0_p)
930 return other;
931 else if (&other==_num0_p)
932 return *this;
933
934 return dynallocate<numeric>(value + other.value);
935}
936
937
942const numeric &numeric::sub_dyn(const numeric &other) const
943{
944 // Efficiency shortcut: trap the neutral exponent (first by pointer). This
945 // hack is supposed to keep the number of distinct numeric objects low.
946 if (&other==_num0_p || cln::zerop(other.value))
947 return *this;
948
949 return dynallocate<numeric>(value - other.value);
950}
951
952
957const numeric &numeric::mul_dyn(const numeric &other) const
958{
959 // Efficiency shortcut: trap the neutral element by pointer. This hack
960 // is supposed to keep the number of distinct numeric objects low.
961 if (this==_num1_p)
962 return other;
963 else if (&other==_num1_p)
964 return *this;
965
966 return dynallocate<numeric>(value * other.value);
967}
968
969
976const numeric &numeric::div_dyn(const numeric &other) const
977{
978 // Efficiency shortcut: trap the neutral element by pointer. This hack
979 // is supposed to keep the number of distinct numeric objects low.
980 if (&other==_num1_p)
981 return *this;
982 if (cln::zerop(cln::the<cln::cl_N>(other.value)))
983 throw std::overflow_error("division by zero");
984
985 return dynallocate<numeric>(value / other.value);
986}
987
988
993const numeric &numeric::power_dyn(const numeric &other) const
994{
995 // Efficiency shortcut: trap the neutral exponent (first try by pointer, then
996 // try harder, since calls to cln::expt() below may return amazing results for
997 // floating point exponent 1.0).
998 if (&other==_num1_p || cln::equal(other.value, _num1_p->value))
999 return *this;
1000
1001 if (cln::zerop(value)) {
1002 if (cln::zerop(other.value))
1003 throw std::domain_error("numeric::eval(): pow(0,0) is undefined");
1004 else if (cln::zerop(cln::realpart(other.value)))
1005 throw std::domain_error("numeric::eval(): pow(0,I) is undefined");
1006 else if (cln::minusp(cln::realpart(other.value)))
1007 throw std::overflow_error("numeric::eval(): division by zero");
1008 else
1009 return *_num0_p;
1010 }
1011
1012 return dynallocate<numeric>(cln::expt(value, other.value));
1013}
1014
1015
1017{
1018 return operator=(numeric(i));
1019}
1020
1021
1022const numeric &numeric::operator=(unsigned int i)
1023{
1024 return operator=(numeric(i));
1025}
1026
1027
1029{
1030 return operator=(numeric(i));
1031}
1032
1033
1034const numeric &numeric::operator=(unsigned long i)
1035{
1036 return operator=(numeric(i));
1037}
1038
1039
1041{
1042 return operator=(numeric(d));
1043}
1044
1045
1046const numeric &numeric::operator=(const char * s)
1047{
1048 return operator=(numeric(s));
1049}
1050
1051
1054{
1055 if (cln::zerop(value))
1056 throw std::overflow_error("numeric::inverse(): division by zero");
1057 return numeric(cln::recip(value));
1058}
1059
1065{ cln::cl_R r = cln::realpart(value);
1066 if(cln::zerop(r))
1067 return numeric(1,2);
1068 if(cln::plusp(r))
1069 return 1;
1070 return 0;
1071}
1072
1078int numeric::csgn() const
1079{
1080 if (cln::zerop(value))
1081 return 0;
1082 cln::cl_R r = cln::realpart(value);
1083 if (!cln::zerop(r)) {
1084 if (cln::plusp(r))
1085 return 1;
1086 else
1087 return -1;
1088 } else {
1089 if (cln::plusp(cln::imagpart(value)))
1090 return 1;
1091 else
1092 return -1;
1093 }
1094}
1095
1096
1104int numeric::compare(const numeric &other) const
1105{
1106 // Comparing two real numbers?
1107 if (cln::instanceof(value, cln::cl_R_ring) &&
1108 cln::instanceof(other.value, cln::cl_R_ring))
1109 // Yes, so just cln::compare them
1110 return cln::compare(cln::the<cln::cl_R>(value), cln::the<cln::cl_R>(other.value));
1111 else {
1112 // No, first cln::compare real parts...
1113 cl_signean real_cmp = cln::compare(cln::realpart(value), cln::realpart(other.value));
1114 if (real_cmp)
1115 return real_cmp;
1116 // ...and then the imaginary parts.
1117 return cln::compare(cln::imagpart(value), cln::imagpart(other.value));
1118 }
1119}
1120
1121
1122bool numeric::is_equal(const numeric &other) const
1123{
1124 return cln::equal(value, other.value);
1125}
1126
1127
1130{
1131 return cln::zerop(value);
1132}
1133
1134
1137{
1138 if (cln::instanceof(value, cln::cl_R_ring)) // real?
1139 return cln::plusp(cln::the<cln::cl_R>(value));
1140 return false;
1141}
1142
1143
1146{
1147 if (cln::instanceof(value, cln::cl_R_ring)) // real?
1148 return cln::minusp(cln::the<cln::cl_R>(value));
1149 return false;
1150}
1151
1152
1155{
1156 return cln::instanceof(value, cln::cl_I_ring);
1157}
1158
1159
1162{
1163 return (cln::instanceof(value, cln::cl_I_ring) && cln::plusp(cln::the<cln::cl_I>(value)));
1164}
1165
1166
1169{
1170 return (cln::instanceof(value, cln::cl_I_ring) && !cln::minusp(cln::the<cln::cl_I>(value)));
1171}
1172
1173
1176{
1177 return (cln::instanceof(value, cln::cl_I_ring) && cln::evenp(cln::the<cln::cl_I>(value)));
1178}
1179
1180
1183{
1184 return (cln::instanceof(value, cln::cl_I_ring) && cln::oddp(cln::the<cln::cl_I>(value)));
1185}
1186
1187
1192{
1193 return (cln::instanceof(value, cln::cl_I_ring) // integer?
1194 && cln::plusp(cln::the<cln::cl_I>(value)) // positive?
1195 && cln::isprobprime(cln::the<cln::cl_I>(value)));
1196}
1197
1198
1202{
1203 return cln::instanceof(value, cln::cl_RA_ring);
1204}
1205
1206
1209{
1210 return cln::instanceof(value, cln::cl_R_ring);
1211}
1212
1213
1214bool numeric::operator==(const numeric &other) const
1215{
1216 return cln::equal(value, other.value);
1217}
1218
1219
1220bool numeric::operator!=(const numeric &other) const
1221{
1222 return !cln::equal(value, other.value);
1223}
1224
1225
1229{
1230 if (cln::instanceof(value, cln::cl_I_ring))
1231 return true;
1232 else if (!this->is_real()) { // complex case, handle n+m*I
1233 if (cln::instanceof(cln::realpart(value), cln::cl_I_ring) &&
1234 cln::instanceof(cln::imagpart(value), cln::cl_I_ring))
1235 return true;
1236 }
1237 return false;
1238}
1239
1240
1244{
1245 if (cln::instanceof(value, cln::cl_RA_ring))
1246 return true;
1247 else if (!this->is_real()) { // complex case, handle Q(i):
1248 if (cln::instanceof(cln::realpart(value), cln::cl_RA_ring) &&
1249 cln::instanceof(cln::imagpart(value), cln::cl_RA_ring))
1250 return true;
1251 }
1252 return false;
1253}
1254
1255
1259bool numeric::operator<(const numeric &other) const
1260{
1261 if (this->is_real() && other.is_real())
1262 return (cln::the<cln::cl_R>(value) < cln::the<cln::cl_R>(other.value));
1263 throw std::invalid_argument("numeric::operator<(): complex inequality");
1264}
1265
1266
1270bool numeric::operator<=(const numeric &other) const
1271{
1272 if (this->is_real() && other.is_real())
1273 return (cln::the<cln::cl_R>(value) <= cln::the<cln::cl_R>(other.value));
1274 throw std::invalid_argument("numeric::operator<=(): complex inequality");
1275}
1276
1277
1281bool numeric::operator>(const numeric &other) const
1282{
1283 if (this->is_real() && other.is_real())
1284 return (cln::the<cln::cl_R>(value) > cln::the<cln::cl_R>(other.value));
1285 throw std::invalid_argument("numeric::operator>(): complex inequality");
1286}
1287
1288
1292bool numeric::operator>=(const numeric &other) const
1293{
1294 if (this->is_real() && other.is_real())
1295 return (cln::the<cln::cl_R>(value) >= cln::the<cln::cl_R>(other.value));
1296 throw std::invalid_argument("numeric::operator>=(): complex inequality");
1297}
1298
1299
1304{
1305 GINAC_ASSERT(this->is_integer());
1306 return cln::cl_I_to_int(cln::the<cln::cl_I>(value));
1307}
1308
1309
1314{
1315 GINAC_ASSERT(this->is_integer());
1316 return cln::cl_I_to_long(cln::the<cln::cl_I>(value));
1317}
1318
1319
1323{
1324 GINAC_ASSERT(this->is_real());
1325 return cln::double_approx(cln::realpart(value));
1326}
1327
1328
1332cln::cl_N numeric::to_cl_N() const
1333{
1334 return value;
1335}
1336
1337
1340{
1341 return numeric(cln::realpart(value));
1342}
1343
1344
1347{
1348 return numeric(cln::imagpart(value));
1349}
1350
1351
1357{
1358 if (cln::instanceof(value, cln::cl_I_ring))
1359 return numeric(*this); // integer case
1360
1361 else if (cln::instanceof(value, cln::cl_RA_ring))
1362 return numeric(cln::numerator(cln::the<cln::cl_RA>(value)));
1363
1364 else if (!this->is_real()) { // complex case, handle Q(i):
1365 const cln::cl_RA r = cln::the<cln::cl_RA>(cln::realpart(value));
1366 const cln::cl_RA i = cln::the<cln::cl_RA>(cln::imagpart(value));
1367 if (cln::instanceof(r, cln::cl_I_ring) && cln::instanceof(i, cln::cl_I_ring))
1368 return numeric(*this);
1369 if (cln::instanceof(r, cln::cl_I_ring) && cln::instanceof(i, cln::cl_RA_ring))
1370 return numeric(cln::complex(r*cln::denominator(i), cln::numerator(i)));
1371 if (cln::instanceof(r, cln::cl_RA_ring) && cln::instanceof(i, cln::cl_I_ring))
1372 return numeric(cln::complex(cln::numerator(r), i*cln::denominator(r)));
1373 if (cln::instanceof(r, cln::cl_RA_ring) && cln::instanceof(i, cln::cl_RA_ring)) {
1374 const cln::cl_I s = cln::lcm(cln::denominator(r), cln::denominator(i));
1375 return numeric(cln::complex(cln::numerator(r)*(cln::exquo(s,cln::denominator(r))),
1376 cln::numerator(i)*(cln::exquo(s,cln::denominator(i)))));
1377 }
1378 }
1379 // at least one float encountered
1380 return numeric(*this);
1381}
1382
1383
1388{
1389 if (cln::instanceof(value, cln::cl_I_ring))
1390 return *_num1_p; // integer case
1391
1392 if (cln::instanceof(value, cln::cl_RA_ring))
1393 return numeric(cln::denominator(cln::the<cln::cl_RA>(value)));
1394
1395 if (!this->is_real()) { // complex case, handle Q(i):
1396 const cln::cl_RA r = cln::the<cln::cl_RA>(cln::realpart(value));
1397 const cln::cl_RA i = cln::the<cln::cl_RA>(cln::imagpart(value));
1398 if (cln::instanceof(r, cln::cl_I_ring) && cln::instanceof(i, cln::cl_I_ring))
1399 return *_num1_p;
1400 if (cln::instanceof(r, cln::cl_I_ring) && cln::instanceof(i, cln::cl_RA_ring))
1401 return numeric(cln::denominator(i));
1402 if (cln::instanceof(r, cln::cl_RA_ring) && cln::instanceof(i, cln::cl_I_ring))
1403 return numeric(cln::denominator(r));
1404 if (cln::instanceof(r, cln::cl_RA_ring) && cln::instanceof(i, cln::cl_RA_ring))
1405 return numeric(cln::lcm(cln::denominator(r), cln::denominator(i)));
1406 }
1407 // at least one float encountered
1408 return *_num1_p;
1409}
1410
1411
1419{
1420 if (cln::instanceof(value, cln::cl_I_ring))
1421 return cln::integer_length(cln::the<cln::cl_I>(value));
1422 else
1423 return 0;
1424}
1425
1427// global constants
1429
1433const numeric I = numeric(cln::complex(cln::cl_I(0),cln::cl_I(1)));
1434
1435
1439const numeric exp(const numeric &x)
1440{
1441 return numeric(cln::exp(x.to_cl_N()));
1442}
1443
1444
1450const numeric log(const numeric &x)
1451{
1452 if (x.is_zero())
1453 throw pole_error("log(): logarithmic pole",0);
1454 return numeric(cln::log(x.to_cl_N()));
1455}
1456
1457
1461const numeric sin(const numeric &x)
1462{
1463 return numeric(cln::sin(x.to_cl_N()));
1464}
1465
1466
1470const numeric cos(const numeric &x)
1471{
1472 return numeric(cln::cos(x.to_cl_N()));
1473}
1474
1475
1479const numeric tan(const numeric &x)
1480{
1481 return numeric(cln::tan(x.to_cl_N()));
1482}
1483
1484
1488const numeric asin(const numeric &x)
1489{
1490 return numeric(cln::asin(x.to_cl_N()));
1491}
1492
1493
1497const numeric acos(const numeric &x)
1498{
1499 return numeric(cln::acos(x.to_cl_N()));
1500}
1501
1502
1508const numeric atan(const numeric &x)
1509{
1510 if (!x.is_real() &&
1511 x.real().is_zero() &&
1512 abs(x.imag()).is_equal(*_num1_p))
1513 throw pole_error("atan(): logarithmic pole",0);
1514 return numeric(cln::atan(x.to_cl_N()));
1515}
1516
1517
1525const numeric atan(const numeric &y, const numeric &x)
1526{
1527 if (x.is_zero() && y.is_zero())
1528 return *_num0_p;
1529 if (x.is_real() && y.is_real())
1530 return numeric(cln::atan(cln::the<cln::cl_R>(x.to_cl_N()),
1531 cln::the<cln::cl_R>(y.to_cl_N())));
1532
1533 // Compute -I*log((x+I*y)/sqrt(x^2+y^2))
1534 // == -I*log((x+I*y)/sqrt((x+I*y)*(x-I*y)))
1535 // Do not "simplify" this to -I/2*log((x+I*y)/(x-I*y))) or likewise.
1536 // The branch cuts are easily messed up.
1537 const cln::cl_N aux_p = x.to_cl_N()+cln::complex(0,1)*y.to_cl_N();
1538 if (cln::zerop(aux_p)) {
1539 // x+I*y==0 => y/x==I, so this is a pole (we have x!=0).
1540 throw pole_error("atan(): logarithmic pole",0);
1541 }
1542 const cln::cl_N aux_m = x.to_cl_N()-cln::complex(0,1)*y.to_cl_N();
1543 if (cln::zerop(aux_m)) {
1544 // x-I*y==0 => y/x==-I, so this is a pole (we have x!=0).
1545 throw pole_error("atan(): logarithmic pole",0);
1546 }
1547 return numeric(cln::complex(0,-1)*cln::log(aux_p/cln::sqrt(aux_p*aux_m)));
1548}
1549
1550
1554const numeric sinh(const numeric &x)
1555{
1556 return numeric(cln::sinh(x.to_cl_N()));
1557}
1558
1559
1563const numeric cosh(const numeric &x)
1564{
1565 return numeric(cln::cosh(x.to_cl_N()));
1566}
1567
1568
1572const numeric tanh(const numeric &x)
1573{
1574 return numeric(cln::tanh(x.to_cl_N()));
1575}
1576
1577
1581const numeric asinh(const numeric &x)
1582{
1583 return numeric(cln::asinh(x.to_cl_N()));
1584}
1585
1586
1590const numeric acosh(const numeric &x)
1591{
1592 return numeric(cln::acosh(x.to_cl_N()));
1593}
1594
1595
1599const numeric atanh(const numeric &x)
1600{
1601 return numeric(cln::atanh(x.to_cl_N()));
1602}
1603
1604
1605/*static cln::cl_N Li2_series(const ::cl_N &x,
1606 const ::float_format_t &prec)
1607{
1608 // Note: argument must be in the unit circle
1609 // This is very inefficient unless we have fast floating point Bernoulli
1610 // numbers implemented!
1611 cln::cl_N c1 = -cln::log(1-x);
1612 cln::cl_N c2 = c1;
1613 // hard-wire the first two Bernoulli numbers
1614 cln::cl_N acc = c1 - cln::square(c1)/4;
1615 cln::cl_N aug;
1616 cln::cl_F pisq = cln::square(cln::cl_pi(prec)); // pi^2
1617 cln::cl_F piac = cln::cl_float(1, prec); // accumulator: pi^(2*i)
1618 unsigned i = 1;
1619 c1 = cln::square(c1);
1620 do {
1621 c2 = c1 * c2;
1622 piac = piac * pisq;
1623 aug = c2 * (*(bernoulli(numeric(2*i)).clnptr())) / cln::factorial(2*i+1);
1624 // aug = c2 * cln::cl_I(i%2 ? 1 : -1) / cln::cl_I(2*i+1) * cln::cl_zeta(2*i, prec) / piac / (cln::cl_I(1)<<(2*i-1));
1625 acc = acc + aug;
1626 ++i;
1627 } while (acc != acc+aug);
1628 return acc;
1629}*/
1630
1633static cln::cl_N Li2_series(const cln::cl_N &x,
1634 const cln::float_format_t &prec)
1635{
1636 // Note: argument must be in the unit circle
1637 cln::cl_N aug, acc;
1638 cln::cl_N num = cln::complex(cln::cl_float(1, prec), 0);
1639 cln::cl_I den = 0;
1640 unsigned i = 1;
1641 do {
1642 num = num * x;
1643 den = den + i; // 1, 4, 9, 16, ...
1644 i += 2;
1645 aug = num / den;
1646 acc = acc + aug;
1647 } while (acc != acc+aug);
1648 return acc;
1649}
1650
1652static cln::cl_N Li2_projection(const cln::cl_N &x,
1653 const cln::float_format_t &prec)
1654{
1655 const cln::cl_R re = cln::realpart(x);
1656 const cln::cl_R im = cln::imagpart(x);
1657 if (re > cln::cl_F(".5"))
1658 // zeta(2) - Li2(1-x) - log(x)*log(1-x)
1659 return(cln::zeta(2)
1660 - Li2_series(1-x, prec)
1661 - cln::log(x)*cln::log(1-x));
1662 if ((re <= 0 && cln::abs(im) > cln::cl_F(".75")) || (re < cln::cl_F("-.5")))
1663 // -log(1-x)^2 / 2 - Li2(x/(x-1))
1664 return(- cln::square(cln::log(1-x))/2
1665 - Li2_series(x/(x-1), prec));
1666 if (re > 0 && cln::abs(im) > cln::cl_LF(".75"))
1667 // Li2(x^2)/2 - Li2(-x)
1668 return(Li2_projection(cln::square(x), prec)/2
1669 - Li2_projection(-x, prec));
1670 return Li2_series(x, prec);
1671}
1672
1673
1679const cln::cl_N Li2_(const cln::cl_N& value)
1680{
1681 if (zerop(value))
1682 return 0;
1683
1684 // what is the desired float format?
1685 // first guess: default format
1686 cln::float_format_t prec = cln::default_float_format;
1687 // second guess: the argument's format
1688 if (!instanceof(realpart(value), cln::cl_RA_ring))
1689 prec = cln::float_format(cln::the<cln::cl_F>(cln::realpart(value)));
1690 else if (!instanceof(imagpart(value), cln::cl_RA_ring))
1691 prec = cln::float_format(cln::the<cln::cl_F>(cln::imagpart(value)));
1692
1693 if (value==1) // may cause trouble with log(1-x)
1694 return cln::zeta(2, prec);
1695
1696 if (cln::abs(value) > 1)
1697 // -log(-x)^2 / 2 - zeta(2) - Li2(1/x)
1698 return(- cln::square(cln::log(-value))/2
1699 - cln::zeta(2, prec)
1700 - Li2_projection(cln::recip(value), prec));
1701 else
1702 return Li2_projection(value, prec);
1703}
1704
1705const numeric Li2(const numeric &x)
1706{
1707 const cln::cl_N x_ = x.to_cl_N();
1708 if (zerop(x_))
1709 return *_num0_p;
1710 const cln::cl_N result = Li2_(x_);
1711 return numeric(result);
1712}
1713
1714
1717const numeric zeta(const numeric &x)
1718{
1719 // A dirty hack to allow for things like zeta(3.0), since CLN currently
1720 // only knows about integer arguments and zeta(3).evalf() automatically
1721 // cascades down to zeta(3.0).evalf(). The trick is to rely on 3.0-3
1722 // being an exact zero for CLN, which can be tested and then we can just
1723 // pass the number casted to an int:
1724 if (x.is_real()) {
1725 const int aux = (int)(cln::double_approx(cln::the<cln::cl_R>(x.to_cl_N())));
1726 if (cln::zerop(x.to_cl_N()-aux))
1727 return numeric(cln::zeta(aux));
1728 }
1729 throw dunno();
1730}
1731
1733{
1734 public:
1736 bool sufficiently_accurate(int digits);
1737 int get_order() const { return current_vector->size(); }
1738 cln::cl_N calc_lanczos_A(const cln::cl_N &) const;
1739 private:
1740 // coeffs[0] is used in case Digits <= 20.
1741 // coeffs[1] is used in case Digits <= 50.
1742 // coeffs[2] is used in case Digits <= 100.
1743 // coeffs[3] is used in case Digits <= 200.
1744 static std::vector<cln::cl_N> *coeffs;
1745 // Pointer to the vector that is currently in use.
1746 std::vector<cln::cl_N> *current_vector;
1747};
1748
1749std::vector<cln::cl_N>* lanczos_coeffs::coeffs = nullptr;
1750
1752{ if (digits<=20) {
1753 current_vector = &(coeffs[0]);
1754 return true;
1755 }
1756 if (digits<=50) {
1757 current_vector = &(coeffs[1]);
1758 return true;
1759 }
1760 if (digits<=100) {
1761 current_vector = &(coeffs[2]);
1762 return true;
1763 }
1764 if (digits<=200) {
1765 current_vector = &(coeffs[3]);
1766 return true;
1767 }
1768 return false;
1769}
1770
1771cln::cl_N lanczos_coeffs::calc_lanczos_A(const cln::cl_N &x) const
1772{
1773 cln::cl_N A = (*current_vector)[0];
1774 int size = current_vector->size();
1775 for (int i=1; i<size; ++i)
1776 A = A + (*current_vector)[i]/(x+cln::cl_I(-1+i));
1777 return A;
1778}
1779
1780// The values in this function have been calculated using the program
1781// lanczos.cpp in the directory doc/examples. If you want to add more
1782// digits, be sure to read the comments in that file.
1784{ if (coeffs)
1785 return;
1786 /* Use four different arrays for different accuracies. */
1787 coeffs = new std::vector<cln::cl_N>[4];
1788 std::vector<cln::cl_N> coeffs_12(12);
1789 /* twelve coefficients follow. */
1790 coeffs_12[0] = "1.000000000000000002194974863102775496587";
1791 coeffs_12[1] = "133550.502942477423232096703994753698903";
1792 coeffs_12[2] = "-492930.93529936026920053070245469905582";
1793 coeffs_12[3] = "741287.473697611642492293025524275986598";
1794 coeffs_12[4] = "-585097.37760399665198416642641725036094";
1795 coeffs_12[5] = "260425.270330385275465083772352301818652";
1796 coeffs_12[6] = "-65413.3533961142651069690504470463782994";
1797 coeffs_12[7] = "8801.45963508441793636152568413199291892";
1798 coeffs_12[8] = "-564.805024129362118607692062642312799553";
1799 coeffs_12[9] = "13.80379833961490898061357227729422691903";
1800 coeffs_12[10] = "-0.0807817619724537563116612761921260762075";
1801 coeffs_12[11] = "3.47974801622326717770813986587340515986E-5";
1802 coeffs[0].swap(coeffs_12);
1803 std::vector<cln::cl_N> coeffs_30(30);
1804 /* thirty coefficients follow. */
1805 coeffs_30[0] = "1.0000000000000000000000000000000000000000000000445658922238202528026977308762";
1806 coeffs_30[1] = "1.40445649204966682962030786915579421135474600150789821268713805046080310901683E13";
1807 coeffs_30[2] = "-1.4473384178280338809560100504713144673757322488310852336205875273000116908753E14";
1808 coeffs_30[3] = "6.9392104219998816400402602197781299548036066538116472480223222192156630720206E14";
1809 coeffs_30[4] = "-2.05552680548452350127164925238339710431333013110755662640014074226849466382297E15";
1810 coeffs_30[5] = "4.21346047774975891986783355395961145235696863271597017695734168781011785582523E15";
1811 coeffs_30[6] = "-6.3439111294220458481092019992445750626799029041090235945435769621790257585491E15";
1812 coeffs_30[7] = "7.2684029986336427327225410026373012514882246322145965580608264703248155838791E15";
1813 coeffs_30[8] = "-6.4784969409198000751978874152931803231807770528527455966624850088042561231024E15";
1814 coeffs_30[9] = "4.5545745239457403086706103662737668418631761744785802123770605916210445083544E15";
1815 coeffs_30[10] = "-2.54592491966737919409139938046543941491145224466411852277136834553178078105403E15";
1816 coeffs_30[11] = "1.1356718195163150156198936885250451780214219874255251444701005988134747787666E15";
1817 coeffs_30[12] = "-4.04275236298036712070700727222520609783336229393218886420197964965371362011123E14";
1818 coeffs_30[13] = "1.14472757259832757229433124273590647229089622322597383276758880048004748372644E14";
1819 coeffs_30[14] = "-2.56166271828342920179612184110684658183432315551120625854181503468327037516717E13";
1820 coeffs_30[15] = "4.4861708254018935131376878973710146069395814469656232761173409397653101421558E12";
1821 coeffs_30[16] = "-6.0657495816705687896607821799338217335976369800808791959096705890743701166037E11";
1822 coeffs_30[17] = "6.21975328147406581536747878587069711930541459818297675578654403265380823122363E10";
1823 coeffs_30[18] = "-4.7255003764027411113501086372508071116675161078057298991208060427341079636661E9";
1824 coeffs_30[19] = "2.5814613908651936680441351265410235295992556406609945442133129515256889464315E8";
1825 coeffs_30[20] = "-9752115.5047412418881417732027953903591189993329461844657371497174389592441887";
1826 coeffs_30[21] = "242056.60372411758318197954509546521913927205056839365620249547101194072057318";
1827 coeffs_30[22] = "-3686.17673045938850138289555088011327333352145765167200561022138925168680049115";
1828 coeffs_30[23] = "31.3494924501834034405048975310989414795238339283146314931357877820190435258517";
1829 coeffs_30[24] = "-0.130254774344853676030752542814176943723937677940441021884132211221409382350105";
1830 coeffs_30[25] = "2.16625679868432886771581352257834967866602495378408740265571976698475288337338E-4";
1831 coeffs_30[26] = "-1.05077239977528252603869373455592388508233760416601143477182890107978206726294E-7";
1832 coeffs_30[27] = "8.5728436055212340846907439451102962820713733082683634385104363203776378266115E-12";
1833 coeffs_30[28] = "-3.9175430218003196379961975369936752665267219444417121562332986822123821080906E-17";
1834 coeffs_30[29] = "1.06841715008998384033789050831892757796251622802680860264598247667384268519263E-24";
1835 coeffs[1].swap(coeffs_30);
1836 std::vector<cln::cl_N> coeffs_60(60);
1837 /* sixty coefficients follow. */
1838 coeffs_60[0] = "1.000000000000000000000000000000000000000000000000000000000000000000000000000000000000007301368866363013444179014835363181183419450549774";
1839 coeffs_60[1] = "2.13152397525281235754468356918725048606852617746577461250754322057711822075135461598274984226013367948201688447853106595646692682568953E26";
1840 coeffs_60[2] = "-4.548529924829267669336610112411669181387790087825260737133755173032543313325682598833009521765336124891170163525664509845740222794717604E27";
1841 coeffs_60[3] = "4.6879437426294973235875133160595324795437824160731608900005486977197800919261614723948577079551305728583507312310069280623018775850412E28";
1842 coeffs_60[4] = "-3.10861265267020467624457768823845414206135580030123228715133927538323570190367768297139526311161786169387040978744732051184844409191231E29";
1843 coeffs_60[5] = "1.490599577483981276717037178787147902256911799467742317379590487947009001487476793680630580522955117318124168494382267800788736334308E30";
1844 coeffs_60[6] = "-5.50755504045738806940255910881807353185463857314393682608295373644157298562106198431098170107741597645409216199785852260920496247655646E30";
1845 coeffs_60[7] = "1.631668518639067070100242032960081591016027803392225476881353619523143028349554534276268195490790113905102273979193269720381236708853746E31";
1846 coeffs_60[8] = "-3.9823057865511431381368541930378720290638930941334849821428293955264049587073723565727061718251925950255036781219414607001763225298119E31";
1847 coeffs_60[9] = "8.16425963140638737297557821827674142140347732117757126331775708561852858085860735359056658172512163756926693444882201094206795155146202E31";
1848 coeffs_60[10] = "-1.426548236351667330492229413193359354309705120770113917370333660827270957172393778178051742077714657388432785747112574456061555034588373E32";
1849 coeffs_60[11] = "2.14821861694536170414714365485614715949416083667308573285807894910742621740039595483105992136915471547998283891842897000924199509164799E32";
1850 coeffs_60[12] = "-2.81233281290021706519566203146379395136352592819625378308636458418501787286411189089807465993150834399778687427813779950602826375635436E32";
1851 coeffs_60[13] = "3.222783358826786224404373038021509245352188734386849874296356404770508945395436142634892645963851510893216093037595555902121365717716154E32";
1852 coeffs_60[14] = "-3.250409075716999887328836263791911196138647661969351655925350981785153422033954649154242209471752219326556302767677017396179477496948985E32";
1853 coeffs_60[15] = "2.897783210826628399578158893643627107049805015801395657097255344786041806868455726759715576609013221857885740543509045196763816109465777E32";
1854 coeffs_60[16] = "-2.29136919195969647663887561122314618826917230275433296293059354280077561407373070937197721317435316121212106870152659174216557412788874E32";
1855 coeffs_60[17] = "1.611288006928200619663496306945576194382628760891807800193737346171844871295031418730500946186238469256168610033434708290528870722514911E32";
1856 coeffs_60[18] = "-1.009632466053186015034182792930705530447465885425278324598880797572411588461783484686932989855033967294215840157892487264656571258327313E32";
1857 coeffs_60[19] = "5.64520651042784179741815642438421132518008517154942873706221206276337451930555926854271086501686252334516011905237101877044320182980053E31";
1858 coeffs_60[20] = "-2.81912877441595327683492797147781153304080114512116755424671954256427789550109614317215500473322621746416096887803928883800132453510579E31";
1859 coeffs_60[21] = "1.257934257434294354026338893625531254891110662111965279263894740714811495074726866375858553579650295684850594211744093582249745250079168E31";
1860 coeffs_60[22] = "-5.01544407232599962845688086323662774702854661522104499328570796808858930542190600193190967249971520736397504227594619670310759235566195E30";
1861 coeffs_60[23] = "1.786035425040937365122699272239542501767986628253845452136132211710520249195280548478081559036323184490150479070929923213045153333111476E30";
1862 coeffs_60[24] = "-5.67605430104368150038863866362066081946938075036837029856903803768657069745962581310398542442108872722631658677177822712376500859930109E29";
1863 coeffs_60[25] = "1.607878222558573982505999018371559631909289246981490321219650132406126936263403946310818841465409950661433241956831540547593847161412447E29";
1864 coeffs_60[26] = "-4.05332042374309456146169816144083508836132423024788116321074411679252452773181941601763924562378611113519038766273534176937279867894066E28";
1865 coeffs_60[27] = "9.07493596543985672039002802030098143847503854224661484396413496012780904911929710460264147600378604646912175235271954302119768907744722E27";
1866 coeffs_60[28] = "-1.800074018924350353143489874038038169034914082090587278672411654146678304871125651069902339241049552886098125667720181441150399048551683E27";
1867 coeffs_60[29] = "3.154250688078046681602499411296013099183808016176992164829953752437167774310360166977972581670851790753785195101324694758021403186162394E26";
1868 coeffs_60[30] = "-4.86629244083379932983782216256143990390210226006560452979433243294026128577640975980482675864760717747936401374948595060083674140963469E25";
1869 coeffs_60[31] = "6.58428611248406176613133080039790689602908099995907522692286902207707012485115422092589779128693214784991500936878932461139361901566087E24";
1870 coeffs_60[32] = "-7.77846893445970039116628280774361378296946997639645747353868461156972352366479641995295874152354776734003001337605345817120316052066992E23";
1871 coeffs_60[33] = "7.98268735994772082084918485121285571015813651374688487489679943603727447378945977989630573952891101472578977333720105112837324185659362E22";
1872 coeffs_60[34] = "-7.07562692971089746095546542541499489835693326760069291570193808615779224025348460132750549389189539682228913778397783434269420284483726E21";
1873 coeffs_60[35] = "5.381346729881846847476909845563262674288431852755093265786345982700437823098162630059919716651136095720390719236493773958116646152386075E20";
1874 coeffs_60[36] = "-3.4856856542678356876484367392130359114150104987588151214926676834365219571876912071608359944324610844909103855562977795837329347647911E19";
1875 coeffs_60[37] = "1.90665542883474657677037950113781854248329048412482665873254624417996252139138481002200079466749149325431679310476862249520001277129217E18";
1876 coeffs_60[38] = "-8.72254994006151131395107200045641306281165826830744222866994799005490857259177347821280095689079457417603257537321939951004603693393316E16";
1877 coeffs_60[39] = "3.30066663941625244322555483012774856710545517350986120571194216206848716066355962922968824538055042855044917677713272771363157100391997E15";
1878 coeffs_60[40] = "-1.020092089391030771746960980075254826475625668908623135552682999358854102567810002206013823466362488147261886160954607897574298699485318E14";
1879 coeffs_60[41] = "2.537518136375035057088980117582986067754938584307761188810498418760131416720976321039509027979006220650166651208980823946300429957067604E12";
1880 coeffs_60[42] = "-4.99523339577986301543863423322168947825482352498610406809585164155176248614834684219539096936869521198401912030883142734471627752449382E10";
1881 coeffs_60[43] = "7.62961024898383965152735310352890448678585029645218309944823403624458716639413808284778269959424212699922000610764015063766429510499158E8";
1882 coeffs_60[44] = "-8834336.1370238009649936481782352367054397712953420330251745022286767420934395739052638862442455545176778475848478708230456099596423988";
1883 coeffs_60[45] = "75445.9196169409678879362111492280315111800786619928588067631801224813888137547544321383450353324917130013984795690223150786036557545929";
1884 coeffs_60[46] = "-459.8458738886001056822131294892698769439281099450630714273592488999986769567563218319365007529495798105783705491469742412340762305916056";
1885 coeffs_60[47] = "1.922366163948404706136462977961544621491268971185908661903800938507393909575693892375103171073678191394626251633433930639174604982075991";
1886 coeffs_60[48] = "-0.00524987734300376305383172698735851896799115189212445098242699916121836353753886238290792298378658233479210271064792489583846726184351881";
1887 coeffs_60[49] = "8.81521840386771771843311455937479573971716020932982441671173279504850522350287085310420429874536637110755391716691475171030099411021337E-6";
1888 coeffs_60[50] = "-8.42883518072336499031504944519862331274440110738275125460829656821173301216150526266773841539372995424665091651911614576906895281293397E-9";
1889 coeffs_60[51] = "4.1559932977982056953309753711587342647729282359841592558743510304569204546713517319749817560490538963802716194154620384631597656968764E-12";
1890 coeffs_60[52] = "-9.26494376646923216540342478135986593801117330292329759013854851055518195892306285985326338987592590319793280515888731024676428929933443E-16";
1891 coeffs_60[53] = "7.80165274836868312019654872701978288745672229459298320116385383568401529728308916875595120085091565550085090877341856355815270191309086E-20";
1892 coeffs_60[54] = "-1.922049272463411538721456378153955404697617250978865956250065913541261535132290272529565880980548519758359440057376306817458561627984943E-24";
1893 coeffs_60[55] = "9.46189821976955264154519811789356895736753858729897267240554901027053652869864043679401817030067356960879571432881603836052222728024736E-30";
1894 coeffs_60[56] = "-5.06814507370603015985813829025522226614719112357562650414521252967497371724973383019436312018485582224796590023220166954083973156538672E-36";
1895 coeffs_60[57] = "1.022249951013180267209479446016461291488484443236553319305574600271584296178678167457933405768832443689762998392188667506451117069946568E-43";
1896 coeffs_60[58] = "-1.158776990252157075591666544736990249102708476419363164106801472497162421792350234416969073422311477683246469337273059290064112071625785E-47";
1897 coeffs_60[59] = "4.27222387142756413870104074160770434521893587460314314301300261552300727494374933435001642531897059406263033431558827297492879960920275E-49";
1898 coeffs[2].swap(coeffs_60);
1899 std::vector<cln::cl_N> coeffs_120(120);
1900 /* 120 coefficients follow. */
1901 coeffs_120[0] = "1.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000060166025676976656004344991957470171590616719251813003320122316373430327091055571";
1902 coeffs_120[1] = "3.4497260317073952007403696383770947678893302981614719279265682622766639811173298171511730607823612517530376844024218507032522459279662180470113961839690189982241536061314319614353993672315096520499373373015802582693649149063603309572777186560148513524E52";
1903 coeffs_120[2] = "-1.4975581565000729527538170857594663742319328925831469933998274880997450758924704742659571258591716460336677591345828722528085692201176737000527729600671680178988361119859420301844184208079614468449296788394212801103162564922199859549237082372776667464E54";
1904 coeffs_120[3] = "3.1957762163065481328529158845807843312720427291703934903666695190945338610786360201875291048323381336567812569171891600400186742244091402566230953251621720778096033490814848238417212345597975915378369497445590090951446115848410773972658485451963575288E55";
1905 coeffs_120[4] = "-4.4689623509319752841609439083871154399631153121231062689347162975834499076693093642474289117173045421812089871506249999929076992135798925381959196225961791389783472385803138226317976820364502651110639008585046458007356178875618627927171581950486124233E56";
1906 coeffs_120[5] = "4.606068718424276543329442566011849623375399823565351941825685310847310447457609082356012685588953435307896055516214072529445026693975872604267789672469025113562486157850515006504573881812473997762948360804814769118883992998548055557441646946685125118E57";
1907 coeffs_120[6] = "-3.7314461146854666499272326592212099391213696621869706562566612605818861385928266960370453310708465394226398321257508947092784006446784523328681347046673172481746936234783770854350210504707173921547794426833735429199925024679815789545465854297845328325E58";
1908 coeffs_120[7] = "2.474425401670711256989398808079221298913654027234786607507813220440957186918973475366048940039541074278444160674001228864321389049663140487504402096319272526201782217412803784224929141788255724630940381342478088455751340159338461174261577243566175687E59";
1909 coeffs_120[8] = "-1.3811875718622847750042362590249762290599823842851465148429257970104907280458901604054390293828410620002370526629527048636126473391278330353375163563724888073254512227198849135923692811222561965740181944727170495185714496890490479692693474125883791901E60";
1910 coeffs_120[9] = "6.623089858532754482582703479109160446021743439335073883710993620625687271109284320721410901325182604938578905712329203551531862389936804947105415805829404869727743706364603519433193421234231031076682156125442577335383798263985569601899041876776866622E60";
1911 coeffs_120[10] = "-2.7709515004299938864490083840820063124223529009388282231525445615826433364331567602934962481829061542793349831611106716261513624279121506887680318284535361848032886450351898892264386237450622827397559067350672965967202437971930333676917000390477963866E61";
1912 coeffs_120[11] = "1.02386112293172223921263435003659366453292875147351461165091656394534393086780717052422266565203902889367201592668259202439166666819852985689989767402099479793087277263747942659943270101657408462079787397068550734516045511611701546009078868077038808757E62";
1913 coeffs_120[12] = "-3.3740197731917655541744976218513993073761175468772389726802124778433432226803314067431898210976006853342921093194297198044021414900546886804610561082663076825192459864843102368108908666053756409152492134638014803233805912009476407113691438596300794146E62";
1914 coeffs_120[13] = "9.996217786487670355655374796561399578704294298563457268841140703036898520360123177193155340144551120260016445533739357030180277693840431766824113840895797510199955331557143980793267795747200088810293047731873192410526786931879590684673414288653913515E62";
1915 coeffs_120[14] = "-2.6804750990199908441443350402311488850281543531194918304012545530803283220192092419107511475988099394746800512008906823331244710178292896561401750166818497729239682879419868799442954945496319510685448344062610897698253544876306888341881254056091234759E63";
1916 coeffs_120[15] = "6.5422964482833531603879610057815197301035372862466791995455246163778529556707384145891730234453157337303612060344197138180893720879196243783337539071284141345021864817147590781643393947019750353147151780290464319306645652085743359080495090595200531486E63";
1917 coeffs_120[16] = "-1.4604487304348366496825146715570516556564950771546738885215899741781982964860978993963272314092830563320794184042847908967120542212316261920409301852237223308467032419706968676861616456179895880956772385510853673982424825597152850339588189159102980666E64";
1918 coeffs_120[17] = "2.9942297466313630831467808691292548682230644559492580161942357031681185068971393754352871129412787966878287389513657398203481589163498279625760316093736277896138061249076616695157422053188087353540756151586375196486093987258640269607104978950906670704E64";
1919 coeffs_120[18] = "-5.658399283776588293772725313973093187743120982052603944865098913586526167668102733207163739469977271584007101869254711458133873627143366757941713180350370056955237604551850024423889291598422917971467957836705917204903687959901869098540153925178692732E64";
1920 coeffs_120[19] = "9.887368951584101622633538892976576123080629424367489037686110916264512731398396560326756128205833849608930564615629875435100785011872254223744155330328477703592008501954532369042429700051416733748454350165515933757314793533786385104271308839525639768E64";
1921 coeffs_120[20] = "-1.6019504550228766725078575508839635707919311420327486864048705642201106895239857903763208049376932672160478820626774934879424498715258948985194011690204294886396827446040036506699933786721588971678753877371518212675519147446728054067639530675249526082E65";
1922 coeffs_120[21] = "2.4124568469636899540706437405441629738413207418758399778576327598069435452295650039157974716832514441625728576753250737726840109004878753294786785674578138926529507088264657400701828947949531197915861820274684954206665488761473274445827472596875582911E65";
1923 coeffs_120[22] = "-3.3841653726400000079488483558717068873181168418395106876260246491163166726612427450773591871178866824643300679819366574162583413250423974373322308130319007820863363304629451933781204964221002853140392226489420463827400812929748772154909106349410663293E65";
1924 coeffs_120[23] = "4.4305670380812288478773282114598811227924298131011853412998479811262358077680067168455361591598296346480072528806092976336961470360354620203822421524751468329936930212919915114854135818230382164555078957880154875221176513434392525189922941290050575762E65";
1925 coeffs_120[24] = "-5.4228176962574428947233160003094662570284359565811627941401342797491445636152854865132166939274138115146035207618348708829039395974942115203986578386666664945394109693178927438991059414217518334491360514633536224841961444935232548483014691997071543828E65";
1926 coeffs_120[25] = "6.214554789078092267222051275213928685756510105900211846145778269883351640710249139978059486185007208670776100912863866582278800642692097830681092656540813877576256048148229340562594504915197956922387464825593941922429396202734006609196778697870436014E65";
1927 coeffs_120[26] = "-6.6773485088895517986512141063848395783979405189075416643094283756118912554557672721632998501682483143868731647507940026369035991063923616298815637819145806214374157182512600196214559297579802178103007615921637577873304407436850546650711237281572008424E65";
1928 coeffs_120[27] = "6.734925330317694704469314845373778479111077864464012553672292377883525864326847400954413754291163739900219432201437895152976917857427306427115230048061308424221525123820493252697918698598513232640014129066982507718245232516657455821629338155744427538E65";
1929 coeffs_120[28] = "-6.383582878496429871173501676061533991181960023885889537277705274319508246322757005217436814481703326467002683699047193244918123789600842413060331898515872574523803039779899326755393070345055586059441271293717500426377884349137309244757708993087958455E65";
1930 coeffs_120[29] = "5.6913529405959275511022780614007027176288843526260372650173869440228336395668389555081751187360483397341349300975285817498083216487282169140596290796279875175764991375447348355187090404486257481827615256024271536396461908482904537799521891879785332367E65";
1931 coeffs_120[30] = "-4.776996734587211249165031248400648409423153869394988746115380756083311986805300361459383722536698540926452976310737678416019979202990255666249869917768885659350216400547190883549730059461513588008706974008085270389354525600694962352952715682056375518E65";
1932 coeffs_120[31] = "3.7775367524287124255443145064623569295746034916892464094281613465046063954544055573214473155479196552309207209647540614474216985097792266203411723082566949978062697757983354600199199984244302856099811940389910544756210676400851240882142140969864764468E65";
1933 coeffs_120[32] = "-2.8161966171919236962021901287860232075259781334554793534017516884995332700369401674058517414969240048359891934178343992080557338603528540157030635217682829894098359736903943078409166055640608627322968554856315650475005062493399450913753277478547118352E65";
1934 coeffs_120[33] = "1.9804651678733327456903212258413521470612733719543558365536494344764973229749132899499862883369665827727506916597326744330471802598610837032598656197205238983585794266213317465548361566327435762497208877015986690267754534342053368396181078097467171858E65";
1935 coeffs_120[34] = "-1.3144275813663231527166312401997093907605894997476799416306355417933431514642211250592825223377757973148122542735038736133300194844844655961425683877005418926364412294006123974642296395931311307760050290069031276972832755406161248577410950671224318855E65";
1936 coeffs_120[35] = "8.2367260670024829522614096155108151082106397954565823313893008773930966293786646885943761866773022391428854862805955553810619924412431932999726399857050871862122529700098570542876369425991842818202826823540112018849926644955200888291063471724203391548E64";
1937 coeffs_120[36] = "-4.8749964750377069822933994525197085013480654713783888755556109773660249389776804499013517227967180500633060271953473316017147397601291325922904139209860429881054757911243087427393920494271315804033914011087815785282473032714919188637172020633929566123E64";
1938 coeffs_120[37] = "2.7259664773094932979328467102942769029907299417406744864696200699394594868759231280169149208728483197299648608091313447896342349454038879581019820193316159535211365363553004387852005780736869678460092714636910972426808304270369152189989142121207224142E64";
1939 coeffs_120[38] = "-1.440426226855027726783521340050349148103881707415523724377763633849488875095817796257895327883428230885349760692732068174527147156893314818583058727424251827006457849321094911262818557954829248070170426870959233263267490276774734065709978749400927185E64";
1940 coeffs_120[39] = "7.193790858249547212173205531149034887209275529426061411129294234841122474820371873361100215884757249851960370114629943083807936135915003201800204713978377250292881453568756354858194614039311248345228434431020394729104593125888325843724239404594830488E63";
1941 coeffs_120[40] = "-3.3960029336234301755324970935705944032408435186630159101426062821929524761770439420961993430248258136340087498829339209014794230274407979103789924433683527009234592433480445831820377517333956042612961562022604325181492952329031432513768020816986814393E63";
1942 coeffs_120[41] = "1.5154618904112106565112797443687014834429200069480460967081898435635890576815349145926430052596468033907024005478559584915319911380449387176530845634833237204659108290330613043367085829373476690728522550189678729181372902816898536141595215616716630939E63";
1943 coeffs_120[42] = "-6.3927647843464050458917092484911245813170740434503951669888756878206365814594631676413018245438405308353724023007754523096143775098898268650326908751515418318201372985246418468844138298345777180517875695389655616000832495210812684049030674085212697428E62";
1944 coeffs_120[43] = "2.5490394366379355452002449693074954071810215414182359403355645652443600688717811337587901850157210686351097591461582890354662732336749618027675479531031836144519267481752770036252137747675754903974915999567019837855523058289177148692481402871253211324E62";
1945 coeffs_120[44] = "-9.606466879185328464666445215840505657671157752044466089989040292763536710311599947887918708456526669882072519263973105599580140713596301388561639705589314111762600854460059589939760935803484446888352368360433606245369171819922425771642408570388554052E61";
1946 coeffs_120[45] = "3.4212392804723358445152430359637323789304939688937873921904941412927756295848104328630952153624979489607759834359194243032109828811134607612715016533909375981353098879969472700079242226099049323998286020341979178782935852542355220403299144612362244738E61";
1947 coeffs_120[46] = "-1.15119134701605919461057899755821946453925102458815313053351247263978303346790555715641513756644038607667203392289423834966320935498856390723555744530789850290369071103208529608463210398231590077340268751005311531529356083188150256469829678612245577446E61";
1948 coeffs_120[47] = "3.6588622583411432033523084711047679684233731128914152509273818448610176621874654252431411048902598388935105893946323641003087000410095802098177375492833543391040706755511234323104846636415419597151008153829618275459606044459923718022154035121167198784E60";
1949 coeffs_120[48] = "-1.0981148514467449476248066871827754422009180048705085132882492434176164929454140182025449006310206725429473330511884213470600326782740663313311256352613244044500057688932314549669435095761340307817735687643806167483576999980691227831561891265615422486E60";
1950 coeffs_120[49] = "3.110998209448997739767747906196101611409160829345058138064861244336130082424927251851805875584197897229644157110035272012393338413235208343342708685139492629786435072305986349067452739209758702078026647999828517440754895711519542954337931090643534216E59";
1951 coeffs_120[50] = "-8.3162485922574890007748232799240657004521608654422032389269811102140449056333167761296051794842882201869698963586030628312914066893199727852512779320175952962772072653493447297721128265231294406156925752496310087025926300388984242024436858845487466277E58";
1952 coeffs_120[51] = "2.0966904516945699848169820408710416999765756367767199815424586610234585829069218729220161654233351574517459523275756901094737085187558904179251813051891939079067686519817858153690134828671544815635956527611986498479411756457222935682849773436423295467E58";
1953 coeffs_120[52] = "-4.983121305881207125553776640558094509942884568949257704810973508397697839859902664482541160531856121365759763455699578413261749913567077796586919935391984240753355552646184306812426079133011894826183873855851966310877619118554510972675999316631346679E57";
1954 coeffs_120[53] = "1.1158023601951707374356047495258406415892974604387009613173591921419195864040428221070481312383179580486787822935456571355463718115785982888531393271665510645725283439572279946304699780331972095822869500426555507626639723865965516308476400920600382357E57";
1955 coeffs_120[54] = "-2.3524850615012075127499506758220926725372558166170912192116695445007095502575329450463479860779122789467638956004572617263549199692255055063165454868102165975951768676031140009643202074220557325155838768661030361538572755082660730808847591840060467064E56";
1956 coeffs_120[55] = "4.6669318431895615057208826641721251136909284138581355667925884903657855204100373961676117747969449100495897986226609480142908763981931305129946569690612924941456739524153327260627627771254850382983581593260532259539447965597396206625726656509884058042E55";
1957 coeffs_120[56] = "-8.7053773217442419007560462613131691749734845382618514999712446313788486289774350240165530159591402631439776213579542026449818009956904779042347595401565525081115611496250192338958392965746523979241969677734430475813057146574920495171984815351708574336E54";
1958 coeffs_120[57] = "1.5256552489620511464542280446639568546874380361953025589702692266626310669215652044048704882910412155084167930513006634430352568411276836880182348033924636960897794333644980768878022821035659978039286230061734024129667272393315169114199838321062607299E54";
1959 coeffs_120[58] = "-2.5099934505534008439782195609383796207770494575364994376922414269548303512602084430128307108303305643530918354709126474742035537827601791192999467996479881350277448357927640707861695639576629921988481117017137420422963638277868648516492581097660522547E53";
1960 coeffs_120[59] = "3.872963359882179682964169603201046384616694634651871844057456079738892419308420856725974686574980381399016464501318163662938118593626674643538005780375691959391996340141057698193381380484420715733863044826589570328349973407598034428591146829028071358E52";
1961 coeffs_120[60] = "-5.59947633823301408044455223877913062308847941596689956112764416031828413291312481723036534632655608672535030921469531903033364444816678754679807809159478411100820014592865068932440734964265842594875758737421026093110624848762070026616564150314951394E51";
1962 coeffs_120[61] = "7.57762861280525531438216991274899157834431478755285945898172885086150762425529113816148806028462888396660067975773261101497666568988246606837690320098870044112671149076084444095163491848634465373822951831018725769263871497616640732007420499659069842E50";
1963 coeffs_120[62] = "-9.587786106526273406187878833167940811862067040706459726637556599860244751467528905534431960251166924163661188573831350928972391892492380823531476387272791432306808700507685765850397294118719242350333451452137838374120658600691461454898577711260078952E49";
1964 coeffs_120[63] = "1.13288726401696728230264357306938076698155303500407071418573081766541065136778223998897791839613776442037036668986628122296219518360439574147622758002647495909592177914657175019781723803408732148262293125845657503039410078589916085532057725749397276232E49";
1965 coeffs_120[64] = "-1.24849787223197441956280303618704887038709792250544105638342097080498907831514597860418910331910245753340059089147824955071899315894649696314820492532126554883819507650973976145456660786429117569053901704116877128391672511345177517877672824534972448216E48";
1966 coeffs_120[65] = "1.2815463720972693091316233381473056495608681859925407504190742949467232967966271661733907550222983737930524555721493736920130260377888287772008209963158064973076933575966719577456540496444474944074979736374259087350416613616719928507635667369740203319E47";
1967 coeffs_120[66] = "-1.2234887340201843394744986892310393596065877342193196880417674427168862926389642850813687099959036354499094230765541977493433449153438766822382486040215211159359175689369230076522107734270943423777076523650345103234411047700646432924770659676420158487E46";
1968 coeffs_120[67] = "1.0847187881607033339631651118075716564835185723270640503055198532318419482330026641941088359447807553514405522074008969583213861070993661224871455023365601323302778638456843760403418046238489404394483720438784739822580385277055304353975028280477740796E45";
1969 coeffs_120[68] = "-8.9160881476675795743767277986448579964735858351472748620623279571408606135698760493224031735408212513500922230670883171668702983221921543376953865813604783695111225412173880768170509738290662806468458720236121755965944855709552219268353813402612336565E43";
1970 coeffs_120[69] = "6.782864920104031936272293608616215844503387641476821968620772153274069873138756405621471099960069602613619793775294358177761533027360002770186566164041138064221354961783144649476276625776241973967317262115970868665380343599565811072109785000646703404E42";
1971 coeffs_120[70] = "-4.7667808452660756441368384708874451089976319738852731080495062883240643961463680300964077232336439626019128672679703771884184482488932861160134911816225569323838390204451496983578077563176966732010513231048738892639707790407292070646798259086924770995E41";
1972 coeffs_120[71] = "3.0885057140860079424719232591765602418793465632939298397987628606701994268384966881159469651774584648643122830739130127593326652998108850492039117928976417052691273804304806596509726701594300563830431015215234640024338277573401498998072908815285293868E40";
1973 coeffs_120[72] = "-1.8410405906573614531857309495652487774337134256805076777639383854080936219680656594060736479739035202182601529001321266214227848431889644620036213870966329509961114940541333851155401637197303308322414678191211465563854205816313387785764908216851396633E39";
1974 coeffs_120[73] = "1.0073694433024942271325653907485159683302928496826793112696958500366488338508211620934892875328717073528902110227362794694820010124321343709182901273795782541866547318841893692957109947576483162095037812781379193423759617638948859880051822460818418552E38";
1975 coeffs_120[74] = "-5.0475051506252944853315611134428802424958512917967945464108691542854207821486654807141339210375899950551724141366521361887864357385178212628348794663127149312605456165451981719848656127310229221238908657530297751682848475855876378576874607521597136906E36";
1976 coeffs_120[75] = "2.3099766115359817610656986443137072041797751710805647712896098246833051023271876304983288225638204962631413469467959017768113430777226924099787875749611560913177631681394153889301715579572842026181746028117354815826836594637709952294015960031772162547E35";
1977 coeffs_120[76] = "-9.629053850440590569772960665435833408449876392175761493622541259322053209458881628458334353756739601360772251654643632187697620334088992038575944303101187678397564511853344433267011583960451100374611538881978045643233876974513962362084978095067025623E33";
1978 coeffs_120[77] = "3.6452126546120530579393646694066971671091434168707822859890104373691687449831950255953317231572802167174179528347370588567969602221261721708890001616085516755796796282628169745443137768549800602834096924025507345446292715781107949529692160434800323E32";
1979 coeffs_120[78] = "-1.2492564030201607643388368733220662634846470405464496879151879822123866671204541555507638492613046717628358162773937737774832271305618491107140304474323049182605167775847584622690299098207979849043605983558768056117581593008210986863088433891075743152E31";
1980 coeffs_120[79] = "3.8627447638297686357472526935538070834588578920414538227245516723308987020816841052950727259618753144711425856434270832495754300189881199851254605718213699755258867641301730599979474865704144160112269948588154919128986989885090481959424806312935273075E29";
1981 coeffs_120[80] = "-1.0736758703963497284148841547397192249226725101007524773889805877171959717011395181953504058502607435217886087332761920207901621377557079619638699346496468750455986591040017334237734940082333954589067611955107878899677189289648293223359861027746438121E28";
1982 coeffs_120[81] = "2.6722714785740082059347577649909834926335247252399259683264830680945466475595847553753509546415283809619181144796536494882020159787371993099998263815645014317923922311421330376008111312767167437401741178863083976628261471599264811824656877164988491393E26";
1983 coeffs_120[82] = "-5.9304047185329750657632568788530498935629656326502947505210292278638825286675833282579834326765999907183142489791905921257123760969245535649745876992946512033156167841406724363867902645010435996961270021857807247440211477908060243655541266857227638988E24";
1984 coeffs_120[83] = "1.16817022089143274700208191285335154155497013626172270535715899131321522799010543339535307798264602677955894930046454353008462671803498794203612585729705145312299224155123919877760274781582850868001155383467754608529345730226972329454404720862870618607E23";
1985 coeffs_120[84] = "-2.03239515657536501213472165328009690017090356606547792466197690386716728380893226886179282271040418637806139515373566132123131620086873213475424131345589653019635327048678766191769576650893957440830876852296666120473866301097954633389040518870395767125E21";
1986 coeffs_120[85] = "3.1065334503269182605978912331263087603258864771943471481540265718169544724355602987297631515907391374943512439350265433478241465606056187134785807375293801936399644663199667496663518171930757047012102683120173353568660795955174938680248863153863947508E19";
1987 coeffs_120[86] = "-4.1476244154347831048636005592317388215032295704489937704602030038303705695463546496640638505584602502764898113504560236629804442607426019604639559875021291459916615723777004493344143132459204229291886967479716413925814352313734234340863490128872380307E17";
1988 coeffs_120[87] = "4.8067293487250079673131214670887682215073707729621636364424152483295071605326220176372385638491275365750175037404843071051780212494354459897540110089573898336327006157766256896984455454193271731091632286742192439925748114360605084629432813597189767538E15";
1989 coeffs_120[88] = "-4.8023544548381246628003457039588616467438691159189277447469028024236284353593054364114519649309416187375157096932150251663679454372678125518452171003992957433311257042292636706448339781439297178835786059318810522834929923770539615271536113963729385909E13";
1990 coeffs_120[89] = "4.1055087514683476865343055835875083237542317413651906253058979029083965525058905726360233143503628224856307545474786181299719957472120906835233967660557875100202077212004953379299507351564181758434304881046845705855303854083493519588411179065109026834E11";
1991 coeffs_120[90] = "-2.9787503393847675871205038539267895335240592213878943742323972872214441728681744433089698206110260166068266926018988659692353298939109421567999207730700359726920482465669373553804927535369930188390246988033893916611435406224816632683980860607732310186E9";
1992 coeffs_120[91] = "1.8178328110729629877907010659834277046059726898311908447099830056830012488194646687474150289147446390570639168063598563291822008033517936194534129929881015025633519502485415000390171249019651579295905194415531994026553693578406432674734610095421683863E7";
1993 coeffs_120[92] = "-92391.136314434380495997449781381513978328604842061708454699991154771188446049720445502194923435235472458378926242100033122111143321209059959788378033220861638093951546784186137626553022963832352496255851690092415165826965388502958309163995296640164754";
1994 coeffs_120[93] = "386.82763074890451546182061419449593717951707520472938425276820204065379182568600735469831672149785863654956632602671563997131280046154927653332261114114005498875447205079045401364007035880825957300757663780618819785476980699579657587509130753204519233";
1995 coeffs_120[94] = "-1.3181204292571874302358432444324779303744749959754136019600954846045028319805074783759764870805734807693739252625657350494147444011046941331047057337345953605042408524072436811726898109072388160378243068564382623631658424851676817690976343859083960324";
1996 coeffs_120[95] = "0.003606538673252695455085947121496196507159591230095595764694813152630524319596509155920374890595867709349176662036024214476302717902368680224618116411588086562230407996267622244422187853090635901906175373997993725355114393033631058067900506212434600015";
1997 coeffs_120[96] = "-7.805244503909439374422205381130511738566245024242591464192744568789876715121004646510755612128565674260161510215430132815223049297785205382643947556846567064565241387424696940674258789227398935846571768027456535982674711768030751512030174841314425949E-6";
1998 coeffs_120[97] = "1.31373705470989377112938364152965446631228819123896570245455699237549295870321627234421140232628798373711221392827979836922621437205363811871692678679625916100572037589291239046725228767017131155814257944742981208252138821140381478767814046301821211856E-8";
1999 coeffs_120[98] = "-1.6872873094408224472617181717534409090015431593544429529131126514352910895332010213914243717484771690790552077128803350550170014347729272790464826195676369023970955260051387240496705602732313607409271794413329062030561818907163134089683283286623809325E-11";
2000 coeffs_120[99] = "1.6183083251905685095057354853863188515437903228178486856957070037813756492593759658405336450433607296873747595037080703825755020175480385843762609522889527239577435110258291566585028919336090916225831079571865425410181260759913688103716786795647286451E-14";
2001 coeffs_120[100] = "-1.13097359411474028225398794102354853670936316496817819635688647804142428962171772690717075128208102537660772310780986623575005236651312181907812813813504742701120603881086064664411899253566047514905519888629604717647221817372977488600336785871295304013E-17";
2002 coeffs_120[101] = "5.599216369109121957949255319730053610385733330502739423509794477602247233276045188197007198417289907263120960704056657544648432653622931077692740961599655386871075693202473992087883485704436336279135221721374640982826144708808646466699352755417123702E-21";
2003 coeffs_120[102] = "-1.9009180102993021108185348502624676395148544369474718879637745630712451378711342634099259114111847962752555305470572286326367888004493816251811794947276966269738750207359305252041104539066278002044545942171476984766923991983055271262414217352967659228E-24";
2004 coeffs_120[103] = "4.261262509940940316499754264112111685174274727656165126333137554124192224955656564229887938745508952447664695831728428607673797269945824475565104978593072684829487175697371245288754204324544164474840153141042852153497051337607734150135978754952561336E-28";
2005 coeffs_120[104] = "-6.033854291373449912236926137860325602686312455380825767485673949251953414778800668020214699151728472172651816317924130614791108454134597377848088327850505473503152696524861086193124979489104732214189466703901268332265826882296309653009237279831825243E-32";
2006 coeffs_120[105] = "5.1208402745272379096703574714836785944518835939702823617280147111145234914591060871138496110227453241036619229980622243972303295470574470937679143516006222494480144845809123492603651773613707216680534850900104861326332900592715684757980394834998321888E-36";
2007 coeffs_120[106] = "-2.4463535717946588550832618025289907099586319384566637643650142186828541109926588999585266911960640972919441499109750654299062004147686492034166034659422424525984094382368955916181276646903453872999065929058429821759475215620044891133652431220664754175E-40";
2008 coeffs_120[107] = "6.0973480699773886324239008989591793773608942051497498591908583910660358857815864266160341286217871697703362816166340947142517661604423899536979689047275448159991318658879804351288744125363072102852651926942302209139318098544348348564409845011546432615E-45";
2009 coeffs_120[108] = "-7.2234185761285078775026471720270426097727212523472472797635230392183067756271499246654638332288950167477129840028892565652782123508855602380279653475510712205780583313834027906297063690370430285856541927759405826980856379432703473274890527421175151858E-50";
2010 coeffs_120[109] = "3.6217112680215791206171182969894344487335819731880124290544082848140757826983885738735436324684863867140575000400288923606439193119990961489053513339202655922248092157737577138929144240507796562250602457839068582279379672722261563501188150876583184441E-55";
2011 coeffs_120[110] = "-6.6329300032795486066608594142675837603786558782159646987663521197523704085781830169369726460621246948945196657495305819768951424025780824076252490918306538895670861455244641773606980519824591785816943621538721352987553804824051051144609050417497894495E-61";
2012 coeffs_120[111] = "3.6664720904335295532012711597888717227860988776477301054518326674835421172405060906940404374163713097964932859351917152390238690399278248344863365606468942320103392909602843987855082225592776850615943708151738327210634139824601616072015258461809772448E-67";
2013 coeffs_120[112] = "-4.7466013179695826928232672846686064011594588664906398407027593213652099998530859940288723349213099851532139911079905393494419637612780994270110734378146177806681489226896952731800026849872070824592339117757940119304241732812925979963178130104280115315E-74";
2014 coeffs_120[113] = "1.0163707785221910939390789816391472677729665860532352695801597334766068288835382195560328979864550624486740471947632369344045378626680607890520366137741785540226552923584183986350590955499329375427326072319268396685478606934920507703868118038891818762E-81";
2015 coeffs_120[114] = "-3.4814151260242800905467399051937942442621710748397374123807284826536707678408888416026868585492229216524609739211131993326633970334388991812593549702868877534701822990946125111761892723042376117665640296993581745994557803052315791392349639065203872505E-90";
2016 coeffs_120[115] = "1.18525924288117432386770939895670573772658621857195305986011196724304231598127227408839423385042572374412446842112646168302015480830234100570192462192015131968307084609177540911503689228342834030959242458698413980031135644018348590823980902427540799814E-91";
2017 coeffs_120[116] = "-8.5714961216566153236700116412888006837408819915951896129362859520462766617634320531162919426026429378433105901035364956643086394331335747930198070611009941831387116980941022864465946989065467218665543814574849964435089931072761832853235509961870476035E-93";
2018 coeffs_120[117] = "4.5681983751743456413033268196376305093509590040595182930261094908859252761697530924655649930852283295534503341542929581967081012867692190108698698006237799801339418962091877730207560007839789937153876806052229193448161273005984514504886230869730232561E-94";
2019 coeffs_120[118] = "-1.5943139155457706045530478744891549581317663177038648406493256399589001327414318955746453934207742828511041930090849236963271943244329753764497401819704943705370596846318480510254313447057477914171472190541408193443142906466279172123681623644325254209E-95";
2020 coeffs_120[119] = "2.7319125666863032595604997603472305262880292377469053594326527505796348018540179196191192420176181194669607935656210005192217186286873953583571180312679155204061051208771126804209623533044988888808754656646355388901404252058383561064953226611421609762E-97";
2021 coeffs[3].swap(coeffs_120);
2022}
2023
2024static cln::float_format_t guess_precision(const cln::cl_N& x)
2025{
2026 cln::float_format_t prec = cln::default_float_format;
2027 if (!instanceof(realpart(x), cln::cl_RA_ring))
2028 prec = cln::float_format(cln::the<cln::cl_F>(realpart(x)));
2029 if (!instanceof(imagpart(x), cln::cl_RA_ring))
2030 prec = cln::float_format(cln::the<cln::cl_F>(imagpart(x)));
2031 return prec;
2032}
2033
2039const cln::cl_N lgamma(const cln::cl_N &x)
2040{
2041 cln::float_format_t prec = guess_precision(x);
2042 lanczos_coeffs lc;
2043 if (lc.sufficiently_accurate(prec)) {
2044 cln::cl_N pi_val = cln::pi(prec);
2045 if (realpart(x) < 0.5)
2046 return cln::log(pi_val) - cln::log(sin(pi_val*x))
2047 - lgamma(1 - x);
2048 cln::cl_N A = lc.calc_lanczos_A(x);
2049 cln::cl_N temp = x + lc.get_order() - cln::cl_N(1)/2;
2050 cln::cl_N result = log(cln::cl_I(2)*pi_val)/2
2051 + (x-cln::cl_N(1)/2)*log(temp)
2052 - temp
2053 + log(A);
2054 return result;
2055 }
2056 else
2057 throw dunno();
2058}
2059
2061{
2062 const cln::cl_N x_ = x.to_cl_N();
2063 const cln::cl_N result = lgamma(x_);
2064 return numeric(result);
2065}
2066
2067const cln::cl_N tgamma(const cln::cl_N &x)
2068{
2069 cln::float_format_t prec = guess_precision(x);
2070 lanczos_coeffs lc;
2071 if (lc.sufficiently_accurate(prec)) {
2072 cln::cl_N pi_val = cln::pi(prec);
2073 if (realpart(x) < 0.5)
2074 return pi_val/(cln::sin(pi_val*x))/tgamma(1 - x);
2075 cln::cl_N A = lc.calc_lanczos_A(x);
2076 cln::cl_N temp = x + lc.get_order() - cln::cl_N(1)/2;
2077 cln::cl_N result = sqrt(cln::cl_I(2)*pi_val)
2078 * expt(temp, x - cln::cl_N(1)/2)
2079 * exp(-temp) * A;
2080 return result;
2081 }
2082 else
2083 throw dunno();
2084}
2085
2087{
2088 const cln::cl_N x_ = x.to_cl_N();
2089 const cln::cl_N result = tgamma(x_);
2090 return numeric(result);
2091}
2092
2095const numeric psi(const numeric &x)
2096{
2097 throw dunno();
2098}
2099
2100
2103const numeric psi(const numeric &n, const numeric &x)
2104{
2105 throw dunno();
2106}
2107
2108
2114{
2115 if (!n.is_nonneg_integer())
2116 throw std::range_error("numeric::factorial(): argument must be integer >= 0");
2117 return numeric(cln::factorial(n.to_int()));
2118}
2119
2120
2128{
2129 if (n.is_equal(*_num_1_p))
2130 return *_num1_p;
2131
2132 if (!n.is_nonneg_integer())
2133 throw std::range_error("numeric::doublefactorial(): argument must be integer >= -1");
2134
2135 return numeric(cln::doublefactorial(n.to_int()));
2136}
2137
2138
2145const numeric binomial(const numeric &n, const numeric &k)
2146{
2147 if (n.is_integer() && k.is_integer()) {
2148 if (n.is_nonneg_integer()) {
2149 if (k.compare(n)!=1 && k.compare(*_num0_p)!=-1)
2150 return numeric(cln::binomial(n.to_int(),k.to_int()));
2151 else
2152 return *_num0_p;
2153 } else {
2154 if (k.is_nonneg_integer())
2155 return _num_1_p->power(k)*binomial(k-n-(*_num1_p), k);
2156 else
2157 return _num_1_p->power(n-k)*binomial(-k-(*_num1_p), n-k);
2158 }
2159 }
2160
2161 // should really be gamma(n+1)/gamma(k+1)/gamma(n-k+1) or a suitable limit
2162 throw std::range_error("numeric::binomial(): don't know how to evaluate that.");
2163}
2164
2165
2171const numeric bernoulli(const numeric &nn)
2172{
2173 if (!nn.is_integer() || nn.is_negative())
2174 throw std::range_error("numeric::bernoulli(): argument must be integer >= 0");
2175
2176 // Method:
2177 //
2178 // The Bernoulli numbers are rational numbers that may be computed using
2179 // the relation
2180 //
2181 // B_n = - 1/(n+1) * sum_{k=0}^{n-1}(binomial(n+1,k)*B_k)
2182 //
2183 // with B(0) = 1. Since the n'th Bernoulli number depends on all the
2184 // previous ones, the computation is necessarily very expensive. There are
2185 // several other ways of computing them, a particularly good one being
2186 // cl_I s = 1;
2187 // cl_I c = n+1;
2188 // cl_RA Bern = 0;
2189 // for (unsigned i=0; i<n; i++) {
2190 // c = exquo(c*(i-n),(i+2));
2191 // Bern = Bern + c*s/(i+2);
2192 // s = s + expt_pos(cl_I(i+2),n);
2193 // }
2194 // return Bern;
2195 //
2196 // But if somebody works with the n'th Bernoulli number she is likely to
2197 // also need all previous Bernoulli numbers. So we need a complete remember
2198 // table and above divide and conquer algorithm is not suited to build one
2199 // up. The formula below accomplishes this. It is a modification of the
2200 // defining formula above but the computation of the binomial coefficients
2201 // is carried along in an inline fashion. It also honors the fact that
2202 // B_n is zero when n is odd and greater than 1.
2203 //
2204 // (There is an interesting relation with the tangent polynomials described
2205 // in `Concrete Mathematics', which leads to a program a little faster as
2206 // our implementation below, but it requires storing one such polynomial in
2207 // addition to the remember table. This doubles the memory footprint so
2208 // we don't use it.)
2209
2210 const unsigned n = nn.to_int();
2211
2212 // the special cases not covered by the algorithm below
2213 if (n & 1)
2214 return (n==1) ? (*_num_1_2_p) : (*_num0_p);
2215 if (!n)
2216 return *_num1_p;
2217
2218 // store nonvanishing Bernoulli numbers here
2219 static std::vector< cln::cl_RA > results;
2220 static unsigned next_r = 0;
2221
2222 // algorithm not applicable to B(2), so just store it
2223 if (!next_r) {
2224 results.push_back(cln::recip(cln::cl_RA(6)));
2225 next_r = 4;
2226 }
2227 if (n<next_r)
2228 return numeric(results[n/2-1]);
2229
2230 results.reserve(n/2);
2231 for (unsigned p=next_r; p<=n; p+=2) {
2232 cln::cl_I c = 1; // seed for binomial coefficients
2233 cln::cl_RA b = cln::cl_RA(p-1)/-2;
2234 // The CLN manual says: "The conversion from `unsigned int' works only
2235 // if the argument is < 2^29" (This is for 32 Bit machines. More
2236 // generally, cl_value_len is the limiting exponent of 2. We must make
2237 // sure that no intermediates are created which exceed this value. The
2238 // largest intermediate is (p+3-2*k)*(p/2-k+1) <= (p^2+p)/2.
2239 if (p < (1UL<<cl_value_len/2)) {
2240 for (unsigned k=1; k<=p/2-1; ++k) {
2241 c = cln::exquo(c * ((p+3-2*k) * (p/2-k+1)), (2*k-1)*k);
2242 b = b + c*results[k-1];
2243 }
2244 } else {
2245 for (unsigned k=1; k<=p/2-1; ++k) {
2246 c = cln::exquo((c * (p+3-2*k)) * (p/2-k+1), cln::cl_I(2*k-1)*k);
2247 b = b + c*results[k-1];
2248 }
2249 }
2250 results.push_back(-b/(p+1));
2251 }
2252 next_r = n+2;
2253 return numeric(results[n/2-1]);
2254}
2255
2256
2264{
2265 if (!n.is_integer())
2266 throw std::range_error("numeric::fibonacci(): argument must be integer");
2267 // Method:
2268 //
2269 // The following addition formula holds:
2270 //
2271 // F(n+m) = F(m-1)*F(n) + F(m)*F(n+1) for m >= 1, n >= 0.
2272 //
2273 // (Proof: For fixed m, the LHS and the RHS satisfy the same recurrence
2274 // w.r.t. n, and the initial values (n=0, n=1) agree. Hence all values
2275 // agree.)
2276 // Replace m by m+1:
2277 // F(n+m+1) = F(m)*F(n) + F(m+1)*F(n+1) for m >= 0, n >= 0
2278 // Now put in m = n, to get
2279 // F(2n) = (F(n+1)-F(n))*F(n) + F(n)*F(n+1) = F(n)*(2*F(n+1) - F(n))
2280 // F(2n+1) = F(n)^2 + F(n+1)^2
2281 // hence
2282 // F(2n+2) = F(n+1)*(2*F(n) + F(n+1))
2283 if (n.is_zero())
2284 return *_num0_p;
2285 if (n.is_negative()) {
2286 if (n.is_even()) {
2287 return -fibonacci(-n);
2288 }
2289 else {
2290 return fibonacci(-n);
2291 }
2292 }
2293
2294 cln::cl_I u(0);
2295 cln::cl_I v(1);
2296 cln::cl_I m = cln::the<cln::cl_I>(n.to_cl_N()) >> 1L; // floor(n/2);
2297 for (uintL bit=cln::integer_length(m); bit>0; --bit) {
2298 // Since a squaring is cheaper than a multiplication, better use
2299 // three squarings instead of one multiplication and two squarings.
2300 cln::cl_I u2 = cln::square(u);
2301 cln::cl_I v2 = cln::square(v);
2302 if (cln::logbitp(bit-1, m)) {
2303 v = cln::square(u + v) - u2;
2304 u = u2 + v2;
2305 } else {
2306 u = v2 - cln::square(v - u);
2307 v = u2 + v2;
2308 }
2309 }
2310 if (n.is_even())
2311 // Here we don't use the squaring formula because one multiplication
2312 // is cheaper than two squarings.
2313 return numeric(u * ((v << 1) - u));
2314 else
2315 return numeric(cln::square(u) + cln::square(v));
2316}
2317
2318
2320const numeric abs(const numeric& x)
2321{
2322 return numeric(cln::abs(x.to_cl_N()));
2323}
2324
2325
2333const numeric mod(const numeric &a, const numeric &b)
2334{
2335 if (a.is_integer() && b.is_integer())
2336 return numeric(cln::mod(cln::the<cln::cl_I>(a.to_cl_N()),
2337 cln::the<cln::cl_I>(b.to_cl_N())));
2338 else
2339 return *_num0_p;
2340}
2341
2342
2346const numeric smod(const numeric &a_, const numeric &b_)
2347{
2348 if (a_.is_integer() && b_.is_integer()) {
2349 const cln::cl_I a = cln::the<cln::cl_I>(a_.to_cl_N());
2350 const cln::cl_I b = cln::the<cln::cl_I>(b_.to_cl_N());
2351 const cln::cl_I b2 = b >> 1;
2352 const cln::cl_I m = cln::mod(a, b);
2353 const cln::cl_I m_b = m - b;
2354 const cln::cl_I ret = m > b2 ? m_b : m;
2355 return numeric(ret);
2356 } else
2357 return *_num0_p;
2358}
2359
2360
2368const numeric irem(const numeric &a, const numeric &b)
2369{
2370 if (b.is_zero())
2371 throw std::overflow_error("numeric::irem(): division by zero");
2372 if (a.is_integer() && b.is_integer())
2373 return numeric(cln::rem(cln::the<cln::cl_I>(a.to_cl_N()),
2374 cln::the<cln::cl_I>(b.to_cl_N())));
2375 else
2376 return *_num0_p;
2377}
2378
2379
2388const numeric irem(const numeric &a, const numeric &b, numeric &q)
2389{
2390 if (b.is_zero())
2391 throw std::overflow_error("numeric::irem(): division by zero");
2392 if (a.is_integer() && b.is_integer()) {
2393 const cln::cl_I_div_t rem_quo = cln::truncate2(cln::the<cln::cl_I>(a.to_cl_N()),
2394 cln::the<cln::cl_I>(b.to_cl_N()));
2395 q = numeric(rem_quo.quotient);
2396 return numeric(rem_quo.remainder);
2397 } else {
2398 q = *_num0_p;
2399 return *_num0_p;
2400 }
2401}
2402
2403
2409const numeric iquo(const numeric &a, const numeric &b)
2410{
2411 if (b.is_zero())
2412 throw std::overflow_error("numeric::iquo(): division by zero");
2413 if (a.is_integer() && b.is_integer())
2414 return numeric(cln::truncate1(cln::the<cln::cl_I>(a.to_cl_N()),
2415 cln::the<cln::cl_I>(b.to_cl_N())));
2416 else
2417 return *_num0_p;
2418}
2419
2420
2428const numeric iquo(const numeric &a, const numeric &b, numeric &r)
2429{
2430 if (b.is_zero())
2431 throw std::overflow_error("numeric::iquo(): division by zero");
2432 if (a.is_integer() && b.is_integer()) {
2433 const cln::cl_I_div_t rem_quo = cln::truncate2(cln::the<cln::cl_I>(a.to_cl_N()),
2434 cln::the<cln::cl_I>(b.to_cl_N()));
2435 r = numeric(rem_quo.remainder);
2436 return numeric(rem_quo.quotient);
2437 } else {
2438 r = *_num0_p;
2439 return *_num0_p;
2440 }
2441}
2442
2443
2448const numeric gcd(const numeric &a, const numeric &b)
2449{
2450 if (a.is_integer() && b.is_integer())
2451 return numeric(cln::gcd(cln::the<cln::cl_I>(a.to_cl_N()),
2452 cln::the<cln::cl_I>(b.to_cl_N())));
2453 else
2454 return *_num1_p;
2455}
2456
2457
2462const numeric lcm(const numeric &a, const numeric &b)
2463{
2464 if (a.is_integer() && b.is_integer())
2465 return numeric(cln::lcm(cln::the<cln::cl_I>(a.to_cl_N()),
2466 cln::the<cln::cl_I>(b.to_cl_N())));
2467 else
2468 return a.mul(b);
2469}
2470
2471
2480const numeric sqrt(const numeric &x)
2481{
2482 return numeric(cln::sqrt(x.to_cl_N()));
2483}
2484
2485
2487const numeric isqrt(const numeric &x)
2488{
2489 if (x.is_integer()) {
2490 cln::cl_I root;
2491 cln::isqrt(cln::the<cln::cl_I>(x.to_cl_N()), &root);
2492 return numeric(root);
2493 } else
2494 return *_num0_p;
2495}
2496
2497
2500{
2501 return numeric(cln::pi(cln::default_float_format));
2502}
2503
2504
2507{
2508 return numeric(cln::eulerconst(cln::default_float_format));
2509}
2510
2511
2514{
2515 return numeric(cln::catalanconst(cln::default_float_format));
2516}
2517
2518
2521 : digits(17)
2522{
2523 // It initializes to 17 digits, because in CLN float_format(17) turns out
2524 // to be 61 (<64) while float_format(18)=65. The reason is we want to
2525 // have a cl_LF instead of cl_SF, cl_FF or cl_DF.
2526 if (too_late)
2527 throw(std::runtime_error("I told you not to do instantiate me!"));
2528 too_late = true;
2529 cln::default_float_format = cln::float_format(17);
2530
2531 // add callbacks for built-in functions
2532 // like ... add_callback(Li_lookuptable);
2533}
2534
2535
2538{
2539 long digitsdiff = prec - digits;
2540 digits = prec;
2541 cln::default_float_format = cln::float_format(prec);
2542
2543 // call registered callbacks
2544 for (auto it : callbacklist) {
2545 (it)(digitsdiff);
2546 }
2547
2548 return *this;
2549}
2550
2551
2553_numeric_digits::operator long()
2554{
2555 // BTW, this is approx. unsigned(cln::default_float_format*0.301)-1
2556 return (long)digits;
2557}
2558
2559
2561void _numeric_digits::print(std::ostream &os) const
2562{
2563 os << digits;
2564}
2565
2566
2569{
2570 callbacklist.push_back(callback);
2571}
2572
2573
2574std::ostream& operator<<(std::ostream &os, const _numeric_digits &e)
2575{
2576 e.print(os);
2577 return os;
2578}
2579
2581// static member variables
2583
2584// private
2585
2586bool _numeric_digits::too_late = false;
2587
2588
2592
2593} // namespace GiNaC
Archiving of GiNaC expressions.
#define GINAC_ASSERT(X)
Assertion macro for checking invariances.
Definition: assertion.h:33
This class is used to instantiate a global singleton object Digits which behaves just like Maple's Di...
Definition: numeric.h:52
_numeric_digits & operator=(long prec)
Assign a native long to global Digits object.
Definition: numeric.cpp:2537
_numeric_digits()
_numeric_digits default ctor, checking for singleton invariance.
Definition: numeric.cpp:2520
std::vector< digits_changed_callback > callbacklist
Definition: numeric.h:65
void print(std::ostream &os) const
Append global Digits object to ostream.
Definition: numeric.cpp:2561
static bool too_late
Already one object present.
Definition: numeric.h:63
void add_callback(digits_changed_callback callback)
Add a new callback function.
Definition: numeric.cpp:2568
long digits
Number of decimal digits.
Definition: numeric.h:62
This class stores all properties needed to record/retrieve the state of one object of class basic (or...
Definition: archive.h:49
This class is the ABC (abstract base class) of GiNaC's class hierarchy.
Definition: basic.h:105
const basic & setflag(unsigned f) const
Set some status_flags.
Definition: basic.h:288
unsigned hashvalue
hash value
Definition: basic.h:303
unsigned flags
of type status_flags
Definition: basic.h:302
const basic & hold() const
Stop further evaluation.
Definition: basic.cpp:887
virtual int compare_same_type(const basic &other) const
Returns order relation between two objects of same type.
Definition: basic.cpp:719
Wrapper template for making GiNaC classes out of STL containers.
Definition: container.h:73
Exception class thrown by functions to signal unimplemented functionality so the expression may just ...
Definition: utils.h:37
Lightweight wrapper for GiNaC's symbolic objects.
Definition: ex.h:72
bool is_zero() const
Definition: ex.h:213
@ integer_polynomial
Definition: flags.h:256
@ cinteger_polynomial
Definition: flags.h:257
@ crational_polynomial
Definition: flags.h:259
@ rational_polynomial
Definition: flags.h:258
bool sufficiently_accurate(int digits)
Definition: numeric.cpp:1751
int get_order() const
Definition: numeric.cpp:1737
cln::cl_N calc_lanczos_A(const cln::cl_N &) const
Definition: numeric.cpp:1771
std::vector< cln::cl_N > * current_vector
Definition: numeric.cpp:1746
static std::vector< cln::cl_N > * coeffs
Definition: numeric.cpp:1744
This class is a wrapper around CLN-numbers within the GiNaC class hierarchy.
Definition: numeric.h:82
void do_print(const print_context &c, unsigned level) const
Definition: numeric.cpp:605
bool is_pos_integer() const
True if object is an exact integer greater than zero.
Definition: numeric.cpp:1161
const numeric & operator=(int i)
Definition: numeric.cpp:1016
void read_archive(const archive_node &n, lst &syms) override
Read (a.k.a.
Definition: numeric.cpp:289
const numeric & sub_dyn(const numeric &other) const
Numerical subtraction method.
Definition: numeric.cpp:942
unsigned calchash() const override
Compute the hash value of an object and if it makes sense to store it in the objects status_flags,...
Definition: numeric.cpp:838
bool is_equal_same_type(const basic &other) const override
Returns true if two objects of same type are equal.
Definition: numeric.cpp:829
const numeric sub(const numeric &other) const
Numerical subtraction method.
Definition: numeric.cpp:872
const numeric & mul_dyn(const numeric &other) const
Numerical multiplication method.
Definition: numeric.cpp:957
void do_print_csrc(const print_csrc &c, unsigned level) const
Definition: numeric.cpp:615
bool is_cinteger() const
True if object is element of the domain of integers extended by I, i.e.
Definition: numeric.cpp:1228
unsigned precedence() const override
Return relative operator precedence (for parenthezing output).
Definition: numeric.h:101
numeric(int i)
Definition: numeric.cpp:85
bool is_polynomial(const ex &var) const override
Check whether this is a polynomial in the given variables.
Definition: numeric.cpp:728
numeric step() const
Return the step function of a numeric.
Definition: numeric.cpp:1064
bool is_crational() const
True if object is an exact rational number, may even be complex (denominator may be unity).
Definition: numeric.cpp:1243
ex imag_part() const override
Definition: numeric.cpp:813
bool is_rational() const
True if object is an exact rational number, may even be complex (denominator may be unity).
Definition: numeric.cpp:1201
bool operator>(const numeric &other) const
Numerical comparison: greater.
Definition: numeric.cpp:1281
void archive(archive_node &n) const override
Save (a.k.a.
Definition: numeric.cpp:344
bool info(unsigned inf) const override
Information about the object.
Definition: numeric.cpp:684
ex real_part() const override
Definition: numeric.cpp:808
int ldegree(const ex &s) const override
Return degree of lowest power in object s.
Definition: numeric.cpp:738
const numeric real() const
Real part of a number.
Definition: numeric.cpp:1339
ex eval() const override
Evaluation of numbers doesn't do anything at all.
Definition: numeric.cpp:783
bool is_prime() const
Probabilistic primality test.
Definition: numeric.cpp:1191
bool has(const ex &other, unsigned options=0) const override
Disassemble real part and imaginary part to scan for the occurrence of a single number.
Definition: numeric.cpp:754
long to_long() const
Converts numeric types to machine's long.
Definition: numeric.cpp:1313
void do_print_latex(const print_latex &c, unsigned level) const
Definition: numeric.cpp:610
void do_print_tree(const print_tree &c, unsigned level) const
Definition: numeric.cpp:669
ex coeff(const ex &s, int n=1) const override
Return coefficient of degree n in object s.
Definition: numeric.cpp:743
const numeric & power_dyn(const numeric &other) const
Numerical exponentiation.
Definition: numeric.cpp:993
void do_print_python_repr(const print_python_repr &c, unsigned level) const
Definition: numeric.cpp:677
int compare(const numeric &other) const
This method establishes a canonical order on all numbers.
Definition: numeric.cpp:1104
bool is_nonneg_integer() const
True if object is an exact integer greater or equal zero.
Definition: numeric.cpp:1168
bool is_positive() const
True if object is not complex and greater than zero.
Definition: numeric.cpp:1136
ex conjugate() const override
Definition: numeric.cpp:800
bool is_real() const
True if object is a real integer, rational or float (but not complex).
Definition: numeric.cpp:1208
const numeric numer() const
Numerator.
Definition: numeric.cpp:1356
cln::cl_N value
Definition: numeric.h:200
bool is_integer() const
True if object is a non-complex integer.
Definition: numeric.cpp:1154
const numeric power(const numeric &other) const
Numerical exponentiation.
Definition: numeric.cpp:900
ex evalf() const override
Cast numeric into a floating-point object.
Definition: numeric.cpp:795
const numeric denom() const
Denominator.
Definition: numeric.cpp:1387
bool is_negative() const
True if object is not complex and less than zero.
Definition: numeric.cpp:1145
bool is_odd() const
True if object is an exact odd integer.
Definition: numeric.cpp:1182
cln::cl_N to_cl_N() const
Returns a new CLN object of type cl_N, representing the value of *this.
Definition: numeric.cpp:1332
const numeric imag() const
Imaginary part of a number.
Definition: numeric.cpp:1346
const numeric mul(const numeric &other) const
Numerical multiplication method.
Definition: numeric.cpp:880
bool is_even() const
True if object is an exact even integer.
Definition: numeric.cpp:1175
const numeric & add_dyn(const numeric &other) const
Numerical addition method.
Definition: numeric.cpp:925
int degree(const ex &s) const override
Return degree of highest power in object s.
Definition: numeric.cpp:733
bool operator<=(const numeric &other) const
Numerical comparison: less or equal.
Definition: numeric.cpp:1270
int csgn() const
Return the complex half-plane (left or right) in which the number lies.
Definition: numeric.cpp:1078
bool operator==(const numeric &other) const
Definition: numeric.cpp:1214
bool is_equal(const numeric &other) const
Definition: numeric.cpp:1122
const numeric & div_dyn(const numeric &other) const
Numerical division method.
Definition: numeric.cpp:976
bool operator>=(const numeric &other) const
Numerical comparison: greater or equal.
Definition: numeric.cpp:1292
bool operator!=(const numeric &other) const
Definition: numeric.cpp:1220
int to_int() const
Converts numeric types to machine's int.
Definition: numeric.cpp:1303
int int_length() const
Size in binary notation.
Definition: numeric.cpp:1418
void print_numeric(const print_context &c, const char *par_open, const char *par_close, const char *imag_sym, const char *mul_sym, unsigned level) const
Definition: numeric.cpp:542
void do_print_csrc_cl_N(const print_csrc_cl_N &c, unsigned level) const
Definition: numeric.cpp:651
double to_double() const
Converts numeric types to machine's double.
Definition: numeric.cpp:1322
const numeric inverse() const
Inverse of a number.
Definition: numeric.cpp:1053
bool operator<(const numeric &other) const
Numerical comparison: less.
Definition: numeric.cpp:1259
const numeric add(const numeric &other) const
Numerical addition method.
Definition: numeric.cpp:864
bool is_zero() const
True if object is zero.
Definition: numeric.cpp:1129
const numeric div(const numeric &other) const
Numerical division method.
Definition: numeric.cpp:890
Exception class thrown when a singularity is encountered.
Definition: numeric.h:70
Base class for print_contexts.
Definition: print.h:103
std::ostream & s
stream to output to
Definition: print.h:109
Context for C source output using CLN numbers.
Definition: print.h:182
Base context for C source output.
Definition: print.h:158
Context for latex-parsable output.
Definition: print.h:123
Context for python-parsable output.
Definition: print.h:139
Context for tree-like output for debugging.
Definition: print.h:147
@ expanded
.expand(0) has already done its job (other expand() options ignore this flag)
Definition: flags.h:204
@ evaluated
.eval() has already done its job
Definition: flags.h:203
@ hash_calculated
.calchash() has already done its job
Definition: flags.h:205
Interface to GiNaC's light-weight expression handles.
static const bool value
Definition: factor.cpp:199
unsigned options
Definition: factor.cpp:2475
vector< int > k
Definition: factor.cpp:1435
size_t n
Definition: factor.cpp:1432
size_t c
Definition: factor.cpp:757
ex x
Definition: factor.cpp:1610
size_t r
Definition: factor.cpp:757
mvec m
Definition: factor.cpp:758
Definition: add.cpp:38
const numeric gcd(const numeric &a, const numeric &b)
Greatest Common Divisor.
Definition: numeric.cpp:2448
const numeric I
Imaginary unit.
Definition: numeric.cpp:1433
const numeric atan(const numeric &x)
Numeric arcustangent.
Definition: numeric.cpp:1508
const numeric * _num_1_p
Definition: utils.cpp:351
std::ostream & operator<<(std::ostream &os, const archive_node &n)
Write archive_node to binary data stream.
Definition: archive.cpp:200
ex denom(const ex &thisex)
Definition: ex.h:763
const numeric atan(const numeric &y, const numeric &x)
Numeric arcustangent of two arguments, analytically continued in a suitable way.
Definition: numeric.cpp:1525
unsigned golden_ratio_hash(uintptr_t n)
Truncated multiplication with golden ratio, for computing hash values.
Definition: utils.h:68
const numeric zeta(const numeric &x)
Numeric evaluation of Riemann's Zeta function.
Definition: numeric.cpp:1717
const numeric bernoulli(const numeric &nn)
Bernoulli number.
Definition: numeric.cpp:2171
const numeric cosh(const numeric &x)
Numeric hyperbolic cosine (trigonometric function).
Definition: numeric.cpp:1563
const numeric mod(const numeric &a, const numeric &b)
Modulus (in positive representation).
Definition: numeric.cpp:2333
const numeric abs(const numeric &x)
Absolute value.
Definition: numeric.cpp:2320
ex EulerEvalf()
Floating point evaluation of Euler's constant gamma.
Definition: numeric.cpp:2506
static void print_real_csrc(const print_context &c, const cln::cl_R &x)
Helper function to print real number in C++ source format.
Definition: numeric.cpp:440
const numeric asin(const numeric &x)
Numeric inverse sine (trigonometric function).
Definition: numeric.cpp:1488
function zeta(const T1 &p1)
Definition: inifcns.h:111
ex PiEvalf()
Floating point evaluation of Archimedes' constant Pi.
Definition: numeric.cpp:2499
const numeric fibonacci(const numeric &n)
Fibonacci number.
Definition: numeric.cpp:2263
const numeric doublefactorial(const numeric &n)
The double factorial combinatorial function.
Definition: numeric.cpp:2127
const numeric tanh(const numeric &x)
Numeric hyperbolic tangent (trigonometric function).
Definition: numeric.cpp:1572
const numeric Li2(const numeric &x)
Definition: numeric.cpp:1705
const numeric acos(const numeric &x)
Numeric inverse cosine (trigonometric function).
Definition: numeric.cpp:1497
ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned options)
Compute GCD (Greatest Common Divisor) of multivariate polynomials a(X) and b(X) in Z[X].
Definition: normal.cpp:1433
function psi(const T1 &p1)
Definition: inifcns.h:165
const numeric sqrt(const numeric &x)
Numeric square root.
Definition: numeric.cpp:2480
const numeric irem(const numeric &a, const numeric &b)
Numeric integer remainder.
Definition: numeric.cpp:2368
ex conjugate(const ex &thisex)
Definition: ex.h:733
const cln::cl_N tgamma(const cln::cl_N &x)
Definition: numeric.cpp:2067
const numeric sinh(const numeric &x)
Numeric hyperbolic sine (trigonometric function).
Definition: numeric.cpp:1554
void(* digits_changed_callback)(long)
Function pointer to implement callbacks in the case 'Digits' gets changed.
Definition: numeric.h:40
static const cln::cl_F make_real_float(const cln::cl_idecoded_float &dec)
Construct a floating point number from sign, mantissa, and exponent.
Definition: numeric.cpp:269
const numeric binomial(const numeric &n, const numeric &k)
The Binomial coefficients.
Definition: numeric.cpp:2145
const numeric exp(const numeric &x)
Exponential function.
Definition: numeric.cpp:1439
static cln::cl_N Li2_projection(const cln::cl_N &x, const cln::float_format_t &prec)
Folds Li2's argument inside a small rectangle to enhance convergence.
Definition: numeric.cpp:1652
const numeric factorial(const numeric &n)
Factorial combinatorial function.
Definition: numeric.cpp:2113
const numeric acosh(const numeric &x)
Numeric inverse hyperbolic cosine (trigonometric function).
Definition: numeric.cpp:1590
const numeric cos(const numeric &x)
Numeric cosine (trigonometric function).
Definition: numeric.cpp:1470
const numeric smod(const numeric &a_, const numeric &b_)
Modulus (in symmetric representation).
Definition: numeric.cpp:2346
const numeric atanh(const numeric &x)
Numeric inverse hyperbolic tangent (trigonometric function).
Definition: numeric.cpp:1599
print_func< print_context >(&varidx::do_print). print_func< print_latex >(&varidx
Definition: idx.cpp:45
const cln::cl_N Li2_(const cln::cl_N &value)
Numeric evaluation of Dilogarithm.
Definition: numeric.cpp:1679
GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(add, expairseq, print_func< print_context >(&add::do_print). print_func< print_latex >(&add::do_print_latex). print_func< print_csrc >(&add::do_print_csrc). print_func< print_tree >(&add::do_print_tree). print_func< print_python_repr >(&add::do_print_python_repr)) add
Definition: add.cpp:40
ex lcm(const ex &a, const ex &b, bool check_args)
Compute LCM (Least Common Multiple) of multivariate polynomials in Z[X].
Definition: normal.cpp:1775
const numeric iquo(const numeric &a, const numeric &b)
Numeric integer quotient.
Definition: numeric.cpp:2409
const numeric isqrt(const numeric &x)
Integer numeric square root.
Definition: numeric.cpp:2487
const numeric log(const numeric &x)
Natural logarithm.
Definition: numeric.cpp:1450
static void print_integer_csrc(const print_context &c, const cln::cl_I &x)
Helper function to print integer number in C++ source format.
Definition: numeric.cpp:426
_numeric_digits Digits
Accuracy in decimal digits.
Definition: numeric.cpp:2591
const numeric sin(const numeric &x)
Numeric sine (trigonometric function).
Definition: numeric.cpp:1461
static void print_real_number(const print_context &c, const cln::cl_R &x)
Helper function to print a real number in a nicer way than is CLN's default.
Definition: numeric.cpp:397
static void print_real_cl_N(const print_context &c, const cln::cl_R &x)
Helper function to print real number in C++ source format using cl_N types.
Definition: numeric.cpp:507
const numeric * _num1_p
Definition: utils.cpp:384
static void write_real_float(std::ostream &s, const cln::cl_R &n)
Definition: numeric.cpp:338
ex numer(const ex &thisex)
Definition: ex.h:760
ex CatalanEvalf()
Floating point evaluation of Catalan's constant.
Definition: numeric.cpp:2513
static cln::float_format_t guess_precision(const cln::cl_N &x)
Definition: numeric.cpp:2024
ex rem(const ex &a, const ex &b, const ex &x, bool check_args)
Remainder r(x) of polynomials a(x) and b(x) in Q[x].
Definition: normal.cpp:423
const numeric asinh(const numeric &x)
Numeric inverse hyperbolic sine (trigonometric function).
Definition: numeric.cpp:1581
const numeric tan(const numeric &x)
Numeric tangent (trigonometric function).
Definition: numeric.cpp:1479
const numeric lcm(const numeric &a, const numeric &b)
Least Common Multiple.
Definition: numeric.cpp:2462
GINAC_IMPLEMENT_REGISTERED_CLASS_OPT_T(lst, basic, print_func< print_context >(&lst::do_print). print_func< print_tree >(&lst::do_print_tree)) template<> bool lst GINAC_BIND_UNARCHIVER(lst)
Specialization of container::info() for lst.
Definition: lst.cpp:42
static const cln::cl_F read_real_float(std::istream &s)
Read serialized floating point number.
Definition: numeric.cpp:281
static bool coerce(T1 &dst, const T2 &arg)
const cln::cl_N lgamma(const cln::cl_N &x)
The Gamma function.
Definition: numeric.cpp:2039
const ex _ex0
Definition: utils.cpp:369
static ex Li2_series(const ex &x, const relational &rel, int order, unsigned options)
Definition: inifcns.cpp:718
const numeric * _num0_p
Definition: utils.cpp:367
Makes the interface to the underlying bignum package available.
Interface to GiNaC's overloaded operators.
Interface to several small and furry utilities needed within GiNaC but not of any interest to the use...

This page is part of the GiNaC developer's reference. It was generated automatically by doxygen. For an introduction, see the tutorial.