- GiNaC-devel - ginac.de

[PATCH 06/10] gcd_pf_mul: get rid of duplicate code.
by Alexei Sheplyakov 25 Aug '08

25 Aug '08

This function (which helps gcd() handle partially factored expressions) contained two copies of the same code. Remove one redundant copy. --- ginac/normal.cpp | 62 ++++++++++++++++++++---------------------------------- 1 files changed, 23 insertions(+), 39 deletions(-) diff --git a/ginac/normal.cpp b/ginac/normal.cpp index 5d044fc..2bb3a43 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -1746,45 +1746,29 @@ static ex gcd_pf_pow(const ex& a, const ex& b, ex* ca, ex* cb, bool check_args) static ex gcd_pf_mul(const ex& a, const ex& b, ex* ca, ex* cb, bool check_args) { - if (is_exactly_a<mul>(a)) { - if (is_exactly_a<mul>(b) && b.nops() > a.nops()) - goto factored_b; -factored_a: - size_t num = a.nops(); - exvector g; g.reserve(num); - exvector acc_ca; acc_ca.reserve(num); - ex part_b = b; - for (size_t i=0; i<num; i++) { - ex part_ca, part_cb; - g.push_back(gcd(a.op(i), part_b, &part_ca, &part_cb, check_args)); - acc_ca.push_back(part_ca); - part_b = part_cb; - } - if (ca) - *ca = (new mul(acc_ca))->setflag(status_flags::dynallocated); - if (cb) - *cb = part_b; - return (new mul(g))->setflag(status_flags::dynallocated); - } else if (is_exactly_a<mul>(b)) { - if (is_exactly_a<mul>(a) && a.nops() > b.nops()) - goto factored_a; -factored_b: - size_t num = b.nops(); - exvector g; g.reserve(num); - exvector acc_cb; acc_cb.reserve(num); - ex part_a = a; - for (size_t i=0; i<num; i++) { - ex part_ca, part_cb; - g.push_back(gcd(part_a, b.op(i), &part_ca, &part_cb, check_args)); - acc_cb.push_back(part_cb); - part_a = part_ca; - } - if (ca) - *ca = part_a; - if (cb) - *cb = (new mul(acc_cb))->setflag(status_flags::dynallocated); - return (new mul(g))->setflag(status_flags::dynallocated); - } + if (is_exactly_a<mul>(a) && is_exactly_a<mul>(b) + && (b.nops() > a.nops())) + return gcd_pf_mul(b, a, cb, ca, check_args); + + if (is_exactly_a<mul>(b) && (!is_exactly_a<mul>(a))) + return gcd_pf_mul(b, a, cb, ca, check_args); + + GINAC_ASSERT(is_exactly_a<mul>(a)); + size_t num = a.nops(); + exvector g; g.reserve(num); + exvector acc_ca; acc_ca.reserve(num); + ex part_b = b; + for (size_t i=0; i<num; i++) { + ex part_ca, part_cb; + g.push_back(gcd(a.op(i), part_b, &part_ca, &part_cb, check_args)); + acc_ca.push_back(part_ca); + part_b = part_cb; + } + if (ca) + *ca = (new mul(acc_ca))->setflag(status_flags::dynallocated); + if (cb) + *cb = part_b; + return (new mul(g))->setflag(status_flags::dynallocated); } /** Compute LCM (Least Common Multiple) of multivariate polynomials in Z[X]. -- 1.5.6 Best regards, Alexei -- All science is either physics or stamp collecting.

1 0

[PATCH 05/10] gcd(): allow user to override (some of) heuristics.
by Alexei Sheplyakov 25 Aug '08

25 Aug '08

GiNaC tries to avoid expanding expressions while computing GCDs and applies a number of heuristics. Usually this improves performance, but in some cases it doesn't. Allow user to switch off heuristics. Part 5: * gcd(): don't use heuristic GCD algorithm if gcd_options::no_heur_gcd flag is set. * gcd(): don't handle partially factored expressions in a special way if gcd_options::no_part_factored flag is set. --- ginac/normal.cpp | 40 ++++++++++++++++++++++------------------ 1 files changed, 22 insertions(+), 18 deletions(-) diff --git a/ginac/normal.cpp b/ginac/normal.cpp index 0227f4e..5d044fc 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -1465,12 +1465,14 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned optio } // Partially factored cases (to avoid expanding large expressions) - if (is_exactly_a<mul>(a) || is_exactly_a<mul>(b)) - return gcd_pf_mul(a, b, ca, cb, check_args); + if (!(options & gcd_options::no_part_factored)) { + if (is_exactly_a<mul>(a) || is_exactly_a<mul>(b)) + return gcd_pf_mul(a, b, ca, cb, check_args); #if FAST_COMPARE - if (is_exactly_a<power>(a) || is_exactly_a<power>(b)) - return gcd_pf_pow(a, b, ca, cb, check_args); + if (is_exactly_a<power>(a) || is_exactly_a<power>(b)) + return gcd_pf_pow(a, b, ca, cb, check_args); #endif + } // Some trivial cases ex aex = a.expand(), bex = b.expand(); @@ -1601,23 +1603,25 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned optio // Try heuristic algorithm first, fall back to PRS if that failed ex g; - bool found = heur_gcd(g, aex, bex, ca, cb, var); - if (found) { - // heur_gcd have already computed cofactors... - if (g.is_equal(_ex1)) { - // ... but we want to keep them factored if possible. - if (ca) - *ca = a; - if (cb) - *cb = b; + if (!(options & gcd_options::no_heur_gcd)) { + bool found = heur_gcd(g, aex, bex, ca, cb, var); + if (found) { + // heur_gcd have already computed cofactors... + if (g.is_equal(_ex1)) { + // ... but we want to keep them factored if possible. + if (ca) + *ca = a; + if (cb) + *cb = b; + } + return g; } - return g; - } #if STATISTICS - else { - heur_gcd_failed++; - } + else { + heur_gcd_failed++; + } #endif + } g = sr_gcd(aex, bex, var); if (g.is_equal(_ex1)) { -- 1.5.6 Best regards, Alexei -- All science is either physics or stamp collecting.

1 0

[PATCH 04/10] refactor gcd() a little bit (no functional changes).
by Alexei Sheplyakov 25 Aug '08

25 Aug '08

GiNaC tries to avoid expanding expressions while computing GCDs and applies a number of heuristics. Usually this improves performance, but in some cases it doesn't. Allow user to switch off heuristics. Part 4: refactor gcd() a little bit, so subsequent patch(es) won't be so big and ugly. --- ginac/normal.cpp | 41 ++++++++++++++++++++++------------------- 1 files changed, 22 insertions(+), 19 deletions(-) diff --git a/ginac/normal.cpp b/ginac/normal.cpp index 0cb9100..0227f4e 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -1602,33 +1602,36 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned optio // Try heuristic algorithm first, fall back to PRS if that failed ex g; bool found = heur_gcd(g, aex, bex, ca, cb, var); - if (!found) { -#if STATISTICS - heur_gcd_failed++; -#endif - g = sr_gcd(aex, bex, var); - if (g.is_equal(_ex1)) { - // Keep cofactors factored if possible - if (ca) - *ca = a; - if (cb) - *cb = b; - } else { - if (ca) - divide(aex, g, *ca, false); - if (cb) - divide(bex, g, *cb, false); - } - } else { + if (found) { + // heur_gcd have already computed cofactors... if (g.is_equal(_ex1)) { - // Keep cofactors factored if possible + // ... but we want to keep them factored if possible. if (ca) *ca = a; if (cb) *cb = b; } + return g; } +#if STATISTICS + else { + heur_gcd_failed++; + } +#endif + g = sr_gcd(aex, bex, var); + if (g.is_equal(_ex1)) { + // Keep cofactors factored if possible + if (ca) + *ca = a; + if (cb) + *cb = b; + } else { + if (ca) + divide(aex, g, *ca, false); + if (cb) + divide(bex, g, *cb, false); + } return g; } -- 1.5.6 Best regards, Alexei -- All science is either physics or stamp collecting.

1 0

[PATCH 03/10] introduce gcd_pf_mul: gcd helper to handle partially factored expressions.
by Alexei Sheplyakov 25 Aug '08

25 Aug '08

GiNaC tries to avoid expanding expressions while computing GCDs and applies a number of heuristics. Usually this improves performance, but in some cases it doesn't. Allow user to switch off heuristics. Part 3: Move the code handling products from gcd() into a separate function. This is *really* only code move, everything else should be considered a bug. --- ginac/normal.cpp | 89 +++++++++++++++++++++++++++++------------------------ 1 files changed, 49 insertions(+), 40 deletions(-) diff --git a/ginac/normal.cpp b/ginac/normal.cpp index 0af5aad..0cb9100 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -1419,6 +1419,10 @@ static bool heur_gcd(ex& res, const ex& a, const ex& b, ex *ca, ex *cb, // large expressions). At least one of the arguments should be a power. static ex gcd_pf_pow(const ex& a, const ex& b, ex* ca, ex* cb, bool check_args); +// gcd helper to handle partially factored polynomials (to avoid expanding +// large expressions). At least one of the arguments should be a product. +static ex gcd_pf_mul(const ex& a, const ex& b, ex* ca, ex* cb, bool check_args); + /** Compute GCD (Greatest Common Divisor) of multivariate polynomials a(X) * and b(X) in Z[X]. Optionally also compute the cofactors of a and b, * defined by a = ca * gcd(a, b) and b = cb * gcd(a, b). @@ -1461,46 +1465,8 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned optio } // Partially factored cases (to avoid expanding large expressions) - if (is_exactly_a<mul>(a)) { - if (is_exactly_a<mul>(b) && b.nops() > a.nops()) - goto factored_b; -factored_a: - size_t num = a.nops(); - exvector g; g.reserve(num); - exvector acc_ca; acc_ca.reserve(num); - ex part_b = b; - for (size_t i=0; i<num; i++) { - ex part_ca, part_cb; - g.push_back(gcd(a.op(i), part_b, &part_ca, &part_cb, check_args)); - acc_ca.push_back(part_ca); - part_b = part_cb; - } - if (ca) - *ca = (new mul(acc_ca))->setflag(status_flags::dynallocated); - if (cb) - *cb = part_b; - return (new mul(g))->setflag(status_flags::dynallocated); - } else if (is_exactly_a<mul>(b)) { - if (is_exactly_a<mul>(a) && a.nops() > b.nops()) - goto factored_a; -factored_b: - size_t num = b.nops(); - exvector g; g.reserve(num); - exvector acc_cb; acc_cb.reserve(num); - ex part_a = a; - for (size_t i=0; i<num; i++) { - ex part_ca, part_cb; - g.push_back(gcd(part_a, b.op(i), &part_ca, &part_cb, check_args)); - acc_cb.push_back(part_cb); - part_a = part_ca; - } - if (ca) - *ca = part_a; - if (cb) - *cb = (new mul(acc_cb))->setflag(status_flags::dynallocated); - return (new mul(g))->setflag(status_flags::dynallocated); - } - + if (is_exactly_a<mul>(a) || is_exactly_a<mul>(b)) + return gcd_pf_mul(a, b, ca, cb, check_args); #if FAST_COMPARE if (is_exactly_a<power>(a) || is_exactly_a<power>(b)) return gcd_pf_pow(a, b, ca, cb, check_args); @@ -1771,6 +1737,49 @@ static ex gcd_pf_pow(const ex& a, const ex& b, ex* ca, ex* cb, bool check_args) } } +static ex gcd_pf_mul(const ex& a, const ex& b, ex* ca, ex* cb, bool check_args) +{ + if (is_exactly_a<mul>(a)) { + if (is_exactly_a<mul>(b) && b.nops() > a.nops()) + goto factored_b; +factored_a: + size_t num = a.nops(); + exvector g; g.reserve(num); + exvector acc_ca; acc_ca.reserve(num); + ex part_b = b; + for (size_t i=0; i<num; i++) { + ex part_ca, part_cb; + g.push_back(gcd(a.op(i), part_b, &part_ca, &part_cb, check_args)); + acc_ca.push_back(part_ca); + part_b = part_cb; + } + if (ca) + *ca = (new mul(acc_ca))->setflag(status_flags::dynallocated); + if (cb) + *cb = part_b; + return (new mul(g))->setflag(status_flags::dynallocated); + } else if (is_exactly_a<mul>(b)) { + if (is_exactly_a<mul>(a) && a.nops() > b.nops()) + goto factored_a; +factored_b: + size_t num = b.nops(); + exvector g; g.reserve(num); + exvector acc_cb; acc_cb.reserve(num); + ex part_a = a; + for (size_t i=0; i<num; i++) { + ex part_ca, part_cb; + g.push_back(gcd(part_a, b.op(i), &part_ca, &part_cb, check_args)); + acc_cb.push_back(part_cb); + part_a = part_ca; + } + if (ca) + *ca = part_a; + if (cb) + *cb = (new mul(acc_cb))->setflag(status_flags::dynallocated); + return (new mul(g))->setflag(status_flags::dynallocated); + } +} + /** Compute LCM (Least Common Multiple) of multivariate polynomials in Z[X]. * * @param a first multivariate polynomial -- 1.5.6 Best regards, Alexei -- All science is either physics or stamp collecting.

1 0

[PATCH 02/10] introduce gcd_pf_pow: gcd helper to handle partially factored expressions.
by Alexei Sheplyakov 25 Aug '08

25 Aug '08

GiNaC tries to avoid expanding expressions while computing GCDs and applies a number of heuristics. Usually this improves performance, but in some cases it doesn't. Allow user to switch off heuristics. Part 2: Move the code handling powers from gcd() into a separate function. This is *really* only code move, everything else should be considered a bug. --- ginac/normal.cpp | 212 ++++++++++++++++++++++++++++-------------------------- 1 files changed, 110 insertions(+), 102 deletions(-) diff --git a/ginac/normal.cpp b/ginac/normal.cpp index 9ec7574..0af5aad 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -1415,6 +1415,10 @@ static bool heur_gcd(ex& res, const ex& a, const ex& b, ex *ca, ex *cb, } +// gcd helper to handle partially factored polynomials (to avoid expanding +// large expressions). At least one of the arguments should be a power. +static ex gcd_pf_pow(const ex& a, const ex& b, ex* ca, ex* cb, bool check_args); + /** Compute GCD (Greatest Common Divisor) of multivariate polynomials a(X) * and b(X) in Z[X]. Optionally also compute the cofactors of a and b, * defined by a = ca * gcd(a, b) and b = cb * gcd(a, b). @@ -1498,108 +1502,8 @@ factored_b: } #if FAST_COMPARE - // Input polynomials of the form poly^n are sometimes also trivial - if (is_exactly_a<power>(a)) { - ex p = a.op(0); - const ex& exp_a = a.op(1); - if (is_exactly_a<power>(b)) { - ex pb = b.op(0); - const ex& exp_b = b.op(1); - if (p.is_equal(pb)) { - // a = p^n, b = p^m, gcd = p^min(n, m) - if (exp_a < exp_b) { - if (ca) - *ca = _ex1; - if (cb) - *cb = power(p, exp_b - exp_a); - return power(p, exp_a); - } else { - if (ca) - *ca = power(p, exp_a - exp_b); - if (cb) - *cb = _ex1; - return power(p, exp_b); - } - } else { - ex p_co, pb_co; - ex p_gcd = gcd(p, pb, &p_co, &pb_co, check_args); - if (p_gcd.is_equal(_ex1)) { - // a(x) = p(x)^n, b(x) = p_b(x)^m, gcd (p, p_b) = 1 ==> - // gcd(a,b) = 1 - if (ca) - *ca = a; - if (cb) - *cb = b; - return _ex1; - // XXX: do I need to check for p_gcd = -1? - } else { - // there are common factors: - // a(x) = g(x)^n A(x)^n, b(x) = g(x)^m B(x)^m ==> - // gcd(a, b) = g(x)^n gcd(A(x)^n, g(x)^(n-m) B(x)^m - if (exp_a < exp_b) { - return power(p_gcd, exp_a)* - gcd(power(p_co, exp_a), power(p_gcd, exp_b-exp_a)*power(pb_co, exp_b), ca, cb, false); - } else { - return power(p_gcd, exp_b)* - gcd(power(p_gcd, exp_a - exp_b)*power(p_co, exp_a), power(pb_co, exp_b), ca, cb, false); - } - } // p_gcd.is_equal(_ex1) - } // p.is_equal(pb) - - } else { - if (p.is_equal(b)) { - // a = p^n, b = p, gcd = p - if (ca) - *ca = power(p, a.op(1) - 1); - if (cb) - *cb = _ex1; - return p; - } - - ex p_co, bpart_co; - ex p_gcd = gcd(p, b, &p_co, &bpart_co, false); - - if (p_gcd.is_equal(_ex1)) { - // a(x) = p(x)^n, gcd(p, b) = 1 ==> gcd(a, b) = 1 - if (ca) - *ca = a; - if (cb) - *cb = b; - return _ex1; - } else { - // 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)) - return p_gcd*gcd(power(p_gcd, exp_a-1)*power(p_co, exp_a), bpart_co, ca, cb, false); - } - } // is_exactly_a<power>(b) - - } else if (is_exactly_a<power>(b)) { - ex p = b.op(0); - if (p.is_equal(a)) { - // a = p, b = p^n, gcd = p - if (ca) - *ca = _ex1; - if (cb) - *cb = power(p, b.op(1) - 1); - return p; - } - - ex p_co, apart_co; - const ex& exp_b(b.op(1)); - ex p_gcd = gcd(a, p, &apart_co, &p_co, false); - if (p_gcd.is_equal(_ex1)) { - // b=p(x)^n, gcd(a, p) = 1 ==> gcd(a, b) == 1 - if (ca) - *ca = a; - if (cb) - *cb = b; - return _ex1; - } else { - // there are common factors: - // a(x) = g(x) A(x), b(x) = g(x)^n B(x)^n ==> gcd = g(x) gcd(g(x)^(n-1) A(x)^n, B(x)) - - return p_gcd*gcd(apart_co, power(p_gcd, exp_b-1)*power(p_co, exp_b), ca, cb, false); - } // p_gcd.is_equal(_ex1) - } + if (is_exactly_a<power>(a) || is_exactly_a<power>(b)) + return gcd_pf_pow(a, b, ca, cb, check_args); #endif // Some trivial cases @@ -1762,6 +1666,110 @@ factored_b: return g; } +static ex gcd_pf_pow(const ex& a, const ex& b, ex* ca, ex* cb, bool check_args) +{ + if (is_exactly_a<power>(a)) { + ex p = a.op(0); + const ex& exp_a = a.op(1); + if (is_exactly_a<power>(b)) { + ex pb = b.op(0); + const ex& exp_b = b.op(1); + if (p.is_equal(pb)) { + // a = p^n, b = p^m, gcd = p^min(n, m) + if (exp_a < exp_b) { + if (ca) + *ca = _ex1; + if (cb) + *cb = power(p, exp_b - exp_a); + return power(p, exp_a); + } else { + if (ca) + *ca = power(p, exp_a - exp_b); + if (cb) + *cb = _ex1; + return power(p, exp_b); + } + } else { + ex p_co, pb_co; + ex p_gcd = gcd(p, pb, &p_co, &pb_co, check_args); + if (p_gcd.is_equal(_ex1)) { + // a(x) = p(x)^n, b(x) = p_b(x)^m, gcd (p, p_b) = 1 ==> + // gcd(a,b) = 1 + if (ca) + *ca = a; + if (cb) + *cb = b; + return _ex1; + // XXX: do I need to check for p_gcd = -1? + } else { + // there are common factors: + // a(x) = g(x)^n A(x)^n, b(x) = g(x)^m B(x)^m ==> + // gcd(a, b) = g(x)^n gcd(A(x)^n, g(x)^(n-m) B(x)^m + if (exp_a < exp_b) { + return power(p_gcd, exp_a)* + gcd(power(p_co, exp_a), power(p_gcd, exp_b-exp_a)*power(pb_co, exp_b), ca, cb, false); + } else { + return power(p_gcd, exp_b)* + gcd(power(p_gcd, exp_a - exp_b)*power(p_co, exp_a), power(pb_co, exp_b), ca, cb, false); + } + } // p_gcd.is_equal(_ex1) + } // p.is_equal(pb) + + } else { + if (p.is_equal(b)) { + // a = p^n, b = p, gcd = p + if (ca) + *ca = power(p, a.op(1) - 1); + if (cb) + *cb = _ex1; + return p; + } + + ex p_co, bpart_co; + ex p_gcd = gcd(p, b, &p_co, &bpart_co, false); + + if (p_gcd.is_equal(_ex1)) { + // a(x) = p(x)^n, gcd(p, b) = 1 ==> gcd(a, b) = 1 + if (ca) + *ca = a; + if (cb) + *cb = b; + return _ex1; + } else { + // 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)) + return p_gcd*gcd(power(p_gcd, exp_a-1)*power(p_co, exp_a), bpart_co, ca, cb, false); + } + } // is_exactly_a<power>(b) + + } else if (is_exactly_a<power>(b)) { + ex p = b.op(0); + if (p.is_equal(a)) { + // a = p, b = p^n, gcd = p + if (ca) + *ca = _ex1; + if (cb) + *cb = power(p, b.op(1) - 1); + return p; + } + + ex p_co, apart_co; + const ex& exp_b(b.op(1)); + ex p_gcd = gcd(a, p, &apart_co, &p_co, false); + if (p_gcd.is_equal(_ex1)) { + // b=p(x)^n, gcd(a, p) = 1 ==> gcd(a, b) == 1 + if (ca) + *ca = a; + if (cb) + *cb = b; + return _ex1; + } else { + // there are common factors: + // a(x) = g(x) A(x), b(x) = g(x)^n B(x)^n ==> gcd = g(x) gcd(g(x)^(n-1) A(x)^n, B(x)) + + return p_gcd*gcd(apart_co, power(p_gcd, exp_b-1)*power(p_co, exp_b), ca, cb, false); + } // p_gcd.is_equal(_ex1) + } +} /** Compute LCM (Least Common Multiple) of multivariate polynomials in Z[X]. * -- 1.5.6 Best regards, Alexei -- All science is either physics or stamp collecting.

1 0

[PATCH 01/10] gcd: allow user to override (some of) heuristics.
by Alexei Sheplyakov 25 Aug '08

25 Aug '08

GiNaC tries to avoid expanding expressions while computing GCDs and applies a number of heuristics. Usually this improves performance, but in some cases it doesn't. Allow user to switch off heuristics. Part 1: * add new (optional) argument to gcd() to control its behaviour. * introduce gcd_options structure. N.B. No actual code changes so far, the actual handling of newly introduced options is the subject of further patches. --- ginac/normal.cpp | 2 +- ginac/normal.h | 25 ++++++++++++++++++++++++- 2 files changed, 25 insertions(+), 2 deletions(-) diff --git a/ginac/normal.cpp b/ginac/normal.cpp index e756a67..9ec7574 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -1426,7 +1426,7 @@ static bool heur_gcd(ex& res, const ex& a, const ex& b, ex *ca, ex *cb, * @param check_args check whether a and b are polynomials with rational * coefficients (defaults to "true") * @return the GCD as a new expression */ -ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args) +ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned options) { #if STATISTICS gcd_called++; diff --git a/ginac/normal.h b/ginac/normal.h index 5f44427..5a735e3 100644 --- a/ginac/normal.h +++ b/ginac/normal.h @@ -30,6 +30,28 @@ namespace GiNaC { +/** + * Flags to control the behaviour of gcd() and friends + */ +struct gcd_options +{ + enum { + /** + * Usually GiNaC tries heuristic GCD algorithm before PRS. + * Some people don't like this, so here's a flag to disable it. + */ + no_heur_gcd = 2, + /** + * GiNaC tries to avoid expanding expressions when computing + * GCDs. This is a good idea, but some people dislike it. + * Hence the flag to disable special handling of partially + * factored polynomials. DON'T SET THIS unless you *really* + * know what are you doing! + */ + no_part_factored = 4 + }; +}; + class ex; class symbol; @@ -52,7 +74,8 @@ extern ex sprem(const ex &a, const ex &b, const ex &x, bool check_args = true); extern bool divide(const ex &a, const ex &b, ex &q, bool check_args = true); // Polynomial GCD in Z[X], cofactors are returned in ca and cb, if desired -extern ex gcd(const ex &a, const ex &b, ex *ca = NULL, ex *cb = NULL, bool check_args = true); +extern ex gcd(const ex &a, const ex &b, ex *ca = NULL, ex *cb = NULL, + bool check_args = true, unsigned options = 0); // Polynomial LCM in Z[X] extern ex lcm(const ex &a, const ex &b, bool check_args = true); -- 1.5.6 Best regards, Alexei -- All science is either physics or stamp collecting.

1 0

[PATCH] inifscn_nstdsums: make functions handling Li/G transformations reentrant.
by Alexei Sheplyakov 24 Aug '08

24 Aug '08

Explicitly pass the dummy symbols instead of using a global variable. As a side effect, the code is more clear now (that's a bit subjective, but anyway). --- ginac/inifcns_nstdsums.cpp | 128 +++++++++++++++++++++++++------------------ 1 files changed, 74 insertions(+), 54 deletions(-) diff --git a/ginac/inifcns_nstdsums.cpp b/ginac/inifcns_nstdsums.cpp index 01a5756..dd8625a 100644 --- a/ginac/inifcns_nstdsums.cpp +++ b/ginac/inifcns_nstdsums.cpp @@ -518,16 +518,12 @@ cln::cl_N mLi_do_summation(const lst& m, const lst& x) lst convert_parameter_Li_to_H(const lst& m, const lst& x, ex& pf); -// holding dummy-symbols for the G/Li transformations -std::vector<ex> gsyms; - - // type used by the transformation functions for G typedef std::vector<int> Gparameter; // G_eval1-function for G transformations -ex G_eval1(int a, int scale) +ex G_eval1(int a, int scale, const exvector& gsyms) { if (a != 0) { const ex& scs = gsyms[std::abs(scale)]; @@ -544,7 +540,7 @@ ex G_eval1(int a, int scale) // G_eval-function for G transformations -ex G_eval(const Gparameter& a, int scale) +ex G_eval(const Gparameter& a, int scale, const exvector& gsyms) { // check for properties of G ex sc = gsyms[std::abs(scale)]; @@ -578,7 +574,7 @@ ex G_eval(const Gparameter& a, int scale) for (; it != a.end(); ++it) { short_a.push_back(*it); } - ex result = G_eval1(a.front(), scale) * G_eval(short_a, scale); + ex result = G_eval1(a.front(), scale, gsyms) * G_eval(short_a, scale, gsyms); it = short_a.begin(); for (int i=1; i<count_ones; ++i) { ++it; @@ -596,14 +592,14 @@ ex G_eval(const Gparameter& a, int scale) for (; it2 != short_a.end(); ++it2) { newa.push_back(*it2); } - result -= G_eval(newa, scale); + result -= G_eval(newa, scale, gsyms); } return result / count_ones; } // G({1,...,1};y) -> G({1};y)^k / k! if (all_ones && a.size() > 1) { - return pow(G_eval1(a.front(),scale), count_ones) / factorial(count_ones); + return pow(G_eval1(a.front(),scale, gsyms), count_ones) / factorial(count_ones); } // G({0,...,0};y) -> log(y)^k / k! @@ -696,7 +692,7 @@ Gparameter prepare_pending_integrals(const Gparameter& pending_integrals, int sc // handles trailing zeroes for an otherwise convergent integral -ex trailing_zeros_G(const Gparameter& a, int scale) +ex trailing_zeros_G(const Gparameter& a, int scale, const exvector& gsyms) { bool convergent; int depth, trailing_zeros; @@ -708,23 +704,23 @@ ex trailing_zeros_G(const Gparameter& a, int scale) if ((trailing_zeros > 0) && (depth > 0)) { ex result; Gparameter new_a(a.begin(), a.end()-1); - result += G_eval1(0, scale) * trailing_zeros_G(new_a, scale); + result += G_eval1(0, scale, gsyms) * trailing_zeros_G(new_a, scale, gsyms); for (Gparameter::const_iterator it = a.begin(); it != last; ++it) { Gparameter new_a(a.begin(), it); new_a.push_back(0); new_a.insert(new_a.end(), it, a.end()-1); - result -= trailing_zeros_G(new_a, scale); + result -= trailing_zeros_G(new_a, scale, gsyms); } return result / trailing_zeros; } else { - return G_eval(a, scale); + return G_eval(a, scale, gsyms); } } // G transformation [VSW] (57),(58) -ex depth_one_trafo_G(const Gparameter& pending_integrals, const Gparameter& a, int scale) +ex depth_one_trafo_G(const Gparameter& pending_integrals, const Gparameter& a, int scale, const exvector& gsyms) { // pendint = ( y1, b1, ..., br ) // a = ( 0, ..., 0, amin ) @@ -755,15 +751,21 @@ ex depth_one_trafo_G(const Gparameter& pending_integrals, const Gparameter& a, i result -= I*Pi; } if (psize) { - result *= trailing_zeros_G(convert_pending_integrals_G(pending_integrals), pending_integrals.front()); + result *= trailing_zeros_G(convert_pending_integrals_G(pending_integrals), + pending_integrals.front(), + gsyms); } // G(y2_{-+}; sr) - result += trailing_zeros_G(convert_pending_integrals_G(new_pending_integrals), new_pending_integrals.front()); + result += trailing_zeros_G(convert_pending_integrals_G(new_pending_integrals), + new_pending_integrals.front(), + gsyms); // G(0; sr) new_pending_integrals.back() = 0; - result -= trailing_zeros_G(convert_pending_integrals_G(new_pending_integrals), new_pending_integrals.front()); + result -= trailing_zeros_G(convert_pending_integrals_G(new_pending_integrals), + new_pending_integrals.front(), + gsyms); return result; } @@ -775,14 +777,16 @@ ex depth_one_trafo_G(const Gparameter& pending_integrals, const Gparameter& a, i //term zeta_m result -= zeta(a.size()); if (psize) { - result *= trailing_zeros_G(convert_pending_integrals_G(pending_integrals), pending_integrals.front()); + result *= trailing_zeros_G(convert_pending_integrals_G(pending_integrals), + pending_integrals.front(), + gsyms); } // term int_0^sr dt/t G_{m-1}( (1/y2)_{+-}; 1/t ) // = int_0^sr dt/t G_{m-1}( t_{+-}; y2 ) Gparameter new_a(a.begin()+1, a.end()); new_pending_integrals.push_back(0); - result -= depth_one_trafo_G(new_pending_integrals, new_a, scale); + result -= depth_one_trafo_G(new_pending_integrals, new_a, scale, gsyms); // term int_0^y2 dt/t G_{m-1}( (1/y2)_{+-}; 1/t ) // = int_0^y2 dt/t G_{m-1}( t_{+-}; y2 ) @@ -790,10 +794,12 @@ ex depth_one_trafo_G(const Gparameter& pending_integrals, const Gparameter& a, i new_pending_integrals_2.push_back(scale); new_pending_integrals_2.push_back(0); if (psize) { - result += trailing_zeros_G(convert_pending_integrals_G(pending_integrals), pending_integrals.front()) - * depth_one_trafo_G(new_pending_integrals_2, new_a, scale); + result += trailing_zeros_G(convert_pending_integrals_G(pending_integrals), + pending_integrals.front(), + gsyms) + * depth_one_trafo_G(new_pending_integrals_2, new_a, scale, gsyms); } else { - result += depth_one_trafo_G(new_pending_integrals_2, new_a, scale); + result += depth_one_trafo_G(new_pending_integrals_2, new_a, scale, gsyms); } return result; @@ -802,11 +808,13 @@ ex depth_one_trafo_G(const Gparameter& pending_integrals, const Gparameter& a, i // forward declaration ex shuffle_G(const Gparameter & a0, const Gparameter & a1, const Gparameter & a2, - const Gparameter& pendint, const Gparameter& a_old, int scale); + const Gparameter& pendint, const Gparameter& a_old, int scale, + const exvector& gsyms); // G transformation [VSW] -ex G_transform(const Gparameter& pendint, const Gparameter& a, int scale) +ex G_transform(const Gparameter& pendint, const Gparameter& a, int scale, + const exvector& gsyms) { // main recursion routine // @@ -832,10 +840,12 @@ ex G_transform(const Gparameter& pendint, const Gparameter& a, int scale) if (a.size() == 0) { result = 1; } else { - result = G_eval(a, scale); + result = G_eval(a, scale, gsyms); } if (pendint.size() > 0) { - result *= trailing_zeros_G(convert_pending_integrals_G(pendint), pendint.front()); + result *= trailing_zeros_G(convert_pending_integrals_G(pendint), + pendint.front(), + gsyms); } return result; } @@ -844,12 +854,12 @@ ex G_transform(const Gparameter& pendint, const Gparameter& a, int scale) if (trailing_zeros > 0) { ex result; Gparameter new_a(a.begin(), a.end()-1); - result += G_eval1(0, scale) * G_transform(pendint, new_a, scale); + result += G_eval1(0, scale, gsyms) * G_transform(pendint, new_a, scale, gsyms); for (Gparameter::const_iterator it = a.begin(); it != firstzero; ++it) { Gparameter new_a(a.begin(), it); new_a.push_back(0); new_a.insert(new_a.end(), it, a.end()-1); - result -= G_transform(pendint, new_a, scale); + result -= G_transform(pendint, new_a, scale, gsyms); } return result / trailing_zeros; } @@ -857,15 +867,17 @@ ex G_transform(const Gparameter& pendint, const Gparameter& a, int scale) // convergence case if (convergent) { if (pendint.size() > 0) { - return G_eval(convert_pending_integrals_G(pendint), pendint.front()) * G_eval(a, scale); + return G_eval(convert_pending_integrals_G(pendint), + pendint.front(), gsyms)* + G_eval(a, scale, gsyms); } else { - return G_eval(a, scale); + return G_eval(a, scale, gsyms); } } // call basic transformation for depth equal one if (depth == 1) { - return depth_one_trafo_G(pendint, a, scale); + return depth_one_trafo_G(pendint, a, scale, gsyms); } // do recursion @@ -880,9 +892,10 @@ ex G_transform(const Gparameter& pendint, const Gparameter& a, int scale) Gparameter a1(a.begin(),min_it+1); Gparameter a2(min_it+1,a.end()); - ex result = G_transform(pendint,a2,scale)*G_transform(empty,a1,scale); + ex result = G_transform(pendint, a2, scale, gsyms)* + G_transform(empty, a1, scale, gsyms); - result -= shuffle_G(empty,a1,a2,pendint,a,scale); + result -= shuffle_G(empty, a1, a2, pendint, a, scale, gsyms); return result; } @@ -893,9 +906,10 @@ ex G_transform(const Gparameter& pendint, const Gparameter& a, int scale) Gparameter new_pendint = prepare_pending_integrals(pendint, a[min_it_pos]); Gparameter new_a = a; new_a[min_it_pos] = 0; - ex result = G_transform(empty, new_a, scale); + ex result = G_transform(empty, new_a, scale, gsyms); if (pendint.size() > 0) { - result *= trailing_zeros_G(convert_pending_integrals_G(pendint), pendint.front()); + result *= trailing_zeros_G(convert_pending_integrals_G(pendint), + pendint.front(), gsyms); } // other terms @@ -904,29 +918,33 @@ ex G_transform(const Gparameter& pendint, const Gparameter& a, int scale) if (changeit != new_a.begin()) { // smallest in the middle new_pendint.push_back(*changeit); - result -= trailing_zeros_G(convert_pending_integrals_G(new_pendint), new_pendint.front()) - * G_transform(empty, new_a, scale); + result -= trailing_zeros_G(convert_pending_integrals_G(new_pendint), + new_pendint.front(), gsyms)* + G_transform(empty, new_a, scale, gsyms); int buffer = *changeit; *changeit = *min_it; - result += G_transform(new_pendint, new_a, scale); + result += G_transform(new_pendint, new_a, scale, gsyms); *changeit = buffer; new_pendint.pop_back(); --changeit; new_pendint.push_back(*changeit); - result += trailing_zeros_G(convert_pending_integrals_G(new_pendint), new_pendint.front()) - * G_transform(empty, new_a, scale); + result += trailing_zeros_G(convert_pending_integrals_G(new_pendint), + new_pendint.front(), gsyms)* + G_transform(empty, new_a, scale, gsyms); *changeit = *min_it; - result -= G_transform(new_pendint, new_a, scale); + result -= G_transform(new_pendint, new_a, scale, gsyms); } else { // smallest at the front new_pendint.push_back(scale); - result += trailing_zeros_G(convert_pending_integrals_G(new_pendint), new_pendint.front()) - * G_transform(empty, new_a, scale); + result += trailing_zeros_G(convert_pending_integrals_G(new_pendint), + new_pendint.front(), gsyms)* + G_transform(empty, new_a, scale, gsyms); new_pendint.back() = *changeit; - result -= trailing_zeros_G(convert_pending_integrals_G(new_pendint), new_pendint.front()) - * G_transform(empty, new_a, scale); + result -= trailing_zeros_G(convert_pending_integrals_G(new_pendint), + new_pendint.front(), gsyms)* + G_transform(empty, new_a, scale, gsyms); *changeit = *min_it; - result += G_transform(new_pendint, new_a, scale); + result += G_transform(new_pendint, new_a, scale, gsyms); } return result; } @@ -935,27 +953,28 @@ ex G_transform(const Gparameter& pendint, const Gparameter& a, int scale) // shuffles the two parameter list a1 and a2 and calls G_transform for every term except // for the one that is equal to a_old ex shuffle_G(const Gparameter & a0, const Gparameter & a1, const Gparameter & a2, - const Gparameter& pendint, const Gparameter& a_old, int scale) + const Gparameter& pendint, const Gparameter& a_old, int scale, + const exvector& gsyms) { if (a1.size()==0 && a2.size()==0) { // veto the one configuration we don't want if ( a0 == a_old ) return 0; - return G_transform(pendint,a0,scale); + return G_transform(pendint, a0, scale, gsyms); } if (a2.size()==0) { Gparameter empty; Gparameter aa0 = a0; aa0.insert(aa0.end(),a1.begin(),a1.end()); - return shuffle_G(aa0,empty,empty,pendint,a_old,scale); + return shuffle_G(aa0, empty, empty, pendint, a_old, scale, gsyms); } if (a1.size()==0) { Gparameter empty; Gparameter aa0 = a0; aa0.insert(aa0.end(),a2.begin(),a2.end()); - return shuffle_G(aa0,empty,empty,pendint,a_old,scale); + return shuffle_G(aa0, empty, empty, pendint, a_old, scale, gsyms); } Gparameter a1_removed(a1.begin()+1,a1.end()); @@ -967,8 +986,8 @@ ex shuffle_G(const Gparameter & a0, const Gparameter & a1, const Gparameter & a2 a01.push_back( a1[0] ); a02.push_back( a2[0] ); - return shuffle_G(a01,a1_removed,a2,pendint,a_old,scale) - + shuffle_G(a02,a1,a2_removed,pendint,a_old,scale); + return shuffle_G(a01, a1_removed, a2, pendint, a_old, scale, gsyms) + + shuffle_G(a02, a1, a2_removed, pendint, a_old, scale, gsyms); } @@ -1067,7 +1086,8 @@ ex G_numeric(const lst& x, const lst& s, const ex& y) // generate missing dummy-symbols int i = 1; - gsyms.clear(); + // holding dummy-symbols for the G/Li transformations + exvector gsyms; gsyms.push_back(symbol("GSYMS_ERROR")); ex lastentry; for (std::multimap<ex,int>::const_iterator it = sortmap.begin(); it != sortmap.end(); ++it) { @@ -1117,7 +1137,7 @@ ex G_numeric(const lst& x, const lst& s, const ex& y) // do transformation Gparameter pendint; - ex result = G_transform(pendint, a, scale); + ex result = G_transform(pendint, a, scale, gsyms); // replace dummy symbols with their values result = result.eval().expand(); result = result.subs(subslst).evalf(); -- 1.5.6 Best regards, Alexei -- All science is either physics or stamp collecting.

1 0

[SCM] GiNaC -- a C++ library for symbolic computations branch, master, updated. release_1-4-0-64-g1222eac
by git＠ginac.de 22 Aug '08

22 Aug '08

This is an automated email from the git hooks/post-receive script. It was generated because a ref change was pushed to the repository containing the project "GiNaC -- a C++ library for symbolic computations". The branch, master has been updated via 1222eac51cee964961d2aad889dc4ceccb144a36 (commit) from 321d5c941f81b1d0b86de91f3f17c5af5b6e4642 (commit) Those revisions listed above that are new to this repository have not appeared on any other notification email; so we list those revisions in full, below. - Log ----------------------------------------------------------------- commit 1222eac51cee964961d2aad889dc4ceccb144a36 Author: Alexei Sheplyakov <varg(a)theor.jinr.ru> Date: Tue Aug 19 23:48:03 2008 +0400 Faster, better (recursive descent) expression parser. Ouch, it fails to compile. Here is a correct version: From: Alexei Sheplyakov <varg(a)theor.jinr.ru> Date: Tue, 19 Aug 2008 21:50:03 +0400 Subject: [PATCH] Faster, better (recursive descent) expression parser. bison generated parser has a number of problems: 1. Bad performance. Parsing a sum seems to be O(N^{2 + a}) (a > 0) [1]. For example, parsing a sum (actually, a univariate polynomial) of 32768 terms takes about 90 sec. on my box, parsing a sum of 10^6 terms takes "eternity". 2. The user is expected to provide list of all symbols in the expression. Often this is very annoying (and useless). 3. Parser is not reentrant (bison *can* produce reentrant parsers, but that won't solve other problems). 4. Parser is difficult to extend. Hence the new parser. Features: 1. Parsing large sums and products is O(N). 2. Parser is reentrant (well, almost). 3. It's possible to insert (shell style) comments inside the expressions. 4. matrices, lists, FAIL are NOT handled. Yes, this *is* a feature :-) Limitations: 1. Error handling is a bit terse: on error exception is thrown, and that's it. 2. Binary, octal, and hexadecimal numbers can not be parsed (yet). 3. Tensors, noncommutative products, etc. can not be parsed. Other notes: 1. ex ctor still uses the old parser. 2. ginsh still uses the old parser. [1] Mesured by this script (requires gnuplot): make_expression_string () { printf "1 + x" local n=2 while test $n -le $1; do printf " + %s*x^%s" $n $n n=$(expr $n + 1) done } single_test () { printf "$1 \t" ( printf "time("; make_expression_string $1; printf ");" ) | \ ginsh | sed -e 's/s$//' } benchmark () { local N=$1 while test $N -le $2; do single_test $N N=$(expr $N \* 2) done } gnuplot_header () { echo "set logscale xy" echo "set xlabel 'degree (# of terms)'" echo "set ylabel 'time, sec.'" echo "set xrange [${1}:${2}]" echo "plot '-' using 1:2 with lines title '1+x+2*x^2+...+n*x^n'" } gnuplot_footer () { echo "e" } benchmark_and_plot () { ( gnuplot_header $1 $2 benchmark $1 $2 | tee $3.txt gnuplot_footer ) | \ gnuplot -persist '-' } N_ini=${1:-1024} N_fin=${2:-32768} out_base=${3:-parser_benchmark} benchmark_and_plot $N_ini $N_fin $out_base ----------------------------------------------------------------------- Summary of changes: INSTALL | 3 +- check/Makefile.am | 8 ++- check/time_parser.cpp | 79 +++++++++++++++++ configure.ac | 4 + ginac/Makefile.am | 31 ++++++- ginac/parser/builtin_fcns.def | 90 +++++++++++++++++++ ginac/parser/debug.hpp | 35 ++++++++ ginac/parser/default_reader.tpl | 50 +++++++++++ ginac/parser/lexer.cpp | 132 ++++++++++++++++++++++++++++ ginac/parser/lexer.hpp | 45 ++++++++++ ginac/parser/parse_binop_rhs.cpp | 176 ++++++++++++++++++++++++++++++++++++++ ginac/parser/parse_context.cpp | 27 ++++++ ginac/parser/parse_context.hpp | 78 +++++++++++++++++ ginac/parser/parser.cpp | 173 +++++++++++++++++++++++++++++++++++++ ginac/parser/parser.hpp | 95 ++++++++++++++++++++ 15 files changed, 1021 insertions(+), 5 deletions(-) create mode 100644 check/time_parser.cpp create mode 100644 ginac/parser/builtin_fcns.def create mode 100644 ginac/parser/debug.hpp create mode 100644 ginac/parser/default_reader.tpl create mode 100644 ginac/parser/lexer.cpp create mode 100644 ginac/parser/lexer.hpp create mode 100644 ginac/parser/parse_binop_rhs.cpp create mode 100644 ginac/parser/parse_context.cpp create mode 100644 ginac/parser/parse_context.hpp create mode 100644 ginac/parser/parser.cpp create mode 100644 ginac/parser/parser.hpp hooks/post-receive -- GiNaC -- a C++ library for symbolic computations

3 3

[SCM] GiNaC -- a C++ library for symbolic computations branch, master, updated. release_1-4-0-63-g321d5c9
by git＠ginac.de 21 Aug '08

21 Aug '08

This is an automated email from the git hooks/post-receive script. It was generated because a ref change was pushed to the repository containing the project "GiNaC -- a C++ library for symbolic computations". The branch, master has been updated via 321d5c941f81b1d0b86de91f3f17c5af5b6e4642 (commit) from 37c418c602485be9b023378a1b60ebdb4a592f2a (commit) Those revisions listed above that are new to this repository have not appeared on any other notification email; so we list those revisions in full, below. - Log ----------------------------------------------------------------- commit 321d5c941f81b1d0b86de91f3f17c5af5b6e4642 Author: Jens Vollinga <jensv(a)nikhef.nl> Date: Thu Aug 21 22:37:09 2008 +0200 Daily bugfix in the polynomial factorization (code didn't catch polynomial "x" and cln::sqrt conversion was buggy). ----------------------------------------------------------------------- Summary of changes: check/exam_factor.cpp | 7 ++++++- ginac/factor.cpp | 5 ++++- 2 files changed, 10 insertions(+), 2 deletions(-) hooks/post-receive -- GiNaC -- a C++ library for symbolic computations

1 0

[PATCH] Faster, better (recursive descent) expression parser.
by Alexei Sheplyakov 21 Aug '08

21 Aug '08

bison generated parser has a number of problems: 1. Bad performance. Parsing a sum seems to be O(N^{2 + a}) (a > 0) [1]. For example, parsing a sum (actually, a univariate polynomial) of 32768 terms takes about 90 sec. on my box, parsing a sum of 10^6 terms takes "eternity". 2. The user is expected to provide list of all symbols in the expression. Often this is very annoying (and useless). 3. Parser is not reentrant (bison *can* produce reentrant parsers, but that won't solve other problems). 4. Parser is difficult to extend. Hence the new parser. Features: 1. Parsing large sums and products is O(N). 2. Parser is reentrant (well, almost). 3. It's possible to insert (shell style) comments inside the expressions. 4. matrices, lists, FAIL are NOT handled. Yes, this *is* a feature :-) Limitations: 1. Error handling is a bit terse: on error exception is thrown, and that's it. 2. Binary, octal, and hexadecimal numbers can not be parsed (yet). 3. Tensors, noncommutative products, etc. can not be parsed. Other notes: 1. ex ctor still uses the old parser. 2. ginsh still uses the old parser. [1] Mesured by this script (requires gnuplot): make_expression_string () { printf "1 + x" local n=2 while test $n -le $1; do printf " + %s*x^%s" $n $n n=$(expr $n + 1) done } single_test () { printf "$1 \t" ( printf "time("; make_expression_string $1; printf ");" ) | \ ginsh | sed -e 's/s$//' } benchmark () { local N=$1 while test $N -le $2; do single_test $N N=$(expr $N \* 2) done } gnuplot_header () { echo "set logscale xy" echo "set xlabel 'degree (# of terms)'" echo "set ylabel 'time, sec.'" echo "set xrange [${1}:${2}]" echo "plot '-' using 1:2 with lines title '1+x+2*x^2+...+n*x^n'" } gnuplot_footer () { echo "e" } benchmark_and_plot () { ( gnuplot_header $1 $2 benchmark $1 $2 | tee $3.txt gnuplot_footer ) | \ gnuplot -persist '-' } N_ini=${1:-1024} N_fin=${2:-32768} out_base=${3:-parser_benchmark} benchmark_and_plot $N_ini $N_fin $out_base --- INSTALL | 3 +- check/Makefile.am | 8 ++- check/time_parser.cpp | 79 +++++++++++++++++ configure.ac | 4 + ginac/Makefile.am | 31 ++++++- ginac/parser/builtin_fcns.def | 90 +++++++++++++++++++ ginac/parser/debug.hpp | 35 ++++++++ ginac/parser/default_reader.tpl | 47 ++++++++++ ginac/parser/lexer.cpp | 132 ++++++++++++++++++++++++++++ ginac/parser/lexer.hpp | 45 ++++++++++ ginac/parser/parse_binop_rhs.cpp | 179 ++++++++++++++++++++++++++++++++++++++ ginac/parser/parse_context.cpp | 27 ++++++ ginac/parser/parse_context.hpp | 82 +++++++++++++++++ ginac/parser/parser.cpp | 172 ++++++++++++++++++++++++++++++++++++ ginac/parser/parser.hpp | 99 +++++++++++++++++++++ 15 files changed, 1028 insertions(+), 5 deletions(-) create mode 100644 check/time_parser.cpp create mode 100644 ginac/parser/builtin_fcns.def create mode 100644 ginac/parser/debug.hpp create mode 100644 ginac/parser/default_reader.tpl create mode 100644 ginac/parser/lexer.cpp create mode 100644 ginac/parser/lexer.hpp create mode 100644 ginac/parser/parse_binop_rhs.cpp create mode 100644 ginac/parser/parse_context.cpp create mode 100644 ginac/parser/parse_context.hpp create mode 100644 ginac/parser/parser.cpp create mode 100644 ginac/parser/parser.hpp diff --git a/INSTALL b/INSTALL index b542d9e..1cad277 100644 --- a/INSTALL +++ b/INSTALL @@ -30,7 +30,8 @@ Known not to work with: is missing there. If you install from CVS, you also need GNU autoconf (>=2.59), automake (>=1.7), -libtool (>= 1.5), bison (>= 2.3), flex (>= 2.5.33) to be installed. +libtool (>= 1.5), bison (>= 2.3), flex (>= 2.5.33), autogen (>= 5.6.0) to be +installed. INSTALLATION diff --git a/check/Makefile.am b/check/Makefile.am index a34bb9f..b656c63 100644 --- a/check/Makefile.am +++ b/check/Makefile.am @@ -50,7 +50,8 @@ TIMES = time_dennyfliegner \ time_lw_Q \ time_lw_Qprime \ time_antipode \ - time_fateman_expand + time_fateman_expand \ + time_parser TESTS = $(CHECKS) $(EXAMS) $(TIMES) check_PROGRAMS = $(CHECKS) $(EXAMS) $(TIMES) @@ -227,6 +228,11 @@ time_fateman_expand_SOURCES = time_fateman_expand.cpp \ randomize_serials.cpp timer.cpp timer.h time_fateman_expand_LDADD = ../ginac/libginac.la +time_parser_SOURCES = time_parser.cpp \ + randomize_serials.cpp timer.cpp timer.h +time_parser_LDADD = ../ginac/libginac.la + + AM_CPPFLAGS = -I$(srcdir)/../ginac -I../ginac CLEANFILES = exam.gar diff --git a/check/time_parser.cpp b/check/time_parser.cpp new file mode 100644 index 0000000..ca35fb6 --- /dev/null +++ b/check/time_parser.cpp @@ -0,0 +1,79 @@ +#include <iostream> +#include <string> +#include <sstream> +#include <cstddef> +#include <cstdlib> +#include <vector> +#include <stdexcept> +#include "ginac.h" +#include "parser/parser.hpp" +#include "timer.h" +extern void randomify_symbol_serials(); +using namespace std; +using namespace GiNaC; + +/// make a string "1+x+2*x^2+...+n*x^n" +static string prepare_str(const unsigned n, const char x = 'x') +{ + ostringstream s; + s << x; + for (unsigned i = 2; i < n; i++) + s << '+' << i << '*' << x << '^' << i; + return s.str(); +} + +void benchmark_and_cmp(const string& srep, double& t_new, double& t_old) +{ + parser the_parser; + timer RSD10; + RSD10.start(); + ex e = the_parser(srep); + t_new = RSD10.read(); + RSD10.stop(); + if (t_new > 2.0) + cout << '.' << flush; + + symtab syms = the_parser.get_syms(); + const symbol x = find_or_insert_symbol("x", syms, true); + lst sl; + sl = x; + RSD10.start(); + ex e2(srep, sl); + t_old = RSD10.read(); + + if (t_old > 2.0) + cout << '.' << flush; + + ex dif = (e - e2).expand(); + if (!dif.is_zero()) { + cerr << "Got a difference: " << dif << endl; + throw std::logic_error("New and old parser give different results"); + } +} + +int main(int argc, char** argv) +{ + cout << "timing GiNaC parser..." << flush; + randomify_symbol_serials(); + unsigned n_min = 1024; + unsigned n_max = 32768; + if (argc > 1) + n_max = atoi(argv[1]); + + vector<double> times_new, times_old; + vector<unsigned> ns; + for (unsigned n = n_min; n <= n_max; n = n << 1) { + double t_new, t_old; + string srep = prepare_str(n); + benchmark_and_cmp(srep, t_new, t_old); + times_new.push_back(t_new); + times_old.push_back(t_old); + ns.push_back(n); + } + + cout << "OK" << endl; + cout << "# terms new parser, s old parser, s" << endl; + for (size_t i = 0; i < times_new.size(); i++) + cout << " " << ns[i] << '\t' << times_new[i] << '\t' << times_old[i] << endl; + return 0; +} diff --git a/configure.ac b/configure.ac index e4a1e60..22f9c96 100644 --- a/configure.ac +++ b/configure.ac @@ -149,6 +149,10 @@ AM_CONDITIONAL(CONFIG_TEX, [test ! \( -z "$LATEX" -o -z $"PDFLATEX" -o -z "$MAKE AC_PATH_PROG(FIG2DEV, fig2dev, "") AM_CONDITIONAL(CONFIG_FIG2DEV, [test ! -z "$FIG2DEV"]) +dnl generate boilerplate code for the (new) parser. +dnl Only developers need this tool. +AC_PATH_PROG(AUTOGEN, autogen, "") + dnl Output makefiles etc. AC_CONFIG_FILES([ Makefile diff --git a/ginac/Makefile.am b/ginac/Makefile.am index d2ef614..0570fde 100644 --- a/ginac/Makefile.am +++ b/ginac/Makefile.am @@ -9,7 +9,15 @@ libginac_la_SOURCES = add.cpp archive.cpp basic.cpp clifford.cpp color.cpp \ operators.cpp power.cpp registrar.cpp relational.cpp remember.cpp \ pseries.cpp print.cpp symbol.cpp symmetry.cpp tensor.cpp \ utils.cpp wildcard.cpp input_parser.yy input_lexer.ll \ - input_lexer.h remember.h tostring.h utils.h compiler.h + input_lexer.h remember.h tostring.h utils.h compiler.h \ + parser/parse_binop_rhs.cpp \ + parser/parser.cpp \ + parser/parse_context.cpp \ + parser/builtin_fcns.cpp \ + parser/lexer.cpp \ + parser/lexer.hpp \ + parser/debug.hpp + libginac_la_LDFLAGS = -version-info $(LT_VERSION_INFO) -release $(LT_RELEASE) libginac_la_LIBADD = $(DL_LIBS) ginacincludedir = $(includedir)/ginac @@ -18,10 +26,27 @@ ginacinclude_HEADERS = ginac.h add.h archive.h assertion.h basic.h class_info.h exprseq.h fail.h factor.h fderivative.h flags.h function.h hash_map.h idx.h indexed.h \ inifcns.h integral.h lst.h matrix.h mul.h ncmul.h normal.h numeric.h operators.h \ power.h print.h pseries.h ptr.h registrar.h relational.h structure.h \ - symbol.h symmetry.h tensor.h version.h wildcard.h + symbol.h symmetry.h tensor.h version.h wildcard.h \ + parser/parser.hpp \ + parser/parse_context.hpp + AM_LFLAGS = -Pginac_yy -olex.yy.c AM_YFLAGS = -p ginac_yy -d -EXTRA_DIST = function.pl input_parser.h version.h.in +EXTRA_DIST = function.pl input_parser.h version.h.in \ +parser/default_reader.tpl parser/builtin_fcns.def + +# Files produced by autogen(1) from templates +$(srcdir)/parser/builtin_fcns.cpp: $(srcdir)/parser/builtin_fcns.def $(srcdir)/parser/default_reader.tpl + set -e; if [ -n "$(AUTOGEN)" ]; then \ + cd $(srcdir)/parser; \ + $(AUTOGEN) -T default_reader.tpl builtin_fcns.def; \ + elif [ -f $@ ]; then \ + echo "WARNING: AutoGen is not available, the \"$@\" file WON'T be re-generated"; \ + else \ + echo "*** ERROR: the \"$@\" file does not exist, and AutoGen is not installed on your system"; \ + echo "*** Please install AutoGen (http://www.gnu.org/software/autogen)"; \ + fi + # Files which are generated by perl scripts $(srcdir)/function.h $(srcdir)/function.cpp: $(srcdir)/function.pl diff --git a/ginac/parser/builtin_fcns.def b/ginac/parser/builtin_fcns.def new file mode 100644 index 0000000..897b1f7 --- /dev/null +++ b/ginac/parser/builtin_fcns.def @@ -0,0 +1,90 @@ +Autogen definitions ginacfcns; + +function = { name = "log"; }; +function = { name = "exp"; }; +function = { name = "sin"; }; +function = { name = "cos"; }; +function = { name = "tan"; }; +function = { name = "asin"; }; +function = { name = "acos"; }; +function = { name = "atan"; }; + +function = { name = "sinh"; }; +function = { name = "cosh"; }; +function = { name = "tanh"; }; +function = { name = "asinh"; }; +function = { name = "acosh"; }; +function = { name = "atanh"; }; + +function = { + name = "atan2"; + args = 2; +}; + +function = { + name = "Li2"; + comment = "Dilogarithm"; +}; + +function = { + name = "Li3"; + comment = "Trilogarithm"; +}; + +function = { + name = "zetaderiv"; + comment = "Derivatives of Riemann's Zeta-function"; + args = 2; +}; + +function = { + name = "Li"; + args = 2; + comment = "Polylogarithm and multiple polylogarithm"; +}; + +function = { + name = "S"; + args = 3; + comment = "Nielsen's generalized polylogarithm"; +}; + +function = { + name = "H"; + args = 2; + comment = "Harmonic polylogarithm"; +}; + +function = { name = "lgamma"; }; +function = { name = "tgamma"; }; + +function = { + name = "beta"; + args = 2; + comment = "Beta-function"; +}; + +function = { name = "factorial"; }; + +function = { + name = "binomial"; + args = 2; +}; + +function = { + name = "Order"; + comment = "Order term function (for truncated power series)"; +}; + +/* Thease are not functions, but anyway ... */ +function = { name = "sqrt"; }; + +function = { + name = "pow"; + args = 2; +}; + +function = { + name = "power"; + args = 2; +}; diff --git a/ginac/parser/debug.hpp b/ginac/parser/debug.hpp new file mode 100644 index 0000000..a43cc7a --- /dev/null +++ b/ginac/parser/debug.hpp @@ -0,0 +1,35 @@ +#ifndef GINAC_PARSER_DEBUG_HPP +#define GINAC_PARSER_DEBUG_HPP +#include <iosfwd> +#include <sstream> +#include <stdexcept> +#include "compiler.h" +#ifndef __GNUC__ +#if __STDC_VERSION__ < 199901L +#define __PRETTY_FUNCTION__ "<unknown>" +#else +#define __PRETTY_FUNCTION__ __func__ +#endif +#endif + +#define bail_out(exception, message) \ +do { \ + std::ostringstream err; \ + err << __PRETTY_FUNCTION__ << "(" << __FILE__ << ':' << __LINE__ << ": "; \ + err << message; \ + throw exception(err.str()); \ +} while (0) + +#define bug(message) bail_out(std::logic_error, message) + +#define dout(condition, message) \ +do { \ + if (unlikely(condition)) { \ + std::cerr << __PRETTY_FUNCTION__ \ + << " (" << __FILE__ << ':' << __LINE__ << "): " \ + << message << std::endl; \ + } \ +} while (0) + +#endif // GINAC_PARSER_DEBUG_HPP + diff --git a/ginac/parser/default_reader.tpl b/ginac/parser/default_reader.tpl new file mode 100644 index 0000000..f3b150d --- /dev/null +++ b/ginac/parser/default_reader.tpl @@ -0,0 +1,47 @@ +[+ AutoGen5 template .cpp +][+ +COMMENT a part of GiNaC parser -- construct functions from a byte stream. ++][+ +(use-modules (ice-9 format)) + +(define (sequence start end . step) + (let ((step (if (null? step) 1 (car step)))) + (let loop ((n start)) + (if (> n end) '() (cons n (loop (+ step n))))))) ++]/* +[+ (dne " * " " * " ) +] + * + * If you want to change this file, edit either `[+ (def-file) +]' or + * `[+ (tpl-file) +]' file, and run the following command: + * + * autogen -T [+ (tpl-file) +] [+ (def-file) +] + */ +#include "parse_context.hpp" + +namespace GiNaC +{ +[+ FOR function +] +static ex [+ (get "name") +]_reader(const exvector& ev) +{ + return GiNaC::[+ (get "name") +]([+ + (let ((nargs (if (exist? "args") + (string->number (get "args")) 1))) + (format '#f "~{ev[~a]~^, ~}" (sequence 0 (- nargs 1)))) +]); +}[+ ENDFOR +] + +const prototype_table& get_default_reader() +{ + using std::make_pair; + static bool initialized = false; + static prototype_table reader; + if (!initialized) { +[+ FOR function +] + reader[make_pair("[+ (get "name") +]", [+ + (if (exist? "args") (get "args") "1") + +])] = [+ (get "name") +]_reader;[+ + ENDFOR +] + initialized = true; + } + return reader; +} +} // namespace GiNaC + diff --git a/ginac/parser/lexer.cpp b/ginac/parser/lexer.cpp new file mode 100644 index 0000000..ed1f894 --- /dev/null +++ b/ginac/parser/lexer.cpp @@ -0,0 +1,132 @@ +#include <string> +#include <iostream> +#include <sstream> +#include "lexer.hpp" +#include "compiler.h" + +namespace GiNaC +{ +/// Skip to the end of line +static int skipline(std::istream* s); +/// Skip to the next non-whitespace character +static int skipspace(std::istream* s, int c, std::size_t& line); +/// Check if the identifier is predefined literal +static bool literal_p(const std::string& name); + +/// gettok - Return the next token from standard input. +int lexer::gettok() +{ + // Skip any whitespace. + c = skipspace(input, c, line_num); + + // identifier: [a-zA-Z][a-zA-Z0-9]* + if (isalpha(c)) { + str = c; + do { + c = input->get(); + if (isalnum(c)) + str += c; + else + break; + } while (true); + if (unlikely(literal_p(str))) + return token_type::literal; + else + return token_type::identifier; + } + + // Number: [0-9.]+ + if (isdigit(c) || c == '.') { + str = ""; + do { + str += c; + c = input->get(); + } while (isdigit(c) || c == '.'); + return token_type::number; + } + + // Comment until end of line. + if (c == '#') { + c = skipline(input); + ++line_num; + if (c != EOF) + return gettok(); + } + + // Check for end of file. Don't eat the EOF. + if (c == EOF) + return token_type::eof; + + // Otherwise, just return the character as its ascii value. + int current = c; + c = input->get(); + return current; +} + +static int skipline(std::istream* s) +{ + int c; + do { + c = s->get(); + } while (c != EOF && c != '\n' && c != '\r'); + return c; +} + +static int skipspace(std::istream* s, int c, std::size_t& line) +{ + while (isspace(c)) { + if (c == '\n') + ++line; + c = s->get(); + } + return c; +} + +static bool literal_p(const std::string& name) +{ + if (name == "I") + return true; + else if (name == "Pi") + return true; + else if (name == "Euler") + return true; + else if (name == "Catalan") + return true; + else + return false; +} + +lexer::lexer(std::istream* in, std::ostream* out, std::ostream* err) +{ + if (in) + input = in; + else + in = &std::cin; + + if (out) + output = out; + else + output = &std::cout; + + if (err) + error = err; + else + error = &std::cerr; + + c = ' '; + str = ""; + line_num = 0; + column = 0; +} + +lexer::~lexer() { } + +void lexer::switch_input(std::istream* in) +{ + input = in; + line_num = 0; + column = 0; +} + +} // namespace GiNaC + diff --git a/ginac/parser/lexer.hpp b/ginac/parser/lexer.hpp new file mode 100644 index 0000000..b53d08f --- /dev/null +++ b/ginac/parser/lexer.hpp @@ -0,0 +1,45 @@ +#ifndef GINAC_LEXER_HPP_ +#define GINAC_LEXER_HPP_ +#include <iosfwd> +#include <string> +#include <cstddef> +namespace GiNaC +{ + +class lexer +{ + std::istream* input; + std::ostream* output; + std::ostream* error; + /// last character read from stream + int c; + /// identifier and number tokens are stored here + std::string str; + std::size_t line_num; + std::size_t column; + friend class parser; +public: + + lexer(std::istream* in = 0, std::ostream* out = 0, std::ostream* err = 0); + ~lexer(); + + int gettok(); + void switch_input(std::istream* in); + + struct token_type + { + enum + { + eof = -1, + identifier = -4, + number = -5, + literal = -6 + + }; + }; +}; + +} // namespace GiNaC + +#endif // GINAC_LEXER_HPP_ + diff --git a/ginac/parser/parse_binop_rhs.cpp b/ginac/parser/parse_binop_rhs.cpp new file mode 100644 index 0000000..b118197 --- /dev/null +++ b/ginac/parser/parse_binop_rhs.cpp @@ -0,0 +1,179 @@ +#ifndef IN_GINAC +#include <ginac/ginac.h> +#else +#include "ex.h" +#include "symbol.h" +#include "mul.h" +#include "add.h" +#include "power.h" +#endif +#include <stdexcept> +#include <sstream> +#include "parser.hpp" +#include "lexer.hpp" +#include "debug.hpp" + +namespace GiNaC +{ + +/// Make a sum or a product. +static ex make_binop_expr(const int binop, const exvector& args); +/// Check if the token is a binary operator. +static inline bool is_binop(const int c); +/// Get the precedence of the pending binary operator. +static int get_tok_prec(const int c); + +/// binoprhs: ([+*/^-] primary)* +ex parser::parse_binop_rhs(int expr_prec, ex& lhs) +{ + exvector args; + args.push_back(lhs); + int binop = -1, orig_binop = -1; + bool need_sign_flip = false; + while (1) { + // check if this is a binop + if (!is_binop(token)) { + if (args.size() > 1) + return make_binop_expr(orig_binop, args); + else + return lhs; + } + + // Okay, we know this is a binop. + if (args.size() == 1) + orig_binop = token; + + binop = token; + + // If this is a binop that binds at least as tightly as + // the current binop, consume it, otherwise we are done. + int tok_prec = get_tok_prec(token); + if (tok_prec < expr_prec) { + if (args.size() > 1) + return make_binop_expr(orig_binop, args); + else + return lhs; + } + + get_next_tok(); // eat binop + + // Parse the primary expression after the binary operator. + ex rhs = parse_primary(); + + // If binop binds less tightly with rhs than the operator after + // rhs, let the pending operator take rhs as its lhs. + int next_prec = get_tok_prec(token); + if (tok_prec < next_prec) + rhs = parse_binop_rhs(tok_prec + 1, rhs); + + // previous operator was '+', and current one is '-' + // (or vice a versa). + if (need_sign_flip) + rhs = - rhs; + + args.push_back(rhs); + + // Minimize the number of eval() and ctor calls. This is + // crucial for a reasonable performance. If the next operator + // is compatible with the pending one (or the same) don't create + // the expression and continue collecting operands instead. + if (binop == token) + continue; + else if (binop == '+' && token == '-') { + need_sign_flip = token != orig_binop; + continue; + } else if (binop == '-' && token == '+') { + need_sign_flip = token != orig_binop; + continue; + } else { + if (args.size() <= 1) + bug("binop has " << args.size() << " arguments, expected >= 2"); + lhs = make_binop_expr(orig_binop, args); + args.clear(); + args.push_back(lhs); + } + } +} + +extern numeric* _num_1_p; + +static ex make_minus_expr(const exvector& args) +{ + exvector rest_args; + rest_args.reserve(args.size() - 1); + std::copy(args.begin() + 1, args.end(), std::back_inserter(rest_args)); + ex rest_base = (new add(rest_args))->setflag(status_flags::dynallocated); + ex rest = (new mul(rest_base, *_num_1_p))->setflag(status_flags::dynallocated); + ex ret = (new add(args[0], rest))->setflag(status_flags::dynallocated); + return ret; +} + +static ex make_divide_expr(const exvector& args) +{ + exvector rest_args; + rest_args.reserve(args.size() - 1); + std::copy(args.begin() + 1, args.end(), std::back_inserter(rest_args)); + ex rest_base = (new mul(rest_args))->setflag(status_flags::dynallocated); + ex rest = pow(rest_base, *_num_1_p); + return (new mul(args[0], rest))->setflag(status_flags::dynallocated); +} + +static ex make_binop_expr(const int binop, const exvector& args) +{ + switch (binop) { + case '+': + return (new add(args))->setflag(status_flags::dynallocated); + case '-': + return make_minus_expr(args); + case '*': + return (new mul(args))->setflag(status_flags::dynallocated); + case '/': + return make_divide_expr(args); + case '^': + if (args.size() != 2) + throw std::invalid_argument( + std::string(__func__) + + ": power should have exactly 2 operands"); + return pow(args[0], args[1]); + default: + throw std::invalid_argument( + std::string(__func__) + + ": invalid binary operation: " + + char(binop)); + } +} + +static inline bool is_binop(const int c) +{ + switch (c) { + case '+': + case '-': + case '*': + case '/': + case '^': + return true; + default: + return false; + } +} + +/// Get the precedence of the pending binary operator. +static int get_tok_prec(const int c) +{ + switch (c) { + case '+': + case '-': + return 20; + case '*': + case '/': + return 40; + case '^': + return 60; + default: + return -1; + // means 'this is not a binary operator' + } +} + +} // namespace GiNaC + diff --git a/ginac/parser/parse_context.cpp b/ginac/parser/parse_context.cpp new file mode 100644 index 0000000..e4956ba --- /dev/null +++ b/ginac/parser/parse_context.cpp @@ -0,0 +1,27 @@ +#include "parse_context.hpp" +#include <sstream> +#include <stdexcept> +namespace GiNaC +{ + +const symbol& +find_or_insert_symbol(const std::string& name, symtab& syms, const bool strict) +{ + symtab::const_iterator p = syms.find(name); + if (p != syms.end()) + return p->second.first; + + if (strict) + throw std::invalid_argument( + std::string("find_or_insert_symbol: symbol \"") + + name + "\" not found"); + + // false means this symbol was created by parser + const std::pair<symbol, bool> tmp = std::make_pair(symbol(name), false); + + symtab::iterator i = syms.insert(symtab::value_type(name, tmp)).first; + return i->second.first; +} + +} + diff --git a/ginac/parser/parse_context.hpp b/ginac/parser/parse_context.hpp new file mode 100644 index 0000000..c7cc19b --- /dev/null +++ b/ginac/parser/parse_context.hpp @@ -0,0 +1,82 @@ +#ifndef _GINAC_PARSE_CONTEXT_HPP +#define _GINAC_PARSE_CONTEXT_HPP +#include <string> +#include <cstddef> // size_t +#ifndef IN_GINAC +#include <ginac/ginac.h> +#else +#include "ex.h" +#include "symbol.h" +#endif +#include <map> +#include <utility> + +namespace GiNaC +{ + +/** + * Establishes correspondence between the strings and symbols. + * The parser will create missing symbols (if not instructed otherwise, + * in which case it fails if the expression contains unknown symbols). + * The .second element of pair helps to distinguish between the user + * supplied symbols and parser generated ones. The .first is the symbol + * itself + */ +typedef std::map<std::string, std::pair<symbol, bool> > symtab; + +/** + * Find the symbol with the @a name in the symbol table @a syms. + * + * If symbol is missing and @a strict = false, insert it, otherwise + * throw an exception. + */ +extern const symbol& +find_or_insert_symbol(const std::string& name, symtab& syms, + const bool strict); + +/** + * Function (or class ctor) prototype + * .first is the name of function(or ctor), + * .second is the number of arguments (each of type ex) + */ +typedef std::pair<std::string, std::size_t> prototype; + +/** + * A (C++) function for reading functions and classes from the stream. + * + * The parser uses (an associative array of) such functions to construct + * (GiNaC) classes and functions from a sequence of characters. + */ +typedef ex (*reader_func)(const exvector& args); + +/** + * Prototype table. + * + * If parser sees an expression which looks like a function call (e.g. + * foo(x+y, z^2, t)), it looks up such a table to find out which + * function (or class) corresponds to the given name and has the given + * number of the arguments. + * + * N.B. + * + * 1. The function don't have to return a (GiNaC) function or class, it + * can return any expression. + * 2. Overloaded functions/ctors are paritally supported, i.e. there might + * be several functions with the same name, but they should take different + * number of arguments. + * 3. User can extend the parser via custom prototype tables. It's possible + * to read user defined classes, create abbreviations, etc. + */ +typedef std::map<prototype, reader_func> prototype_table; + +/** + * Default prototype table. + * + * It supports most of builtin GiNaC functions. + */ +extern const prototype_table& get_default_reader(); + +} + +#endif // _GINAC_PARSE_CONTEXT_HPP + diff --git a/ginac/parser/parser.cpp b/ginac/parser/parser.cpp new file mode 100644 index 0000000..bd44137 --- /dev/null +++ b/ginac/parser/parser.cpp @@ -0,0 +1,172 @@ +#include <ginac/ginac.h> +#include <stdexcept> +#include <sstream> +#include "parser.hpp" +#include "lexer.hpp" +#include "debug.hpp" + +namespace GiNaC +{ + +/// identifier_expr: identifier | identifier '(' expression* ')' +ex parser::parse_identifier_expr() +{ + std::string name = scanner->str; + get_next_tok(); // eat identifier. + + if (token != '(') // symbol + return find_or_insert_symbol(name, syms, strict); + + // function/ctor call. + get_next_tok(); // eat ( + exvector args; + if (token != ')') { + while (true) { + ex e = parse_expression(); + args.push_back(e); + + if (token == ')') + break; + + if (token != ',') + throw std::invalid_argument("Expected ')' or ',' in argument list"); + + get_next_tok(); + } + } + // Eat the ')'. + get_next_tok(); + prototype the_prototype = make_pair(name, args.size()); + prototype_table::const_iterator reader = funcs.find(the_prototype); + if (reader == funcs.end()) { + bail_out(std::invalid_argument, + "no function \"" << name << "\" with " << args.size() + << " arguments"); + } + ex ret = reader->second(args); + return ret; +} + +/// paren_expr: '(' expression ')' +ex parser::parse_paren_expr() +{ + get_next_tok(); // eat (. + ex e = parse_expression(); + + if (token != ')') + throw std::invalid_argument("parser::parse_paren_expr: expected ')'"); + get_next_tok(); // eat ). + return e; +} + +extern numeric* _num_1_p; + +/// unary_expr: [+-] expression +ex parser::parse_unary_expr(const int s) +{ + // consume '-' (or '+') + get_next_tok(); + ex v = parse_expression(); + switch (s) { + case '-': + return (new mul(v, *_num_1_p))->setflag(status_flags::dynallocated); + case '+': + return v; + default: + throw std::invalid_argument( + std::string(__func__) + + ": invalid unary operator \"" + + char(s) + "\""); + } +} + +/// primary: identifier_expr | number_expr | paren_expr | unary_expr +ex parser::parse_primary() +{ + switch (token) { + case lexer::token_type::identifier: + return parse_identifier_expr(); + case lexer::token_type::number: + return parse_number_expr(); + case '(': + return parse_paren_expr(); + case '-': + return parse_unary_expr('-'); + case '+': + return parse_unary_expr('+'); + case lexer::token_type::literal: + return parse_literal_expr(); + case lexer::token_type::eof: + bail_out(std::invalid_argument, "got EOF while parsing the expression"); + default: + bail_out(std::invalid_argument, "unknown token " << + token << " (\"" << + (token ? std::string("") + char(token) : "NULL") + << "\")"); + } +} + +/// expression ::= primary binoprhs +ex parser::parse_expression() +{ + ex lhs = parse_primary(); + ex res = parse_binop_rhs(0, lhs); + return res; +} + +/// number_expr: number +ex parser::parse_number_expr() +{ + ex n = numeric(scanner->str.c_str()); + get_next_tok(); // consume the number + return n; +} + +/// literal_expr: 'I' | 'Pi' | 'Euler' | 'Catalan' +ex parser::parse_literal_expr() +{ + if (scanner->str == "I") + return I; + else if (scanner->str == "Pi") + return Pi; + else if (scanner->str == "Euler") + return Euler; + else if (scanner->str == "Catalan") + return Catalan; + bug("unknown literal: \"" << scanner->str << "\""); +} + +ex parser::operator()(std::istream& input) +{ + scanner->switch_input(&input); + get_next_tok(); + ex ret = parse_expression(); + return ret; +} + +ex parser::operator()(const std::string& input) +{ + std::istringstream is(input); + ex ret = operator()(is); + return ret; +} + +int parser::get_next_tok() +{ + token = scanner->gettok(); + return token; +} + +parser::parser(const symtab& syms_, const prototype_table& funcs_, + const bool strict_) : strict(strict_), funcs(funcs_), + syms(syms_) +{ + scanner = new lexer(); +} + +parser::~parser() +{ + delete scanner; +} + +} // namespace GiNaC diff --git a/ginac/parser/parser.hpp b/ginac/parser/parser.hpp new file mode 100644 index 0000000..40e7e07 --- /dev/null +++ b/ginac/parser/parser.hpp @@ -0,0 +1,99 @@ +#ifndef GINAC_PARSER_HPP_ + +#include "parse_context.hpp" +#include <stdexcept> +#ifndef IN_GINAC +#include <ginac/ginac.h> +#else +#include "ex.h" +#endif + +namespace GiNaC +{ + +class lexer; + +/** + * Recursive descent parser for GiNaC expressions. + */ +class parser +{ + // The actual parser rules (in EBNF-alike notation): + + /// expression: primary binoprhs + ex parse_expression(); + + /// primary: indentifier_expr | number_expr | paren_expr | unary_expr + ex parse_primary(); + + /// binoprhs: ([+*/^-] primary)* + ex parse_binop_rhs(int, ex&); + + /// identifier_expr: identifier | + /// identifier '(' expression (',' expression)* ')' + ex parse_identifier_expr(); + + /// paren_expr: '(' expression ')' + ex parse_paren_expr(); + + /// number_expr: number + ex parse_number_expr(); + + /// unary_expr: [+-] expression + ex parse_unary_expr(const int c); + + /// literal_expr: 'I' | 'Pi' | 'Euler' | 'Catalan' + ex parse_literal_expr(); + +public: + /** + * @param syms_ symbol table. + * @param funcs_ function/ctors table. + * @param strict_ if true, throw an exception if unknown + * symbol is encountered. + */ + parser(const symtab& syms_ = symtab(), + const prototype_table& funcs_ = get_default_reader(), + const bool strict_ = false); + ~parser(); + + /// parse the stream @a input + ex operator()(std::istream& input); + /// parse the string @a input + ex operator()(const std::string& input); + + /// report the symbol table used by parser + symtab get_syms() const + { + return syms; + } + +private: + /// If true, throw an exception if an unknown symbol is encountered. + const bool strict; + /** + * Function/ctor table, maps a prototype (which is a name and number + * arguments) to a C++ function. Used for parsing identifier_expr's + * (see parse_identifier_expr). If expression contains unknown + * prototype, an exception will be thrown. + */ + const prototype_table funcs; + /** + * Symbol (variable) table. Used for parsing identifier_expr's + * (see parse_identifier_expr). If parser is strict, exception is + * thrown if an unknown symbol is encountered. Non-strict parser + * appends unknown symbols to the symbol table. + */ + symtab syms; + /// token scanner + lexer* scanner; + /// current token the parser is looking at + int token; + /// read the next token from the scanner + int get_next_tok(); +}; + +} // namespace GiNaC + +#endif // GINAC_PARSER_HPP_ + -- 1.5.6 Best regards, Alexei -- All science is either physics or stamp collecting.

3 8