GiNaC  1.8.0
normal.cpp
Go to the documentation of this file.
1 
8 /*
9  * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
10  *
11  * This program is free software; you can redistribute it and/or modify
12  * it under the terms of the GNU General Public License as published by
13  * the Free Software Foundation; either version 2 of the License, or
14  * (at your option) any later version.
15  *
16  * This program is distributed in the hope that it will be useful,
17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19  * GNU General Public License for more details.
20  *
21  * You should have received a copy of the GNU General Public License
22  * along with this program; if not, write to the Free Software
23  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
24  */
25 
26 #include "normal.h"
27 #include "basic.h"
28 #include "ex.h"
29 #include "add.h"
30 #include "constant.h"
31 #include "expairseq.h"
32 #include "fail.h"
33 #include "inifcns.h"
34 #include "lst.h"
35 #include "mul.h"
36 #include "numeric.h"
37 #include "power.h"
38 #include "relational.h"
39 #include "operators.h"
40 #include "matrix.h"
41 #include "pseries.h"
42 #include "symbol.h"
43 #include "utils.h"
44 #include "polynomial/chinrem_gcd.h"
45 
46 #include <algorithm>
47 #include <map>
48 
49 namespace GiNaC {
50 
51 // If comparing expressions (ex::compare()) is fast, you can set this to 1.
52 // Some routines like quo(), rem() and gcd() will then return a quick answer
53 // when they are called with two identical arguments.
54 #define FAST_COMPARE 1
55 
56 // Set this if you want divide_in_z() to use remembering
57 #define USE_REMEMBER 0
58 
59 // Set this if you want divide_in_z() to use trial division followed by
60 // polynomial interpolation (always slower except for completely dense
61 // polynomials)
62 #define USE_TRIAL_DIVISION 0
63 
64 // Set this to enable some statistical output for the GCD routines
65 #define STATISTICS 0
66 
67 
68 #if STATISTICS
69 // Statistics variables
70 static int gcd_called = 0;
71 static int sr_gcd_called = 0;
72 static int heur_gcd_called = 0;
73 static int heur_gcd_failed = 0;
74 
75 // Print statistics at end of program
76 static struct _stat_print {
77  _stat_print() {}
78  ~_stat_print() {
79  std::cout << "gcd() called " << gcd_called << " times\n";
80  std::cout << "sr_gcd() called " << sr_gcd_called << " times\n";
81  std::cout << "heur_gcd() called " << heur_gcd_called << " times\n";
82  std::cout << "heur_gcd() failed " << heur_gcd_failed << " times\n";
83  }
84 } stat_print;
85 #endif
86 
87 
95 static bool get_first_symbol(const ex &e, ex &x)
96 {
97  if (is_a<symbol>(e)) {
98  x = e;
99  return true;
100  } else if (is_exactly_a<add>(e) || is_exactly_a<mul>(e)) {
101  for (size_t i=0; i<e.nops(); i++)
102  if (get_first_symbol(e.op(i), x))
103  return true;
104  } else if (is_exactly_a<power>(e)) {
105  if (get_first_symbol(e.op(0), x))
106  return true;
107  }
108  return false;
109 }
110 
111 
112 /*
113  * Statistical information about symbols in polynomials
114  */
115 
122 struct sym_desc {
124  sym_desc(const ex& s)
125  : sym(s), deg_a(0), deg_b(0), ldeg_a(0), ldeg_b(0), max_deg(0), max_lcnops(0)
126  { }
127 
130 
132  int deg_a;
133 
135  int deg_b;
136 
138  int ldeg_a;
139 
141  int ldeg_b;
142 
144  int max_deg;
145 
147  size_t max_lcnops;
148 
150  bool operator<(const sym_desc &x) const
151  {
152  if (max_deg == x.max_deg)
153  return max_lcnops < x.max_lcnops;
154  else
155  return max_deg < x.max_deg;
156  }
157 };
158 
159 // Vector of sym_desc structures
160 typedef std::vector<sym_desc> sym_desc_vec;
161 
162 // Add symbol the sym_desc_vec (used internally by get_symbol_stats())
163 static void add_symbol(const ex &s, sym_desc_vec &v)
164 {
165  for (auto & it : v)
166  if (it.sym.is_equal(s)) // If it's already in there, don't add it a second time
167  return;
168 
169  v.push_back(sym_desc(s));
170 }
171 
172 // Collect all symbols of an expression (used internally by get_symbol_stats())
173 static void collect_symbols(const ex &e, sym_desc_vec &v)
174 {
175  if (is_a<symbol>(e)) {
176  add_symbol(e, v);
177  } else if (is_exactly_a<add>(e) || is_exactly_a<mul>(e)) {
178  for (size_t i=0; i<e.nops(); i++)
179  collect_symbols(e.op(i), v);
180  } else if (is_exactly_a<power>(e)) {
181  collect_symbols(e.op(0), v);
182  }
183 }
184 
197 static void get_symbol_stats(const ex &a, const ex &b, sym_desc_vec &v)
198 {
199  collect_symbols(a, v);
200  collect_symbols(b, v);
201  for (auto & it : v) {
202  int deg_a = a.degree(it.sym);
203  int deg_b = b.degree(it.sym);
204  it.deg_a = deg_a;
205  it.deg_b = deg_b;
206  it.max_deg = std::max(deg_a, deg_b);
207  it.max_lcnops = std::max(a.lcoeff(it.sym).nops(), b.lcoeff(it.sym).nops());
208  it.ldeg_a = a.ldegree(it.sym);
209  it.ldeg_b = b.ldegree(it.sym);
210  }
211  std::sort(v.begin(), v.end());
212 
213 #if 0
214  std::clog << "Symbols:\n";
215  auto it = v.begin(), itend = v.end();
216  while (it != itend) {
217  std::clog << " " << it->sym << ": deg_a=" << it->deg_a << ", deg_b=" << it->deg_b << ", ldeg_a=" << it->ldeg_a << ", ldeg_b=" << it->ldeg_b << ", max_deg=" << it->max_deg << ", max_lcnops=" << it->max_lcnops << std::endl;
218  std::clog << " lcoeff_a=" << a.lcoeff(it->sym) << ", lcoeff_b=" << b.lcoeff(it->sym) << std::endl;
219  ++it;
220  }
221 #endif
222 }
223 
224 
225 /*
226  * Computation of LCM of denominators of coefficients of a polynomial
227  */
228 
229 // Compute LCM of denominators of coefficients by going through the
230 // expression recursively (used internally by lcm_of_coefficients_denominators())
231 static numeric lcmcoeff(const ex &e, const numeric &l)
232 {
233  if (e.info(info_flags::rational))
234  return lcm(ex_to<numeric>(e).denom(), l);
235  else if (is_exactly_a<add>(e)) {
236  numeric c = *_num1_p;
237  for (size_t i=0; i<e.nops(); i++)
238  c = lcmcoeff(e.op(i), c);
239  return lcm(c, l);
240  } else if (is_exactly_a<mul>(e)) {
241  numeric c = *_num1_p;
242  for (size_t i=0; i<e.nops(); i++)
243  c *= lcmcoeff(e.op(i), *_num1_p);
244  return lcm(c, l);
245  } else if (is_exactly_a<power>(e)) {
246  if (is_a<symbol>(e.op(0)))
247  return l;
248  else
249  return pow(lcmcoeff(e.op(0), l), ex_to<numeric>(e.op(1)));
250  }
251  return l;
252 }
253 
262 {
263  return lcmcoeff(e, *_num1_p);
264 }
265 
271 static ex multiply_lcm(const ex &e, const numeric &lcm)
272 {
273  if (lcm.is_equal(*_num1_p))
274  // e * 1 -> e;
275  return e;
276 
277  if (is_exactly_a<mul>(e)) {
278  // (a*b*...)*lcm -> (a*lcma)*(b*lcmb)*...*(lcm/(lcma*lcmb*...))
279  size_t num = e.nops();
280  exvector v;
281  v.reserve(num + 1);
282  numeric lcm_accum = *_num1_p;
283  for (size_t i=0; i<num; i++) {
284  numeric op_lcm = lcmcoeff(e.op(i), *_num1_p);
285  v.push_back(multiply_lcm(e.op(i), op_lcm));
286  lcm_accum *= op_lcm;
287  }
288  v.push_back(lcm / lcm_accum);
289  return dynallocate<mul>(v);
290  } else if (is_exactly_a<add>(e)) {
291  // (a+b+...)*lcm -> a*lcm+b*lcm+...
292  size_t num = e.nops();
293  exvector v;
294  v.reserve(num);
295  for (size_t i=0; i<num; i++)
296  v.push_back(multiply_lcm(e.op(i), lcm));
297  return dynallocate<add>(v);
298  } else if (is_exactly_a<power>(e)) {
299  if (!is_a<symbol>(e.op(0))) {
300  // (b^e)*lcm -> (b*lcm^(1/e))^e if lcm^(1/e) ∈ ℚ (i.e. not a float)
301  // but not for symbolic b, as evaluation would undo this again
302  numeric root_of_lcm = lcm.power(ex_to<numeric>(e.op(1)).inverse());
303  if (root_of_lcm.is_rational())
304  return pow(multiply_lcm(e.op(0), root_of_lcm), e.op(1));
305  }
306  }
307  // can't recurse down into e
308  return dynallocate<mul>(e, lcm);
309 }
310 
311 
319 {
320  return bp->integer_content();
321 }
322 
324 {
325  return *_num1_p;
326 }
327 
329 {
330  return abs(*this);
331 }
332 
334 {
335  numeric c = *_num0_p, l = *_num1_p;
336  for (auto & it : seq) {
337  GINAC_ASSERT(!is_exactly_a<numeric>(it.rest));
338  GINAC_ASSERT(is_exactly_a<numeric>(it.coeff));
339  c = gcd(ex_to<numeric>(it.coeff).numer(), c);
340  l = lcm(ex_to<numeric>(it.coeff).denom(), l);
341  }
342  GINAC_ASSERT(is_exactly_a<numeric>(overall_coeff));
343  c = gcd(ex_to<numeric>(overall_coeff).numer(), c);
344  l = lcm(ex_to<numeric>(overall_coeff).denom(), l);
345  return c/l;
346 }
347 
349 {
350 #ifdef DO_GINAC_ASSERT
351  for (auto & it : seq) {
352  GINAC_ASSERT(!is_exactly_a<numeric>(recombine_pair_to_ex(it)));
353  }
354 #endif // def DO_GINAC_ASSERT
355  GINAC_ASSERT(is_exactly_a<numeric>(overall_coeff));
356  return abs(ex_to<numeric>(overall_coeff));
357 }
358 
359 
360 /*
361  * Polynomial quotients and remainders
362  */
363 
373 ex quo(const ex &a, const ex &b, const ex &x, bool check_args)
374 {
375  if (b.is_zero())
376  throw(std::overflow_error("quo: division by zero"));
377  if (is_exactly_a<numeric>(a) && is_exactly_a<numeric>(b))
378  return a / b;
379 #if FAST_COMPARE
380  if (a.is_equal(b))
381  return _ex1;
382 #endif
384  throw(std::invalid_argument("quo: arguments must be polynomials over the rationals"));
385 
386  // Polynomial long division
387  ex r = a.expand();
388  if (r.is_zero())
389  return r;
390  int bdeg = b.degree(x);
391  int rdeg = r.degree(x);
392  ex blcoeff = b.expand().coeff(x, bdeg);
393  bool blcoeff_is_numeric = is_exactly_a<numeric>(blcoeff);
394  exvector v; v.reserve(std::max(rdeg - bdeg + 1, 0));
395  while (rdeg >= bdeg) {
396  ex term, rcoeff = r.coeff(x, rdeg);
397  if (blcoeff_is_numeric)
398  term = rcoeff / blcoeff;
399  else {
400  if (!divide(rcoeff, blcoeff, term, false))
401  return dynallocate<fail>();
402  }
403  term *= pow(x, rdeg - bdeg);
404  v.push_back(term);
405  r -= (term * b).expand();
406  if (r.is_zero())
407  break;
408  rdeg = r.degree(x);
409  }
410  return dynallocate<add>(v);
411 }
412 
413 
423 ex rem(const ex &a, const ex &b, const ex &x, bool check_args)
424 {
425  if (b.is_zero())
426  throw(std::overflow_error("rem: division by zero"));
427  if (is_exactly_a<numeric>(a)) {
428  if (is_exactly_a<numeric>(b))
429  return _ex0;
430  else
431  return a;
432  }
433 #if FAST_COMPARE
434  if (a.is_equal(b))
435  return _ex0;
436 #endif
438  throw(std::invalid_argument("rem: arguments must be polynomials over the rationals"));
439 
440  // Polynomial long division
441  ex r = a.expand();
442  if (r.is_zero())
443  return r;
444  int bdeg = b.degree(x);
445  int rdeg = r.degree(x);
446  ex blcoeff = b.expand().coeff(x, bdeg);
447  bool blcoeff_is_numeric = is_exactly_a<numeric>(blcoeff);
448  while (rdeg >= bdeg) {
449  ex term, rcoeff = r.coeff(x, rdeg);
450  if (blcoeff_is_numeric)
451  term = rcoeff / blcoeff;
452  else {
453  if (!divide(rcoeff, blcoeff, term, false))
454  return dynallocate<fail>();
455  }
456  term *= pow(x, rdeg - bdeg);
457  r -= (term * b).expand();
458  if (r.is_zero())
459  break;
460  rdeg = r.degree(x);
461  }
462  return r;
463 }
464 
465 
472 ex decomp_rational(const ex &a, const ex &x)
473 {
474  ex nd = numer_denom(a);
475  ex numer = nd.op(0), denom = nd.op(1);
476  ex q = quo(numer, denom, x);
477  if (is_exactly_a<fail>(q))
478  return a;
479  else
480  return q + rem(numer, denom, x) / denom;
481 }
482 
483 
492 ex prem(const ex &a, const ex &b, const ex &x, bool check_args)
493 {
494  if (b.is_zero())
495  throw(std::overflow_error("prem: division by zero"));
496  if (is_exactly_a<numeric>(a)) {
497  if (is_exactly_a<numeric>(b))
498  return _ex0;
499  else
500  return b;
501  }
503  throw(std::invalid_argument("prem: arguments must be polynomials over the rationals"));
504 
505  // Polynomial long division
506  ex r = a.expand();
507  ex eb = b.expand();
508  int rdeg = r.degree(x);
509  int bdeg = eb.degree(x);
510  ex blcoeff;
511  if (bdeg <= rdeg) {
512  blcoeff = eb.coeff(x, bdeg);
513  if (bdeg == 0)
514  eb = _ex0;
515  else
516  eb -= blcoeff * pow(x, bdeg);
517  } else
518  blcoeff = _ex1;
519 
520  int delta = rdeg - bdeg + 1, i = 0;
521  while (rdeg >= bdeg && !r.is_zero()) {
522  ex rlcoeff = r.coeff(x, rdeg);
523  ex term = (pow(x, rdeg - bdeg) * eb * rlcoeff).expand();
524  if (rdeg == 0)
525  r = _ex0;
526  else
527  r -= rlcoeff * pow(x, rdeg);
528  r = (blcoeff * r).expand() - term;
529  rdeg = r.degree(x);
530  i++;
531  }
532  return pow(blcoeff, delta - i) * r;
533 }
534 
535 
544 ex sprem(const ex &a, const ex &b, const ex &x, bool check_args)
545 {
546  if (b.is_zero())
547  throw(std::overflow_error("prem: division by zero"));
548  if (is_exactly_a<numeric>(a)) {
549  if (is_exactly_a<numeric>(b))
550  return _ex0;
551  else
552  return b;
553  }
555  throw(std::invalid_argument("prem: arguments must be polynomials over the rationals"));
556 
557  // Polynomial long division
558  ex r = a.expand();
559  ex eb = b.expand();
560  int rdeg = r.degree(x);
561  int bdeg = eb.degree(x);
562  ex blcoeff;
563  if (bdeg <= rdeg) {
564  blcoeff = eb.coeff(x, bdeg);
565  if (bdeg == 0)
566  eb = _ex0;
567  else
568  eb -= blcoeff * pow(x, bdeg);
569  } else
570  blcoeff = _ex1;
571 
572  while (rdeg >= bdeg && !r.is_zero()) {
573  ex rlcoeff = r.coeff(x, rdeg);
574  ex term = (pow(x, rdeg - bdeg) * eb * rlcoeff).expand();
575  if (rdeg == 0)
576  r = _ex0;
577  else
578  r -= rlcoeff * pow(x, rdeg);
579  r = (blcoeff * r).expand() - term;
580  rdeg = r.degree(x);
581  }
582  return r;
583 }
584 
585 
595 bool divide(const ex &a, const ex &b, ex &q, bool check_args)
596 {
597  if (b.is_zero())
598  throw(std::overflow_error("divide: division by zero"));
599  if (a.is_zero()) {
600  q = _ex0;
601  return true;
602  }
603  if (is_exactly_a<numeric>(b)) {
604  q = a / b;
605  return true;
606  } else if (is_exactly_a<numeric>(a))
607  return false;
608 #if FAST_COMPARE
609  if (a.is_equal(b)) {
610  q = _ex1;
611  return true;
612  }
613 #endif
614  if (check_args && (!a.info(info_flags::rational_polynomial) ||
616  throw(std::invalid_argument("divide: arguments must be polynomials over the rationals"));
617 
618  // Find first symbol
619  ex x;
620  if (!get_first_symbol(a, x) && !get_first_symbol(b, x))
621  throw(std::invalid_argument("invalid expression in divide()"));
622 
623  // Try to avoid expanding partially factored expressions.
624  if (is_exactly_a<mul>(b)) {
625  // Divide sequentially by each term
626  ex rem_new, rem_old = a;
627  for (size_t i=0; i < b.nops(); i++) {
628  if (! divide(rem_old, b.op(i), rem_new, false))
629  return false;
630  rem_old = rem_new;
631  }
632  q = rem_new;
633  return true;
634  } else if (is_exactly_a<power>(b)) {
635  const ex& bb(b.op(0));
636  int exp_b = ex_to<numeric>(b.op(1)).to_int();
637  ex rem_new, rem_old = a;
638  for (int i=exp_b; i>0; i--) {
639  if (! divide(rem_old, bb, rem_new, false))
640  return false;
641  rem_old = rem_new;
642  }
643  q = rem_new;
644  return true;
645  }
646 
647  if (is_exactly_a<mul>(a)) {
648  // Divide sequentially each term. If some term in a is divisible
649  // by b we are done... and if not, we can't really say anything.
650  size_t i;
651  ex rem_i;
652  bool divisible_p = false;
653  for (i=0; i < a.nops(); ++i) {
654  if (divide(a.op(i), b, rem_i, false)) {
655  divisible_p = true;
656  break;
657  }
658  }
659  if (divisible_p) {
660  exvector resv;
661  resv.reserve(a.nops());
662  for (size_t j=0; j < a.nops(); j++) {
663  if (j==i)
664  resv.push_back(rem_i);
665  else
666  resv.push_back(a.op(j));
667  }
668  q = dynallocate<mul>(resv);
669  return true;
670  }
671  } else if (is_exactly_a<power>(a)) {
672  // The base itself might be divisible by b, in that case we don't
673  // need to expand a
674  const ex& ab(a.op(0));
675  int a_exp = ex_to<numeric>(a.op(1)).to_int();
676  ex rem_i;
677  if (divide(ab, b, rem_i, false)) {
678  q = rem_i * pow(ab, a_exp - 1);
679  return true;
680  }
681 // code below is commented-out because it leads to a significant slowdown
682 // for (int i=2; i < a_exp; i++) {
683 // if (divide(power(ab, i), b, rem_i, false)) {
684 // q = rem_i*power(ab, a_exp - i);
685 // return true;
686 // }
687 // } // ... so we *really* need to expand expression.
688  }
689 
690  // Polynomial long division (recursive)
691  ex r = a.expand();
692  if (r.is_zero()) {
693  q = _ex0;
694  return true;
695  }
696  int bdeg = b.degree(x);
697  int rdeg = r.degree(x);
698  ex blcoeff = b.expand().coeff(x, bdeg);
699  bool blcoeff_is_numeric = is_exactly_a<numeric>(blcoeff);
700  exvector v; v.reserve(std::max(rdeg - bdeg + 1, 0));
701  while (rdeg >= bdeg) {
702  ex term, rcoeff = r.coeff(x, rdeg);
703  if (blcoeff_is_numeric)
704  term = rcoeff / blcoeff;
705  else
706  if (!divide(rcoeff, blcoeff, term, false))
707  return false;
708  term *= pow(x, rdeg - bdeg);
709  v.push_back(term);
710  r -= (term * b).expand();
711  if (r.is_zero()) {
712  q = dynallocate<add>(v);
713  return true;
714  }
715  rdeg = r.degree(x);
716  }
717  return false;
718 }
719 
720 
721 #if USE_REMEMBER
722 /*
723  * Remembering
724  */
725 
726 typedef std::pair<ex, ex> ex2;
727 typedef std::pair<ex, bool> exbool;
728 
729 struct ex2_less {
730  bool operator() (const ex2 &p, const ex2 &q) const
731  {
732  int cmp = p.first.compare(q.first);
733  return ((cmp<0) || (!(cmp>0) && p.second.compare(q.second)<0));
734  }
735 };
736 
737 typedef std::map<ex2, exbool, ex2_less> ex2_exbool_remember;
738 #endif
739 
740 
757 static bool divide_in_z(const ex &a, const ex &b, ex &q, sym_desc_vec::const_iterator var)
758 {
759  q = _ex0;
760  if (b.is_zero())
761  throw(std::overflow_error("divide_in_z: division by zero"));
762  if (b.is_equal(_ex1)) {
763  q = a;
764  return true;
765  }
766  if (is_exactly_a<numeric>(a)) {
767  if (is_exactly_a<numeric>(b)) {
768  q = a / b;
769  return q.info(info_flags::integer);
770  } else
771  return false;
772  }
773 #if FAST_COMPARE
774  if (a.is_equal(b)) {
775  q = _ex1;
776  return true;
777  }
778 #endif
779 
780 #if USE_REMEMBER
781  // Remembering
782  static ex2_exbool_remember dr_remember;
783  ex2_exbool_remember::const_iterator remembered = dr_remember.find(ex2(a, b));
784  if (remembered != dr_remember.end()) {
785  q = remembered->second.first;
786  return remembered->second.second;
787  }
788 #endif
789 
790  if (is_exactly_a<power>(b)) {
791  const ex& bb(b.op(0));
792  ex qbar = a;
793  int exp_b = ex_to<numeric>(b.op(1)).to_int();
794  for (int i=exp_b; i>0; i--) {
795  if (!divide_in_z(qbar, bb, q, var))
796  return false;
797  qbar = q;
798  }
799  return true;
800  }
801 
802  if (is_exactly_a<mul>(b)) {
803  ex qbar = a;
804  for (const auto & it : b) {
805  sym_desc_vec sym_stats;
806  get_symbol_stats(a, it, sym_stats);
807  if (!divide_in_z(qbar, it, q, sym_stats.begin()))
808  return false;
809 
810  qbar = q;
811  }
812  return true;
813  }
814 
815  // Main symbol
816  const ex &x = var->sym;
817 
818  // Compare degrees
819  int adeg = a.degree(x), bdeg = b.degree(x);
820  if (bdeg > adeg)
821  return false;
822 
823 #if USE_TRIAL_DIVISION
824 
825  // Trial division with polynomial interpolation
826  int i, k;
827 
828  // Compute values at evaluation points 0..adeg
829  vector<numeric> alpha; alpha.reserve(adeg + 1);
830  exvector u; u.reserve(adeg + 1);
831  numeric point = *_num0_p;
832  ex c;
833  for (i=0; i<=adeg; i++) {
834  ex bs = b.subs(x == point, subs_options::no_pattern);
835  while (bs.is_zero()) {
836  point += *_num1_p;
837  bs = b.subs(x == point, subs_options::no_pattern);
838  }
839  if (!divide_in_z(a.subs(x == point, subs_options::no_pattern), bs, c, var+1))
840  return false;
841  alpha.push_back(point);
842  u.push_back(c);
843  point += *_num1_p;
844  }
845 
846  // Compute inverses
847  vector<numeric> rcp; rcp.reserve(adeg + 1);
848  rcp.push_back(*_num0_p);
849  for (k=1; k<=adeg; k++) {
850  numeric product = alpha[k] - alpha[0];
851  for (i=1; i<k; i++)
852  product *= alpha[k] - alpha[i];
853  rcp.push_back(product.inverse());
854  }
855 
856  // Compute Newton coefficients
857  exvector v; v.reserve(adeg + 1);
858  v.push_back(u[0]);
859  for (k=1; k<=adeg; k++) {
860  ex temp = v[k - 1];
861  for (i=k-2; i>=0; i--)
862  temp = temp * (alpha[k] - alpha[i]) + v[i];
863  v.push_back((u[k] - temp) * rcp[k]);
864  }
865 
866  // Convert from Newton form to standard form
867  c = v[adeg];
868  for (k=adeg-1; k>=0; k--)
869  c = c * (x - alpha[k]) + v[k];
870 
871  if (c.degree(x) == (adeg - bdeg)) {
872  q = c.expand();
873  return true;
874  } else
875  return false;
876 
877 #else
878 
879  // Polynomial long division (recursive)
880  ex r = a.expand();
881  if (r.is_zero())
882  return true;
883  int rdeg = adeg;
884  ex eb = b.expand();
885  ex blcoeff = eb.coeff(x, bdeg);
886  exvector v; v.reserve(std::max(rdeg - bdeg + 1, 0));
887  while (rdeg >= bdeg) {
888  ex term, rcoeff = r.coeff(x, rdeg);
889  if (!divide_in_z(rcoeff, blcoeff, term, var+1))
890  break;
891  term = (term * pow(x, rdeg - bdeg)).expand();
892  v.push_back(term);
893  r -= (term * eb).expand();
894  if (r.is_zero()) {
895  q = dynallocate<add>(v);
896 #if USE_REMEMBER
897  dr_remember[ex2(a, b)] = exbool(q, true);
898 #endif
899  return true;
900  }
901  rdeg = r.degree(x);
902  }
903 #if USE_REMEMBER
904  dr_remember[ex2(a, b)] = exbool(q, false);
905 #endif
906  return false;
907 
908 #endif
909 }
910 
911 
912 /*
913  * Separation of unit part, content part and primitive part of polynomials
914  */
915 
923 ex ex::unit(const ex &x) const
924 {
925  ex c = expand().lcoeff(x);
926  if (is_exactly_a<numeric>(c))
927  return c.info(info_flags::negative) ?_ex_1 : _ex1;
928  else {
929  ex y;
930  if (get_first_symbol(c, y))
931  return c.unit(y);
932  else
933  throw(std::invalid_argument("invalid expression in unit()"));
934  }
935 }
936 
937 
945 ex ex::content(const ex &x) const
946 {
947  if (is_exactly_a<numeric>(*this))
948  return info(info_flags::negative) ? -*this : *this;
949 
950  ex e = expand();
951  if (e.is_zero())
952  return _ex0;
953 
954  // First, divide out the integer content (which we can calculate very efficiently).
955  // If the leading coefficient of the quotient is an integer, we are done.
956  ex c = e.integer_content();
957  ex r = e / c;
958  int deg = r.degree(x);
959  ex lcoeff = r.coeff(x, deg);
961  return c;
962 
963  // GCD of all coefficients
964  int ldeg = r.ldegree(x);
965  if (deg == ldeg)
966  return lcoeff * c / lcoeff.unit(x);
967  ex cont = _ex0;
968  for (int i=ldeg; i<=deg; i++)
969  cont = gcd(r.coeff(x, i), cont, nullptr, nullptr, false);
970  return cont * c;
971 }
972 
973 
981 ex ex::primpart(const ex &x) const
982 {
983  // We need to compute the unit and content anyway, so call unitcontprim()
984  ex u, c, p;
985  unitcontprim(x, u, c, p);
986  return p;
987 }
988 
989 
997 ex ex::primpart(const ex &x, const ex &c) const
998 {
999  if (is_zero() || c.is_zero())
1000  return _ex0;
1001  if (is_exactly_a<numeric>(*this))
1002  return _ex1;
1003 
1004  // Divide by unit and content to get primitive part
1005  ex u = unit(x);
1006  if (is_exactly_a<numeric>(c))
1007  return *this / (c * u);
1008  else
1009  return quo(*this, c * u, x, false);
1010 }
1011 
1012 
1022 void ex::unitcontprim(const ex &x, ex &u, ex &c, ex &p) const
1023 {
1024  // Quick check for zero (avoid expanding)
1025  if (is_zero()) {
1026  u = _ex1;
1027  c = p = _ex0;
1028  return;
1029  }
1030 
1031  // Special case: input is a number
1032  if (is_exactly_a<numeric>(*this)) {
1033  if (info(info_flags::negative)) {
1034  u = _ex_1;
1035  c = abs(ex_to<numeric>(*this));
1036  } else {
1037  u = _ex1;
1038  c = *this;
1039  }
1040  p = _ex1;
1041  return;
1042  }
1043 
1044  // Expand input polynomial
1045  ex e = expand();
1046  if (e.is_zero()) {
1047  u = _ex1;
1048  c = p = _ex0;
1049  return;
1050  }
1051 
1052  // Compute unit and content
1053  u = unit(x);
1054  c = content(x);
1055 
1056  // Divide by unit and content to get primitive part
1057  if (c.is_zero()) {
1058  p = _ex0;
1059  return;
1060  }
1061  if (is_exactly_a<numeric>(c))
1062  p = *this / (c * u);
1063  else
1064  p = quo(e, c * u, x, false);
1065 }
1066 
1067 
1068 /*
1069  * GCD of multivariate polynomials
1070  */
1071 
1081 static ex sr_gcd(const ex &a, const ex &b, sym_desc_vec::const_iterator var)
1082 {
1083 #if STATISTICS
1084  sr_gcd_called++;
1085 #endif
1086 
1087  // The first symbol is our main variable
1088  const ex &x = var->sym;
1089 
1090  // Sort c and d so that c has higher degree
1091  ex c, d;
1092  int adeg = a.degree(x), bdeg = b.degree(x);
1093  int cdeg, ddeg;
1094  if (adeg >= bdeg) {
1095  c = a;
1096  d = b;
1097  cdeg = adeg;
1098  ddeg = bdeg;
1099  } else {
1100  c = b;
1101  d = a;
1102  cdeg = bdeg;
1103  ddeg = adeg;
1104  }
1105 
1106  // Remove content from c and d, to be attached to GCD later
1107  ex cont_c = c.content(x);
1108  ex cont_d = d.content(x);
1109  ex gamma = gcd(cont_c, cont_d, nullptr, nullptr, false);
1110  if (ddeg == 0)
1111  return gamma;
1112  c = c.primpart(x, cont_c);
1113  d = d.primpart(x, cont_d);
1114 
1115  // First element of subresultant sequence
1116  ex r = _ex0, ri = _ex1, psi = _ex1;
1117  int delta = cdeg - ddeg;
1118 
1119  for (;;) {
1120 
1121  // Calculate polynomial pseudo-remainder
1122  r = prem(c, d, x, false);
1123  if (r.is_zero())
1124  return gamma * d.primpart(x);
1125 
1126  c = d;
1127  cdeg = ddeg;
1128  if (!divide_in_z(r, ri * pow(psi, delta), d, var))
1129  throw(std::runtime_error("invalid expression in sr_gcd(), division failed"));
1130  ddeg = d.degree(x);
1131  if (ddeg == 0) {
1132  if (is_exactly_a<numeric>(r))
1133  return gamma;
1134  else
1135  return gamma * r.primpart(x);
1136  }
1137 
1138  // Next element of subresultant sequence
1139  ri = c.expand().lcoeff(x);
1140  if (delta == 1)
1141  psi = ri;
1142  else if (delta)
1143  divide_in_z(pow(ri, delta), pow(psi, delta-1), psi, var+1);
1144  delta = cdeg - ddeg;
1145  }
1146 }
1147 
1148 
1155 {
1156  return bp->max_coefficient();
1157 }
1158 
1162 {
1163  return *_num1_p;
1164 }
1165 
1167 {
1168  return abs(*this);
1169 }
1170 
1172 {
1173  GINAC_ASSERT(is_exactly_a<numeric>(overall_coeff));
1174  numeric cur_max = abs(ex_to<numeric>(overall_coeff));
1175  for (auto & it : seq) {
1176  numeric a;
1177  GINAC_ASSERT(!is_exactly_a<numeric>(it.rest));
1178  a = abs(ex_to<numeric>(it.coeff));
1179  if (a > cur_max)
1180  cur_max = a;
1181  }
1182  return cur_max;
1183 }
1184 
1186 {
1187 #ifdef DO_GINAC_ASSERT
1188  for (auto & it : seq) {
1189  GINAC_ASSERT(!is_exactly_a<numeric>(recombine_pair_to_ex(it)));
1190  }
1191 #endif // def DO_GINAC_ASSERT
1192  GINAC_ASSERT(is_exactly_a<numeric>(overall_coeff));
1193  return abs(ex_to<numeric>(overall_coeff));
1194 }
1195 
1196 
1203 ex basic::smod(const numeric &xi) const
1204 {
1205  return *this;
1206 }
1207 
1208 ex numeric::smod(const numeric &xi) const
1209 {
1210  return GiNaC::smod(*this, xi);
1211 }
1212 
1213 ex add::smod(const numeric &xi) const
1214 {
1215  epvector newseq;
1216  newseq.reserve(seq.size()+1);
1217  for (auto & it : seq) {
1218  GINAC_ASSERT(!is_exactly_a<numeric>(it.rest));
1219  numeric coeff = GiNaC::smod(ex_to<numeric>(it.coeff), xi);
1220  if (!coeff.is_zero())
1221  newseq.push_back(expair(it.rest, coeff));
1222  }
1223  GINAC_ASSERT(is_exactly_a<numeric>(overall_coeff));
1224  numeric coeff = GiNaC::smod(ex_to<numeric>(overall_coeff), xi);
1225  return dynallocate<add>(std::move(newseq), coeff);
1226 }
1227 
1228 ex mul::smod(const numeric &xi) const
1229 {
1230 #ifdef DO_GINAC_ASSERT
1231  for (auto & it : seq) {
1232  GINAC_ASSERT(!is_exactly_a<numeric>(recombine_pair_to_ex(it)));
1233  }
1234 #endif // def DO_GINAC_ASSERT
1235  mul & mulcopy = dynallocate<mul>(*this);
1236  GINAC_ASSERT(is_exactly_a<numeric>(overall_coeff));
1237  mulcopy.overall_coeff = GiNaC::smod(ex_to<numeric>(overall_coeff),xi);
1240  return mulcopy;
1241 }
1242 
1243 
1245 static ex interpolate(const ex &gamma, const numeric &xi, const ex &x, int degree_hint = 1)
1246 {
1247  exvector g; g.reserve(degree_hint);
1248  ex e = gamma;
1249  numeric rxi = xi.inverse();
1250  for (int i=0; !e.is_zero(); i++) {
1251  ex gi = e.smod(xi);
1252  g.push_back(gi * pow(x, i));
1253  e = (e - gi) * rxi;
1254  }
1255  return dynallocate<add>(g);
1256 }
1257 
1259 class gcdheu_failed {};
1260 
1277 static bool heur_gcd_z(ex& res, const ex &a, const ex &b, ex *ca, ex *cb,
1278  sym_desc_vec::const_iterator var)
1279 {
1280 #if STATISTICS
1281  heur_gcd_called++;
1282 #endif
1283 
1284  // Algorithm only works for non-vanishing input polynomials
1285  if (a.is_zero() || b.is_zero())
1286  return false;
1287 
1288  // GCD of two numeric values -> CLN
1289  if (is_exactly_a<numeric>(a) && is_exactly_a<numeric>(b)) {
1290  numeric g = gcd(ex_to<numeric>(a), ex_to<numeric>(b));
1291  if (ca)
1292  *ca = ex_to<numeric>(a) / g;
1293  if (cb)
1294  *cb = ex_to<numeric>(b) / g;
1295  res = g;
1296  return true;
1297  }
1298 
1299  // The first symbol is our main variable
1300  const ex &x = var->sym;
1301 
1302  // Remove integer content
1303  numeric gc = gcd(a.integer_content(), b.integer_content());
1304  numeric rgc = gc.inverse();
1305  ex p = a * rgc;
1306  ex q = b * rgc;
1307  int maxdeg = std::max(p.degree(x), q.degree(x));
1308 
1309  // Find evaluation point
1310  numeric mp = p.max_coefficient();
1311  numeric mq = q.max_coefficient();
1312  numeric xi;
1313  if (mp > mq)
1314  xi = mq * (*_num2_p) + (*_num2_p);
1315  else
1316  xi = mp * (*_num2_p) + (*_num2_p);
1317 
1318  // 6 tries maximum
1319  for (int t=0; t<6; t++) {
1320  if (xi.int_length() * maxdeg > 100000) {
1321  throw gcdheu_failed();
1322  }
1323 
1324  // Apply evaluation homomorphism and calculate GCD
1325  ex cp, cq;
1326  ex gamma;
1327  bool found = heur_gcd_z(gamma,
1328  p.subs(x == xi, subs_options::no_pattern),
1329  q.subs(x == xi, subs_options::no_pattern),
1330  &cp, &cq, var+1);
1331  if (found) {
1332  gamma = gamma.expand();
1333  // Reconstruct polynomial from GCD of mapped polynomials
1334  ex g = interpolate(gamma, xi, x, maxdeg);
1335 
1336  // Remove integer content
1337  g /= g.integer_content();
1338 
1339  // If the calculated polynomial divides both p and q, this is the GCD
1340  ex dummy;
1341  if (divide_in_z(p, g, ca ? *ca : dummy, var) && divide_in_z(q, g, cb ? *cb : dummy, var)) {
1342  g *= gc;
1343  res = g;
1344  return true;
1345  }
1346  }
1347 
1348  // Next evaluation point
1349  xi = iquo(xi * isqrt(isqrt(xi)) * numeric(73794), numeric(27011));
1350  }
1351  return false;
1352 }
1353 
1371 static bool heur_gcd(ex& res, const ex& a, const ex& b, ex *ca, ex *cb,
1372  sym_desc_vec::const_iterator var)
1373 {
1376  try {
1377  return heur_gcd_z(res, a, b, ca, cb, var);
1378  } catch (gcdheu_failed) {
1379  return false;
1380  }
1381  }
1382 
1383  // convert polynomials to Z[X]
1384  const numeric a_lcm = lcm_of_coefficients_denominators(a);
1385  const numeric ab_lcm = lcmcoeff(b, a_lcm);
1386 
1387  const ex ai = a*ab_lcm;
1388  const ex bi = b*ab_lcm;
1390  throw std::logic_error("heur_gcd: not an integer polynomial [1]");
1391 
1393  throw std::logic_error("heur_gcd: not an integer polynomial [2]");
1394 
1395  bool found = false;
1396  try {
1397  found = heur_gcd_z(res, ai, bi, ca, cb, var);
1398  } catch (gcdheu_failed) {
1399  return false;
1400  }
1401 
1402  // GCD is not unique, it's defined up to a unit (i.e. invertible
1403  // element). If the coefficient ring is a field, every its element is
1404  // invertible, so one can multiply the polynomial GCD with any element
1405  // of the coefficient field. We use this ambiguity to make cofactors
1406  // integer polynomials.
1407  if (found)
1408  res /= ab_lcm;
1409  return found;
1410 }
1411 
1412 
1413 // gcd helper to handle partially factored polynomials (to avoid expanding
1414 // large expressions). At least one of the arguments should be a power.
1415 static ex gcd_pf_pow(const ex& a, const ex& b, ex* ca, ex* cb);
1416 
1417 // gcd helper to handle partially factored polynomials (to avoid expanding
1418 // large expressions). At least one of the arguments should be a product.
1419 static ex gcd_pf_mul(const ex& a, const ex& b, ex* ca, ex* cb);
1420 
1432 ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned options)
1433 {
1434 #if STATISTICS
1435  gcd_called++;
1436 #endif
1437 
1438  // GCD of numerics -> CLN
1439  if (is_exactly_a<numeric>(a) && is_exactly_a<numeric>(b)) {
1440  numeric g = gcd(ex_to<numeric>(a), ex_to<numeric>(b));
1441  if (ca || cb) {
1442  if (g.is_zero()) {
1443  if (ca)
1444  *ca = _ex0;
1445  if (cb)
1446  *cb = _ex0;
1447  } else {
1448  if (ca)
1449  *ca = ex_to<numeric>(a) / g;
1450  if (cb)
1451  *cb = ex_to<numeric>(b) / g;
1452  }
1453  }
1454  return g;
1455  }
1456 
1457  // Check arguments
1459  throw(std::invalid_argument("gcd: arguments must be polynomials over the rationals"));
1460  }
1461 
1462  // Partially factored cases (to avoid expanding large expressions)
1464  if (is_exactly_a<mul>(a) || is_exactly_a<mul>(b))
1465  return gcd_pf_mul(a, b, ca, cb);
1466 #if FAST_COMPARE
1467  if (is_exactly_a<power>(a) || is_exactly_a<power>(b))
1468  return gcd_pf_pow(a, b, ca, cb);
1469 #endif
1470  }
1471 
1472  // Some trivial cases
1473  ex aex = a.expand();
1474  if (aex.is_zero()) {
1475  if (ca)
1476  *ca = _ex0;
1477  if (cb)
1478  *cb = _ex1;
1479  return b;
1480  }
1481  ex bex = b.expand();
1482  if (bex.is_zero()) {
1483  if (ca)
1484  *ca = _ex1;
1485  if (cb)
1486  *cb = _ex0;
1487  return a;
1488  }
1489  if (aex.is_equal(_ex1) || bex.is_equal(_ex1)) {
1490  if (ca)
1491  *ca = a;
1492  if (cb)
1493  *cb = b;
1494  return _ex1;
1495  }
1496 #if FAST_COMPARE
1497  if (a.is_equal(b)) {
1498  if (ca)
1499  *ca = _ex1;
1500  if (cb)
1501  *cb = _ex1;
1502  return a;
1503  }
1504 #endif
1505 
1506  if (is_a<symbol>(aex)) {
1507  if (! bex.subs(aex==_ex0, subs_options::no_pattern).is_zero()) {
1508  if (ca)
1509  *ca = a;
1510  if (cb)
1511  *cb = b;
1512  return _ex1;
1513  }
1514  }
1515 
1516  if (is_a<symbol>(bex)) {
1517  if (! aex.subs(bex==_ex0, subs_options::no_pattern).is_zero()) {
1518  if (ca)
1519  *ca = a;
1520  if (cb)
1521  *cb = b;
1522  return _ex1;
1523  }
1524  }
1525 
1526  if (is_exactly_a<numeric>(aex)) {
1527  numeric bcont = bex.integer_content();
1528  numeric g = gcd(ex_to<numeric>(aex), bcont);
1529  if (ca)
1530  *ca = ex_to<numeric>(aex)/g;
1531  if (cb)
1532  *cb = bex/g;
1533  return g;
1534  }
1535 
1536  if (is_exactly_a<numeric>(bex)) {
1537  numeric acont = aex.integer_content();
1538  numeric g = gcd(ex_to<numeric>(bex), acont);
1539  if (ca)
1540  *ca = aex/g;
1541  if (cb)
1542  *cb = ex_to<numeric>(bex)/g;
1543  return g;
1544  }
1545 
1546  // Gather symbol statistics
1547  sym_desc_vec sym_stats;
1548  get_symbol_stats(a, b, sym_stats);
1549 
1550  // The symbol with least degree which is contained in both polynomials
1551  // is our main variable
1552  auto vari = sym_stats.begin();
1553  while ((vari != sym_stats.end()) &&
1554  (((vari->ldeg_b == 0) && (vari->deg_b == 0)) ||
1555  ((vari->ldeg_a == 0) && (vari->deg_a == 0))))
1556  vari++;
1557 
1558  // No common symbols at all, just return 1:
1559  if (vari == sym_stats.end()) {
1560  // N.B: keep cofactors factored
1561  if (ca)
1562  *ca = a;
1563  if (cb)
1564  *cb = b;
1565  return _ex1;
1566  }
1567  // move symbol contained only in one of the polynomials to the end:
1568  rotate(sym_stats.begin(), vari, sym_stats.end());
1569 
1570  sym_desc_vec::const_iterator var = sym_stats.begin();
1571  const ex &x = var->sym;
1572 
1573  // Cancel trivial common factor
1574  int ldeg_a = var->ldeg_a;
1575  int ldeg_b = var->ldeg_b;
1576  int min_ldeg = std::min(ldeg_a,ldeg_b);
1577  if (min_ldeg > 0) {
1578  ex common = pow(x, min_ldeg);
1579  return gcd((aex / common).expand(), (bex / common).expand(), ca, cb, false) * common;
1580  }
1581 
1582  // Try to eliminate variables
1583  if (var->deg_a == 0 && var->deg_b != 0 ) {
1584  ex bex_u, bex_c, bex_p;
1585  bex.unitcontprim(x, bex_u, bex_c, bex_p);
1586  ex g = gcd(aex, bex_c, ca, cb, false);
1587  if (cb)
1588  *cb *= bex_u * bex_p;
1589  return g;
1590  } else if (var->deg_b == 0 && var->deg_a != 0) {
1591  ex aex_u, aex_c, aex_p;
1592  aex.unitcontprim(x, aex_u, aex_c, aex_p);
1593  ex g = gcd(aex_c, bex, ca, cb, false);
1594  if (ca)
1595  *ca *= aex_u * aex_p;
1596  return g;
1597  }
1598 
1599  // Try heuristic algorithm first, fall back to PRS if that failed
1600  ex g;
1601  if (!(options & gcd_options::no_heur_gcd)) {
1602  bool found = heur_gcd(g, aex, bex, ca, cb, var);
1603  if (found) {
1604  // heur_gcd have already computed cofactors...
1605  if (g.is_equal(_ex1)) {
1606  // ... but we want to keep them factored if possible.
1607  if (ca)
1608  *ca = a;
1609  if (cb)
1610  *cb = b;
1611  }
1612  return g;
1613  }
1614 #if STATISTICS
1615  else {
1616  heur_gcd_failed++;
1617  }
1618 #endif
1619  }
1621  g = sr_gcd(aex, bex, var);
1622  } else {
1623  exvector vars;
1624  for (std::size_t n = sym_stats.size(); n-- != 0; )
1625  vars.push_back(sym_stats[n].sym);
1626  g = chinrem_gcd(aex, bex, vars);
1627  }
1628 
1629  if (g.is_equal(_ex1)) {
1630  // Keep cofactors factored if possible
1631  if (ca)
1632  *ca = a;
1633  if (cb)
1634  *cb = b;
1635  } else {
1636  if (ca)
1637  divide(aex, g, *ca, false);
1638  if (cb)
1639  divide(bex, g, *cb, false);
1640  }
1641  return g;
1642 }
1643 
1644 // gcd helper to handle partially factored polynomials (to avoid expanding
1645 // large expressions). Both arguments should be powers.
1646 static ex gcd_pf_pow_pow(const ex& a, const ex& b, ex* ca, ex* cb)
1647 {
1648  ex p = a.op(0);
1649  const ex& exp_a = a.op(1);
1650  ex pb = b.op(0);
1651  const ex& exp_b = b.op(1);
1652 
1653  // a = p^n, b = p^m, gcd = p^min(n, m)
1654  if (p.is_equal(pb)) {
1655  if (exp_a < exp_b) {
1656  if (ca)
1657  *ca = _ex1;
1658  if (cb)
1659  *cb = pow(p, exp_b - exp_a);
1660  return pow(p, exp_a);
1661  } else {
1662  if (ca)
1663  *ca = pow(p, exp_a - exp_b);
1664  if (cb)
1665  *cb = _ex1;
1666  return pow(p, exp_b);
1667  }
1668  }
1669 
1670  ex p_co, pb_co;
1671  ex p_gcd = gcd(p, pb, &p_co, &pb_co, false);
1672  // a(x) = p(x)^n, b(x) = p_b(x)^m, gcd (p, p_b) = 1 ==> gcd(a,b) = 1
1673  if (p_gcd.is_equal(_ex1)) {
1674  if (ca)
1675  *ca = a;
1676  if (cb)
1677  *cb = b;
1678  return _ex1;
1679  }
1680 
1681  // there are common factors:
1682  // a(x) = g(x)^n A(x)^n, b(x) = g(x)^m B(x)^m ==>
1683  // gcd(a, b) = g(x)^n gcd(A(x)^n, g(x)^(n-m) B(x)^m
1684  if (exp_a < exp_b) {
1685  ex pg = gcd(pow(p_co, exp_a), pow(p_gcd, exp_b-exp_a)*pow(pb_co, exp_b), ca, cb, false);
1686  return pow(p_gcd, exp_a)*pg;
1687  } else {
1688  ex pg = gcd(pow(p_gcd, exp_a - exp_b)*pow(p_co, exp_a), pow(pb_co, exp_b), ca, cb, false);
1689  return pow(p_gcd, exp_b)*pg;
1690  }
1691 }
1692 
1693 static ex gcd_pf_pow(const ex& a, const ex& b, ex* ca, ex* cb)
1694 {
1695  if (is_exactly_a<power>(a) && is_exactly_a<power>(b))
1696  return gcd_pf_pow_pow(a, b, ca, cb);
1697 
1698  if (is_exactly_a<power>(b) && (! is_exactly_a<power>(a)))
1699  return gcd_pf_pow(b, a, cb, ca);
1700 
1701  GINAC_ASSERT(is_exactly_a<power>(a));
1702 
1703  ex p = a.op(0);
1704  const ex& exp_a = a.op(1);
1705  if (p.is_equal(b)) {
1706  // a = p^n, b = p, gcd = p
1707  if (ca)
1708  *ca = pow(p, exp_a - 1);
1709  if (cb)
1710  *cb = _ex1;
1711  return p;
1712  }
1713  if (is_a<symbol>(p)) {
1714  // Cancel trivial common factor
1715  int ldeg_a = ex_to<numeric>(exp_a).to_int();
1716  int ldeg_b = b.ldegree(p);
1717  int min_ldeg = std::min(ldeg_a, ldeg_b);
1718  if (min_ldeg > 0) {
1719  ex common = pow(p, min_ldeg);
1720  return gcd(pow(p, ldeg_a - min_ldeg), (b / common).expand(), ca, cb, false) * common;
1721  }
1722  }
1723 
1724  ex p_co, bpart_co;
1725  ex p_gcd = gcd(p, b, &p_co, &bpart_co, false);
1726 
1727  if (p_gcd.is_equal(_ex1)) {
1728  // a(x) = p(x)^n, gcd(p, b) = 1 ==> gcd(a, b) = 1
1729  if (ca)
1730  *ca = a;
1731  if (cb)
1732  *cb = b;
1733  return _ex1;
1734  }
1735  // a(x) = g(x)^n A(x)^n, b(x) = g(x) B(x) ==> gcd(a, b) = g(x) gcd(g(x)^(n-1) A(x)^n, B(x))
1736  ex rg = gcd(pow(p_gcd, exp_a-1)*pow(p_co, exp_a), bpart_co, ca, cb, false);
1737  return p_gcd*rg;
1738 }
1739 
1740 static ex gcd_pf_mul(const ex& a, const ex& b, ex* ca, ex* cb)
1741 {
1742  if (is_exactly_a<mul>(a) && is_exactly_a<mul>(b)
1743  && (b.nops() > a.nops()))
1744  return gcd_pf_mul(b, a, cb, ca);
1745 
1746  if (is_exactly_a<mul>(b) && (!is_exactly_a<mul>(a)))
1747  return gcd_pf_mul(b, a, cb, ca);
1748 
1749  GINAC_ASSERT(is_exactly_a<mul>(a));
1750  size_t num = a.nops();
1751  exvector g; g.reserve(num);
1752  exvector acc_ca; acc_ca.reserve(num);
1753  ex part_b = b;
1754  for (size_t i=0; i<num; i++) {
1755  ex part_ca, part_cb;
1756  g.push_back(gcd(a.op(i), part_b, &part_ca, &part_cb, false));
1757  acc_ca.push_back(part_ca);
1758  part_b = part_cb;
1759  }
1760  if (ca)
1761  *ca = dynallocate<mul>(acc_ca);
1762  if (cb)
1763  *cb = part_b;
1764  return dynallocate<mul>(g);
1765 }
1766 
1774 ex lcm(const ex &a, const ex &b, bool check_args)
1775 {
1776  if (is_exactly_a<numeric>(a) && is_exactly_a<numeric>(b))
1777  return lcm(ex_to<numeric>(a), ex_to<numeric>(b));
1779  throw(std::invalid_argument("lcm: arguments must be polynomials over the rationals"));
1780 
1781  ex ca, cb;
1782  ex g = gcd(a, b, &ca, &cb, false);
1783  return ca * cb * g;
1784 }
1785 
1786 
1787 /*
1788  * Square-free factorization
1789  */
1790 
1798 static epvector sqrfree_yun(const ex &a, const symbol &x)
1799 {
1800  ex w = a;
1801  ex z = w.diff(x);
1802  ex g = gcd(w, z);
1803  if (g.is_zero()) {
1804  // manifest zero or hidden zero
1805  return {};
1806  }
1807  if (g.is_equal(_ex1)) {
1808  // w(x) and w'(x) share no factors: w(x) is square-free
1809  return {expair(a, _ex1)};
1810  }
1811 
1812  epvector factors;
1813  ex i = 0; // exponent
1814  do {
1815  w = quo(w, g, x);
1816  if (w.is_zero()) {
1817  // hidden zero
1818  break;
1819  }
1820  z = quo(z, g, x) - w.diff(x);
1821  i += 1;
1822  if (w.is_equal(x)) {
1823  // shortcut for x^n with n ∈ ℕ
1824  i += quo(z, w.diff(x), x);
1825  factors.push_back(expair(w, i));
1826  break;
1827  }
1828  g = gcd(w, z);
1829  if (!g.is_equal(_ex1)) {
1830  factors.push_back(expair(g, i));
1831  }
1832  } while (!z.is_zero());
1833 
1834  // correct for lost factor
1835  // (being based on GCDs, Yun's algorithm only finds factors up to a unit)
1836  const ex lost_factor = quo(a, mul{factors}, x);
1837  if (lost_factor.is_equal(_ex1)) {
1838  // trivial lost factor
1839  return factors;
1840  }
1841  if (!factors.empty() && factors[0].coeff.is_equal(1)) {
1842  // multiply factor^1 with lost_factor
1843  factors[0].rest *= lost_factor;
1844  return factors;
1845  }
1846  // no factor^1: prepend lost_factor^1 to the results
1847  epvector results = {expair(lost_factor, 1)};
1848  std::move(factors.begin(), factors.end(), std::back_inserter(results));
1849  return results;
1850 }
1851 
1852 
1888 ex sqrfree(const ex &a, const lst &l)
1889 {
1890  if (is_exactly_a<numeric>(a) ||
1891  is_a<symbol>(a)) // shortcuts
1892  return a;
1893 
1894  // If no lst of variables to factorize in was specified we have to
1895  // invent one now. Maybe one can optimize here by reversing the order
1896  // or so, I don't know.
1897  lst args;
1898  if (l.nops()==0) {
1899  sym_desc_vec sdv;
1900  get_symbol_stats(a, _ex0, sdv);
1901  for (auto & it : sdv)
1902  args.append(it.sym);
1903  } else {
1904  args = l;
1905  }
1906 
1907  // Find the symbol to factor in at this stage
1908  if (!is_a<symbol>(args.op(0)))
1909  throw (std::runtime_error("sqrfree(): invalid factorization variable"));
1910  const symbol &x = ex_to<symbol>(args.op(0));
1911 
1912  // convert the argument from something in Q[X] to something in Z[X]
1914  const ex tmp = multiply_lcm(a, lcm);
1915 
1916  // find the factors
1917  epvector factors = sqrfree_yun(tmp, x);
1918  if (factors.empty()) {
1919  // the polynomial was a hidden zero
1920  return _ex0;
1921  }
1922 
1923  // remove symbol x and proceed recursively with the remaining symbols
1924  args.remove_first();
1925 
1926  // recurse down the factors in remaining variables
1927  if (args.nops()>0) {
1928  for (auto & it : factors)
1929  it.rest = sqrfree(it.rest, args);
1930  }
1931 
1932  // Done with recursion, now construct the final result
1933  ex result = mul(factors);
1934 
1935  // Put in the rational overall factor again and return
1936  return result * lcm.inverse();
1937 }
1938 
1939 
1947 ex sqrfree_parfrac(const ex & a, const symbol & x)
1948 {
1949  // Find numerator and denominator
1950  ex nd = numer_denom(a);
1951  ex numer = nd.op(0), denom = nd.op(1);
1952 //std::clog << "numer = " << numer << ", denom = " << denom << std::endl;
1953 
1954  // Convert N(x)/D(x) -> Q(x) + R(x)/D(x), so degree(R) < degree(D)
1955  ex red_poly = quo(numer, denom, x), red_numer = rem(numer, denom, x).expand();
1956 //std::clog << "red_poly = " << red_poly << ", red_numer = " << red_numer << std::endl;
1957 
1958  // Factorize denominator and compute cofactors
1959  epvector yun = sqrfree_yun(denom, x);
1960  size_t yun_max_exponent = yun.empty() ? 0 : ex_to<numeric>(yun.back().coeff).to_int();
1961  exvector factor, cofac;
1962  for (size_t i=0; i<yun.size(); i++) {
1963  numeric i_exponent = ex_to<numeric>(yun[i].coeff);
1964  for (size_t j=0; j<i_exponent; j++) {
1965  factor.push_back(pow(yun[i].rest, j+1));
1966  ex prod = _ex1;
1967  for (size_t k=0; k<yun.size(); k++) {
1968  if (yun[k].coeff == i_exponent)
1969  prod *= pow(yun[k].rest, i_exponent-1-j);
1970  else
1971  prod *= pow(yun[k].rest, yun[k].coeff);
1972  }
1973  cofac.push_back(prod.expand());
1974  }
1975  }
1976  size_t num_factors = factor.size();
1977 //std::clog << "factors : " << exprseq(factor) << std::endl;
1978 //std::clog << "cofactors: " << exprseq(cofac) << std::endl;
1979 
1980  // Construct coefficient matrix for decomposition
1981  int max_denom_deg = denom.degree(x);
1982  matrix sys(max_denom_deg + 1, num_factors);
1983  matrix rhs(max_denom_deg + 1, 1);
1984  for (int i=0; i<=max_denom_deg; i++) {
1985  for (size_t j=0; j<num_factors; j++)
1986  sys(i, j) = cofac[j].coeff(x, i);
1987  rhs(i, 0) = red_numer.coeff(x, i);
1988  }
1989 //std::clog << "coeffs: " << sys << std::endl;
1990 //std::clog << "rhs : " << rhs << std::endl;
1991 
1992  // Solve resulting linear system
1993  matrix vars(num_factors, 1);
1994  for (size_t i=0; i<num_factors; i++)
1995  vars(i, 0) = symbol();
1996  matrix sol = sys.solve(vars, rhs);
1997 
1998  // Sum up decomposed fractions
1999  ex sum = 0;
2000  for (size_t i=0; i<num_factors; i++)
2001  sum += sol(i, 0) / factor[i];
2002 
2003  return red_poly + sum;
2004 }
2005 
2006 
2007 /*
2008  * Normal form of rational functions
2009  */
2010 
2011 /*
2012  * Note: The internal normal() functions (= basic::normal() and overloaded
2013  * functions) all return lists of the form {numerator, denominator}. This
2014  * is to get around mul::eval()'s automatic expansion of numeric coefficients.
2015  * E.g. (a+b)/3 is automatically converted to a/3+b/3 but we want to keep
2016  * the information that (a+b) is the numerator and 3 is the denominator.
2017  */
2018 
2019 
2052 static ex replace_with_symbol(const ex & e, exmap & repl, exmap & rev_lookup, lst & modifier)
2053 {
2054  // Since the repl contains replaced expressions we should search for them
2055  ex e_replaced = e.subs(repl, subs_options::no_pattern);
2056 
2057  // Expression already replaced? Then return the assigned symbol
2058  auto it = rev_lookup.find(e_replaced);
2059  if (it != rev_lookup.end())
2060  return it->second;
2061 
2062  // The expression can be the base of substituted power, which requires a more careful search
2063  if (! is_a<numeric>(e_replaced))
2064  for (auto & it : repl)
2065  if (is_a<power>(it.second) && e_replaced.is_equal(it.second.op(0))) {
2066  ex degree = pow(it.second.op(1), _ex_1);
2067  if (is_a<numeric>(degree) && ex_to<numeric>(degree).is_integer())
2068  return pow(it.first, degree);
2069  }
2070 
2071  // We treat powers and the exponent functions differently because
2072  // they can be rationalised more efficiently
2073  if (is_a<function>(e_replaced) && is_ex_the_function(e_replaced, exp)) {
2074  for (auto & it : repl) {
2075  if (is_a<function>(it.second) && is_ex_the_function(e_replaced, exp)) {
2076  ex ratio = normal(e_replaced.op(0) / it.second.op(0));
2077  if (is_a<numeric>(ratio) && ex_to<numeric>(ratio).is_rational()) {
2078  // Different exponents can be treated as powers of the same basic equation
2079  if (ex_to<numeric>(ratio).is_integer()) {
2080  // If ratio is an integer then this is simply the power of the existing symbol.
2081  // std::clog << e_replaced << " is a " << ratio << " power of " << it.first << std::endl;
2082  return dynallocate<power>(it.first, ratio);
2083  } else {
2084  // otherwise we need to give the replacement pattern to change
2085  // the previous expression...
2086  ex es = dynallocate<symbol>();
2087  ex Num = numer(ratio);
2088  modifier.append(it.first == power(es, denom(ratio)));
2089  // std::clog << e_replaced << " is power " << Num << " and "
2090  // << it.first << " is power " << denom(ratio) << " of the common base "
2091  // << exp(e_replaced.op(0)/Num) << std::endl;
2092  // ... and modify the replacement tables
2093  rev_lookup.erase(it.second);
2094  rev_lookup.insert({exp(e_replaced.op(0)/Num), es});
2095  repl.erase(it.first);
2096  repl.insert({es, exp(e_replaced.op(0)/Num)});
2097  return dynallocate<power>(es, Num);
2098  }
2099  }
2100  }
2101  }
2102  } else if (is_a<power>(e_replaced) && !is_a<numeric>(e_replaced.op(0)) // We do not replace simple monomials like x^3 or sqrt(2)
2103  && ! (is_a<symbol>(e_replaced.op(0))
2104  && is_a<numeric>(e_replaced.op(1)) && ex_to<numeric>(e_replaced.op(1)).is_integer())) {
2105  for (auto & it : repl) {
2106  if (e_replaced.op(0).is_equal(it.second) // The base is an allocated symbol or base of power
2107  || (is_a<power>(it.second) && e_replaced.op(0).is_equal(it.second.op(0)))) {
2108  ex ratio; // We bind together two above cases
2109  if (is_a<power>(it.second))
2110  ratio = normal(e_replaced.op(1) / it.second.op(1));
2111  else
2112  ratio = e_replaced.op(1);
2113  if (is_a<numeric>(ratio) && ex_to<numeric>(ratio).is_rational()) {
2114  // Different powers can be treated as powers of the same basic equation
2115  if (ex_to<numeric>(ratio).is_integer()) {
2116  // If ratio is an integer then this is simply the power of the existing symbol.
2117  //std::clog << e_replaced << " is a " << ratio << " power of " << it.first << std::endl;
2118  return dynallocate<power>(it.first, ratio);
2119  } else {
2120  // otherwise we need to give the replacement pattern to change
2121  // the previous expression...
2122  ex es = dynallocate<symbol>();
2123  ex Num = numer(ratio);
2124  modifier.append(it.first == power(es, denom(ratio)));
2125  //std::clog << e_replaced << " is power " << Num << " and "
2126  // << it.first << " is power " << denom(ratio) << " of the common base "
2127  // << pow(e_replaced.op(0), e_replaced.op(1)/Num) << std::endl;
2128  // ... and modify the replacement tables
2129  rev_lookup.erase(it.second);
2130  rev_lookup.insert({pow(e_replaced.op(0), e_replaced.op(1)/Num), es});
2131  repl.erase(it.first);
2132  repl.insert({es, pow(e_replaced.op(0), e_replaced.op(1)/Num)});
2133  return dynallocate<power>(es, Num);
2134  }
2135  }
2136  }
2137  }
2138  // There is no existing substitution, thus we are creating a new one.
2139  // This needs to be done separately to treat possible occurrences of
2140  // b = e_replaced.op(0) elsewhere in the expression as pow(b, 1).
2141  ex degree = pow(e_replaced.op(1), _ex_1);
2142  if (is_a<numeric>(degree) && ex_to<numeric>(degree).is_integer()) {
2143  ex es = dynallocate<symbol>();
2144  modifier.append(e_replaced.op(0) == power(es, degree));
2145  repl.insert({es, e_replaced});
2146  rev_lookup.insert({e_replaced, es});
2147  return es;
2148  }
2149  }
2150 
2151  // Otherwise create new symbol and add to list, taking care that the
2152  // replacement expression doesn't itself contain symbols from repl,
2153  // because subs() is not recursive
2154  ex es = dynallocate<symbol>();
2155  repl.insert(std::make_pair(es, e_replaced));
2156  rev_lookup.insert(std::make_pair(e_replaced, es));
2157  return es;
2158 }
2159 
2165 static ex replace_with_symbol(const ex & e, exmap & repl)
2166 {
2167  // Since the repl contains replaced expressions we should search for them
2168  ex e_replaced = e.subs(repl, subs_options::no_pattern);
2169 
2170  // Expression already replaced? Then return the assigned symbol
2171  for (auto & it : repl)
2172  if (it.second.is_equal(e_replaced))
2173  return it.first;
2174 
2175  // Otherwise create new symbol and add to list, taking care that the
2176  // replacement expression doesn't itself contain symbols from repl,
2177  // because subs() is not recursive
2178  ex es = dynallocate<symbol>();
2179  repl.insert(std::make_pair(es, e_replaced));
2180  return es;
2181 }
2182 
2183 
2186  ex operator()(const ex & e) override { return normal(e); }
2187 };
2188 
2192 ex basic::normal(exmap & repl, exmap & rev_lookup, lst & modifier) const
2193 {
2194  if (nops() == 0)
2195  return dynallocate<lst>({replace_with_symbol(*this, repl, rev_lookup, modifier), _ex1});
2196 
2197  normal_map_function map_normal;
2198  int nmod = modifier.nops(); // To watch new modifiers to the replacement list
2199  lst result = dynallocate<lst>({replace_with_symbol(map(map_normal), repl, rev_lookup, modifier), _ex1});
2200  for (int imod = nmod; imod < modifier.nops(); ++imod) {
2201  exmap this_repl;
2202  this_repl.insert(std::make_pair(modifier.op(imod).op(0), modifier.op(imod).op(1)));
2203  result = ex_to<lst>(result.subs(this_repl, subs_options::no_pattern));
2204  }
2205 
2206  return result;
2207 }
2208 
2209 
2212 ex symbol::normal(exmap & repl, exmap & rev_lookup, lst & modifier) const
2213 {
2214  return dynallocate<lst>({*this, _ex1});
2215 }
2216 
2217 
2222 ex numeric::normal(exmap & repl, exmap & rev_lookup, lst & modifier) const
2223 {
2224  numeric num = numer();
2225  ex numex = num;
2226 
2227  if (num.is_real()) {
2228  if (!num.is_integer())
2229  numex = replace_with_symbol(numex, repl, rev_lookup, modifier);
2230  } else { // complex
2231  numeric re = num.real(), im = num.imag();
2232  ex re_ex = re.is_rational() ? re : replace_with_symbol(re, repl, rev_lookup, modifier);
2233  ex im_ex = im.is_rational() ? im : replace_with_symbol(im, repl, rev_lookup, modifier);
2234  numex = re_ex + im_ex * replace_with_symbol(I, repl, rev_lookup, modifier);
2235  }
2236 
2237  // Denominator is always a real integer (see numeric::denom())
2238  return dynallocate<lst>({numex, denom()});
2239 }
2240 
2241 
2246 static ex frac_cancel(const ex &n, const ex &d)
2247 {
2248  ex num = n;
2249  ex den = d;
2250  numeric pre_factor = *_num1_p;
2251 
2252 //std::clog << "frac_cancel num = " << num << ", den = " << den << std::endl;
2253 
2254  // Handle trivial case where denominator is 1
2255  if (den.is_equal(_ex1))
2256  return dynallocate<lst>({num, den});
2257 
2258  // Handle special cases where numerator or denominator is 0
2259  if (num.is_zero())
2260  return dynallocate<lst>({num, _ex1});
2261  if (den.expand().is_zero())
2262  throw(std::overflow_error("frac_cancel: division by zero in frac_cancel"));
2263 
2264  // Bring numerator and denominator to Z[X] by multiplying with
2265  // LCM of all coefficients' denominators
2268  num = multiply_lcm(num, num_lcm);
2269  den = multiply_lcm(den, den_lcm);
2270  pre_factor = den_lcm / num_lcm;
2271 
2272  // Cancel GCD from numerator and denominator
2273  ex cnum, cden;
2274  if (gcd(num, den, &cnum, &cden, false) != _ex1) {
2275  num = cnum;
2276  den = cden;
2277  }
2278 
2279  // Make denominator unit normal (i.e. coefficient of first symbol
2280  // as defined by get_first_symbol() is made positive)
2281  if (is_exactly_a<numeric>(den)) {
2282  if (ex_to<numeric>(den).is_negative()) {
2283  num *= _ex_1;
2284  den *= _ex_1;
2285  }
2286  } else {
2287  ex x;
2288  if (get_first_symbol(den, x)) {
2289  GINAC_ASSERT(is_exactly_a<numeric>(den.unit(x)));
2290  if (ex_to<numeric>(den.unit(x)).is_negative()) {
2291  num *= _ex_1;
2292  den *= _ex_1;
2293  }
2294  }
2295  }
2296 
2297  // Return result as list
2298 //std::clog << " returns num = " << num << ", den = " << den << ", pre_factor = " << pre_factor << std::endl;
2299  return dynallocate<lst>({num * pre_factor.numer(), den * pre_factor.denom()});
2300 }
2301 
2302 
2306 ex add::normal(exmap & repl, exmap & rev_lookup, lst & modifier) const
2307 {
2308  // Normalize children and split each one into numerator and denominator
2309  exvector nums, dens;
2310  nums.reserve(seq.size()+1);
2311  dens.reserve(seq.size()+1);
2312  int nmod = modifier.nops(); // To watch new modifiers to the replacement list
2313  for (auto & it : seq) {
2314  ex n = ex_to<basic>(recombine_pair_to_ex(it)).normal(repl, rev_lookup, modifier);
2315  nums.push_back(n.op(0));
2316  dens.push_back(n.op(1));
2317  }
2318  ex n = ex_to<numeric>(overall_coeff).normal(repl, rev_lookup, modifier);
2319  nums.push_back(n.op(0));
2320  dens.push_back(n.op(1));
2321  GINAC_ASSERT(nums.size() == dens.size());
2322 
2323  // Now, nums is a vector of all numerators and dens is a vector of
2324  // all denominators
2325 //std::clog << "add::normal uses " << nums.size() << " summands:\n";
2326 
2327  // Add fractions sequentially
2328  auto num_it = nums.begin(), num_itend = nums.end();
2329  auto den_it = dens.begin(), den_itend = dens.end();
2330 //std::clog << " num = " << *num_it << ", den = " << *den_it << std::endl;
2331  for (int imod = nmod; imod < modifier.nops(); ++imod) {
2332  while (num_it != num_itend) {
2333  *num_it = num_it->subs(modifier.op(imod), subs_options::no_pattern);
2334  ++num_it;
2335  *den_it = den_it->subs(modifier.op(imod), subs_options::no_pattern);
2336  ++den_it;
2337  }
2338  // Reset iterators for the next round
2339  num_it = nums.begin();
2340  den_it = dens.begin();
2341  }
2342 
2343  ex num = *num_it++, den = *den_it++;
2344  while (num_it != num_itend) {
2345 //std::clog << " num = " << *num_it << ", den = " << *den_it << std::endl;
2346  ex next_num = *num_it++, next_den = *den_it++;
2347 
2348  // Trivially add sequences of fractions with identical denominators
2349  while ((den_it != den_itend) && next_den.is_equal(*den_it)) {
2350  next_num += *num_it;
2351  num_it++; den_it++;
2352  }
2353 
2354  // Addition of two fractions, taking advantage of the fact that
2355  // the heuristic GCD algorithm computes the cofactors at no extra cost
2356  ex co_den1, co_den2;
2357  ex g = gcd(den, next_den, &co_den1, &co_den2, false);
2358  num = ((num * co_den2) + (next_num * co_den1)).expand();
2359  den *= co_den2; // this is the lcm(den, next_den)
2360  }
2361 //std::clog << " common denominator = " << den << std::endl;
2362 
2363  // Cancel common factors from num/den
2364  return frac_cancel(num, den);
2365 }
2366 
2367 
2371 ex mul::normal(exmap & repl, exmap & rev_lookup, lst & modifier) const
2372 {
2373  // Normalize children, separate into numerator and denominator
2374  exvector num; num.reserve(seq.size());
2375  exvector den; den.reserve(seq.size());
2376  ex n;
2377  int nmod = modifier.nops(); // To watch new modifiers to the replacement list
2378  for (auto & it : seq) {
2379  n = ex_to<basic>(recombine_pair_to_ex(it)).normal(repl, rev_lookup, modifier);
2380  num.push_back(n.op(0));
2381  den.push_back(n.op(1));
2382  }
2383  n = ex_to<numeric>(overall_coeff).normal(repl, rev_lookup, modifier);
2384  num.push_back(n.op(0));
2385  den.push_back(n.op(1));
2386  auto num_it = num.begin(), num_itend = num.end();
2387  auto den_it = den.begin(), den_itend = den.end();
2388  for (int imod = nmod; imod < modifier.nops(); ++imod) {
2389  while (num_it != num_itend) {
2390  *num_it = num_it->subs(modifier.op(imod), subs_options::no_pattern);
2391  ++num_it;
2392  *den_it = den_it->subs(modifier.op(imod), subs_options::no_pattern);
2393  ++den_it;
2394  }
2395  num_it = num.begin();
2396  den_it = den.begin();
2397  }
2398 
2399  // Perform fraction cancellation
2400  return frac_cancel(dynallocate<mul>(num), dynallocate<mul>(den));
2401 }
2402 
2403 
2408 ex power::normal(exmap & repl, exmap & rev_lookup, lst & modifier) const
2409 {
2410  // Normalize basis and exponent (exponent gets reassembled)
2411  int nmod = modifier.nops(); // To watch new modifiers to the replacement list
2412  ex n_basis = ex_to<basic>(basis).normal(repl, rev_lookup, modifier);
2413  for (int imod = nmod; imod < modifier.nops(); ++imod)
2414  n_basis = n_basis.subs(modifier.op(imod), subs_options::no_pattern);
2415 
2416  nmod = modifier.nops();
2417  ex n_exponent = ex_to<basic>(exponent).normal(repl, rev_lookup, modifier);
2418  for (int imod = nmod; imod < modifier.nops(); ++imod)
2419  n_exponent = n_exponent.subs(modifier.op(imod), subs_options::no_pattern);
2420  n_exponent = n_exponent.op(0) / n_exponent.op(1);
2421 
2422  if (n_exponent.info(info_flags::integer)) {
2423 
2424  if (n_exponent.info(info_flags::positive)) {
2425 
2426  // (a/b)^n -> {a^n, b^n}
2427  return dynallocate<lst>({pow(n_basis.op(0), n_exponent), pow(n_basis.op(1), n_exponent)});
2428 
2429  } else if (n_exponent.info(info_flags::negative)) {
2430 
2431  // (a/b)^-n -> {b^n, a^n}
2432  return dynallocate<lst>({pow(n_basis.op(1), -n_exponent), pow(n_basis.op(0), -n_exponent)});
2433  }
2434 
2435  } else {
2436 
2437  if (n_exponent.info(info_flags::positive)) {
2438 
2439  // (a/b)^x -> {sym((a/b)^x), 1}
2440  return dynallocate<lst>({replace_with_symbol(pow(n_basis.op(0) / n_basis.op(1), n_exponent), repl, rev_lookup, modifier), _ex1});
2441 
2442  } else if (n_exponent.info(info_flags::negative)) {
2443 
2444  if (n_basis.op(1).is_equal(_ex1)) {
2445 
2446  // a^-x -> {1, sym(a^x)}
2447  return dynallocate<lst>({_ex1, replace_with_symbol(pow(n_basis.op(0), -n_exponent), repl, rev_lookup, modifier)});
2448 
2449  } else {
2450 
2451  // (a/b)^-x -> {sym((b/a)^x), 1}
2452  return dynallocate<lst>({replace_with_symbol(pow(n_basis.op(1) / n_basis.op(0), -n_exponent), repl, rev_lookup, modifier), _ex1});
2453  }
2454  }
2455  }
2456 
2457  // (a/b)^x -> {sym((a/b)^x, 1}
2458  return dynallocate<lst>({replace_with_symbol(pow(n_basis.op(0) / n_basis.op(1), n_exponent), repl, rev_lookup, modifier), _ex1});
2459 }
2460 
2461 
2465 ex pseries::normal(exmap & repl, exmap & rev_lookup, lst & modifier) const
2466 {
2467  epvector newseq;
2468  for (auto & it : seq) {
2469  ex restexp = it.rest.normal();
2470  if (!restexp.is_zero())
2471  newseq.push_back(expair(restexp, it.coeff));
2472  }
2473  ex n = pseries(relational(var,point), std::move(newseq));
2474  return dynallocate<lst>({replace_with_symbol(n, repl, rev_lookup, modifier), _ex1});
2475 }
2476 
2477 
2490 {
2491  exmap repl, rev_lookup;
2492  lst modifier;
2493 
2494  ex e = bp->normal(repl, rev_lookup, modifier);
2495  GINAC_ASSERT(is_a<lst>(e));
2496 
2497  // Re-insert replaced symbols
2498  if (!repl.empty()) {
2499  for(int i=0; i < modifier.nops(); ++i)
2500  e = e.subs(modifier.op(i), subs_options::no_pattern);
2501  e = e.subs(repl, subs_options::no_pattern);
2502  }
2503 
2504  // Convert {numerator, denominator} form back to fraction
2505  return e.op(0) / e.op(1);
2506 }
2507 
2514 ex ex::numer() const
2515 {
2516  exmap repl, rev_lookup;
2517  lst modifier;
2518 
2519  ex e = bp->normal(repl, rev_lookup, modifier);
2520  GINAC_ASSERT(is_a<lst>(e));
2521 
2522  // Re-insert replaced symbols
2523  if (repl.empty())
2524  return e.op(0);
2525  else {
2526  for(int i=0; i < modifier.nops(); ++i)
2527  e = e.subs(modifier.op(i), subs_options::no_pattern);
2528 
2529  return e.op(0).subs(repl, subs_options::no_pattern);
2530  }
2531 }
2532 
2539 ex ex::denom() const
2540 {
2541  exmap repl, rev_lookup;
2542  lst modifier;
2543 
2544  ex e = bp->normal(repl, rev_lookup, modifier);
2545  GINAC_ASSERT(is_a<lst>(e));
2546 
2547  // Re-insert replaced symbols
2548  if (repl.empty())
2549  return e.op(1);
2550  else {
2551  for(int i=0; i < modifier.nops(); ++i)
2552  e = e.subs(modifier.op(i), subs_options::no_pattern);
2553 
2554  return e.op(1).subs(repl, subs_options::no_pattern);
2555  }
2556 }
2557 
2565 {
2566  exmap repl, rev_lookup;
2567  lst modifier;
2568 
2569  ex e = bp->normal(repl, rev_lookup, modifier);
2570  GINAC_ASSERT(is_a<lst>(e));
2571 
2572  // Re-insert replaced symbols
2573  if (repl.empty())
2574  return e;
2575  else {
2576  for(int i=0; i < modifier.nops(); ++i)
2577  e = e.subs(modifier.op(i), subs_options::no_pattern);
2578 
2579  return e.subs(repl, subs_options::no_pattern);
2580  }
2581 }
2582 
2583 
2597 ex ex::to_rational(exmap & repl) const
2598 {
2599  return bp->to_rational(repl);
2600 }
2601 
2603 {
2604  return bp->to_polynomial(repl);
2605 }
2606 
2610 {
2611  return replace_with_symbol(*this, repl);
2612 }
2613 
2615 {
2616  return replace_with_symbol(*this, repl);
2617 }
2618 
2619 
2623 {
2624  return *this;
2625 }
2626 
2630 {
2631  return *this;
2632 }
2633 
2634 
2639 {
2640  if (is_real()) {
2641  if (!is_rational())
2642  return replace_with_symbol(*this, repl);
2643  } else { // complex
2644  numeric re = real();
2645  numeric im = imag();
2646  ex re_ex = re.is_rational() ? re : replace_with_symbol(re, repl);
2647  ex im_ex = im.is_rational() ? im : replace_with_symbol(im, repl);
2648  return re_ex + im_ex * replace_with_symbol(I, repl);
2649  }
2650  return *this;
2651 }
2652 
2657 {
2658  if (is_real()) {
2659  if (!is_integer())
2660  return replace_with_symbol(*this, repl);
2661  } else { // complex
2662  numeric re = real();
2663  numeric im = imag();
2664  ex re_ex = re.is_integer() ? re : replace_with_symbol(re, repl);
2665  ex im_ex = im.is_integer() ? im : replace_with_symbol(im, repl);
2666  return re_ex + im_ex * replace_with_symbol(I, repl);
2667  }
2668  return *this;
2669 }
2670 
2671 
2675 {
2677  return pow(basis.to_rational(repl), exponent);
2678  else
2679  return replace_with_symbol(*this, repl);
2680 }
2681 
2685 {
2687  return pow(basis.to_rational(repl), exponent);
2688  else if (exponent.info(info_flags::negint))
2689  {
2690  ex basis_pref = collect_common_factors(basis);
2691  if (is_exactly_a<mul>(basis_pref) || is_exactly_a<power>(basis_pref)) {
2692  // (A*B)^n will be automagically transformed to A^n*B^n
2693  ex t = pow(basis_pref, exponent);
2694  return t.to_polynomial(repl);
2695  }
2696  else
2697  return pow(replace_with_symbol(pow(basis, _ex_1), repl), -exponent);
2698  }
2699  else
2700  return replace_with_symbol(*this, repl);
2701 }
2702 
2703 
2706 {
2707  epvector s;
2708  s.reserve(seq.size());
2709  for (auto & it : seq)
2710  s.push_back(split_ex_to_pair(recombine_pair_to_ex(it).to_rational(repl)));
2711 
2712  ex oc = overall_coeff.to_rational(repl);
2713  if (oc.info(info_flags::numeric))
2714  return thisexpairseq(std::move(s), overall_coeff);
2715  else
2716  s.push_back(expair(oc, _ex1));
2717  return thisexpairseq(std::move(s), default_overall_coeff());
2718 }
2719 
2722 {
2723  epvector s;
2724  s.reserve(seq.size());
2725  for (auto & it : seq)
2726  s.push_back(split_ex_to_pair(recombine_pair_to_ex(it).to_polynomial(repl)));
2727 
2728  ex oc = overall_coeff.to_polynomial(repl);
2729  if (oc.info(info_flags::numeric))
2730  return thisexpairseq(std::move(s), overall_coeff);
2731  else
2732  s.push_back(expair(oc, _ex1));
2733  return thisexpairseq(std::move(s), default_overall_coeff());
2734 }
2735 
2736 
2740 static ex find_common_factor(const ex & e, ex & factor, exmap & repl)
2741 {
2742  if (is_exactly_a<add>(e)) {
2743 
2744  size_t num = e.nops();
2745  exvector terms; terms.reserve(num);
2746  ex gc;
2747 
2748  // Find the common GCD
2749  for (size_t i=0; i<num; i++) {
2750  ex x = e.op(i).to_polynomial(repl);
2751 
2752  if (is_exactly_a<add>(x) || is_exactly_a<mul>(x) || is_a<power>(x)) {
2753  ex f = 1;
2754  x = find_common_factor(x, f, repl);
2755  x *= f;
2756  }
2757 
2758  if (gc.is_zero())
2759  gc = x;
2760  else
2761  gc = gcd(gc, x);
2762 
2763  terms.push_back(x);
2764  }
2765 
2766  if (gc.is_equal(_ex1))
2767  return e;
2768 
2769  if (gc.is_zero())
2770  return _ex0;
2771 
2772  // The GCD is the factor we pull out
2773  factor *= gc;
2774 
2775  // Now divide all terms by the GCD
2776  for (size_t i=0; i<num; i++) {
2777  ex x;
2778 
2779  // Try to avoid divide() because it expands the polynomial
2780  ex &t = terms[i];
2781  if (is_exactly_a<mul>(t)) {
2782  for (size_t j=0; j<t.nops(); j++) {
2783  if (t.op(j).is_equal(gc)) {
2784  exvector v; v.reserve(t.nops());
2785  for (size_t k=0; k<t.nops(); k++) {
2786  if (k == j)
2787  v.push_back(_ex1);
2788  else
2789  v.push_back(t.op(k));
2790  }
2791  t = dynallocate<mul>(v);
2792  goto term_done;
2793  }
2794  }
2795  }
2796 
2797  divide(t, gc, x);
2798  t = x;
2799 term_done: ;
2800  }
2801  return dynallocate<add>(terms);
2802 
2803  } else if (is_exactly_a<mul>(e)) {
2804 
2805  size_t num = e.nops();
2806  exvector v; v.reserve(num);
2807 
2808  for (size_t i=0; i<num; i++)
2809  v.push_back(find_common_factor(e.op(i), factor, repl));
2810 
2811  return dynallocate<mul>(v);
2812 
2813  } else if (is_exactly_a<power>(e)) {
2814  const ex e_exp(e.op(1));
2815  if (e_exp.info(info_flags::integer)) {
2816  ex eb = e.op(0).to_polynomial(repl);
2817  ex factor_local(_ex1);
2818  ex pre_res = find_common_factor(eb, factor_local, repl);
2819  factor *= pow(factor_local, e_exp);
2820  return pow(pre_res, e_exp);
2821 
2822  } else
2823  return e.to_polynomial(repl);
2824 
2825  } else
2826  return e;
2827 }
2828 
2829 
2833 {
2834  if (is_exactly_a<add>(e) || is_exactly_a<mul>(e) || is_exactly_a<power>(e)) {
2835 
2836  exmap repl;
2837  ex factor = 1;
2838  ex r = find_common_factor(e, factor, repl);
2839  return factor.subs(repl, subs_options::no_pattern) * r.subs(repl, subs_options::no_pattern);
2840 
2841  } else
2842  return e;
2843 }
2844 
2845 
2848 ex resultant(const ex & e1, const ex & e2, const ex & s)
2849 {
2850  const ex ee1 = e1.expand();
2851  const ex ee2 = e2.expand();
2852  if (!ee1.info(info_flags::polynomial) ||
2854  throw(std::runtime_error("resultant(): arguments must be polynomials"));
2855 
2856  const int h1 = ee1.degree(s);
2857  const int l1 = ee1.ldegree(s);
2858  const int h2 = ee2.degree(s);
2859  const int l2 = ee2.ldegree(s);
2860 
2861  const int msize = h1 + h2;
2862  matrix m(msize, msize);
2863 
2864  for (int l = h1; l >= l1; --l) {
2865  const ex e = ee1.coeff(s, l);
2866  for (int k = 0; k < h2; ++k)
2867  m(k, k+h1-l) = e;
2868  }
2869  for (int l = h2; l >= l2; --l) {
2870  const ex e = ee2.coeff(s, l);
2871  for (int k = 0; k < h1; ++k)
2872  m(k+h2, k+h2-l) = e;
2873  }
2874 
2875  return m.determinant();
2876 }
2877 
2878 
2879 } // namespace GiNaC
inifcns.h
Interface to GiNaC's initially known functions.
GiNaC::container::append
container & append(const ex &b)
Add element at back.
Definition: container.h:391
GiNaC::epvector
std::vector< expair > epvector
expair-vector
Definition: expairseq.h:33
GiNaC::ex::expand
ex expand(unsigned options=0) const
Definition: ex.cpp:73
fail.h
Interface to class signaling failure of operation.
GiNaC::heur_gcd
static bool heur_gcd(ex &res, const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const_iterator var)
Compute GCD of multivariate polynomials using the heuristic GCD algorithm.
Definition: normal.cpp:1371
constant.h
Interface to GiNaC's constant types and some special constants.
GiNaC::sr_gcd
static ex sr_gcd(const ex &a, const ex &b, sym_desc_vec::const_iterator var)
Compute GCD of multivariate polynomials using the subresultant PRS algorithm.
Definition: normal.cpp:1081
GiNaC::rhs
ex rhs(const ex &thisex)
Definition: ex.h:817
GiNaC::info_flags::integer
@ integer
Definition: flags.h:223
GiNaC::numeric::is_integer
bool is_integer() const
True if object is a non-complex integer.
Definition: numeric.cpp:1154
GiNaC::add::max_coefficient
numeric max_coefficient() const override
Implementation ex::max_coefficient().
Definition: normal.cpp:1171
r
size_t r
Definition: factor.cpp:770
GiNaC::multiply_lcm
static ex multiply_lcm(const ex &e, const numeric &lcm)
Bring polynomial from Q[X] to Z[X] by multiplying in the previously determined LCM of the coefficient...
Definition: normal.cpp:271
GiNaC::info_flags::negative
@ negative
Definition: flags.h:227
mul.h
Interface to GiNaC's products of expressions.
GiNaC::ex::to_rational
ex to_rational(exmap &repl) const
Rationalization of non-rational functions.
Definition: normal.cpp:2597
GiNaC::ex::normal
ex normal() const
Normalization of rational functions.
Definition: normal.cpp:2489
GiNaC::quo
ex quo(const ex &a, const ex &b, const ex &x, bool check_args)
Quotient q(x) of polynomials a(x) and b(x) in Q[x].
Definition: normal.cpp:373
GiNaC::pseries::normal
ex normal(exmap &repl, exmap &rev_lookup, lst &modifier) const override
Implementation of ex::normal() for pseries.
Definition: normal.cpp:2465
GiNaC::symbol::normal
ex normal(exmap &repl, exmap &rev_lookup, lst &modifier) const override
Implementation of ex::normal() for symbols.
Definition: normal.cpp:2212
GiNaC::pseries::pseries
pseries(const ex &rel_, const epvector &ops_)
Construct pseries from a vector of coefficients and powers.
Definition: pseries.cpp:71
GiNaC::ex::coeff
ex coeff(const ex &s, int n=1) const
Definition: ex.h:175
GiNaC::_ex0
const ex _ex0
Definition: utils.cpp:177
GiNaC::iquo
const numeric iquo(const numeric &a, const numeric &b)
Numeric integer quotient.
Definition: numeric.cpp:2404
GiNaC::mul::normal
ex normal(exmap &repl, exmap &rev_lookup, lst &modifier) const override
Implementation of ex::normal() for a product.
Definition: normal.cpp:2371
GiNaC::relational
This class holds a relation consisting of two expressions and a logical relation between them.
Definition: relational.h:35
numeric.h
Makes the interface to the underlying bignum package available.
GiNaC::map_function
Function object for map().
Definition: basic.h:85
GiNaC::mul::max_coefficient
numeric max_coefficient() const override
Implementation ex::max_coefficient().
Definition: normal.cpp:1185
GiNaC::frac_cancel
static ex frac_cancel(const ex &n, const ex &d)
Fraction cancellation.
Definition: normal.cpp:2246
GiNaC::numeric::to_rational
ex to_rational(exmap &repl) const override
Implementation of ex::to_rational() for a numeric.
Definition: normal.cpp:2638
GiNaC::ex::subs
ex subs(const exmap &m, unsigned options=0) const
Definition: ex.h:826
GiNaC::numeric::is_rational
bool is_rational() const
True if object is an exact rational number, may even be complex (denominator may be unity).
Definition: numeric.cpp:1201
add.h
Interface to GiNaC's sums of expressions.
k
vector< int > k
Definition: factor.cpp:1466
GiNaC::psi
function psi(const T1 &p1)
Definition: inifcns.h:165
GiNaC::numeric::imag
const numeric imag() const
Imaginary part of a number.
Definition: numeric.cpp:1346
GiNaC::ex::denom
ex denom() const
Get denominator of an expression.
Definition: normal.cpp:2539
GiNaC::sym_desc::sym_desc
sym_desc(const ex &s)
Initialize symbol, leave other variables uninitialized.
Definition: normal.cpp:124
power.h
Interface to GiNaC's symbolic exponentiation (basis^exponent).
GiNaC::exvector
std::vector< ex > exvector
Definition: basic.h:46
GiNaC::ex::nops
size_t nops() const
Definition: ex.h:135
GiNaC::info_flags::rational_polynomial
@ rational_polynomial
Definition: flags.h:258
GiNaC::power
This class holds a two-component object, a basis and and exponent representing exponentiation.
Definition: power.h:39
GiNaC::sqrfree
ex sqrfree(const ex &a, const lst &l)
Compute a square-free factorization of a multivariate polynomial in Q[X].
Definition: normal.cpp:1888
GiNaC::status_flags::evaluated
@ evaluated
.eval() has already done its job
Definition: flags.h:203
GiNaC::sym_desc::max_lcnops
size_t max_lcnops
Maximum number of terms of leading coefficient of symbol in both polynomials.
Definition: normal.cpp:147
GiNaC::is_rational
bool is_rational(const numeric &x)
Definition: numeric.h:290
GiNaC::numeric::is_real
bool is_real() const
True if object is a real integer, rational or float (but not complex).
Definition: numeric.cpp:1208
GiNaC::ex::degree
int degree(const ex &s) const
Definition: ex.h:173
GiNaC::expairseq::split_ex_to_pair
virtual expair split_ex_to_pair(const ex &e) const
GiNaC::power::exponent
ex exponent
Definition: power.h:106
GiNaC::heur_gcd_z
static bool heur_gcd_z(ex &res, const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const_iterator var)
Compute GCD of multivariate polynomials using the heuristic GCD algorithm.
Definition: normal.cpp:1277
GiNaC::basic::normal
virtual ex normal(exmap &repl, exmap &rev_lookup, lst &modifier) const
Default implementation of ex::normal().
Definition: normal.cpp:2192
GiNaC::resultant
ex resultant(const ex &e1, const ex &e2, const ex &s)
Resultant of two expressions e1,e2 with respect to symbol s.
Definition: normal.cpp:2848
GiNaC::sym_desc::deg_b
int deg_b
Highest degree of symbol in polynomial "b".
Definition: normal.cpp:135
GiNaC::_ex1
const ex _ex1
Definition: utils.cpp:193
GiNaC::add_symbol
static void add_symbol(const ex &s, sym_desc_vec &v)
Definition: normal.cpp:163
GiNaC::basic::nops
virtual size_t nops() const
Number of operands/members.
Definition: basic.cpp:229
GiNaC::sym_desc_vec
std::vector< sym_desc > sym_desc_vec
Definition: normal.cpp:160
relational.h
Interface to relations between expressions.
GiNaC::info_flags::positive
@ positive
Definition: flags.h:226
GiNaC::numeric::smod
ex smod(const numeric &xi) const override
Apply symmetric modular homomorphism to an expanded multivariate polynomial.
Definition: normal.cpp:1208
qbar
ex qbar
Definition: integration_kernel.cpp:247
options
unsigned options
Definition: factor.cpp:2480
GiNaC::ex::is_equal
bool is_equal(const ex &other) const
Definition: ex.h:345
GiNaC::normal_map_function::operator()
ex operator()(const ex &e) override
Definition: normal.cpp:2186
factors
upvec factors
Definition: factor.cpp:1461
m
mvec m
Definition: factor.cpp:771
GiNaC::matrix::solve
matrix solve(const matrix &vars, const matrix &rhs, unsigned algo=solve_algo::automatic) const
Solve a linear system consisting of a m x n matrix and a m x p right hand side by applying an elimina...
Definition: matrix.cpp:995
GiNaC::symbol::to_rational
ex to_rational(exmap &repl) const override
Implementation of ex::to_rational() for symbols.
Definition: normal.cpp:2622
GiNaC::ex::unit
ex unit(const ex &x) const
Compute unit part (= sign of leading coefficient) of a multivariate polynomial in Q[x].
Definition: normal.cpp:923
GiNaC::ex::primpart
ex primpart(const ex &x) const
Compute primitive part of a multivariate polynomial in Q[x].
Definition: normal.cpp:981
GiNaC
Definition: add.cpp:38
GiNaC::sym_desc
This structure holds information about the highest and lowest degrees in which a symbol appears in tw...
Definition: normal.cpp:122
GiNaC::power::to_rational
ex to_rational(exmap &repl) const override
Implementation of ex::to_rational() for powers.
Definition: normal.cpp:2674
GiNaC::power::normal
ex normal(exmap &repl, exmap &rev_lookup, lst &modifier) const override
Implementation of ex::normal() for powers.
Definition: normal.cpp:2408
GiNaC::sym_desc::operator<
bool operator<(const sym_desc &x) const
Comparison operator for sorting.
Definition: normal.cpp:150
GiNaC::gcd_options::no_heur_gcd
@ no_heur_gcd
Usually GiNaC tries heuristic GCD first, because typically it's much faster than anything else.
Definition: normal.h:47
GiNaC::ex::numer_denom
ex numer_denom() const
Get numerator and denominator of an expression.
Definition: normal.cpp:2564
GiNaC::basic::clearflag
const basic & clearflag(unsigned f) const
Clear some status_flags.
Definition: basic.h:291
matrix.h
Interface to symbolic matrices.
GiNaC::find_common_factor
static ex find_common_factor(const ex &e, ex &factor, exmap &repl)
Remove the common factor in the terms of a sum 'e' by calculating the GCD, and multiply it into the e...
Definition: normal.cpp:2740
GiNaC::add::normal
ex normal(exmap &repl, exmap &rev_lookup, lst &modifier) const override
Implementation of ex::normal() for a sum.
Definition: normal.cpp:2306
GiNaC::gcd_options::no_part_factored
@ no_part_factored
GiNaC tries to avoid expanding expressions when computing GCDs.
Definition: normal.h:55
x
ex x
Definition: factor.cpp:1641
GiNaC::gcd_pf_pow_pow
static ex gcd_pf_pow_pow(const ex &a, const ex &b, ex *ca, ex *cb)
Definition: normal.cpp:1646
GiNaC::ex::smod
ex smod(const numeric &xi) const
Definition: ex.h:202
GiNaC::ex::op
ex op(size_t i) const
Definition: ex.h:136
GiNaC::ex::info
bool info(unsigned inf) const
Definition: ex.h:132
utils.h
Interface to several small and furry utilities needed within GiNaC but not of any interest to the use...
GiNaC::expairseq::default_overall_coeff
virtual ex default_overall_coeff() const
GiNaC::numeric::normal
ex normal(exmap &repl, exmap &rev_lookup, lst &modifier) const override
Implementation of ex::normal() for a numeric.
Definition: normal.cpp:2222
GiNaC::collect_symbols
static void collect_symbols(const ex &e, sym_desc_vec &v)
Definition: normal.cpp:173
GiNaC::gcd
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
GiNaC::ex
Lightweight wrapper for GiNaC's symbolic objects.
Definition: ex.h:72
GiNaC::basic::smod
virtual ex smod(const numeric &xi) const
Apply symmetric modular homomorphism to an expanded multivariate polynomial.
Definition: normal.cpp:1203
lst.h
Definition of GiNaC's lst.
GiNaC::ex::begin
const_iterator begin() const noexcept
Definition: ex.h:647
GiNaC::numeric::is_zero
bool is_zero() const
True if object is zero.
Definition: numeric.cpp:1129
GiNaC::isqrt
const numeric isqrt(const numeric &x)
Integer numeric square root.
Definition: numeric.cpp:2482
GiNaC::pseries::point
ex point
Expansion point.
Definition: pseries.h:121
GiNaC::normal
ex normal(const ex &thisex)
Definition: ex.h:754
GiNaC::numeric::integer_content
numeric integer_content() const override
Definition: normal.cpp:328
GiNaC::ex::max_coefficient
numeric max_coefficient() const
Return maximum (absolute value) coefficient of a polynomial.
Definition: normal.cpp:1154
GiNaC::sqrfree_yun
static epvector sqrfree_yun(const ex &a, const symbol &x)
Compute square-free factorization of multivariate polynomial a(x) using Yun's algorithm.
Definition: normal.cpp:1798
GiNaC::numer
ex numer(const ex &thisex)
Definition: ex.h:745
GiNaC::ex::find
bool find(const ex &pattern, exset &found) const
Find all occurrences of a pattern.
Definition: ex.cpp:105
GiNaC::exp
const numeric exp(const numeric &x)
Exponential function.
Definition: numeric.cpp:1439
GiNaC::status_flags::hash_calculated
@ hash_calculated
.calchash() has already done its job
Definition: flags.h:205
GiNaC::ex::diff
ex diff(const symbol &s, unsigned nth=1) const
Compute partial derivative of an expression.
Definition: ex.cpp:86
GiNaC::rem
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
GiNaC::ex::ldegree
int ldegree(const ex &s) const
Definition: ex.h:174
GiNaC::ex::to_polynomial
ex to_polynomial(exmap &repl) const
Definition: normal.cpp:2602
GiNaC::numeric::inverse
const numeric inverse() const
Inverse of a number.
Definition: numeric.cpp:1053
GiNaC::expand
ex expand(const ex &thisex, unsigned options=0)
Definition: ex.h:715
GiNaC::info_flags::rational
@ rational
Definition: flags.h:222
GiNaC::add::coeff
ex coeff(const ex &s, int n=1) const override
Return coefficient of degree n in object s.
Definition: add.cpp:294
ex.h
Interface to GiNaC's light-weight expression handles.
GiNaC::expairseq::seq
epvector seq
Definition: expairseq.h:127
GiNaC::power::to_polynomial
ex to_polynomial(exmap &repl) const override
Implementation of ex::to_polynomial() for powers.
Definition: normal.cpp:2684
GiNaC::I
const numeric I
Imaginary unit.
Definition: numeric.cpp:1433
GiNaC::basic::coeff
virtual ex coeff(const ex &s, int n=1) const
Return coefficient of degree n in object s.
Definition: basic.cpp:337
symbol.h
Interface to GiNaC's symbolic objects.
GiNaC::container::op
ex op(size_t i) const override
Return operand/member at position i.
Definition: container.h:295
GiNaC::lcm_of_coefficients_denominators
static numeric lcm_of_coefficients_denominators(const ex &e)
Compute LCM of denominators of coefficients of a polynomial.
Definition: normal.cpp:261
GiNaC::container::remove_first
container & remove_first()
Remove first element.
Definition: container.h:400
GiNaC::interpolate
static ex interpolate(const ex &gamma, const numeric &xi, const ex &x, int degree_hint=1)
xi-adic polynomial interpolation
Definition: normal.cpp:1245
GiNaC::container::subs
ex subs(const exmap &m, unsigned options=0) const override
Substitute a set of objects by arbitrary expressions.
Definition: container.h:316
GiNaC::add::recombine_pair_to_ex
ex recombine_pair_to_ex(const expair &p) const override
Definition: add.cpp:564
GiNaC::collect_common_factors
ex collect_common_factors(const ex &e)
Collect common factors in sums.
Definition: normal.cpp:2832
GiNaC::ex::unitcontprim
void unitcontprim(const ex &x, ex &u, ex &c, ex &p) const
Compute unit part, content part, and primitive part of a multivariate polynomial in Q[x].
Definition: normal.cpp:1022
GiNaC::sqrfree_parfrac
ex sqrfree_parfrac(const ex &a, const symbol &x)
Compute square-free partial fraction decomposition of rational function a(x).
Definition: normal.cpp:1947
GiNaC::denom
ex denom(const ex &thisex)
Definition: ex.h:748
GiNaC::add::smod
ex smod(const numeric &xi) const override
Apply symmetric modular homomorphism to an expanded multivariate polynomial.
Definition: normal.cpp:1213
GiNaC::info_flags::numeric
@ numeric
Definition: flags.h:220
GiNaC::info_flags::polynomial
@ polynomial
Definition: flags.h:255
GiNaC::numeric::int_length
int int_length() const
Size in binary notation.
Definition: numeric.cpp:1418
GiNaC::basic::integer_content
virtual numeric integer_content() const
Definition: normal.cpp:323
GiNaC::pseries::seq
epvector seq
Vector of {coefficient, power} pairs.
Definition: pseries.h:115
GiNaC::expairseq::thisexpairseq
virtual ex thisexpairseq(const epvector &v, const ex &oc, bool do_index_renaming=false) const
GiNaC::container
Wrapper template for making GiNaC classes out of STL containers.
Definition: container.h:73
GiNaC::get_first_symbol
static bool get_first_symbol(const ex &e, ex &x)
Return pointer to first symbol found in expression.
Definition: normal.cpp:95
GiNaC::subs_options::no_pattern
@ no_pattern
disable pattern matching
Definition: flags.h:51
GiNaC::divide_in_z
static bool divide_in_z(const ex &a, const ex &b, ex &q, sym_desc_vec::const_iterator var)
Exact polynomial division of a(X) by b(X) in Z[X].
Definition: normal.cpp:757
GiNaC::basic::to_polynomial
virtual ex to_polynomial(exmap &repl) const
Definition: normal.cpp:2614
c
size_t c
Definition: factor.cpp:770
GiNaC::exmap
std::map< ex, ex, ex_is_less > exmap
Definition: basic.h:50
GiNaC::gcdheu_failed
Exception thrown by heur_gcd() to signal failure.
Definition: normal.cpp:1259
GiNaC::ex::lcoeff
ex lcoeff(const ex &s) const
Definition: ex.h:176
GiNaC::sym_desc::ldeg_b
int ldeg_b
Lowest degree of symbol in polynomial "b".
Definition: normal.cpp:141
GiNaC::mul::smod
ex smod(const numeric &xi) const override
Apply symmetric modular homomorphism to an expanded multivariate polynomial.
Definition: normal.cpp:1228
GiNaC::coeff
ex coeff(const ex &thisex, const ex &s, int n=1)
Definition: ex.h:742
GiNaC::factor
ex factor(const ex &poly, unsigned options)
Interface function to the outside world.
Definition: factor.cpp:2581
GiNaC::expairseq::overall_coeff
ex overall_coeff
Definition: expairseq.h:128
GiNaC::symbol::to_polynomial
ex to_polynomial(exmap &repl) const override
Implementation of ex::to_polynomial() for symbols.
Definition: normal.cpp:2629
GiNaC::pseries::var
ex var
Series variable (holds a symbol)
Definition: pseries.h:118
GiNaC::degree
int degree(const ex &thisex, const ex &s)
Definition: ex.h:736
GiNaC::numer_denom
ex numer_denom(const ex &thisex)
Definition: ex.h:751
basic.h
Interface to GiNaC's ABC.
GiNaC::symbol
Basic CAS symbol.
Definition: symbol.h:39
expairseq.h
Interface to sequences of expression pairs.
n
size_t n
Definition: factor.cpp:1463
GiNaC::normal_map_function
Function object to be applied by basic::normal().
Definition: normal.cpp:2185
GiNaC::pow
const numeric pow(const numeric &x, const numeric &y)
Definition: numeric.h:251
pseries.h
Interface to class for extended truncated power series.
GiNaC::sym_desc::deg_a
int deg_a
Highest degree of symbol in polynomial "a".
Definition: normal.cpp:132
GiNaC::matrix
Symbolic matrices.
Definition: matrix.h:38
GiNaC::add::integer_content
numeric integer_content() const override
Definition: normal.cpp:333
GiNaC::gcd_pf_pow
static ex gcd_pf_pow(const ex &a, const ex &b, ex *ca, ex *cb)
Definition: normal.cpp:1693
GiNaC::numeric::max_coefficient
numeric max_coefficient() const override
Implementation ex::max_coefficient().
Definition: normal.cpp:1166
GiNaC::numeric::to_polynomial
ex to_polynomial(exmap &repl) const override
Implementation of ex::to_polynomial() for a numeric.
Definition: normal.cpp:2656
GiNaC::is_negative
bool is_negative(const numeric &x)
Definition: numeric.h:269
GiNaC::expair
A pair of expressions.
Definition: expair.h:38
GiNaC::_num0_p
const numeric * _num0_p
Definition: utils.cpp:175
GiNaC::info_flags::negint
@ negint
Definition: flags.h:230
GiNaC::power::basis
ex basis
Definition: power.h:105
GiNaC::ex::is_zero
bool is_zero() const
Definition: ex.h:213
is_ex_the_function
#define is_ex_the_function(OBJ, FUNCNAME)
Definition: function.h:765
GiNaC::gcd_pf_mul
static ex gcd_pf_mul(const ex &a, const ex &b, ex *ca, ex *cb)
Definition: normal.cpp:1740
GiNaC::expairseq::to_polynomial
ex to_polynomial(exmap &repl) const override
Implementation of ex::to_polynomial() for expairseqs.
Definition: normal.cpp:2721
GiNaC::lcmcoeff
static numeric lcmcoeff(const ex &e, const numeric &l)
Definition: normal.cpp:231
GiNaC::mul::recombine_pair_to_ex
ex recombine_pair_to_ex(const expair &p) const override
Definition: mul.cpp:999
GiNaC::info_flags::posint
@ posint
Definition: flags.h:229
GiNaC::lcm
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
GiNaC::sym_desc::ldeg_a
int ldeg_a
Lowest degree of symbol in polynomial "a".
Definition: normal.cpp:138
GiNaC::add::expand
ex expand(unsigned options=0) const override
Expand expression, i.e.
Definition: add.cpp:572
GiNaC::ex::numer
ex numer() const
Get numerator of an expression.
Definition: normal.cpp:2514
GiNaC::decomp_rational
ex decomp_rational(const ex &a, const ex &x)
Decompose rational function a(x)=N(x)/D(x) into P(x)+n(x)/D(x) with degree(n, x) < degree(D,...
Definition: normal.cpp:472
GiNaC::sym_desc::max_deg
int max_deg
Maximum of deg_a and deg_b (Used for sorting)
Definition: normal.cpp:144
GiNaC::info_flags::integer_polynomial
@ integer_polynomial
Definition: flags.h:256
GiNaC::basic::max_coefficient
virtual numeric max_coefficient() const
Implementation ex::max_coefficient().
Definition: normal.cpp:1161
GiNaC::abs
const numeric abs(const numeric &x)
Absolute value.
Definition: numeric.cpp:2315
GiNaC::smod
const numeric smod(const numeric &a_, const numeric &b_)
Modulus (in symmetric representation).
Definition: numeric.cpp:2341
GiNaC::numeric::real
const numeric real() const
Real part of a number.
Definition: numeric.cpp:1339
GiNaC::is_integer
bool is_integer(const numeric &x)
Definition: numeric.h:272
GiNaC::container::nops
size_t nops() const override
Number of operands/members.
Definition: container.h:118
GiNaC::_num1_p
const numeric * _num1_p
Definition: utils.cpp:192
GiNaC::numeric::denom
const numeric denom() const
Denominator.
Definition: numeric.cpp:1387
GiNaC::basic::map
virtual ex map(map_function &f) const
Construct new expression by applying the specified function to all sub-expressions (one level only,...
Definition: basic.cpp:294
GiNaC::numeric::numer
const numeric numer() const
Numerator.
Definition: numeric.cpp:1356
GiNaC::replace_with_symbol
static ex replace_with_symbol(const ex &e, exmap &repl, exmap &rev_lookup, lst &modifier)
Create a symbol for replacing the expression "e" (or return a previously assigned symbol).
Definition: normal.cpp:2052
GiNaC::get_symbol_stats
static void get_symbol_stats(const ex &a, const ex &b, sym_desc_vec &v)
Collect statistical information about symbols in polynomials.
Definition: normal.cpp:197
GiNaC::ex::bp
ptr< basic > bp
pointer to basic object managed by this
Definition: ex.h:251
GiNaC::divide
bool divide(const ex &a, const ex &b, ex &q, bool check_args)
Exact polynomial division of a(X) by b(X) in Q[X].
Definition: normal.cpp:595
GiNaC::sprem
ex sprem(const ex &a, const ex &b, const ex &x, bool check_args)
Sparse pseudo-remainder of polynomials a(x) and b(x) in Q[x].
Definition: normal.cpp:544
operators.h
Interface to GiNaC's overloaded operators.
GiNaC::ex::integer_content
numeric integer_content() const
Compute the integer content (= GCD of all numeric coefficients) of an expanded polynomial.
Definition: normal.cpp:318
GiNaC::mul::integer_content
numeric integer_content() const override
Definition: normal.cpp:348
GiNaC::_ex_1
const ex _ex_1
Definition: utils.cpp:160
GiNaC::expairseq::to_rational
ex to_rational(exmap &repl) const override
Implementation of ex::to_rational() for expairseqs.
Definition: normal.cpp:2705
GiNaC::ex::content
ex content(const ex &x) const
Compute content part (= unit normal GCD of all coefficients) of a multivariate polynomial in Q[x].
Definition: normal.cpp:945
GiNaC::expairseq::recombine_pair_to_ex
virtual ex recombine_pair_to_ex(const expair &p) const
GiNaC::numeric
This class is a wrapper around CLN-numbers within the GiNaC class hierarchy.
Definition: numeric.h:82
GiNaC::gcd_options::use_sr_gcd
@ use_sr_gcd
By default GiNaC uses modular GCD algorithm.
Definition: normal.h:61
normal.h
This file defines several functions that work on univariate and multivariate polynomials and rational...
GINAC_ASSERT
#define GINAC_ASSERT(X)
Assertion macro for checking invariances.
Definition: assertion.h:33
GiNaC::basic::to_rational
virtual ex to_rational(exmap &repl) const
Default implementation of ex::to_rational().
Definition: normal.cpp:2609
GiNaC::sym_desc::sym
ex sym
Reference to symbol.
Definition: normal.cpp:129
GiNaC::mul
Product of expressions.
Definition: mul.h:32
GiNaC::prem
ex prem(const ex &a, const ex &b, const ex &x, bool check_args)
Pseudo-remainder of polynomials a(x) and b(x) in Q[x].
Definition: normal.cpp:492

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