GiNaC 1.8.5
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
2143const numeric binomial(const numeric &n, const numeric &k)
2144{
2145 if (n.is_integer() && k.is_integer()) {
2146 if (n.is_nonneg_integer()) {
2147 if (k.compare(n)!=1 && k.compare(*_num0_p)!=-1)
2148 return numeric(cln::binomial(n.to_int(),k.to_int()));
2149 else
2150 return *_num0_p;
2151 } else {
2152 return _num_1_p->power(k)*binomial(k-n-(*_num1_p),k);
2153 }
2154 }
2155
2156 // should really be gamma(n+1)/gamma(k+1)/gamma(n-k+1) or a suitable limit
2157 throw std::range_error("numeric::binomial(): don't know how to evaluate that.");
2158}
2159
2160
2166const numeric bernoulli(const numeric &nn)
2167{
2168 if (!nn.is_integer() || nn.is_negative())
2169 throw std::range_error("numeric::bernoulli(): argument must be integer >= 0");
2170
2171 // Method:
2172 //
2173 // The Bernoulli numbers are rational numbers that may be computed using
2174 // the relation
2175 //
2176 // B_n = - 1/(n+1) * sum_{k=0}^{n-1}(binomial(n+1,k)*B_k)
2177 //
2178 // with B(0) = 1. Since the n'th Bernoulli number depends on all the
2179 // previous ones, the computation is necessarily very expensive. There are
2180 // several other ways of computing them, a particularly good one being
2181 // cl_I s = 1;
2182 // cl_I c = n+1;
2183 // cl_RA Bern = 0;
2184 // for (unsigned i=0; i<n; i++) {
2185 // c = exquo(c*(i-n),(i+2));
2186 // Bern = Bern + c*s/(i+2);
2187 // s = s + expt_pos(cl_I(i+2),n);
2188 // }
2189 // return Bern;
2190 //
2191 // But if somebody works with the n'th Bernoulli number she is likely to
2192 // also need all previous Bernoulli numbers. So we need a complete remember
2193 // table and above divide and conquer algorithm is not suited to build one
2194 // up. The formula below accomplishes this. It is a modification of the
2195 // defining formula above but the computation of the binomial coefficients
2196 // is carried along in an inline fashion. It also honors the fact that
2197 // B_n is zero when n is odd and greater than 1.
2198 //
2199 // (There is an interesting relation with the tangent polynomials described
2200 // in `Concrete Mathematics', which leads to a program a little faster as
2201 // our implementation below, but it requires storing one such polynomial in
2202 // addition to the remember table. This doubles the memory footprint so
2203 // we don't use it.)
2204
2205 const unsigned n = nn.to_int();
2206
2207 // the special cases not covered by the algorithm below
2208 if (n & 1)
2209 return (n==1) ? (*_num_1_2_p) : (*_num0_p);
2210 if (!n)
2211 return *_num1_p;
2212
2213 // store nonvanishing Bernoulli numbers here
2214 static std::vector< cln::cl_RA > results;
2215 static unsigned next_r = 0;
2216
2217 // algorithm not applicable to B(2), so just store it
2218 if (!next_r) {
2219 results.push_back(cln::recip(cln::cl_RA(6)));
2220 next_r = 4;
2221 }
2222 if (n<next_r)
2223 return numeric(results[n/2-1]);
2224
2225 results.reserve(n/2);
2226 for (unsigned p=next_r; p<=n; p+=2) {
2227 cln::cl_I c = 1; // seed for binomial coefficients
2228 cln::cl_RA b = cln::cl_RA(p-1)/-2;
2229 // The CLN manual says: "The conversion from `unsigned int' works only
2230 // if the argument is < 2^29" (This is for 32 Bit machines. More
2231 // generally, cl_value_len is the limiting exponent of 2. We must make
2232 // sure that no intermediates are created which exceed this value. The
2233 // largest intermediate is (p+3-2*k)*(p/2-k+1) <= (p^2+p)/2.
2234 if (p < (1UL<<cl_value_len/2)) {
2235 for (unsigned k=1; k<=p/2-1; ++k) {
2236 c = cln::exquo(c * ((p+3-2*k) * (p/2-k+1)), (2*k-1)*k);
2237 b = b + c*results[k-1];
2238 }
2239 } else {
2240 for (unsigned k=1; k<=p/2-1; ++k) {
2241 c = cln::exquo((c * (p+3-2*k)) * (p/2-k+1), cln::cl_I(2*k-1)*k);
2242 b = b + c*results[k-1];
2243 }
2244 }
2245 results.push_back(-b/(p+1));
2246 }
2247 next_r = n+2;
2248 return numeric(results[n/2-1]);
2249}
2250
2251
2259{
2260 if (!n.is_integer())
2261 throw std::range_error("numeric::fibonacci(): argument must be integer");
2262 // Method:
2263 //
2264 // The following addition formula holds:
2265 //
2266 // F(n+m) = F(m-1)*F(n) + F(m)*F(n+1) for m >= 1, n >= 0.
2267 //
2268 // (Proof: For fixed m, the LHS and the RHS satisfy the same recurrence
2269 // w.r.t. n, and the initial values (n=0, n=1) agree. Hence all values
2270 // agree.)
2271 // Replace m by m+1:
2272 // F(n+m+1) = F(m)*F(n) + F(m+1)*F(n+1) for m >= 0, n >= 0
2273 // Now put in m = n, to get
2274 // F(2n) = (F(n+1)-F(n))*F(n) + F(n)*F(n+1) = F(n)*(2*F(n+1) - F(n))
2275 // F(2n+1) = F(n)^2 + F(n+1)^2
2276 // hence
2277 // F(2n+2) = F(n+1)*(2*F(n) + F(n+1))
2278 if (n.is_zero())
2279 return *_num0_p;
2280 if (n.is_negative()) {
2281 if (n.is_even()) {
2282 return -fibonacci(-n);
2283 }
2284 else {
2285 return fibonacci(-n);
2286 }
2287 }
2288
2289 cln::cl_I u(0);
2290 cln::cl_I v(1);
2291 cln::cl_I m = cln::the<cln::cl_I>(n.to_cl_N()) >> 1L; // floor(n/2);
2292 for (uintL bit=cln::integer_length(m); bit>0; --bit) {
2293 // Since a squaring is cheaper than a multiplication, better use
2294 // three squarings instead of one multiplication and two squarings.
2295 cln::cl_I u2 = cln::square(u);
2296 cln::cl_I v2 = cln::square(v);
2297 if (cln::logbitp(bit-1, m)) {
2298 v = cln::square(u + v) - u2;
2299 u = u2 + v2;
2300 } else {
2301 u = v2 - cln::square(v - u);
2302 v = u2 + v2;
2303 }
2304 }
2305 if (n.is_even())
2306 // Here we don't use the squaring formula because one multiplication
2307 // is cheaper than two squarings.
2308 return numeric(u * ((v << 1) - u));
2309 else
2310 return numeric(cln::square(u) + cln::square(v));
2311}
2312
2313
2315const numeric abs(const numeric& x)
2316{
2317 return numeric(cln::abs(x.to_cl_N()));
2318}
2319
2320
2328const numeric mod(const numeric &a, const numeric &b)
2329{
2330 if (a.is_integer() && b.is_integer())
2331 return numeric(cln::mod(cln::the<cln::cl_I>(a.to_cl_N()),
2332 cln::the<cln::cl_I>(b.to_cl_N())));
2333 else
2334 return *_num0_p;
2335}
2336
2337
2341const numeric smod(const numeric &a_, const numeric &b_)
2342{
2343 if (a_.is_integer() && b_.is_integer()) {
2344 const cln::cl_I a = cln::the<cln::cl_I>(a_.to_cl_N());
2345 const cln::cl_I b = cln::the<cln::cl_I>(b_.to_cl_N());
2346 const cln::cl_I b2 = b >> 1;
2347 const cln::cl_I m = cln::mod(a, b);
2348 const cln::cl_I m_b = m - b;
2349 const cln::cl_I ret = m > b2 ? m_b : m;
2350 return numeric(ret);
2351 } else
2352 return *_num0_p;
2353}
2354
2355
2363const numeric irem(const numeric &a, const numeric &b)
2364{
2365 if (b.is_zero())
2366 throw std::overflow_error("numeric::irem(): division by zero");
2367 if (a.is_integer() && b.is_integer())
2368 return numeric(cln::rem(cln::the<cln::cl_I>(a.to_cl_N()),
2369 cln::the<cln::cl_I>(b.to_cl_N())));
2370 else
2371 return *_num0_p;
2372}
2373
2374
2383const numeric irem(const numeric &a, const numeric &b, numeric &q)
2384{
2385 if (b.is_zero())
2386 throw std::overflow_error("numeric::irem(): division by zero");
2387 if (a.is_integer() && b.is_integer()) {
2388 const cln::cl_I_div_t rem_quo = cln::truncate2(cln::the<cln::cl_I>(a.to_cl_N()),
2389 cln::the<cln::cl_I>(b.to_cl_N()));
2390 q = numeric(rem_quo.quotient);
2391 return numeric(rem_quo.remainder);
2392 } else {
2393 q = *_num0_p;
2394 return *_num0_p;
2395 }
2396}
2397
2398
2404const numeric iquo(const numeric &a, const numeric &b)
2405{
2406 if (b.is_zero())
2407 throw std::overflow_error("numeric::iquo(): division by zero");
2408 if (a.is_integer() && b.is_integer())
2409 return numeric(cln::truncate1(cln::the<cln::cl_I>(a.to_cl_N()),
2410 cln::the<cln::cl_I>(b.to_cl_N())));
2411 else
2412 return *_num0_p;
2413}
2414
2415
2423const numeric iquo(const numeric &a, const numeric &b, numeric &r)
2424{
2425 if (b.is_zero())
2426 throw std::overflow_error("numeric::iquo(): division by zero");
2427 if (a.is_integer() && b.is_integer()) {
2428 const cln::cl_I_div_t rem_quo = cln::truncate2(cln::the<cln::cl_I>(a.to_cl_N()),
2429 cln::the<cln::cl_I>(b.to_cl_N()));
2430 r = numeric(rem_quo.remainder);
2431 return numeric(rem_quo.quotient);
2432 } else {
2433 r = *_num0_p;
2434 return *_num0_p;
2435 }
2436}
2437
2438
2443const numeric gcd(const numeric &a, const numeric &b)
2444{
2445 if (a.is_integer() && b.is_integer())
2446 return numeric(cln::gcd(cln::the<cln::cl_I>(a.to_cl_N()),
2447 cln::the<cln::cl_I>(b.to_cl_N())));
2448 else
2449 return *_num1_p;
2450}
2451
2452
2457const numeric lcm(const numeric &a, const numeric &b)
2458{
2459 if (a.is_integer() && b.is_integer())
2460 return numeric(cln::lcm(cln::the<cln::cl_I>(a.to_cl_N()),
2461 cln::the<cln::cl_I>(b.to_cl_N())));
2462 else
2463 return a.mul(b);
2464}
2465
2466
2475const numeric sqrt(const numeric &x)
2476{
2477 return numeric(cln::sqrt(x.to_cl_N()));
2478}
2479
2480
2482const numeric isqrt(const numeric &x)
2483{
2484 if (x.is_integer()) {
2485 cln::cl_I root;
2486 cln::isqrt(cln::the<cln::cl_I>(x.to_cl_N()), &root);
2487 return numeric(root);
2488 } else
2489 return *_num0_p;
2490}
2491
2492
2495{
2496 return numeric(cln::pi(cln::default_float_format));
2497}
2498
2499
2502{
2503 return numeric(cln::eulerconst(cln::default_float_format));
2504}
2505
2506
2509{
2510 return numeric(cln::catalanconst(cln::default_float_format));
2511}
2512
2513
2516 : digits(17)
2517{
2518 // It initializes to 17 digits, because in CLN float_format(17) turns out
2519 // to be 61 (<64) while float_format(18)=65. The reason is we want to
2520 // have a cl_LF instead of cl_SF, cl_FF or cl_DF.
2521 if (too_late)
2522 throw(std::runtime_error("I told you not to do instantiate me!"));
2523 too_late = true;
2524 cln::default_float_format = cln::float_format(17);
2525
2526 // add callbacks for built-in functions
2527 // like ... add_callback(Li_lookuptable);
2528}
2529
2530
2533{
2534 long digitsdiff = prec - digits;
2535 digits = prec;
2536 cln::default_float_format = cln::float_format(prec);
2537
2538 // call registered callbacks
2539 for (auto it : callbacklist) {
2540 (it)(digitsdiff);
2541 }
2542
2543 return *this;
2544}
2545
2546
2548_numeric_digits::operator long()
2549{
2550 // BTW, this is approx. unsigned(cln::default_float_format*0.301)-1
2551 return (long)digits;
2552}
2553
2554
2556void _numeric_digits::print(std::ostream &os) const
2557{
2558 os << digits;
2559}
2560
2561
2564{
2565 callbacklist.push_back(callback);
2566}
2567
2568
2569std::ostream& operator<<(std::ostream &os, const _numeric_digits &e)
2570{
2571 e.print(os);
2572 return os;
2573}
2574
2576// static member variables
2578
2579// private
2580
2581bool _numeric_digits::too_late = false;
2582
2583
2587
2588} // 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:2532
_numeric_digits()
_numeric_digits default ctor, checking for singleton invariance.
Definition: numeric.cpp:2515
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:2556
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:2563
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:2443
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:2166
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:2328
const numeric abs(const numeric &x)
Absolute value.
Definition: numeric.cpp:2315
ex EulerEvalf()
Floating point evaluation of Euler's constant gamma.
Definition: numeric.cpp:2501
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:2494
const numeric fibonacci(const numeric &n)
Fibonacci number.
Definition: numeric.cpp:2258
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:1432
function psi(const T1 &p1)
Definition: inifcns.h:165
const numeric sqrt(const numeric &x)
Numeric square root.
Definition: numeric.cpp:2475
const numeric irem(const numeric &a, const numeric &b)
Numeric integer remainder.
Definition: numeric.cpp:2363
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:2143
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:2341
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:1774
const numeric iquo(const numeric &a, const numeric &b)
Numeric integer quotient.
Definition: numeric.cpp:2404
const numeric isqrt(const numeric &x)
Integer numeric square root.
Definition: numeric.cpp:2482
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:2586
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:2508
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:2457
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.