GiNaC  1.8.3
normal.cpp
Go to the documentation of this file.
1 
8 /*
9  * GiNaC Copyright (C) 1999-2022 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  exvector factor, cofac;
1961  for (size_t i=0; i<yun.size(); i++) {
1962  numeric i_exponent = ex_to<numeric>(yun[i].coeff);
1963  for (size_t j=0; j<i_exponent; j++) {
1964  factor.push_back(pow(yun[i].rest, j+1));
1965  ex prod = _ex1;
1966  for (size_t k=0; k<yun.size(); k++) {
1967  if (yun[k].coeff == i_exponent)
1968  prod *= pow(yun[k].rest, i_exponent-1-j);
1969  else
1970  prod *= pow(yun[k].rest, yun[k].coeff);
1971  }
1972  cofac.push_back(prod.expand());
1973  }
1974  }
1975  size_t num_factors = factor.size();
1976 //std::clog << "factors : " << exprseq(factor) << std::endl;
1977 //std::clog << "cofactors: " << exprseq(cofac) << std::endl;
1978 
1979  // Construct coefficient matrix for decomposition
1980  int max_denom_deg = denom.degree(x);
1981  matrix sys(max_denom_deg + 1, num_factors);
1982  matrix rhs(max_denom_deg + 1, 1);
1983  for (int i=0; i<=max_denom_deg; i++) {
1984  for (size_t j=0; j<num_factors; j++)
1985  sys(i, j) = cofac[j].coeff(x, i);
1986  rhs(i, 0) = red_numer.coeff(x, i);
1987  }
1988 //std::clog << "coeffs: " << sys << std::endl;
1989 //std::clog << "rhs : " << rhs << std::endl;
1990 
1991  // Solve resulting linear system
1992  matrix vars(num_factors, 1);
1993  for (size_t i=0; i<num_factors; i++)
1994  vars(i, 0) = symbol();
1995  matrix sol = sys.solve(vars, rhs);
1996 
1997  // Sum up decomposed fractions
1998  ex sum = 0;
1999  for (size_t i=0; i<num_factors; i++)
2000  sum += sol(i, 0) / factor[i];
2001 
2002  return red_poly + sum;
2003 }
2004 
2005 
2006 /*
2007  * Normal form of rational functions
2008  */
2009 
2010 /*
2011  * Note: The internal normal() functions (= basic::normal() and overloaded
2012  * functions) all return lists of the form {numerator, denominator}. This
2013  * is to get around mul::eval()'s automatic expansion of numeric coefficients.
2014  * E.g. (a+b)/3 is automatically converted to a/3+b/3 but we want to keep
2015  * the information that (a+b) is the numerator and 3 is the denominator.
2016  */
2017 
2018 
2051 static ex replace_with_symbol(const ex & e, exmap & repl, exmap & rev_lookup, lst & modifier)
2052 {
2053  // Since the repl contains replaced expressions we should search for them
2054  ex e_replaced = e.subs(repl, subs_options::no_pattern);
2055 
2056  // Expression already replaced? Then return the assigned symbol
2057  auto it = rev_lookup.find(e_replaced);
2058  if (it != rev_lookup.end())
2059  return it->second;
2060 
2061  // The expression can be the base of substituted power, which requires a more careful search
2062  if (! is_a<numeric>(e_replaced))
2063  for (auto & it : repl)
2064  if (is_a<power>(it.second) && e_replaced.is_equal(it.second.op(0))) {
2065  ex degree = pow(it.second.op(1), _ex_1);
2066  if (is_a<numeric>(degree) && ex_to<numeric>(degree).is_integer())
2067  return pow(it.first, degree);
2068  }
2069 
2070  // We treat powers and the exponent functions differently because
2071  // they can be rationalised more efficiently
2072  if (is_a<function>(e_replaced) && is_ex_the_function(e_replaced, exp)) {
2073  for (auto & it : repl) {
2074  if (is_a<function>(it.second) && is_ex_the_function(it.second, exp)) {
2075  ex ratio = normal(e_replaced.op(0) / it.second.op(0));
2076  if (is_a<numeric>(ratio) && ex_to<numeric>(ratio).is_rational()) {
2077  // Different exponents can be treated as powers of the same basic equation
2078  if (ex_to<numeric>(ratio).is_integer()) {
2079  // If ratio is an integer then this is simply the power of the existing symbol.
2080  // std::clog << e_replaced << " is a " << ratio << " power of " << it.first << std::endl;
2081  return dynallocate<power>(it.first, ratio);
2082  } else {
2083  // otherwise we need to give the replacement pattern to change
2084  // the previous expression...
2085  ex es = dynallocate<symbol>();
2086  ex Num = numer(ratio);
2087  modifier.append(it.first == power(es, denom(ratio)));
2088  // std::clog << e_replaced << " is power " << Num << " and "
2089  // << it.first << " is power " << denom(ratio) << " of the common base "
2090  // << exp(e_replaced.op(0)/Num) << std::endl;
2091  // ... and modify the replacement tables
2092  rev_lookup.erase(it.second);
2093  rev_lookup.insert({exp(e_replaced.op(0)/Num), es});
2094  repl.erase(it.first);
2095  repl.insert({es, exp(e_replaced.op(0)/Num)});
2096  return dynallocate<power>(es, Num);
2097  }
2098  }
2099  }
2100  }
2101  } 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)
2102  && ! (is_a<symbol>(e_replaced.op(0))
2103  && is_a<numeric>(e_replaced.op(1)) && ex_to<numeric>(e_replaced.op(1)).is_integer())) {
2104  for (auto & it : repl) {
2105  if (e_replaced.op(0).is_equal(it.second) // The base is an allocated symbol or base of power
2106  || (is_a<power>(it.second) && e_replaced.op(0).is_equal(it.second.op(0)))) {
2107  ex ratio; // We bind together two above cases
2108  if (is_a<power>(it.second))
2109  ratio = normal(e_replaced.op(1) / it.second.op(1));
2110  else
2111  ratio = e_replaced.op(1);
2112  if (is_a<numeric>(ratio) && ex_to<numeric>(ratio).is_rational()) {
2113  // Different powers can be treated as powers of the same basic equation
2114  if (ex_to<numeric>(ratio).is_integer()) {
2115  // If ratio is an integer then this is simply the power of the existing symbol.
2116  //std::clog << e_replaced << " is a " << ratio << " power of " << it.first << std::endl;
2117  return dynallocate<power>(it.first, ratio);
2118  } else {
2119  // otherwise we need to give the replacement pattern to change
2120  // the previous expression...
2121  ex es = dynallocate<symbol>();
2122  ex Num = numer(ratio);
2123  modifier.append(it.first == power(es, denom(ratio)));
2124  //std::clog << e_replaced << " is power " << Num << " and "
2125  // << it.first << " is power " << denom(ratio) << " of the common base "
2126  // << pow(e_replaced.op(0), e_replaced.op(1)/Num) << std::endl;
2127  // ... and modify the replacement tables
2128  rev_lookup.erase(it.second);
2129  rev_lookup.insert({pow(e_replaced.op(0), e_replaced.op(1)/Num), es});
2130  repl.erase(it.first);
2131  repl.insert({es, pow(e_replaced.op(0), e_replaced.op(1)/Num)});
2132  return dynallocate<power>(es, Num);
2133  }
2134  }
2135  }
2136  }
2137  // There is no existing substitution, thus we are creating a new one.
2138  // This needs to be done separately to treat possible occurrences of
2139  // b = e_replaced.op(0) elsewhere in the expression as pow(b, 1).
2140  ex degree = pow(e_replaced.op(1), _ex_1);
2141  if (is_a<numeric>(degree) && ex_to<numeric>(degree).is_integer()) {
2142  ex es = dynallocate<symbol>();
2143  modifier.append(e_replaced.op(0) == power(es, degree));
2144  repl.insert({es, e_replaced});
2145  rev_lookup.insert({e_replaced, es});
2146  return es;
2147  }
2148  }
2149 
2150  // Otherwise create new symbol and add to list, taking care that the
2151  // replacement expression doesn't itself contain symbols from repl,
2152  // because subs() is not recursive
2153  ex es = dynallocate<symbol>();
2154  repl.insert(std::make_pair(es, e_replaced));
2155  rev_lookup.insert(std::make_pair(e_replaced, es));
2156  return es;
2157 }
2158 
2164 static ex replace_with_symbol(const ex & e, exmap & repl)
2165 {
2166  // Since the repl contains replaced expressions we should search for them
2167  ex e_replaced = e.subs(repl, subs_options::no_pattern);
2168 
2169  // Expression already replaced? Then return the assigned symbol
2170  for (auto & it : repl)
2171  if (it.second.is_equal(e_replaced))
2172  return it.first;
2173 
2174  // Otherwise create new symbol and add to list, taking care that the
2175  // replacement expression doesn't itself contain symbols from repl,
2176  // because subs() is not recursive
2177  ex es = dynallocate<symbol>();
2178  repl.insert(std::make_pair(es, e_replaced));
2179  return es;
2180 }
2181 
2182 
2185  ex operator()(const ex & e) override { return normal(e); }
2186 };
2187 
2191 ex basic::normal(exmap & repl, exmap & rev_lookup, lst & modifier) const
2192 {
2193  if (nops() == 0)
2194  return dynallocate<lst>({replace_with_symbol(*this, repl, rev_lookup, modifier), _ex1});
2195 
2196  normal_map_function map_normal;
2197  size_t nmod = modifier.nops(); // To watch new modifiers to the replacement list
2198  ex result = replace_with_symbol(map(map_normal), repl, rev_lookup, modifier);
2199  for (size_t imod = nmod; imod < modifier.nops(); ++imod) {
2200  exmap this_repl;
2201  this_repl.insert(std::make_pair(modifier.op(imod).op(0), modifier.op(imod).op(1)));
2202  result = result.subs(this_repl, subs_options::no_pattern);
2203  }
2204 
2205  // Sometimes we may obtain negative powers, they need to be placed to denominator
2206  if (is_a<power>(result) && result.op(1).info(info_flags::negative))
2207  return dynallocate<lst>({_ex1, power(result.op(0), -result.op(1))});
2208  else
2209  return dynallocate<lst>({result, _ex1});
2210 }
2211 
2212 
2215 ex symbol::normal(exmap & repl, exmap & rev_lookup, lst & modifier) const
2216 {
2217  return dynallocate<lst>({*this, _ex1});
2218 }
2219 
2220 
2225 ex numeric::normal(exmap & repl, exmap & rev_lookup, lst & modifier) const
2226 {
2227  numeric num = numer();
2228  ex numex = num;
2229 
2230  if (num.is_real()) {
2231  if (!num.is_integer())
2232  numex = replace_with_symbol(numex, repl, rev_lookup, modifier);
2233  } else { // complex
2234  numeric re = num.real(), im = num.imag();
2235  ex re_ex = re.is_rational() ? re : replace_with_symbol(re, repl, rev_lookup, modifier);
2236  ex im_ex = im.is_rational() ? im : replace_with_symbol(im, repl, rev_lookup, modifier);
2237  numex = re_ex + im_ex * replace_with_symbol(I, repl, rev_lookup, modifier);
2238  }
2239 
2240  // Denominator is always a real integer (see numeric::denom())
2241  return dynallocate<lst>({numex, denom()});
2242 }
2243 
2244 
2249 static ex frac_cancel(const ex &n, const ex &d)
2250 {
2251  ex num = n;
2252  ex den = d;
2253  numeric pre_factor = *_num1_p;
2254 
2255 //std::clog << "frac_cancel num = " << num << ", den = " << den << std::endl;
2256 
2257  // Handle trivial case where denominator is 1
2258  if (den.is_equal(_ex1))
2259  return dynallocate<lst>({num, den});
2260 
2261  // Handle special cases where numerator or denominator is 0
2262  if (num.is_zero())
2263  return dynallocate<lst>({num, _ex1});
2264  if (den.expand().is_zero())
2265  throw(std::overflow_error("frac_cancel: division by zero in frac_cancel"));
2266 
2267  // Bring numerator and denominator to Z[X] by multiplying with
2268  // LCM of all coefficients' denominators
2271  num = multiply_lcm(num, num_lcm);
2272  den = multiply_lcm(den, den_lcm);
2273  pre_factor = den_lcm / num_lcm;
2274 
2275  // Cancel GCD from numerator and denominator
2276  ex cnum, cden;
2277  if (gcd(num, den, &cnum, &cden, false) != _ex1) {
2278  num = cnum;
2279  den = cden;
2280  }
2281 
2282  // Make denominator unit normal (i.e. coefficient of first symbol
2283  // as defined by get_first_symbol() is made positive)
2284  if (is_exactly_a<numeric>(den)) {
2285  if (ex_to<numeric>(den).is_negative()) {
2286  num *= _ex_1;
2287  den *= _ex_1;
2288  }
2289  } else {
2290  ex x;
2291  if (get_first_symbol(den, x)) {
2292  GINAC_ASSERT(is_exactly_a<numeric>(den.unit(x)));
2293  if (ex_to<numeric>(den.unit(x)).is_negative()) {
2294  num *= _ex_1;
2295  den *= _ex_1;
2296  }
2297  }
2298  }
2299 
2300  // Return result as list
2301 //std::clog << " returns num = " << num << ", den = " << den << ", pre_factor = " << pre_factor << std::endl;
2302  return dynallocate<lst>({num * pre_factor.numer(), den * pre_factor.denom()});
2303 }
2304 
2305 
2309 ex add::normal(exmap & repl, exmap & rev_lookup, lst & modifier) const
2310 {
2311  // Normalize children and split each one into numerator and denominator
2312  exvector nums, dens;
2313  nums.reserve(seq.size()+1);
2314  dens.reserve(seq.size()+1);
2315  size_t nmod = modifier.nops(); // To watch new modifiers to the replacement list
2316  for (auto & it : seq) {
2317  ex n = ex_to<basic>(recombine_pair_to_ex(it)).normal(repl, rev_lookup, modifier);
2318  nums.push_back(n.op(0));
2319  dens.push_back(n.op(1));
2320  }
2321  ex n = ex_to<numeric>(overall_coeff).normal(repl, rev_lookup, modifier);
2322  nums.push_back(n.op(0));
2323  dens.push_back(n.op(1));
2324  GINAC_ASSERT(nums.size() == dens.size());
2325 
2326  // Now, nums is a vector of all numerators and dens is a vector of
2327  // all denominators
2328 //std::clog << "add::normal uses " << nums.size() << " summands:\n";
2329 
2330  // Add fractions sequentially
2331  auto num_it = nums.begin(), num_itend = nums.end();
2332  auto den_it = dens.begin(), den_itend = dens.end();
2333 //std::clog << " num = " << *num_it << ", den = " << *den_it << std::endl;
2334  for (size_t imod = nmod; imod < modifier.nops(); ++imod) {
2335  while (num_it != num_itend) {
2336  *num_it = num_it->subs(modifier.op(imod), subs_options::no_pattern);
2337  ++num_it;
2338  *den_it = den_it->subs(modifier.op(imod), subs_options::no_pattern);
2339  ++den_it;
2340  }
2341  // Reset iterators for the next round
2342  num_it = nums.begin();
2343  den_it = dens.begin();
2344  }
2345 
2346  ex num = *num_it++, den = *den_it++;
2347  while (num_it != num_itend) {
2348 //std::clog << " num = " << *num_it << ", den = " << *den_it << std::endl;
2349  ex next_num = *num_it++, next_den = *den_it++;
2350 
2351  // Trivially add sequences of fractions with identical denominators
2352  while ((den_it != den_itend) && next_den.is_equal(*den_it)) {
2353  next_num += *num_it;
2354  num_it++; den_it++;
2355  }
2356 
2357  // Addition of two fractions, taking advantage of the fact that
2358  // the heuristic GCD algorithm computes the cofactors at no extra cost
2359  ex co_den1, co_den2;
2360  ex g = gcd(den, next_den, &co_den1, &co_den2, false);
2361  num = ((num * co_den2) + (next_num * co_den1)).expand();
2362  den *= co_den2; // this is the lcm(den, next_den)
2363  }
2364 //std::clog << " common denominator = " << den << std::endl;
2365 
2366  // Cancel common factors from num/den
2367  return frac_cancel(num, den);
2368 }
2369 
2370 
2374 ex mul::normal(exmap & repl, exmap & rev_lookup, lst & modifier) const
2375 {
2376  // Normalize children, separate into numerator and denominator
2377  exvector num; num.reserve(seq.size());
2378  exvector den; den.reserve(seq.size());
2379  ex n;
2380  size_t nmod = modifier.nops(); // To watch new modifiers to the replacement list
2381  for (auto & it : seq) {
2382  n = ex_to<basic>(recombine_pair_to_ex(it)).normal(repl, rev_lookup, modifier);
2383  num.push_back(n.op(0));
2384  den.push_back(n.op(1));
2385  }
2386  n = ex_to<numeric>(overall_coeff).normal(repl, rev_lookup, modifier);
2387  num.push_back(n.op(0));
2388  den.push_back(n.op(1));
2389  auto num_it = num.begin(), num_itend = num.end();
2390  auto den_it = den.begin();
2391  for (size_t imod = nmod; imod < modifier.nops(); ++imod) {
2392  while (num_it != num_itend) {
2393  *num_it = num_it->subs(modifier.op(imod), subs_options::no_pattern);
2394  ++num_it;
2395  *den_it = den_it->subs(modifier.op(imod), subs_options::no_pattern);
2396  ++den_it;
2397  }
2398  num_it = num.begin();
2399  den_it = den.begin();
2400  }
2401 
2402  // Perform fraction cancellation
2403  return frac_cancel(dynallocate<mul>(num), dynallocate<mul>(den));
2404 }
2405 
2406 
2411 ex power::normal(exmap & repl, exmap & rev_lookup, lst & modifier) const
2412 {
2413  // Normalize basis and exponent (exponent gets reassembled)
2414  size_t nmod = modifier.nops(); // To watch new modifiers to the replacement list
2415  ex n_basis = ex_to<basic>(basis).normal(repl, rev_lookup, modifier);
2416  for (size_t imod = nmod; imod < modifier.nops(); ++imod)
2417  n_basis = n_basis.subs(modifier.op(imod), subs_options::no_pattern);
2418 
2419  nmod = modifier.nops();
2420  ex n_exponent = ex_to<basic>(exponent).normal(repl, rev_lookup, modifier);
2421  for (size_t imod = nmod; imod < modifier.nops(); ++imod)
2422  n_exponent = n_exponent.subs(modifier.op(imod), subs_options::no_pattern);
2423  n_exponent = n_exponent.op(0) / n_exponent.op(1);
2424 
2425  if (n_exponent.info(info_flags::integer)) {
2426 
2427  if (n_exponent.info(info_flags::positive)) {
2428 
2429  // (a/b)^n -> {a^n, b^n}
2430  return dynallocate<lst>({pow(n_basis.op(0), n_exponent), pow(n_basis.op(1), n_exponent)});
2431 
2432  } else if (n_exponent.info(info_flags::negative)) {
2433 
2434  // (a/b)^-n -> {b^n, a^n}
2435  return dynallocate<lst>({pow(n_basis.op(1), -n_exponent), pow(n_basis.op(0), -n_exponent)});
2436  }
2437 
2438  } else {
2439 
2440  if (n_exponent.info(info_flags::positive)) {
2441 
2442  // (a/b)^x -> {sym((a/b)^x), 1}
2443  return dynallocate<lst>({replace_with_symbol(pow(n_basis.op(0) / n_basis.op(1), n_exponent), repl, rev_lookup, modifier), _ex1});
2444 
2445  } else if (n_exponent.info(info_flags::negative)) {
2446 
2447  if (n_basis.op(1).is_equal(_ex1)) {
2448 
2449  // a^-x -> {1, sym(a^x)}
2450  return dynallocate<lst>({_ex1, replace_with_symbol(pow(n_basis.op(0), -n_exponent), repl, rev_lookup, modifier)});
2451 
2452  } else {
2453 
2454  // (a/b)^-x -> {sym((b/a)^x), 1}
2455  return dynallocate<lst>({replace_with_symbol(pow(n_basis.op(1) / n_basis.op(0), -n_exponent), repl, rev_lookup, modifier), _ex1});
2456  }
2457  }
2458  }
2459 
2460  // (a/b)^x -> {sym((a/b)^x, 1}
2461  return dynallocate<lst>({replace_with_symbol(pow(n_basis.op(0) / n_basis.op(1), n_exponent), repl, rev_lookup, modifier), _ex1});
2462 }
2463 
2464 
2468 ex pseries::normal(exmap & repl, exmap & rev_lookup, lst & modifier) const
2469 {
2470  epvector newseq;
2471  for (auto & it : seq) {
2472  ex restexp = it.rest.normal();
2473  if (!restexp.is_zero())
2474  newseq.push_back(expair(restexp, it.coeff));
2475  }
2476  ex n = pseries(relational(var,point), std::move(newseq));
2477  return dynallocate<lst>({replace_with_symbol(n, repl, rev_lookup, modifier), _ex1});
2478 }
2479 
2480 
2493 {
2494  exmap repl, rev_lookup;
2495  lst modifier;
2496 
2497  ex e = bp->normal(repl, rev_lookup, modifier);
2498  GINAC_ASSERT(is_a<lst>(e));
2499 
2500  // Re-insert replaced symbols
2501  if (!repl.empty()) {
2502  for(size_t i=0; i < modifier.nops(); ++i)
2503  e = e.subs(modifier.op(i), subs_options::no_pattern);
2504  e = e.subs(repl, subs_options::no_pattern);
2505  }
2506 
2507  // Convert {numerator, denominator} form back to fraction
2508  return e.op(0) / e.op(1);
2509 }
2510 
2517 ex ex::numer() const
2518 {
2519  exmap repl, rev_lookup;
2520  lst modifier;
2521 
2522  ex e = bp->normal(repl, rev_lookup, modifier);
2523  GINAC_ASSERT(is_a<lst>(e));
2524 
2525  // Re-insert replaced symbols
2526  if (repl.empty())
2527  return e.op(0);
2528  else {
2529  for(size_t i=0; i < modifier.nops(); ++i)
2530  e = e.subs(modifier.op(i), subs_options::no_pattern);
2531 
2532  return e.op(0).subs(repl, subs_options::no_pattern);
2533  }
2534 }
2535 
2542 ex ex::denom() const
2543 {
2544  exmap repl, rev_lookup;
2545  lst modifier;
2546 
2547  ex e = bp->normal(repl, rev_lookup, modifier);
2548  GINAC_ASSERT(is_a<lst>(e));
2549 
2550  // Re-insert replaced symbols
2551  if (repl.empty())
2552  return e.op(1);
2553  else {
2554  for(size_t i=0; i < modifier.nops(); ++i)
2555  e = e.subs(modifier.op(i), subs_options::no_pattern);
2556 
2557  return e.op(1).subs(repl, subs_options::no_pattern);
2558  }
2559 }
2560 
2568 {
2569  exmap repl, rev_lookup;
2570  lst modifier;
2571 
2572  ex e = bp->normal(repl, rev_lookup, modifier);
2573  GINAC_ASSERT(is_a<lst>(e));
2574 
2575  // Re-insert replaced symbols
2576  if (repl.empty())
2577  return e;
2578  else {
2579  for(size_t i=0; i < modifier.nops(); ++i)
2580  e = e.subs(modifier.op(i), subs_options::no_pattern);
2581 
2582  return e.subs(repl, subs_options::no_pattern);
2583  }
2584 }
2585 
2586 
2600 ex ex::to_rational(exmap & repl) const
2601 {
2602  return bp->to_rational(repl);
2603 }
2604 
2606 {
2607  return bp->to_polynomial(repl);
2608 }
2609 
2613 {
2614  return replace_with_symbol(*this, repl);
2615 }
2616 
2618 {
2619  return replace_with_symbol(*this, repl);
2620 }
2621 
2622 
2626 {
2627  return *this;
2628 }
2629 
2633 {
2634  return *this;
2635 }
2636 
2637 
2642 {
2643  if (is_real()) {
2644  if (!is_rational())
2645  return replace_with_symbol(*this, repl);
2646  } else { // complex
2647  numeric re = real();
2648  numeric im = imag();
2649  ex re_ex = re.is_rational() ? re : replace_with_symbol(re, repl);
2650  ex im_ex = im.is_rational() ? im : replace_with_symbol(im, repl);
2651  return re_ex + im_ex * replace_with_symbol(I, repl);
2652  }
2653  return *this;
2654 }
2655 
2660 {
2661  if (is_real()) {
2662  if (!is_integer())
2663  return replace_with_symbol(*this, repl);
2664  } else { // complex
2665  numeric re = real();
2666  numeric im = imag();
2667  ex re_ex = re.is_integer() ? re : replace_with_symbol(re, repl);
2668  ex im_ex = im.is_integer() ? im : replace_with_symbol(im, repl);
2669  return re_ex + im_ex * replace_with_symbol(I, repl);
2670  }
2671  return *this;
2672 }
2673 
2674 
2678 {
2680  return pow(basis.to_rational(repl), exponent);
2681  else
2682  return replace_with_symbol(*this, repl);
2683 }
2684 
2688 {
2690  return pow(basis.to_polynomial(repl), exponent);
2691  else if (exponent.info(info_flags::negint))
2692  {
2693  ex basis_pref = collect_common_factors(basis);
2694  if (is_exactly_a<mul>(basis_pref) || is_exactly_a<power>(basis_pref)) {
2695  // (A*B)^n will be automagically transformed to A^n*B^n
2696  ex t = pow(basis_pref, exponent);
2697  return t.to_polynomial(repl);
2698  }
2699  else
2700  return pow(replace_with_symbol(pow(basis, _ex_1), repl), -exponent);
2701  }
2702  else
2703  return replace_with_symbol(*this, repl);
2704 }
2705 
2706 
2709 {
2710  epvector s;
2711  s.reserve(seq.size());
2712  for (auto & it : seq)
2713  s.push_back(split_ex_to_pair(recombine_pair_to_ex(it).to_rational(repl)));
2714 
2715  ex oc = overall_coeff.to_rational(repl);
2716  if (oc.info(info_flags::numeric))
2717  return thisexpairseq(std::move(s), overall_coeff);
2718  else
2719  s.push_back(expair(oc, _ex1));
2720  return thisexpairseq(std::move(s), default_overall_coeff());
2721 }
2722 
2725 {
2726  epvector s;
2727  s.reserve(seq.size());
2728  for (auto & it : seq)
2729  s.push_back(split_ex_to_pair(recombine_pair_to_ex(it).to_polynomial(repl)));
2730 
2731  ex oc = overall_coeff.to_polynomial(repl);
2732  if (oc.info(info_flags::numeric))
2733  return thisexpairseq(std::move(s), overall_coeff);
2734  else
2735  s.push_back(expair(oc, _ex1));
2736  return thisexpairseq(std::move(s), default_overall_coeff());
2737 }
2738 
2739 
2743 static ex find_common_factor(const ex & e, ex & factor, exmap & repl)
2744 {
2745  if (is_exactly_a<add>(e)) {
2746 
2747  size_t num = e.nops();
2748  exvector terms; terms.reserve(num);
2749  ex gc;
2750 
2751  // Find the common GCD
2752  for (size_t i=0; i<num; i++) {
2753  ex x = e.op(i).to_polynomial(repl);
2754 
2755  if (is_exactly_a<add>(x) || is_exactly_a<mul>(x) || is_a<power>(x)) {
2756  ex f = 1;
2757  x = find_common_factor(x, f, repl);
2758  x *= f;
2759  }
2760 
2761  if (gc.is_zero())
2762  gc = x;
2763  else
2764  gc = gcd(gc, x);
2765 
2766  terms.push_back(x);
2767  }
2768 
2769  if (gc.is_equal(_ex1))
2770  return e;
2771 
2772  if (gc.is_zero())
2773  return _ex0;
2774 
2775  // The GCD is the factor we pull out
2776  factor *= gc;
2777 
2778  // Now divide all terms by the GCD
2779  for (size_t i=0; i<num; i++) {
2780  ex x;
2781 
2782  // Try to avoid divide() because it expands the polynomial
2783  ex &t = terms[i];
2784  if (is_exactly_a<mul>(t)) {
2785  for (size_t j=0; j<t.nops(); j++) {
2786  if (t.op(j).is_equal(gc)) {
2787  exvector v; v.reserve(t.nops());
2788  for (size_t k=0; k<t.nops(); k++) {
2789  if (k == j)
2790  v.push_back(_ex1);
2791  else
2792  v.push_back(t.op(k));
2793  }
2794  t = dynallocate<mul>(v);
2795  goto term_done;
2796  }
2797  }
2798  }
2799 
2800  divide(t, gc, x);
2801  t = x;
2802 term_done: ;
2803  }
2804  return dynallocate<add>(terms);
2805 
2806  } else if (is_exactly_a<mul>(e)) {
2807 
2808  size_t num = e.nops();
2809  exvector v; v.reserve(num);
2810 
2811  for (size_t i=0; i<num; i++)
2812  v.push_back(find_common_factor(e.op(i), factor, repl));
2813 
2814  return dynallocate<mul>(v);
2815 
2816  } else if (is_exactly_a<power>(e)) {
2817  const ex e_exp(e.op(1));
2818  if (e_exp.info(info_flags::integer)) {
2819  ex eb = e.op(0).to_polynomial(repl);
2820  ex factor_local(_ex1);
2821  ex pre_res = find_common_factor(eb, factor_local, repl);
2822  factor *= pow(factor_local, e_exp);
2823  return pow(pre_res, e_exp);
2824 
2825  } else
2826  return e.to_polynomial(repl);
2827 
2828  } else
2829  return e;
2830 }
2831 
2832 
2836 {
2837  if (is_exactly_a<add>(e) || is_exactly_a<mul>(e) || is_exactly_a<power>(e)) {
2838 
2839  exmap repl;
2840  ex factor = 1;
2841  ex r = find_common_factor(e, factor, repl);
2842  return factor.subs(repl, subs_options::no_pattern) * r.subs(repl, subs_options::no_pattern);
2843 
2844  } else
2845  return e;
2846 }
2847 
2848 
2851 ex resultant(const ex & e1, const ex & e2, const ex & s)
2852 {
2853  const ex ee1 = e1.expand();
2854  const ex ee2 = e2.expand();
2855  if (!ee1.info(info_flags::polynomial) ||
2857  throw(std::runtime_error("resultant(): arguments must be polynomials"));
2858 
2859  const int h1 = ee1.degree(s);
2860  const int l1 = ee1.ldegree(s);
2861  const int h2 = ee2.degree(s);
2862  const int l2 = ee2.ldegree(s);
2863 
2864  const int msize = h1 + h2;
2865  matrix m(msize, msize);
2866 
2867  for (int l = h1; l >= l1; --l) {
2868  const ex e = ee1.coeff(s, l);
2869  for (int k = 0; k < h2; ++k)
2870  m(k, k+h1-l) = e;
2871  }
2872  for (int l = h2; l >= l2; --l) {
2873  const ex e = ee2.coeff(s, l);
2874  for (int k = 0; k < h1; ++k)
2875  m(k+h2, k+h2-l) = e;
2876  }
2877 
2878  return m.determinant();
2879 }
2880 
2881 
2882 } // namespace GiNaC
Interface to GiNaC's sums of expressions.
#define GINAC_ASSERT(X)
Assertion macro for checking invariances.
Definition: assertion.h:33
Interface to GiNaC's ABC.
ex coeff(const ex &s, int n=1) const override
Return coefficient of degree n in object s.
Definition: add.cpp:294
numeric integer_content() const override
Definition: normal.cpp:333
numeric max_coefficient() const override
Implementation ex::max_coefficient().
Definition: normal.cpp:1171
ex expand(unsigned options=0) const override
Expand expression, i.e.
Definition: add.cpp:572
ex normal(exmap &repl, exmap &rev_lookup, lst &modifier) const override
Implementation of ex::normal() for a sum.
Definition: normal.cpp:2309
ex recombine_pair_to_ex(const expair &p) const override
Definition: add.cpp:564
ex smod(const numeric &xi) const override
Apply symmetric modular homomorphism to an expanded multivariate polynomial.
Definition: normal.cpp:1213
virtual size_t nops() const
Number of operands/members.
Definition: basic.cpp:229
virtual numeric integer_content() const
Definition: normal.cpp:323
virtual numeric max_coefficient() const
Implementation ex::max_coefficient().
Definition: normal.cpp:1161
virtual ex to_polynomial(exmap &repl) const
Definition: normal.cpp:2617
const basic & clearflag(unsigned f) const
Clear some status_flags.
Definition: basic.h:291
virtual ex to_rational(exmap &repl) const
Default implementation of ex::to_rational().
Definition: normal.cpp:2612
virtual ex coeff(const ex &s, int n=1) const
Return coefficient of degree n in object s.
Definition: basic.cpp:337
virtual ex smod(const numeric &xi) const
Apply symmetric modular homomorphism to an expanded multivariate polynomial.
Definition: normal.cpp:1203
virtual ex normal(exmap &repl, exmap &rev_lookup, lst &modifier) const
Default implementation of ex::normal().
Definition: normal.cpp:2191
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
Wrapper template for making GiNaC classes out of STL containers.
Definition: container.h:73
size_t nops() const override
Number of operands/members.
Definition: container.h:118
ex op(size_t i) const override
Return operand/member at position i.
Definition: container.h:295
container & remove_first()
Remove first element.
Definition: container.h:400
container & append(const ex &b)
Add element at back.
Definition: container.h:391
Lightweight wrapper for GiNaC's symbolic objects.
Definition: ex.h:72
ex to_rational(exmap &repl) const
Rationalization of non-rational functions.
Definition: normal.cpp:2600
ex unit(const ex &x) const
Compute unit part (= sign of leading coefficient) of a multivariate polynomial in Q[x].
Definition: normal.cpp:923
numeric integer_content() const
Compute the integer content (= GCD of all numeric coefficients) of an expanded polynomial.
Definition: normal.cpp:318
ex primpart(const ex &x) const
Compute primitive part of a multivariate polynomial in Q[x].
Definition: normal.cpp:981
ex lcoeff(const ex &s) const
Definition: ex.h:176
const_iterator begin() const noexcept
Definition: ex.h:647
bool find(const ex &pattern, exset &found) const
Find all occurrences of a pattern.
Definition: ex.cpp:105
ex diff(const symbol &s, unsigned nth=1) const
Compute partial derivative of an expression.
Definition: ex.cpp:86
ex expand(unsigned options=0) const
Definition: ex.cpp:73
ex numer_denom() const
Get numerator and denominator of an expression.
Definition: normal.cpp:2567
bool is_equal(const ex &other) const
Definition: ex.h:345
ex normal() const
Normalization of rational functions.
Definition: normal.cpp:2492
int degree(const ex &s) const
Definition: ex.h:173
ptr< basic > bp
pointer to basic object managed by this
Definition: ex.h:251
numeric max_coefficient() const
Return maximum (absolute value) coefficient of a polynomial.
Definition: normal.cpp:1154
ex to_polynomial(exmap &repl) const
Definition: normal.cpp:2605
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
size_t nops() const
Definition: ex.h:135
ex smod(const numeric &xi) const
Definition: ex.h:202
ex subs(const exmap &m, unsigned options=0) const
Definition: ex.h:826
bool info(unsigned inf) const
Definition: ex.h:132
ex denom() const
Get denominator of an expression.
Definition: normal.cpp:2542
bool is_zero() const
Definition: ex.h:213
ex op(size_t i) const
Definition: ex.h:136
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
int ldegree(const ex &s) const
Definition: ex.h:174
ex numer() const
Get numerator of an expression.
Definition: normal.cpp:2517
ex coeff(const ex &s, int n=1) const
Definition: ex.h:175
A pair of expressions.
Definition: expair.h:38
virtual expair split_ex_to_pair(const ex &e) const
epvector seq
Definition: expairseq.h:127
ex to_polynomial(exmap &repl) const override
Implementation of ex::to_polynomial() for expairseqs.
Definition: normal.cpp:2724
virtual ex thisexpairseq(const epvector &v, const ex &oc, bool do_index_renaming=false) const
virtual ex default_overall_coeff() const
virtual ex recombine_pair_to_ex(const expair &p) const
ex to_rational(exmap &repl) const override
Implementation of ex::to_rational() for expairseqs.
Definition: normal.cpp:2708
Exception thrown by heur_gcd() to signal failure.
Definition: normal.cpp:1259
@ integer_polynomial
Definition: flags.h:256
@ rational_polynomial
Definition: flags.h:258
Symbolic matrices.
Definition: matrix.h:38
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
Product of expressions.
Definition: mul.h:32
ex recombine_pair_to_ex(const expair &p) const override
Definition: mul.cpp:999
ex smod(const numeric &xi) const override
Apply symmetric modular homomorphism to an expanded multivariate polynomial.
Definition: normal.cpp:1228
ex normal(exmap &repl, exmap &rev_lookup, lst &modifier) const override
Implementation of ex::normal() for a product.
Definition: normal.cpp:2374
numeric integer_content() const override
Definition: normal.cpp:348
numeric max_coefficient() const override
Implementation ex::max_coefficient().
Definition: normal.cpp:1185
This class is a wrapper around CLN-numbers within the GiNaC class hierarchy.
Definition: numeric.h:82
ex to_rational(exmap &repl) const override
Implementation of ex::to_rational() for a numeric.
Definition: normal.cpp:2641
numeric integer_content() const override
Definition: normal.cpp:328
bool is_rational() const
True if object is an exact rational number, may even be complex (denominator may be unity).
Definition: numeric.cpp:1201
const numeric real() const
Real part of a number.
Definition: numeric.cpp:1339
ex to_polynomial(exmap &repl) const override
Implementation of ex::to_polynomial() for a numeric.
Definition: normal.cpp:2659
bool is_real() const
True if object is a real integer, rational or float (but not complex).
Definition: numeric.cpp:1208
const numeric numer() const
Numerator.
Definition: numeric.cpp:1356
bool is_integer() const
True if object is a non-complex integer.
Definition: numeric.cpp:1154
ex normal(exmap &repl, exmap &rev_lookup, lst &modifier) const override
Implementation of ex::normal() for a numeric.
Definition: normal.cpp:2225
const numeric denom() const
Denominator.
Definition: numeric.cpp:1387
const numeric imag() const
Imaginary part of a number.
Definition: numeric.cpp:1346
numeric max_coefficient() const override
Implementation ex::max_coefficient().
Definition: normal.cpp:1166
int int_length() const
Size in binary notation.
Definition: numeric.cpp:1418
ex smod(const numeric &xi) const override
Apply symmetric modular homomorphism to an expanded multivariate polynomial.
Definition: normal.cpp:1208
const numeric inverse() const
Inverse of a number.
Definition: numeric.cpp:1053
bool is_zero() const
True if object is zero.
Definition: numeric.cpp:1129
This class holds a two-component object, a basis and and exponent representing exponentiation.
Definition: power.h:39
ex normal(exmap &repl, exmap &rev_lookup, lst &modifier) const override
Implementation of ex::normal() for powers.
Definition: normal.cpp:2411
ex basis
Definition: power.h:105
ex to_rational(exmap &repl) const override
Implementation of ex::to_rational() for powers.
Definition: normal.cpp:2677
ex exponent
Definition: power.h:106
ex to_polynomial(exmap &repl) const override
Implementation of ex::to_polynomial() for powers.
Definition: normal.cpp:2687
ex var
Series variable (holds a symbol)
Definition: pseries.h:118
pseries(const ex &rel_, const epvector &ops_)
Construct pseries from a vector of coefficients and powers.
Definition: pseries.cpp:71
epvector seq
Vector of {coefficient, power} pairs.
Definition: pseries.h:115
ex point
Expansion point.
Definition: pseries.h:121
ex normal(exmap &repl, exmap &rev_lookup, lst &modifier) const override
Implementation of ex::normal() for pseries.
Definition: normal.cpp:2468
This class holds a relation consisting of two expressions and a logical relation between them.
Definition: relational.h:35
@ evaluated
.eval() has already done its job
Definition: flags.h:203
@ hash_calculated
.calchash() has already done its job
Definition: flags.h:205
@ no_pattern
disable pattern matching
Definition: flags.h:51
Basic CAS symbol.
Definition: symbol.h:39
ex to_rational(exmap &repl) const override
Implementation of ex::to_rational() for symbols.
Definition: normal.cpp:2625
ex normal(exmap &repl, exmap &rev_lookup, lst &modifier) const override
Implementation of ex::normal() for symbols.
Definition: normal.cpp:2215
ex to_polynomial(exmap &repl) const override
Implementation of ex::to_polynomial() for symbols.
Definition: normal.cpp:2632
Interface to GiNaC's constant types and some special constants.
Interface to GiNaC's light-weight expression handles.
Interface to sequences of expression pairs.
upvec factors
Definition: factor.cpp:1461
unsigned options
Definition: factor.cpp:2480
vector< int > k
Definition: factor.cpp:1466
ex x
Definition: factor.cpp:1641
size_t n
Definition: factor.cpp:1463
size_t c
Definition: factor.cpp:770
size_t r
Definition: factor.cpp:770
mvec m
Definition: factor.cpp:771
Interface to class signaling failure of operation.
#define is_ex_the_function(OBJ, FUNCNAME)
Definition: function.h:765
Interface to GiNaC's initially known functions.
Definition of GiNaC's lst.
Interface to symbolic matrices.
Interface to GiNaC's products of expressions.
Definition: add.cpp:38
const numeric I
Imaginary unit.
Definition: numeric.cpp:1433
ex denom(const ex &thisex)
Definition: ex.h:748
const numeric pow(const numeric &x, const numeric &y)
Definition: numeric.h:251
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:2051
std::map< ex, ex, ex_is_less > exmap
Definition: basic.h:50
ex sqrfree(const ex &a, const lst &l)
Compute a square-free factorization of a multivariate polynomial in Q[X].
Definition: normal.cpp:1888
std::vector< expair > epvector
expair-vector
Definition: expairseq.h:33
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:2851
const numeric abs(const numeric &x)
Absolute value.
Definition: numeric.cpp:2315
static ex interpolate(const ex &gamma, const numeric &xi, const ex &x, int degree_hint=1)
xi-adic polynomial interpolation
Definition: normal.cpp:1245
static void add_symbol(const ex &s, sym_desc_vec &v)
Definition: normal.cpp:163
bool is_negative(const numeric &x)
Definition: numeric.h:269
const ex _ex1
Definition: utils.cpp:385
bool is_rational(const numeric &x)
Definition: numeric.h:290
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
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:2743
function psi(const T1 &p1)
Definition: inifcns.h:165
ex rhs(const ex &thisex)
Definition: ex.h:817
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
static ex gcd_pf_pow_pow(const ex &a, const ex &b, ex *ca, ex *cb)
Definition: normal.cpp:1646
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
const numeric exp(const numeric &x)
Exponential function.
Definition: numeric.cpp:1439
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
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
static bool get_first_symbol(const ex &e, ex &x)
Return pointer to first symbol found in expression.
Definition: normal.cpp:95
const ex _ex_1
Definition: utils.cpp:352
const numeric smod(const numeric &a_, const numeric &b_)
Modulus (in symmetric representation).
Definition: numeric.cpp:2341
int degree(const ex &thisex, const ex &s)
Definition: ex.h:736
static ex gcd_pf_mul(const ex &a, const ex &b, ex *ca, ex *cb)
Definition: normal.cpp:1740
ex sqrfree_parfrac(const ex &a, const symbol &x)
Compute square-free partial fraction decomposition of rational function a(x).
Definition: normal.cpp:1947
ex lcm(const ex &a, const ex &b, bool check_args)
Compute LCM (Least Common Multiple) of multivariate polynomials in Z[X].
Definition: normal.cpp:1774
const numeric iquo(const numeric &a, const numeric &b)
Numeric integer quotient.
Definition: numeric.cpp:2404
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
const numeric isqrt(const numeric &x)
Integer numeric square root.
Definition: numeric.cpp:2482
ex normal(const ex &thisex)
Definition: ex.h:754
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
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
const numeric * _num1_p
Definition: utils.cpp:384
ex coeff(const ex &thisex, const ex &s, int n=1)
Definition: ex.h:742
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
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
ex numer(const ex &thisex)
Definition: ex.h:745
ex factor(const ex &poly, unsigned options)
Interface function to the outside world.
Definition: factor.cpp:2581
ex collect_common_factors(const ex &e)
Collect common factors in sums.
Definition: normal.cpp:2835
static numeric lcm_of_coefficients_denominators(const ex &e)
Compute LCM of denominators of coefficients of a polynomial.
Definition: normal.cpp:261
bool is_integer(const numeric &x)
Definition: numeric.h:272
static numeric lcmcoeff(const ex &e, const numeric &l)
Definition: normal.cpp:231
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
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
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
std::vector< sym_desc > sym_desc_vec
Definition: normal.cpp:160
static ex frac_cancel(const ex &n, const ex &d)
Fraction cancellation.
Definition: normal.cpp:2249
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
static void collect_symbols(const ex &e, sym_desc_vec &v)
Definition: normal.cpp:173
const ex _ex0
Definition: utils.cpp:369
std::vector< ex > exvector
Definition: basic.h:46
ex numer_denom(const ex &thisex)
Definition: ex.h:751
static ex gcd_pf_pow(const ex &a, const ex &b, ex *ca, ex *cb)
Definition: normal.cpp:1693
const numeric * _num0_p
Definition: utils.cpp:367
ex expand(const ex &thisex, unsigned options=0)
Definition: ex.h:715
This file defines several functions that work on univariate and multivariate polynomials and rational...
Makes the interface to the underlying bignum package available.
Interface to GiNaC's overloaded operators.
Interface to GiNaC's symbolic exponentiation (basis^exponent).
Interface to class for extended truncated power series.
Interface to relations between expressions.
@ use_sr_gcd
By default GiNaC uses modular GCD algorithm.
Definition: normal.h:61
@ no_part_factored
GiNaC tries to avoid expanding expressions when computing GCDs.
Definition: normal.h:55
@ no_heur_gcd
Usually GiNaC tries heuristic GCD first, because typically it's much faster than anything else.
Definition: normal.h:47
Function object for map().
Definition: basic.h:85
Function object to be applied by basic::normal().
Definition: normal.cpp:2184
ex operator()(const ex &e) override
Definition: normal.cpp:2185
This structure holds information about the highest and lowest degrees in which a symbol appears in tw...
Definition: normal.cpp:122
int max_deg
Maximum of deg_a and deg_b (Used for sorting)
Definition: normal.cpp:144
int ldeg_a
Lowest degree of symbol in polynomial "a".
Definition: normal.cpp:138
ex sym
Reference to symbol.
Definition: normal.cpp:129
int ldeg_b
Lowest degree of symbol in polynomial "b".
Definition: normal.cpp:141
bool operator<(const sym_desc &x) const
Comparison operator for sorting.
Definition: normal.cpp:150
int deg_b
Highest degree of symbol in polynomial "b".
Definition: normal.cpp:135
int deg_a
Highest degree of symbol in polynomial "a".
Definition: normal.cpp:132
size_t max_lcnops
Maximum number of terms of leading coefficient of symbol in both polynomials.
Definition: normal.cpp:147
sym_desc(const ex &s)
Initialize symbol, leave other variables uninitialized.
Definition: normal.cpp:124
Interface to GiNaC's symbolic objects.
Interface to several small and furry utilities needed within GiNaC but not of any interest to the use...

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