GiNaC  1.8.3
numeric.cpp
Go to the documentation of this file.
1 
9 /*
10  * GiNaC Copyright (C) 1999-2022 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 
58 namespace 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 // 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 
106 numeric::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 
134 numeric::numeric(unsigned long i)
135 {
136  value = cln::cl_I(i);
138 }
139 
140 numeric::numeric(long long i)
141 {
142  value = cln::cl_I(i);
144 }
145 
146 numeric::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 
176 numeric::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 
255 numeric::numeric(const cln::cl_N &z)
256 {
257  value = z;
259 }
260 
261 
263 // archiving
265 
269 static 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 
281 static 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 
289 void numeric::read_archive(const archive_node &n, lst &sym_lst)
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 
338 static 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 
397 static 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 
426 static 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 
440 static 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 
470 template<typename T1, typename T2>
471 static inline bool coerce(T1& dst, const T2& arg);
472 
478 template<>
479 inline 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 
492 template<>
493 inline 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 
507 static 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 
542 void 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 {
570  print_real_number(c, i);
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 {
587  print_real_number(c, i);
588  c.s << mul_sym << imag_sym;
589  }
590  } else {
591  if (i == 1) {
592  c.s << "+" << imag_sym;
593  } else {
594  c.s << "+";
595  print_real_number(c, i);
596  c.s << mul_sym << imag_sym;
597  }
598  }
599  if (precedence() <= level)
600  c.s << par_close;
601  }
602  }
603 }
604 
605 void numeric::do_print(const print_context & c, unsigned level) const
606 {
607  print_numeric(c, "(", ")", "I", "*", level);
608 }
609 
610 void numeric::do_print_latex(const print_latex & c, unsigned level) const
611 {
612  print_numeric(c, "{(", ")}", "i", " ", level);
613 }
614 
615 void 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 
651 void 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 
669 void 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 
677 void 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 
684 bool numeric::info(unsigned inf) const
685 {
686  switch (inf) {
687  case info_flags::numeric:
691  return true;
692  case info_flags::real:
693  return is_real();
696  return is_rational();
699  return is_crational();
700  case info_flags::integer:
702  return is_integer();
705  return is_cinteger();
707  return is_positive();
709  return is_negative();
711  return is_zero() || is_positive();
712  case info_flags::posint:
713  return is_pos_integer();
714  case info_flags::negint:
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();
722  case info_flags::prime:
723  return is_prime();
724  }
725  return false;
726 }
727 
728 bool numeric::is_polynomial(const ex & var) const
729 {
730  return true;
731 }
732 
733 int numeric::degree(const ex & s) const
734 {
735  return 0;
736 }
737 
738 int numeric::ldegree(const ex & s) const
739 {
740  return 0;
741 }
742 
743 ex numeric::coeff(const ex & s, int n) const
744 {
745  return n==0 ? *this : _ex0;
746 }
747 
754 bool 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 
820 int 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 
829 bool 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 
838 unsigned 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 
864 const numeric numeric::add(const numeric &other) const
865 {
866  return numeric(value + other.value);
867 }
868 
869 
872 const numeric numeric::sub(const numeric &other) const
873 {
874  return numeric(value - other.value);
875 }
876 
877 
880 const numeric numeric::mul(const numeric &other) const
881 {
882  return numeric(value * other.value);
883 }
884 
885 
890 const 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 
900 const 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 
925 const 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 
942 const 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 
957 const 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 
976 const 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 
993 const 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 
1022 const numeric &numeric::operator=(unsigned int i)
1023 {
1024  return operator=(numeric(i));
1025 }
1026 
1027 
1029 {
1030  return operator=(numeric(i));
1031 }
1032 
1033 
1034 const numeric &numeric::operator=(unsigned long i)
1035 {
1036  return operator=(numeric(i));
1037 }
1038 
1039 
1040 const numeric &numeric::operator=(double d)
1041 {
1042  return operator=(numeric(d));
1043 }
1044 
1045 
1046 const 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 
1078 int 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 
1104 int 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 
1122 bool numeric::is_equal(const numeric &other) const
1123 {
1124  return cln::equal(value, other.value);
1125 }
1126 
1127 
1129 bool numeric::is_zero() const
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 
1175 bool numeric::is_even() const
1176 {
1177  return (cln::instanceof(value, cln::cl_I_ring) && cln::evenp(cln::the<cln::cl_I>(value)));
1178 }
1179 
1180 
1182 bool numeric::is_odd() const
1183 {
1184  return (cln::instanceof(value, cln::cl_I_ring) && cln::oddp(cln::the<cln::cl_I>(value)));
1185 }
1186 
1187 
1191 bool numeric::is_prime() const
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 
1208 bool numeric::is_real() const
1209 {
1210  return cln::instanceof(value, cln::cl_R_ring);
1211 }
1212 
1213 
1214 bool numeric::operator==(const numeric &other) const
1215 {
1216  return cln::equal(value, other.value);
1217 }
1218 
1219 
1220 bool 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 
1259 bool 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 
1270 bool 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 
1281 bool 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 
1292 bool 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 
1303 int numeric::to_int() const
1304 {
1305  GINAC_ASSERT(this->is_integer());
1306  return cln::cl_I_to_int(cln::the<cln::cl_I>(value));
1307 }
1308 
1309 
1313 long numeric::to_long() const
1314 {
1315  GINAC_ASSERT(this->is_integer());
1316  return cln::cl_I_to_long(cln::the<cln::cl_I>(value));
1317 }
1318 
1319 
1322 double numeric::to_double() const
1323 {
1324  GINAC_ASSERT(this->is_real());
1325  return cln::double_approx(cln::realpart(value));
1326 }
1327 
1328 
1332 cln::cl_N numeric::to_cl_N() const
1333 {
1334  return value;
1335 }
1336 
1337 
1339 const numeric numeric::real() const
1340 {
1341  return numeric(cln::realpart(value));
1342 }
1343 
1344 
1346 const numeric numeric::imag() const
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 
1433 const numeric I = numeric(cln::complex(cln::cl_I(0),cln::cl_I(1)));
1434 
1435 
1439 const numeric exp(const numeric &x)
1440 {
1441  return numeric(cln::exp(x.to_cl_N()));
1442 }
1443 
1444 
1450 const 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 
1461 const numeric sin(const numeric &x)
1462 {
1463  return numeric(cln::sin(x.to_cl_N()));
1464 }
1465 
1466 
1470 const numeric cos(const numeric &x)
1471 {
1472  return numeric(cln::cos(x.to_cl_N()));
1473 }
1474 
1475 
1479 const numeric tan(const numeric &x)
1480 {
1481  return numeric(cln::tan(x.to_cl_N()));
1482 }
1483 
1484 
1488 const numeric asin(const numeric &x)
1489 {
1490  return numeric(cln::asin(x.to_cl_N()));
1491 }
1492 
1493 
1497 const numeric acos(const numeric &x)
1498 {
1499  return numeric(cln::acos(x.to_cl_N()));
1500 }
1501 
1502 
1508 const 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 
1525 const 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 
1554 const numeric sinh(const numeric &x)
1555 {
1556  return numeric(cln::sinh(x.to_cl_N()));
1557 }
1558 
1559 
1563 const numeric cosh(const numeric &x)
1564 {
1565  return numeric(cln::cosh(x.to_cl_N()));
1566 }
1567 
1568 
1572 const numeric tanh(const numeric &x)
1573 {
1574  return numeric(cln::tanh(x.to_cl_N()));
1575 }
1576 
1577 
1581 const numeric asinh(const numeric &x)
1582 {
1583  return numeric(cln::asinh(x.to_cl_N()));
1584 }
1585 
1586 
1590 const numeric acosh(const numeric &x)
1591 {
1592  return numeric(cln::acosh(x.to_cl_N()));
1593 }
1594 
1595 
1599 const 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 
1633 static 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 
1652 static 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 
1679 const 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 
1705 const 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 
1717 const 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:
1735  lanczos_coeffs();
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 
1749 std::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 
1771 cln::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 
2024 static 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 
2039 const 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 
2060 const numeric lgamma(const numeric &x)
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 
2067 const 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 
2086 const numeric tgamma(const numeric &x)
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 
2095 const numeric psi(const numeric &x)
2096 {
2097  throw dunno();
2098 }
2099 
2100 
2103 const 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 
2143 const 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 
2166 const 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 
2315 const numeric abs(const numeric& x)
2316 {
2317  return numeric(cln::abs(x.to_cl_N()));
2318 }
2319 
2320 
2328 const 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 
2341 const 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 
2363 const 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 
2383 const 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 
2404 const 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 
2423 const 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 
2443 const 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 
2457 const 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 
2475 const numeric sqrt(const numeric &x)
2476 {
2477  return numeric(cln::sqrt(x.to_cl_N()));
2478 }
2479 
2480 
2482 const 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 
2556 void _numeric_digits::print(std::ostream &os) const
2557 {
2558  os << digits;
2559 }
2560 
2561 
2564 {
2565  callbacklist.push_back(callback);
2566 }
2567 
2568 
2569 std::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 
2581 bool _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
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
const basic & setflag(unsigned f) const
Set some status_flags.
Definition: basic.h:288
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:231
unsigned options
Definition: factor.cpp:2480
vector< int > k
Definition: factor.cpp:1466
ex x
Definition: factor.cpp:1641
size_t n
Definition: factor.cpp:1463
size_t c
Definition: factor.cpp:770
size_t r
Definition: factor.cpp:770
mvec m
Definition: factor.cpp:771
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
ex denom(const ex &thisex)
Definition: ex.h:748
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:718
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
std::ostream & operator<<(std::ostream &os, const archive_node &n)
Write archive_node to binary data stream.
Definition: archive.cpp:200
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
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
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:745
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
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.