@for file in $(DISTFILES); do \
d=$(srcdir); \
if test -d $$d/$$file; then \
- cp -pr $$d/$$file $(distdir)/$$file; \
+ cp -pr $$/$$file $(distdir)/$$file; \
else \
test -f $(distdir)/$$file \
|| ln $$d/$$file $(distdir)/$$file 2> /dev/null \
General Information
===================
-GiNaC (which stands for "GiNaC is Not a CAS (computer algebra system)) is a
-C++ library for symbolic mathematical calculations. It is designed to allow
+GiNaC (which stands for "GiNaC is Not a CAS" (computer algebra system)) is a
+C++ library for symbolic mathematical calculations. It is designed to allow
the creation of integrated systems that embed symbolic manipulations together
with more established areas of computer science (like computation-intense
numeric applications, graphical interfaces, etc.) under one roof.
==================
If you have identified a bug in GiNaC you are welcome to send a detailed
-bug report to <ginac-bugs@ginac.de>. Please think about your bug! This
+bug report to <ginac-bugs@ginac.de>. Please think about your bug! This
means that you should include
* Information about your system
together with the output you get and the output you expect will
help us to reproduce it quickly.
-Patches are most welcome. If possible please make them with diff -c and
+Patches are most welcome. If possible please make them with diff -c and
include ChangeLog entries.
@for file in $(DISTFILES); do \
d=$(srcdir); \
if test -d $$d/$$file; then \
- cp -pr $$d/$$file $(distdir)/$$file; \
+ cp -pr $$/$$file $(distdir)/$$file; \
else \
test -f $(distdir)/$$file \
|| ln $$d/$$file $(distdir)/$$file 2> /dev/null \
result += inifcns_consist_psi();
result += inifcns_consist_zeta();
- if ( !result ) {
+ if (!result) {
cout << " passed ";
clog << "(no output)" << endl;
} else {
matrix m_i = m.inverse();
ex det = m.determinant().expand();
- if ( (normal(m_i(0,0)*det) != d) ||
- (normal(m_i(0,1)*det) != -b) ||
- (normal(m_i(1,0)*det) != -c) ||
- (normal(m_i(1,1)*det) != a) ) {
+ if ((normal(m_i(0,0)*det) != d) ||
+ (normal(m_i(0,1)*det) != -b) ||
+ (normal(m_i(1,0)*det) != -c) ||
+ (normal(m_i(1,1)*det) != a)) {
clog << "inversion of 2x2 matrix " << m
<< " erroneously returned " << m_i << endl;
return 1;
matrix m_i = m.inverse();
ex det = m.determinant().normal().expand();
- if ( (normal(m_i(0,0)*det) != (e*i-f*h)) ||
- (normal(m_i(0,1)*det) != (c*h-b*i)) ||
- (normal(m_i(0,2)*det) != (b*f-c*e)) ||
- (normal(m_i(1,0)*det) != (f*g-d*i)) ||
- (normal(m_i(1,1)*det) != (a*i-c*g)) ||
- (normal(m_i(1,2)*det) != (c*d-a*f)) ||
- (normal(m_i(2,0)*det) != (d*h-e*g)) ||
- (normal(m_i(2,1)*det) != (b*g-a*h)) ||
- (normal(m_i(2,2)*det) != (a*e-b*d)) ) {
+ if ((normal(m_i(0,0)*det) != (e*i-f*h)) ||
+ (normal(m_i(0,1)*det) != (c*h-b*i)) ||
+ (normal(m_i(0,2)*det) != (b*f-c*e)) ||
+ (normal(m_i(1,0)*det) != (f*g-d*i)) ||
+ (normal(m_i(1,1)*det) != (a*i-c*g)) ||
+ (normal(m_i(1,2)*det) != (c*d-a*f)) ||
+ (normal(m_i(2,0)*det) != (d*h-e*g)) ||
+ (normal(m_i(2,1)*det) != (b*g-a*h)) ||
+ (normal(m_i(2,2)*det) != (a*e-b*d))) {
clog << "inversion of 3x3 matrix " << m
<< " erroneously returned " << m_i << endl;
return 1;
// square roots of squares of integers:
passed = true;
for (int i=0; i<42; ++i) {
- if ( !sqrt(numeric(i*i)).is_integer() ) {
+ if (!sqrt(numeric(i*i)).is_integer()) {
passed = false;
}
}
passed = true;
for (int num=0; num<41; ++num) {
for (int den=1; den<42; ++den) {
- if ( !sqrt(numeric(num*num)/numeric(den*den)).is_rational() ) {
+ if (!sqrt(numeric(num*num)/numeric(den*den)).is_rational()) {
passed = false;
}
}
clog << "---------output of numeric types:" << endl;
unsigned long Digits_before = Digits;
- Digits = 200;
+ Digits = 222;
clog << "Using " << Digits << " digits" << endl;
clog << Pi << " evalfs to: " << Pi.evalf() << endl;
clog << Catalan << " evalfs to: " << Catalan.evalf() << endl;
/** @file paranoia_check.cpp
*
* This set of tests checks for some of GiNaC's oopses which showed up during
- * development. Things were evaluated wrongly and so. It should not find such
- * a sick behaviour again. But since we are paranoic and we want to exclude
- * that behaviour for good... */
+ * development. Things were evaluated wrongly and so. Such a sick behaviour
+ * shouldn't occur any more. But we are paranoic and we want to exclude these
+ * these oopses for good, so we run those stupid tests... */
/*
* GiNaC Copyright (C) 1999 Johannes Gutenberg University Mainz, Germany
g = e / f;
// In the first one expand did not do any job at all:
- if ( !g.expand().is_equal(x) ) {
+ if (!g.expand().is_equal(x)) {
clog << "e = x*y*z; f = y*z; expand(e/f) erroneously returned "
<< g.expand() << endl;
++result;
---------several ex-bugs just out of pure paranoia:
(no output)
---------output of numeric types:
-Using 200 digits
-Pi evalfs to: 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644L0
-Catalan evalfs to: 0.9159655941772190150546035149323841107741493742816721342664981196217630197762547694793565129261151062485744226191961995790358988033258590594315947374811584069953320287733194605190387274781640878659090247L0
-EulerGamma evalfs to: 0.577215664901532860606512090082402431042159335939923598805767234884867726777664670936947063291746749514631447249807082480960504014486542836224173997644923536253500333742937337737673942792595258247094916L0
+Using 222 digits
+Pi evalfs to: 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648235L0
+Catalan evalfs to: 0.915965594177219015054603514932384110774149374281672134266498119621763019776254769479356512926115106248574422619196199579035898803325859059431594737481158406995332028773319460519038727478164087865909024706484152163000228727640942388L0
+EulerGamma evalfs to: 0.577215664901532860606512090082402431042159335939923598805767234884867726777664670936947063291746749514631447249807082480960504014486542836224173997644923536253500333742937337737673942792595258247094916008735203948165670853233151777L0
Complex integers: {(0,0)=0} {(1,0)=1} {(1,1)=1+I} {(0,1)=I} {(-1,1)=-1+I} {(-1,0)=-1} {(-1,-1)=-1-I} {(0,-1)=-I} {(1,-1)=1-I}
---------consistency of numeric types:
(no output)
static symbol x("x");
-static unsigned check_series(const ex &e, const ex &point, const ex &d)
+static unsigned check_series(const ex &e, const ex &point, const ex &d, int order = 8)
{
- ex es = e.series(x, point, 8);
- ex ep = static_cast<series *>(es.bp)->convert_to_poly();
- if ((ep - d).compare(exZERO()) != 0) {
- clog << "series expansion of " << e << " at " << point
+ ex es = e.series(x, point, order);
+ ex ep = static_cast<series *>(es.bp)->convert_to_poly();
+ if ((ep - d).compare(exZERO()) != 0) {
+ clog << "series expansion of " << e << " at " << point
<< " erroneously returned " << ep << " (instead of " << d
<< ")" << endl;
- (ep-d).printtree(clog);
- return 1;
- }
- return 0;
+ (ep-d).printtree(clog);
+ return 1;
+ }
+ return 0;
}
// Series expansion
static unsigned series1(void)
{
- unsigned result = 0;
- ex e, d;
-
- e = sin(x);
- d = x - pow(x, 3) / 6 + pow(x, 5) / 120 - pow(x, 7) / 5040 + Order(pow(x, 8));
- result += check_series(e, exZERO(), d);
-
- e = cos(x);
- d = 1 - pow(x, 2) / 2 + pow(x, 4) / 24 - pow(x, 6) / 720 + Order(pow(x, 8));
- result += check_series(e, exZERO(), d);
-
- e = exp(x);
- d = 1 + x + pow(x, 2) / 2 + pow(x, 3) / 6 + pow(x, 4) / 24 + pow(x, 5) / 120 + pow(x, 6) / 720 + pow(x, 7) / 5040 + Order(pow(x, 8));
- result += check_series(e, exZERO(), d);
-
- e = pow(1 - x, -1);
- d = 1 + x + pow(x, 2) + pow(x, 3) + pow(x, 4) + pow(x, 5) + pow(x, 6) + pow(x, 7) + Order(pow(x, 8));
- result += check_series(e, exZERO(), d);
-
- e = x + pow(x, -1);
- d = x + pow(x, -1);
- result += check_series(e, exZERO(), d);
-
- e = x + pow(x, -1);
- d = 2 + pow(x-1, 2) - pow(x-1, 3) + pow(x-1, 4) - pow(x-1, 5) + pow(x-1, 6) - pow(x-1, 7) + Order(pow(x-1, 8));
- result += check_series(e, exONE(), d);
-
- e = pow(x + pow(x, 3), -1);
- d = pow(x, -1) - x + pow(x, 3) - pow(x, 5) + Order(pow(x, 7));
- result += check_series(e, exZERO(), d);
-
- e = pow(pow(x, 2) + pow(x, 4), -1);
- d = pow(x, -2) - 1 + pow(x, 2) - pow(x, 4) + Order(pow(x, 6));
- result += check_series(e, exZERO(), d);
-
- e = pow(sin(x), -2);
- d = pow(x, -2) + numeric(1,3) + pow(x, 2) / 15 + pow(x, 4) * 2/189 + Order(pow(x, 5));
- result += check_series(e, exZERO(), d);
-
- e = sin(x) / cos(x);
- d = x + pow(x, 3) / 3 + pow(x, 5) * 2/15 + pow(x, 7) * 17/315 + Order(pow(x, 8));
- result += check_series(e, exZERO(), d);
-
- e = cos(x) / sin(x);
- d = pow(x, -1) - x / 3 - pow(x, 3) / 45 - pow(x, 5) * 2/945 + Order(pow(x, 6));
- result += check_series(e, exZERO(), d);
-
- e = pow(numeric(2), x);
- ex t = log(ex(2)) * x;
- d = 1 + t + pow(t, 2) / 2 + pow(t, 3) / 6 + pow(t, 4) / 24 + pow(t, 5) / 120 + pow(t, 6) / 720 + pow(t, 7) / 5040 + Order(pow(x, 8));
- result += check_series(e, exZERO(), d.expand());
-
- e = pow(Pi, x);
- t = log(Pi) * x;
- d = 1 + t + pow(t, 2) / 2 + pow(t, 3) / 6 + pow(t, 4) / 24 + pow(t, 5) / 120 + pow(t, 6) / 720 + pow(t, 7) / 5040 + Order(pow(x, 8));
- result += check_series(e, exZERO(), d.expand());
-
- return result;
+ unsigned result = 0;
+ ex e, d;
+
+ e = sin(x);
+ d = x - pow(x, 3) / 6 + pow(x, 5) / 120 - pow(x, 7) / 5040 + Order(pow(x, 8));
+ result += check_series(e, exZERO(), d);
+
+ e = cos(x);
+ d = 1 - pow(x, 2) / 2 + pow(x, 4) / 24 - pow(x, 6) / 720 + Order(pow(x, 8));
+ result += check_series(e, exZERO(), d);
+
+ e = exp(x);
+ d = 1 + x + pow(x, 2) / 2 + pow(x, 3) / 6 + pow(x, 4) / 24 + pow(x, 5) / 120 + pow(x, 6) / 720 + pow(x, 7) / 5040 + Order(pow(x, 8));
+ result += check_series(e, exZERO(), d);
+
+ e = pow(1 - x, -1);
+ d = 1 + x + pow(x, 2) + pow(x, 3) + pow(x, 4) + pow(x, 5) + pow(x, 6) + pow(x, 7) + Order(pow(x, 8));
+ result += check_series(e, exZERO(), d);
+
+ e = x + pow(x, -1);
+ d = x + pow(x, -1);
+ result += check_series(e, exZERO(), d);
+
+ e = x + pow(x, -1);
+ d = 2 + pow(x-1, 2) - pow(x-1, 3) + pow(x-1, 4) - pow(x-1, 5) + pow(x-1, 6) - pow(x-1, 7) + Order(pow(x-1, 8));
+ result += check_series(e, exONE(), d);
+
+ e = pow(x + pow(x, 3), -1);
+ d = pow(x, -1) - x + pow(x, 3) - pow(x, 5) + Order(pow(x, 7));
+ result += check_series(e, exZERO(), d);
+
+ e = pow(pow(x, 2) + pow(x, 4), -1);
+ d = pow(x, -2) - 1 + pow(x, 2) - pow(x, 4) + Order(pow(x, 6));
+ result += check_series(e, exZERO(), d);
+
+ e = pow(sin(x), -2);
+ d = pow(x, -2) + numeric(1,3) + pow(x, 2) / 15 + pow(x, 4) * 2/189 + Order(pow(x, 5));
+ result += check_series(e, exZERO(), d);
+
+ e = sin(x) / cos(x);
+ d = x + pow(x, 3) / 3 + pow(x, 5) * 2/15 + pow(x, 7) * 17/315 + Order(pow(x, 8));
+ result += check_series(e, exZERO(), d);
+
+ e = cos(x) / sin(x);
+ d = pow(x, -1) - x / 3 - pow(x, 3) / 45 - pow(x, 5) * 2/945 + Order(pow(x, 6));
+ result += check_series(e, exZERO(), d);
+
+ e = pow(numeric(2), x);
+ ex t = log(ex(2)) * x;
+ d = 1 + t + pow(t, 2) / 2 + pow(t, 3) / 6 + pow(t, 4) / 24 + pow(t, 5) / 120 + pow(t, 6) / 720 + pow(t, 7) / 5040 + Order(pow(x, 8));
+ result += check_series(e, exZERO(), d.expand());
+
+ e = pow(Pi, x);
+ t = log(Pi) * x;
+ d = 1 + t + pow(t, 2) / 2 + pow(t, 3) / 6 + pow(t, 4) / 24 + pow(t, 5) / 120 + pow(t, 6) / 720 + pow(t, 7) / 5040 + Order(pow(x, 8));
+ result += check_series(e, exZERO(), d.expand());
+
+ return result;
}
// Series addition
static unsigned series2(void)
{
- unsigned result = 0;
- ex e, d;
-
- e = pow(sin(x), -1).series(x, exZERO(), 8) + pow(sin(-x), -1).series(x, exZERO(), 12);
- d = Order(pow(x, 6));
- result += check_series(e, exZERO(), d);
-
- return result;
+ unsigned result = 0;
+ ex e, d;
+
+ e = pow(sin(x), -1).series(x, exZERO(), 8) + pow(sin(-x), -1).series(x, exZERO(), 12);
+ d = Order(pow(x, 6));
+ result += check_series(e, exZERO(), d);
+
+ return result;
}
// Series multiplication
static unsigned series3(void)
{
- unsigned result = 0;
- ex e, d;
-
- e = sin(x).series(x, exZERO(), 8) * pow(sin(x), -1).series(x, exZERO(), 12);
- d = 1 + Order(pow(x, 7));
- result += check_series(e, exZERO(), d);
+ unsigned result = 0;
+ ex e, d;
+
+ e = sin(x).series(x, exZERO(), 8) * pow(sin(x), -1).series(x, exZERO(), 12);
+ d = 1 + Order(pow(x, 7));
+ result += check_series(e, exZERO(), d);
+
+ return result;
+}
- return result;
+// Series of special functions
+static unsigned series4(void)
+{
+ unsigned result = 0;
+ ex e, d;
+
+ e = gamma(2*x);
+ d = pow(x+1,-1)*numeric(1,4) +
+ pow(x+1,0)*(numeric(3,4) -
+ numeric(1,2)*EulerGamma) +
+ pow(x+1,1)*(numeric(7,4) -
+ numeric(3,2)*EulerGamma +
+ numeric(1,2)*pow(EulerGamma,2) +
+ numeric(1,12)*pow(Pi,2)) +
+ pow(x+1,2)*(numeric(15,4) -
+ numeric(7,2)*EulerGamma -
+ numeric(1,3)*pow(EulerGamma,3) +
+ numeric(1,4)*pow(Pi,2) +
+ numeric(3,2)*pow(EulerGamma,2) -
+ numeric(1,6)*pow(Pi,2)*EulerGamma -
+ numeric(2,3)*zeta(3)) +
+ pow(x+1,3)*(numeric(31,4) - pow(EulerGamma,3) -
+ numeric(15,2)*EulerGamma +
+ numeric(1,6)*pow(EulerGamma,4) +
+ numeric(7,2)*pow(EulerGamma,2) +
+ numeric(7,12)*pow(Pi,2) -
+ numeric(1,2)*pow(Pi,2)*EulerGamma -
+ numeric(2)*zeta(3) +
+ numeric(1,6)*pow(EulerGamma,2)*pow(Pi,2) +
+ numeric(1,40)*pow(Pi,4) +
+ numeric(4,3)*zeta(3)*EulerGamma) +
+ Order(pow(x+1,4));
+ result += check_series(e, -1, d, 4);
+
+ e = tan(x*Pi/2);
+ d = pow(x-1,-1)/Pi*(-2) +
+ pow(x-1,1)*Pi/6 +
+ pow(x-1,3)*pow(Pi,3)/360 +
+ pow(x-1,5)*pow(Pi,5)/15120 +
+ pow(x-1,7)*pow(Pi,7)/604800 +
+ Order(pow(x-1,8));
+ result += check_series(e,1,d,8);
+
+ return result;
}
unsigned series_expansion(void)
{
- unsigned result = 0;
-
- cout << "checking series expansion..." << flush;
- clog << "---------series expansion:" << endl;
-
- result += series1();
- result += series2();
- result += series3();
-
- if (!result) {
- cout << " passed ";
- clog << "(no output)" << endl;
- } else {
- cout << " failed ";
- }
- return result;
+ unsigned result = 0;
+
+ cout << "checking series expansion..." << flush;
+ clog << "---------series expansion:" << endl;
+
+ result += series1();
+ result += series2();
+ result += series3();
+ result += series4();
+
+ if (!result) {
+ cout << " passed ";
+ clog << "(no output)" << endl;
+ } else {
+ cout << " failed ";
+ }
+ return result;
}
@for file in $(DISTFILES); do \
d=$(srcdir); \
if test -d $$d/$$file; then \
- cp -pr $$d/$$file $(distdir)/$$file; \
+ cp -pr $$/$$file $(distdir)/$$file; \
else \
test -f $(distdir)/$$file \
|| ln $$d/$$file $(distdir)/$$file 2> /dev/null \
@for file in $(DISTFILES); do \
d=$(srcdir); \
if test -d $$d/$$file; then \
- cp -pr $$d/$$file $(distdir)/$$file; \
+ cp -pr $$/$$file $(distdir)/$$file; \
else \
test -f $(distdir)/$$file \
|| ln $$d/$$file $(distdir)/$$file 2> /dev/null \
@for file in $(DISTFILES); do \
d=$(srcdir); \
if test -d $$d/$$file; then \
- cp -pr $$d/$$file $(distdir)/$$file; \
+ cp -pr $$/$$file $(distdir)/$$file; \
else \
test -f $(distdir)/$$file \
|| ln $$d/$$file $(distdir)/$$file 2> /dev/null \
1200 2
5 1 1 2 0 7 100 0 -1 4.000 0 0 1 0 1896.500 2807.128 791 579 1973 321 3002 579
1 1 1.00 68.57 137.14
-6 996 1144 2025 1401
-2 2 0 1 0 7 100 0 20 0.000 0 0 -1 0 0 5
- 2025 1144 996 1144 996 1401 2025 1401 2025 1144
-4 1 0 99 0 0 14 0.0000 4 195 840 1511 1350 expairseq\001
--6
6 3053 630 3876 887
2 2 0 1 0 7 100 0 20 0.000 0 0 -1 0 0 5
3053 630 3876 630 3876 887 3053 887 3053 630
4 1 0 99 0 0 14 0.0000 4 150 450 3465 836 basic\001
-6
-6 2340 1935 3195 2205
-2 2 0 1 0 7 100 0 20 0.000 0 0 -1 0 0 5
- 2340 1948 3195 1948 3195 2205 2340 2205 2340 1948
-4 1 0 99 0 0 14 0.0000 4 150 690 2764 2154 exprseq\001
+6 225 450 765 720
+2 2 0 1 0 30 100 0 20 0.000 0 0 -1 0 0 5
+ 225 463 739 463 739 720 225 720 225 463
+4 1 0 99 0 0 14 0.0000 4 105 210 482 668 ex\001
-6
-6 2565 2925 3510 3195
-2 2 0 1 0 7 100 0 20 0.000 0 0 -1 0 0 5
- 2565 2925 3510 2925 3510 3195 2565 3195 2565 2925
-4 1 0 99 0 0 14 0.0000 4 150 675 3034 3131 indexed\001
+6 4770 1980 5760 2340
+1 2 0 1 0 30 100 0 20 0.000 1 0.0000 5287 2182 473 158 4814 2182 5760 2182
+4 1 0 99 0 0 14 0.0000 4 150 705 5310 2250 numeric\001
-6
-6 3060 3645 3420 3735
-1 4 0 1 0 0 100 0 20 0.000 1 0.0000 3127 3690 23 23 3104 3690 3150 3690
-1 4 0 1 0 0 100 0 20 0.000 1 0.0000 3263 3690 23 23 3240 3690 3286 3690
-1 4 0 1 0 0 100 0 20 0.000 1 0.0000 3397 3690 23 23 3374 3690 3420 3690
+6 4770 1530 5760 1890
+1 2 0 1 0 30 100 0 20 0.000 1 0.0000 5287 1732 473 158 4814 1732 5760 1732
+4 1 0 99 0 0 14 0.0000 4 135 735 5310 1800 constant\001
-6
-6 4410 3600 4770 3690
-1 4 0 1 0 0 100 0 20 0.000 1 0.0000 4477 3645 23 23 4454 3645 4500 3645
-1 4 0 1 0 0 100 0 20 0.000 1 0.0000 4613 3645 23 23 4590 3645 4636 3645
-1 4 0 1 0 0 100 0 20 0.000 1 0.0000 4747 3645 23 23 4724 3645 4770 3645
+6 4770 1080 5760 1440
+1 2 0 1 0 30 100 0 20 0.000 1 0.0000 5287 1282 473 158 4814 1282 5760 1282
+4 1 0 99 0 0 14 0.0000 4 195 615 5310 1350 symbol\001
-6
-6 585 1980 1125 2295
+6 996 1080 2025 1337
+2 2 0 1 0 7 100 0 20 0.000 0 0 -1 0 0 5
+ 2025 1080 996 1080 996 1337 2025 1337 2025 1080
+4 1 0 99 0 0 14 0.0000 4 195 840 1511 1286 expairseq\001
+-6
+6 495 1800 1035 2115
2 2 0 1 0 30 100 0 20 0.000 0 0 -1 0 0 5
- 585 2019 1099 2019 1099 2276 585 2276 585 2019
-4 1 0 99 0 0 14 0.0000 4 150 315 842 2224 add\001
+ 495 1839 1009 1839 1009 2096 495 2096 495 1839
+4 1 0 99 0 0 14 0.0000 4 150 315 752 2044 add\001
-6
-6 225 450 765 720
+6 1260 1800 1800 2115
2 2 0 1 0 30 100 0 20 0.000 0 0 -1 0 0 5
- 225 463 739 463 739 720 225 720 225 463
-4 1 0 99 0 0 14 0.0000 4 105 210 482 668 ex\001
+ 1260 1839 1774 1839 1774 2096 1260 2096 1260 1839
+4 1 0 99 0 0 14 0.0000 4 150 315 1517 2044 mul\001
-6
-6 1350 1980 1890 2295
+6 2295 2295 3150 2565
+2 2 0 1 0 7 100 0 20 0.000 0 0 -1 0 0 5
+ 2295 2308 3150 2308 3150 2565 2295 2565 2295 2308
+4 1 0 99 0 0 14 0.0000 4 150 690 2719 2514 exprseq\001
+-6
+6 2520 3285 3465 3555
+2 2 0 1 0 7 100 0 20 0.000 0 0 -1 0 0 5
+ 2520 3285 3465 3285 3465 3555 2520 3555 2520 3285
+4 1 0 99 0 0 14 0.0000 4 150 675 2989 3491 indexed\001
+-6
+6 3015 4005 3375 4095
+1 4 0 1 0 0 100 0 20 0.000 1 0.0000 3082 4050 23 23 3059 4050 3105 4050
+1 4 0 1 0 0 100 0 20 0.000 1 0.0000 3218 4050 23 23 3195 4050 3241 4050
+1 4 0 1 0 0 100 0 20 0.000 1 0.0000 3352 4050 23 23 3329 4050 3375 4050
+-6
+6 1440 3240 2385 3510
2 2 0 1 0 30 100 0 20 0.000 0 0 -1 0 0 5
- 1350 2019 1864 2019 1864 2276 1350 2276 1350 2019
-4 1 0 99 0 0 14 0.0000 4 150 315 1607 2224 mul\001
+ 1440 3240 2385 3240 2385 3510 1440 3510 1440 3240
+4 1 0 99 0 0 14 0.0000 4 150 690 1909 3446 function\001
-6
-6 1305 2520 2025 2790
+6 1800 3870 2475 4140
2 2 0 1 0 30 100 0 20 0.000 0 0 -1 0 0 5
- 1305 2533 2025 2533 2025 2790 1305 2790 1305 2533
-4 1 0 99 0 0 14 0.0000 4 150 525 1672 2739 ncmul\001
+ 1800 3870 2475 3870 2475 4140 1800 4140 1800 3870
+4 1 0 99 0 0 14 0.0000 4 195 585 2134 4076 isospin\001
-6
-6 1485 2880 2430 3150
+6 2520 3870 3015 4140
2 2 0 1 0 30 100 0 20 0.000 0 0 -1 0 0 5
- 1485 2880 2430 2880 2430 3150 1485 3150 1485 2880
-4 1 0 99 0 0 14 0.0000 4 150 690 1954 3086 function\001
+ 2520 3870 3015 3870 3015 4140 2520 4140 2520 3870
+4 1 0 99 0 0 14 0.0000 4 150 435 2764 4076 color\001
-6
-6 1845 3510 2520 3780
+6 1350 2250 2070 2520
2 2 0 1 0 30 100 0 20 0.000 0 0 -1 0 0 5
- 1845 3510 2520 3510 2520 3780 1845 3780 1845 3510
-4 1 0 99 0 0 14 0.0000 4 195 585 2179 3716 isospin\001
+ 1350 2250 2070 2250 2070 2520 1350 2520 1350 2250
+4 1 0 99 0 0 14 0.0000 4 150 510 1684 2456 series\001
-6
-6 2565 3510 3060 3780
+6 1755 1440 2475 1710
2 2 0 1 0 30 100 0 20 0.000 0 0 -1 0 0 5
- 2565 3510 3060 3510 3060 3780 2565 3780 2565 3510
-4 1 0 99 0 0 14 0.0000 4 150 435 2809 3716 color\001
+ 1755 1453 2475 1453 2475 1710 1755 1710 1755 1453
+4 1 0 99 0 0 14 0.0000 4 150 555 2115 1659 power\001
-6
-6 3555 3465 4365 3780
-1 2 0 1 0 30 100 0 20 0.000 1 0.0000 3960 3622 405 157 3555 3622 4365 3622
-4 1 0 99 0 0 14 0.0000 4 150 690 3960 3690 coloridx\001
+6 1170 2835 1890 3105
+2 2 0 1 0 30 100 0 20 0.000 0 0 -1 0 0 5
+ 1170 2848 1890 2848 1890 3105 1170 3105 1170 2848
+4 1 0 99 0 0 14 0.0000 4 150 525 1537 3054 ncmul\001
-6
-6 3600 2880 4275 3195
-1 2 0 1 0 30 100 0 20 0.000 1 0.0000 3937 3037 337 157 3600 3037 4274 3037
-4 1 0 99 0 0 14 0.0000 4 150 255 3960 3105 idx\001
+6 4455 3960 4815 4050
+1 4 0 1 0 0 100 0 20 0.000 1 0.0000 4522 4005 23 23 4499 4005 4545 4005
+1 4 0 1 0 0 100 0 20 0.000 1 0.0000 4658 4005 23 23 4635 4005 4681 4005
+1 4 0 1 0 0 100 0 20 0.000 1 0.0000 4792 4005 23 23 4769 4005 4815 4005
-6
-6 3780 2250 4320 2520
-2 2 0 1 0 30 100 0 20 0.000 0 0 -1 0 0 5
- 3806 2263 4320 2263 4320 2520 3806 2520 3806 2263
-4 1 0 99 0 0 14 0.0000 4 150 195 4063 2468 lst\001
+6 3600 3825 4410 4140
+1 2 0 1 0 30 100 0 20 0.000 1 0.0000 4005 3982 405 157 3600 3982 4410 3982
+4 1 0 99 0 0 14 0.0000 4 150 690 4005 4050 coloridx\001
-6
-6 2160 1485 2880 1755
-2 2 0 1 0 30 100 0 20 0.000 0 0 -1 0 0 5
- 2160 1498 2880 1498 2880 1755 2160 1755 2160 1498
-4 1 0 99 0 0 14 0.0000 4 150 555 2520 1704 power\001
+6 3645 3240 4320 3555
+1 2 0 1 0 30 100 0 20 0.000 1 0.0000 3982 3397 337 157 3645 3397 4319 3397
+4 1 0 99 0 0 14 0.0000 4 150 255 4005 3465 idx\001
-6
-6 4725 2970 5445 3240
+6 4770 2565 5760 2835
2 2 0 1 0 30 100 0 20 0.000 0 0 -1 0 0 5
- 4725 2970 5445 2970 5445 3240 4725 3240 4725 2970
-4 1 0 99 0 0 14 0.0000 4 150 510 5059 3176 series\001
+ 4770 2578 5760 2578 5760 2835 4770 2835 4770 2578
+4 1 0 99 0 0 14 0.0000 4 150 795 5271 2784 relational\001
-6
-6 4815 2520 5805 2790
+6 4455 2970 5265 3240
2 2 0 1 0 30 100 0 20 0.000 0 0 -1 0 0 5
- 4853 2533 5779 2533 5779 2790 4853 2790 4853 2533
-4 1 0 99 0 0 14 0.0000 4 150 555 5316 2739 matrix\001
--6
-6 4770 1980 5760 2340
-1 2 0 1 0 30 100 0 20 0.000 1 0.0000 5287 2182 473 158 4814 2182 5760 2182
-4 1 0 99 0 0 14 0.0000 4 150 705 5310 2250 numeric\001
--6
-6 4770 1530 5760 1890
-1 2 0 1 0 30 100 0 20 0.000 1 0.0000 5287 1732 473 158 4814 1732 5760 1732
-4 1 0 99 0 0 14 0.0000 4 135 735 5310 1800 constant\001
+ 4455 2983 5265 2983 5265 3240 4455 3240 4455 2983
+4 1 0 99 0 0 14 0.0000 4 150 555 4866 3189 matrix\001
-6
-6 4770 1080 5760 1440
-1 2 0 1 0 30 100 0 20 0.000 1 0.0000 5287 1282 473 158 4814 1282 5760 1282
-4 1 0 99 0 0 14 0.0000 4 195 615 5310 1350 symbol\001
+6 3825 2385 4365 2655
+2 2 0 1 0 30 100 0 20 0.000 0 0 -1 0 0 5
+ 3851 2398 4365 2398 4365 2655 3851 2655 3851 2398
+4 1 0 99 0 0 14 0.0000 4 150 195 4108 2603 lst\001
-6
2 1 0 1 0 7 100 0 -1 0.000 0 0 -1 0 1 2
1 1 1.00 60.00 120.00
3825 945 4802 2121
2 1 0 1 0 7 100 0 -1 0.000 0 0 -1 0 1 2
1 1 1.00 60.00 120.00
- 3053 939 2076 1144
+ 3053 939 2070 1080
2 1 0 1 0 7 100 0 -1 0.000 0 0 -1 0 1 2
1 1 1.00 60.00 120.00
- 3156 939 2693 1453
+ 3156 939 2430 1395
2 1 0 1 0 7 100 0 -1 0.000 0 0 -1 0 1 2
1 1 1.00 60.00 120.00
- 3722 939 4802 2584
+ 3722 939 4815 2565
2 1 0 1 0 7 100 0 -1 0.000 0 0 -1 0 1 2
1 1 1.00 60.00 120.00
3915 945 4815 1665
2 1 0 1 0 7 100 0 -1 0.000 0 0 -1 0 1 2
1 1 1.00 60.00 120.00
- 3208 939 2925 1935
+ 3285 945 2790 2250
2 1 0 1 0 7 100 0 -1 0.000 0 0 -1 0 1 2
1 1 1.00 60.00 120.00
- 1459 1453 894 1967
+ 1305 1395 810 1800
2 1 0 1 0 7 100 0 -1 0.000 0 0 -1 0 1 2
1 1 1.00 60.00 120.00
- 1562 1453 1575 1980
+ 1395 1395 1530 1800
2 1 0 1 0 7 100 0 -1 0.000 0 0 -1 0 1 2
1 1 1.00 60.00 120.00
3915 881 4815 1215
2 1 0 1 0 7 100 0 -1 0.000 0 0 -1 0 1 2
1 1 1.00 60.00 120.00
- 3645 945 4770 2925
+ 3645 945 4725 2925
+2 1 0 1 0 7 100 0 -1 0.000 0 0 -1 0 1 2
+ 1 1 1.00 60.00 120.00
+ 3555 945 4050 2340
2 1 0 1 0 7 100 0 -1 0.000 0 0 -1 0 1 2
1 1 1.00 60.00 120.00
- 2430 2250 1935 2475
+ 3420 945 3825 3240
2 1 0 1 0 7 100 0 -1 0.000 0 0 -1 0 1 2
1 1 1.00 60.00 120.00
- 3599 937 4005 2205
+ 2385 2610 1890 2835
2 1 0 1 0 7 100 0 -1 0.000 0 0 -1 0 1 2
1 1 1.00 60.00 120.00
- 2610 2295 2160 2835
+ 2565 2610 2115 3195
2 1 0 1 0 7 100 0 -1 0.000 0 0 -1 0 1 2
1 1 1.00 60.00 120.00
- 2835 2295 2970 2880
+ 2790 2610 2925 3240
2 1 0 1 0 7 100 0 -1 0.000 0 0 -1 0 1 2
1 1 1.00 60.00 120.00
- 3330 945 3780 2880
+ 2610 3600 2250 3825
2 1 0 1 0 7 100 0 -1 0.000 0 0 -1 0 1 2
1 1 1.00 60.00 120.00
- 2655 3240 2295 3465
+ 2790 3600 2790 3825
2 1 0 1 0 7 100 0 -1 0.000 0 0 -1 0 1 2
1 1 1.00 60.00 120.00
- 2835 3240 2835 3465
+ 3240 945 2115 2205
2 1 0 1 0 7 100 0 -1 0.000 0 0 -1 0 1 2
1 1 1.00 60.00 120.00
- 3870 3240 3870 3465
+ 3915 3600 3915 3825
@section How to use it from within C++
The GiNaC open framework for symbolic computation within the C++ programming
-language does not try to define a language of it's own as conventional
+language does not try to define a language of its own as conventional
CAS do. Instead, it extends the capabilities of C++ by symbolic
manipulations. Here is how to generate and print a simple (and rather
pointless) bivariate polynomial with some large coefficients:
@c node-name, next, previous, up
@section What it can do for you
-@cindex @code{ginsh}
+@cindex @command{ginsh}
After invoking @command{ginsh} one can test and experiment with GiNaC's
features much like in other Computer Algebra Systems except that it does
not provide programming constructs like loops or conditionals. For a
series (the third argument of @code{series} is the evaluation point, the
fourth defines the order):
+@cindex Zeta function
@example
> diff(tan(x),x);
tan(x)^2+1
x-1/6*x^3+Order(x^4)
> series(1/tan(x),x,0,4);
x^(-1)-1/3*x+Order(x^2)
+> series(gamma(x),x,0,3);
+x^(-1)-EulerGamma+(1/12*Pi^2+1/2*EulerGamma^2)*x
++(-1/3*zeta(3)-1/12*Pi^2*EulerGamma-1/6*EulerGamma^3)*x^2+Order(x^3)
+> evalf(");
+x^(-1.0)-0.5772156649015328606+(0.98905599532797255544)*x
+-(0.90747907608088628905)*x^2+Order(x^(3.0))
> series(gamma(2*sin(x)-2),x,Pi/2,6);
-(x-1/2*Pi)^(-2)+(-1/12*Pi^2-1/2*EulerGamma^2-1/240)*(x-1/2*Pi)^2
-EulerGamma-1/12+Order((x-1/2*Pi)^3)
@end example
-If you ever wanted to convert units in C or C++ and found this
-is cumbersome, here is the solution. Symbolic types can always be
-used as tags for different types of objects. Converting from wrong
-units to the metric system is therefore easy:
+Here we have made use of the @command{ginsh}-command @code{"} to pop the
+previously evaluated element from @command{ginsh}'s internal stack.
+
+If you ever wanted to convert units in C or C++ and found this is
+cumbersome, here is the solution. Symbolic types can always be used as
+tags for different types of objects. Converting from wrong units to the
+metric system is therefore easy:
@example
> in=.0254*m;
@option{--prefix=@var{PREFIX}}: The directory where the compiled library
and headers are installed. It defaults to @file{/usr/local} which means
that the library is installed in the directory @file{/usr/local/lib},
-the header files in @file{/usr/local/include/GiNaC} and the documentation
+the header files in @file{/usr/local/include/ginac} and the documentation
(like this one) into @file{/usr/local/share/doc/GiNaC}.
@item
@file{@var{PREFIX}/include/ginac/}. For instance, if you specify
@option{--includedir=/usr/include} you will end up with the header files
sitting in the directory @file{/usr/include/ginac/}. Note that the
-subdirectory @file{GiNaC} is enforced by this process in order to
+subdirectory @file{ginac} is enforced by this process in order to
keep the header files separated from others. This avoids some
clashes and allows for an easier deinstallation of GiNaC. This ought
to be considered A Good Thing (tm).
And here is a configuration for a private static GiNaC library with
several components sitting in custom places (site-wide @acronym{GCC} and
private @acronym{CLN}). The compiler is pursuaded to be picky and full
-assertions and debugging are switched on:
+assertions and debugging information are switched on:
@example
$ export CXX=/usr/local/gnu/bin/c++
* Symbols:: Symbolic objects.
* Numbers:: Numerical objects.
* Constants:: Pre-defined constants.
-* Fundamental operations:: The power, add and mul classes.
+* Fundamental containers:: The power, add and mul classes.
* Built-in functions:: Mathematical functions.
+* Relations:: Equality, Inequality and all that.
@end menu
mathematical objects, all of which (except for @code{ex} and some
helpers) are internally derived from one abstract base class called
@code{basic}. You do not have to deal with objects of class
-@code{basic}, instead you'll be dealing with symbols and functions of
-symbols. You'll soon learn in this chapter how many of the functions on
-symbols are really classes. This is because simple symbolic arithmetic
-is not supported by languages like C++ so in a certain way GiNaC has to
-implement its own arithmetic.
-
-To give an idea about what kinds of symbolic composits may be built we
+@code{basic}, instead you'll be dealing with symbols, numbers,
+containers of expressions and so on. You'll soon learn in this chapter
+how many of the functions on symbols are really classes. This is
+because simple symbolic arithmetic is not supported by languages like
+C++ so in a certain way GiNaC has to implement its own arithmetic.
+
+@cindex container
+@cindex atom
+To get an idea about what kinds of symbolic composits may be built we
have a look at the most important classes in the class hierarchy. The
oval classes are atomic ones and the squared classes are containers.
-The dashed line symbolizes a "points to" or "handles" relationship while
-the solid lines stand for "inherits from" relationship in the class
+The dashed line symbolizes a `points to' or `handles' relationship while
+the solid lines stand for `inherits from' relationship in the class
hierarchy:
@image{classhierarchy}
consisting of one expression and a number (@code{numeric}). What
@emph{is} visible to the user are the derived classes @code{add} and
@code{mul}, representing sums of terms and products, respectively.
-We'll come back later to some more details about these two classes and
-motivate the use of pairs in sums and products here.
+@xref{Internal Structures}, where these two classes are described in
+more detail.
@node Symbols, Numbers, The Class Hierarchy, Basic Concepts
@c node-name, next, previous, up
@section Symbols
@cindex Symbols (class @code{symbol})
+@cindex hierarchy of classes
+@cindex atom
Symbols are for symbolic manipulation what atoms are for chemistry. You
can declare objects of class @code{symbol} as any other object simply by
saying @code{symbol x,y;}. There is, however, a catch in here having to
in the calling function. Again, comparing them (using @code{operator==}
for instance) will always reveal their difference. Watch out, please.
+@cindex @code{subs()}
Although symbols can be assigned expressions for internal reasons, you
should not do it (and we are not going to tell you how it is done). If
you want to replace a symbol with something else in an expression, you
@section Numbers
@cindex numbers (class @code{numeric})
+@cindex GMP
@cindex CLN
+@cindex rational
+@cindex fraction
For storing numerical things, GiNaC uses Bruno Haible's library
@acronym{CLN}. The classes therein serve as foundation classes for
GiNaC. @acronym{CLN} stands for Class Library for Numbers or
numerics at work in most of our examples. This design really becomes
convenient when one declares own functions having more than one
parameter but it forbids using implicit constructors because that would
-lead to ambiguities.
+lead to compile-time ambiguities.
It may be tempting to construct numbers writing @code{numeric r(3/2)}.
This would, however, call C's built-in operator @code{/} for integers
first and result in a numeric holding a plain integer 1. @strong{Never
-use @code{/} on integers!} Use the constructor from two integers
-instead, as shown in the example above. Writing @code{numeric(1)/2} may
-look funny but works also.
+use the operator @code{/} on integers} unless you know exactly what you
+are doing! Use the constructor from two integers instead, as shown in
+the example above. Writing @code{numeric(1)/2} may look funny but works
+also.
@cindex @code{Digits}
@cindex accuracy
We have seen now the distinction between exact numbers and floating
point numbers. Clearly, the user should never have to worry about
-dynamically created exact numbers, since their "exactness" always
-determines how they ought to be handled, i.e. how "long" they are. The
+dynamically created exact numbers, since their `exactness' always
+determines how they ought to be handled, i.e. how `long' they are. The
situation is different for floating point numbers. Their accuracy is
controlled by one @emph{global} variable, called @code{Digits}. (For
those readers who know about Maple: it behaves very much like Maple's
following table.
@cartouche
-@multitable @columnfractions .33 .67
-@item Method @tab Returns true if@dots{}
+@multitable @columnfractions .30 .70
+@item @strong{Method} @tab @strong{Returns true if the object is@dots{}}
@item @code{.is_zero()}
-@tab object is equal to zero
+@tab @dots{}equal to zero
@item @code{.is_positive()}
-@tab object is not complex and greater than 0
+@tab @dots{}not complex and greater than 0
@item @code{.is_integer()}
-@tab object is a (non-complex) integer
+@tab @dots{}a (non-complex) integer
@item @code{.is_pos_integer()}
-@tab object is an integer and greater than 0
+@tab @dots{}an integer and greater than 0
@item @code{.is_nonneg_integer()}
-@tab object is an integer and greater equal 0
+@tab @dots{}an integer and greater equal 0
@item @code{.is_even()}
-@tab object is an even integer
+@tab @dots{}an even integer
@item @code{.is_odd()}
-@tab object is an odd integer
+@tab @dots{}an odd integer
@item @code{.is_prime()}
-@tab object is a prime integer (probabilistic primality test)
+@tab @dots{}a prime integer (probabilistic primality test)
@item @code{.is_rational()}
-@tab object is an exact rational number (integers are rational, too)
+@tab @dots{}an exact rational number (integers are rational, too)
@item @code{.is_real()}
-@tab object is a real integer, rational or float (i.e. is not complex)
+@tab @dots{}a real integer, rational or float (i.e. is not complex)
@item @code{.is_cinteger()}
-@tab object is a (complex) integer, such as @math{2-3*I}
+@tab @dots{}a (complex) integer, such as @math{2-3*I}
@item @code{.is_crational()}
-@tab object is an exact (complex) rational number (such as @math{2/3+7/2*I})
+@tab @dots{}an exact (complex) rational number (such as @math{2/3+7/2*I})
@end multitable
@end cartouche
-@node Constants, Fundamental operations, Numbers, Basic Concepts
+@node Constants, Fundamental containers, Numbers, Basic Concepts
@c node-name, next, previous, up
@section Constants
@cindex constants (class @code{constant})
-@cindex @code{evalf()}
+@cindex @code{Pi}
+@cindex @code{Catalan}
+@cindex @code{EulerGamma}
+@cindex @code{evalf()}
Constants behave pretty much like symbols except that they return some
specific number when the method @code{.evalf()} is called.
The predefined known constants are:
@cartouche
-@multitable @columnfractions .14 .29 .57
-@item Name @tab Common Name @tab Numerical Value (35 digits)
+@multitable @columnfractions .14 .30 .56
+@item @strong{Name} @tab @strong{Common Name} @tab @strong{Numerical Value (to 35 digits)}
@item @code{Pi}
@tab Archimedes' constant
@tab 3.14159265358979323846264338327950288
@end cartouche
-@node Fundamental operations, Built-in functions, Constants, Basic Concepts
+@node Fundamental containers, Built-in functions, Constants, Basic Concepts
@c node-name, next, previous, up
-@section Fundamental operations: the @code{power}, @code{add} and @code{mul} classes
-@cindex polynomials
+@section Fundamental containers: the @code{power}, @code{add} and @code{mul} classes
+@cindex polynomial
@cindex @code{add}
@cindex @code{mul}
@cindex @code{power}
@}
@end example
+@cindex @code{pow()}
For exponentiation, you have already seen the somewhat clumsy (though C-ish)
statement @code{pow(x,2);} to represent @code{x} squared. This direct
construction is necessary since we cannot safely overload the constructor
canonical form.
-@node Built-in functions, Important Algorithms, Fundamental operations, Basic Concepts
+@node Built-in functions, Relations, Fundamental containers, Basic Concepts
@c node-name, next, previous, up
@section Built-in functions
+@cindex functions (class @code{function})
+@cindex trigonometric function
+@cindex hyperbolic function
-There are quite a number of useful functions built into GiNaC. They are
-all objects of class @code{function}. They accept one or more
+There are quite a number of useful functions hard-wired into GiNaC. For
+instance, all trigonometric and hyperbolic functions are implemented.
+They are all objects of class @code{function}. They accept one or more
expressions as arguments and return one expression. If the arguments
are not numerical, the evaluation of the function may be halted, as it
does in the next example:
+@cindex Gamma function
+@cindex @code{subs()}
@example
#include <ginac/ginac.h>
using namespace GiNaC;
on. Read the next chapter in order to learn more about this.
-@node Important Algorithms, Polynomial Expansion, Built-in functions, Top
+@node Relations, Important Algorithms, Built-in functions, Basic Concepts
+@c node-name, next, previous, up
+@section Relations
+@cindex relations (class @code{relational})
+
+Sometimes, a relation holding between two expressions must be stored
+somehow. The class @code{relational} is a convenient container for such
+purposes. A relation is by definition a container for two @code{ex} and
+a relation between them that signals equality, inequality and so on.
+They are created by simply using the C++ operators @code{==}, @code{!=},
+@code{<}, @code{<=}, @code{>} and @code{>=} between two expressions.
+
+@xref{Built-in functions}, for examples where various applications of
+the @code{.subs()} method show how objects of class relational are used
+as arguments. There they provide an intuitive syntax for substitutions.
+
+
+@node Important Algorithms, Polynomial Expansion, Relations, Top
@c node-name, next, previous, up
@chapter Important Algorithms
-@cindex polynomials
+@cindex polynomial
In this chapter the most important algorithms provided by GiNaC will be
described. Some of them are implemented as functions on expressions,
@}
@end example
+@cindex @code{subs()}
The general rule is that wherever methods accept one or more parameters
(@var{arg1}, @var{arg2}, @dots{}) the order of arguments the function
wrapper accepts is the same but preceded by the object to act on
with GiNaC's convention. All function wrappers are always implemented
as simple inline functions which just call the corresponding method and
are only provided for users uncomfortable with OO who are dead set to
-avoid method invocations. Generally, a chain of function wrappers is
-much harder to read than a chain of methods and should therefore be
+avoid method invocations. Generally, nested function wrappers are much
+harder to read than a sequence of methods and should therefore be
avoided if possible. On the other hand, not everything in GiNaC is a
method on class @code{ex} and sometimes calling a function cannot be
avoided.
to be able to find the coefficients properly. The range of occuring
coefficients can be checked using the two methods
+@cindex @code{degree()}
+@cindex @code{ldegree()}
@example
int ex::degree(symbol const & s);
int ex::ldegree(symbol const & s);
where @code{degree()} returns the highest coefficient and
@code{ldegree()} the lowest one. (These two methods work also reliably
on non-expanded input polynomials). An application is illustrated in
-the next example, where a multivariate polynomial is analysed:
+the next example, where a multivariate polynomial is analyzed:
@example
#include <ginac/ginac.h>
@c node-name, next, previous, up
@section Symbolic Differentiation
@cindex differentiation
+@cindex @code{diff()}
@cindex chain rule
@cindex product rule
(using @code{Digits}). Some may be evaluated at certain points, but not
generally. This ought to be fixed. However, doing numerical
computations with GiNaC's quite abstract classes is doomed to be
-inefficient. For this purpose, the underlying bignum-package
-@acronym{CLN} is much better suited.
+inefficient. For this purpose, the underlying foundation classes
+provided by @acronym{CLN} are much better suited.
@node Symbolic functions, A Comparison With Other CAS, What does not belong into GiNaC, Extending GiNaC
The easiest and most instructive way to start with is probably to
implement your own function. Objects of class @code{function} are
-inserted into the system via a kind of "registry". They get a serial
+inserted into the system via a kind of `registry'. They get a serial
number that is used internally to identify them but you usually need not
worry about this. What you have to care for are functions that are
called when the user invokes certain methods. These are usual
@example
static ex cos_eval_method(ex const & x)
@{
- // if x%2*Pi return 1
- // if x%Pi return -1
- // if x%Pi/2 return 0
+ // if (!x%(2*Pi)) return 1
+ // if (!x%Pi) return -1
+ // if (!x%Pi/2) return 0
// care for other cases...
return cos(x).hold();
@}
@end example
+@cindex @code{hold()}
+@cindex evaluation
The last line returns @code{cos(x)} if we don't know what else to do and
stops a potential recursive evaluation by saying @code{.hold()}. We
should also implement a method for numerical evaluation and since we are
That's it. May the source be with you!
-@node A Comparison With Other CAS, Internal Structures, Symbolic functions, Top
+@node A Comparison With Other CAS, Advantages, Symbolic functions, Top
@c node-name, next, previous, up
@chapter A Comparison With Other CAS
@cindex advocacy
@emph{Mathematica} or @emph{Reduce}, where it has advantages and
disadvantages over these systems.
+@menu
+* Advantages:: Stengths of the GiNaC approach.
+* Disadvantages:: Weaknesses of the GiNaC approach.
+* Why C++?:: Attractiveness of C++.
+@end menu
-@heading Advantages
+@node Advantages, Disadvantages, A Comparison With Other CAS, A Comparison With Other CAS
+@c node-name, next, previous, up
+@section Advantages
GiNaC has several advantages over traditional Computer
Algebra Systems, like
@item
familiar language: all common CAS implement their own proprietary
grammar which you have to learn first (and maybe learn again when your
-vendor chooses to "enhance" it). With GiNaC you can write your program
+vendor decides to `enhance' it). With GiNaC you can write your program
in common C++, which is standardized.
+@cindex STL
@item
structured data types: you can build up structured data types using
@code{struct}s or @code{class}es together with STL features instead of
are part of your program. You can easily call third party libraries,
e.g. for numerical evaluation or graphical interaction. All other
approaches are much more cumbersome: they range from simply ignoring the
-problem (i.e. @emph{Maple}) to providing a method for "embedding" the
+problem (i.e. @emph{Maple}) to providing a method for `embedding' the
system (i.e. @emph{Yacas}).
@item
@end itemize
-@heading Disadvantages
+@node Disadvantages, Why C++?, Advantages, A Comparison With Other CAS
+@c node-name, next, previous, up
+@section Disadvantages
Of course it also has some disadvantages:
@item
not interactive: GiNaC programs have to be written in an editor,
-compiled and executed. You cannot play with expressions interactively.
+compiled and executed. You cannot play with expressions interactively.
However, such an extension is not inherently forbidden by design. In
-fact, two interactive interfaces are possible: First, a simple shell
-that exposes GiNaC's types to a command line can readily be written (and
-has been written) and second, as a more consistent approach we plan an
-integration with the @acronym{CINT} C++ interpreter.
+fact, two interactive interfaces are possible: First, a shell that
+exposes GiNaC's types to a command line can readily be written (the tiny
+@command{ginsh} that is part of the distribution being an example) and
+second, as a more consistent approach we are working on an integration
+with the @acronym{CINT} C++ interpreter.
@item
advanced features: GiNaC cannot compete with a program like
@emph{Reduce} which exists for more than 30 years now or @emph{Maple}
which grows since 1981 by the work of dozens of programmers, with
-respect to mathematical features. Integration, factorization,
+respect to mathematical features. Integration, factorization,
non-trivial simplifications, limits etc. are missing in GiNaC (and are
not planned for the near future).
compiler), the currently used version of the CLN library (fast large
integer and arbitrary precision arithmetics) can be compiled only on
systems with a recently new C++ compiler from the GNU Compiler
-Collection (@acronym{GCC}). GiNaC uses recent language features like
-explicit constructors, mutable members, RTTI, @code{dynamic_cast}s and
-STL, so ANSI compliance is meant literally. Recent @acronym{GCC}
-versions starting at 2.95, although itself not yet ANSI compliant,
-support all needed features.
+Collection (@acronym{GCC}).@footnote{This is because CLN uses
+PROVIDE/REQUIRE like macros to let the compiler gather all static
+initializations, which works for GNU C++ only.} GiNaC uses recent
+language features like explicit constructors, mutable members, RTTI,
+@code{dynamic_cast}s and STL, so ANSI compliance is meant literally.
+Recent @acronym{GCC} versions starting at 2.95, although itself not yet
+ANSI compliant, support all needed features.
@end itemize
-@heading Why C++?
+@node Why C++?, Internal Structures, Disadvantages, A Comparison With Other CAS
+@c node-name, next, previous, up
+@section Why C++?
Why did we choose to implement GiNaC in C++ instead of Java or any other
-language? C++ is not perfect: type checking is not strict (casting is
+language? C++ is not perfect: type checking is not strict (casting is
possible), separation between interface and implementation is not
complete, object oriented design is not enforced. The main reason is
-the often scolded feature of operator overloading in C++. While it may
+the often scolded feature of operator overloading in C++. While it may
be true that operating on classes with a @code{+} operator is rarely
-meaningful, it is perfectly suited for algebraic expressions. Writing
+meaningful, it is perfectly suited for algebraic expressions. Writing
@math{3x+5y} as @code{3*x+5*y} instead of
-@code{x.times(3).plus(y.times(5))} looks much more natural. Furthermore,
-the main developers are more familiar with C++ than with any other
-programming language.
+@code{x.times(3).plus(y.times(5))} looks much more natural.
+Furthermore, the main developers are more familiar with C++ than with
+any other programming language.
-@node Internal Structures, Expressions are reference counted, A Comparison With Other CAS, Top
+@node Internal Structures, Expressions are reference counted, Why C++? , Top
@c node-name, next, previous, up
@appendix Internal Structures
@cindex pair-wise representation
However, doing so results in a rather deeply nested tree which will
quickly become inefficient to manipulate. We can improve on this by
-representing the sum instead as a sequence of terms, each one being a
-pair of a purely numeric multiplicative coefficient and its rest. In
-the same spirit we can store the multiplication as a sequence of terms,
-each having a numeric exponent and a possibly complicated base, the tree
+representing the sum as a sequence of terms, each one being a pair of a
+purely numeric multiplicative coefficient and its rest. In the same
+spirit we can store the multiplication as a sequence of terms, each
+having a numeric exponent and a possibly complicated base, the tree
becomes much more flat:
@image{reppair}
coefficient. This results in the realistic picture of internal
representation for
@tex
-$2d^3 \left( 4a + 5b - 3 \right)$
+$2d^3 \left( 4a + 5b - 3 \right)$:
@end tex
@ifnottex
-@math{2*d^3*(4*a+5*b-3)}
+@math{2*d^3*(4*a+5*b-3)}:
@end ifnottex
-:
@image{repreal}
+@cindex radical
This also allows for a better handling of numeric radicals, since
@code{sqrt(2)} can now be carried along calculations. Now it should be
clear, why both classes @code{add} and @code{mul} are derived from the
Prints out the linker flags necessary to link a program against GiNaC.
@item --prefix[=@var{PREFIX}]
If @var{PREFIX} is specified, overrides the configured value of @env{$prefix}.
-(And of exec-prefix, unless --exec-prefix is also specified)
+(And of exec-prefix, unless @code{--exec-prefix} is also specified)
Otherwise, prints out the configured value of @env{$prefix}.
@item --exec-prefix[=@var{PREFIX}]
If @var{PREFIX} is specified, overrides the configured value of @env{$exec_prefix}.
Otherwise, prints out the configured value of @env{$exec_prefix}.
@end table
-Typically, @command{ginac-config} will be used within a configure script,
-as described below. It, however, can also be used directly
-from the command line to compile a simple program. For example:
+Typically, @command{ginac-config} will be used within a configure
+script, as described below. It, however, can also be used directly from
+the command line using backquotes to compile a simple program. For
+example:
@example
c++ -o simple `ginac-config --cppflags` simple.cpp `ginac-config --libs`
@env{GINACLIB_CONFIG}.
@item
-Tests the installed libraries to make sure that there version
+Tests the installed libraries to make sure that their version
is later than @var{MINIMUM-VERSION}. (A default version will be used
if not specified)
@display
If a GiNaC version greater than 0.4.0 is found, adds @env{$GINACLIB_LIBS} to
@env{$LIBS} and @env{$GINACLIB_CPPFLAGS} to @env{$CPPFLAGS}. Otherwise, dies
-with the error message "need to have GiNaC installed"
+with the error message `need to have GiNaC installed'
@end display
And the @file{Makefile.am}, which will be used to build the Makefile.
-@set UPDATED 10 December 1999
+@set UPDATED 12 December 1999
@set EDITION 0.4.1
@set VERSION 0.4.1
libginac_la_SOURCES = add.cpp basic.cpp constant.cpp diff.cpp ex.cpp \
expairseq.cpp exprseq.cpp fail.cpp function.cpp inifcns.cpp \
inifcns_trans.cpp inifcns_zeta.cpp inifcns_gamma.cpp matrix.cpp mul.cpp \
- normal.cpp numeric.cpp operators.cpp power.cpp print.cpp printraw.cpp \
- printtree.cpp printcsrc.cpp relational.cpp symbol.cpp utils.cpp series.cpp \
- ncmul.cpp clifford.cpp structure.cpp color.cpp indexed.cpp idx.cpp \
- isospin.cpp exprseq_suppl.cpp lst.cpp lst_suppl.cpp simp_lor.cpp \
- coloridx.cpp lorentzidx.cpp debugmsg.h utils.h
+ normal.cpp numeric.cpp operators.cpp power.cpp relational.cpp symbol.cpp \
+ utils.cpp series.cpp ncmul.cpp clifford.cpp structure.cpp color.cpp \
+ indexed.cpp idx.cpp isospin.cpp exprseq_suppl.cpp lst.cpp lst_suppl.cpp \
+ simp_lor.cpp coloridx.cpp lorentzidx.cpp debugmsg.h utils.h
libginac_la_LDFLAGS = -version-info $(LT_CURRENT):$(LT_REVISION):$(LT_AGE) \
-release $(LT_RELEASE)
ginacincludedir = $(includedir)/ginac
YACC = @YACC@
lib_LTLIBRARIES = libginac.la
-libginac_la_SOURCES = add.cpp basic.cpp constant.cpp diff.cpp ex.cpp expairseq.cpp exprseq.cpp fail.cpp function.cpp inifcns.cpp inifcns_trans.cpp inifcns_zeta.cpp inifcns_gamma.cpp matrix.cpp mul.cpp normal.cpp numeric.cpp operators.cpp power.cpp print.cpp printraw.cpp printtree.cpp printcsrc.cpp relational.cpp symbol.cpp utils.cpp series.cpp ncmul.cpp clifford.cpp structure.cpp color.cpp indexed.cpp idx.cpp isospin.cpp exprseq_suppl.cpp lst.cpp lst_suppl.cpp simp_lor.cpp coloridx.cpp lorentzidx.cpp debugmsg.h utils.h
+libginac_la_SOURCES = add.cpp basic.cpp constant.cpp diff.cpp ex.cpp expairseq.cpp exprseq.cpp fail.cpp function.cpp inifcns.cpp inifcns_trans.cpp inifcns_zeta.cpp inifcns_gamma.cpp matrix.cpp mul.cpp normal.cpp numeric.cpp operators.cpp power.cpp relational.cpp symbol.cpp utils.cpp series.cpp ncmul.cpp clifford.cpp structure.cpp color.cpp indexed.cpp idx.cpp isospin.cpp exprseq_suppl.cpp lst.cpp lst_suppl.cpp simp_lor.cpp coloridx.cpp lorentzidx.cpp debugmsg.h utils.h
libginac_la_LDFLAGS = -version-info $(LT_CURRENT):$(LT_REVISION):$(LT_AGE) -release $(LT_RELEASE)
libginac_la_OBJECTS = add.lo basic.lo constant.lo diff.lo ex.lo \
expairseq.lo exprseq.lo fail.lo function.lo inifcns.lo inifcns_trans.lo \
inifcns_zeta.lo inifcns_gamma.lo matrix.lo mul.lo normal.lo numeric.lo \
-operators.lo power.lo print.lo printraw.lo printtree.lo printcsrc.lo \
-relational.lo symbol.lo utils.lo series.lo ncmul.lo clifford.lo \
-structure.lo color.lo indexed.lo idx.lo isospin.lo exprseq_suppl.lo \
-lst.lo lst_suppl.lo simp_lor.lo coloridx.lo lorentzidx.lo
+operators.lo power.lo relational.lo symbol.lo utils.lo series.lo \
+ncmul.lo clifford.lo structure.lo color.lo indexed.lo idx.lo isospin.lo \
+exprseq_suppl.lo lst.lo lst_suppl.lo simp_lor.lo coloridx.lo \
+lorentzidx.lo
CXXFLAGS = @CXXFLAGS@
CXXCOMPILE = $(CXX) $(DEFS) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS)
LTCXXCOMPILE = $(LIBTOOL) --mode=compile $(CXX) $(DEFS) $(INCLUDES) $(AM_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS)
.deps/inifcns_gamma.P .deps/inifcns_trans.P .deps/inifcns_zeta.P \
.deps/isospin.P .deps/lorentzidx.P .deps/lst.P .deps/lst_suppl.P \
.deps/matrix.P .deps/mul.P .deps/ncmul.P .deps/normal.P .deps/numeric.P \
-.deps/operators.P .deps/power.P .deps/print.P .deps/printcsrc.P \
-.deps/printraw.P .deps/printtree.P .deps/relational.P .deps/series.P \
+.deps/operators.P .deps/power.P .deps/relational.P .deps/series.P \
.deps/simp_lor.P .deps/structure.P .deps/symbol.P .deps/utils.P
SOURCES = $(libginac_la_SOURCES)
OBJECTS = $(libginac_la_OBJECTS)
@for file in $(DISTFILES); do \
d=$(srcdir); \
if test -d $$d/$$file; then \
- cp -pr $$d/$$file $(distdir)/$$file; \
+ cp -pr $$/$$file $(distdir)/$$file; \
else \
test -f $(distdir)/$$file \
|| ln $$d/$$file $(distdir)/$$file 2> /dev/null \
return new add(*this);
}
+void add::print(ostream & os, unsigned upper_precedence) const
+{
+ debugmsg("add print",LOGLEVEL_PRINT);
+ if (precedence<=upper_precedence) os << "(";
+ numeric coeff;
+ bool first=true;
+ for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
+ coeff = ex_to_numeric(cit->coeff);
+ if (!first) {
+ if (coeff.csgn()==-1) os << '-'; else os << '+';
+ } else {
+ if (coeff.csgn()==-1) os << '-';
+ first=false;
+ }
+ if (!coeff.is_equal(numONE()) &&
+ !coeff.is_equal(numMINUSONE())) {
+ if (coeff.csgn()==-1)
+ (numMINUSONE()*coeff).print(os, precedence);
+ else
+ coeff.print(os, precedence);
+ os << '*';
+ }
+ os << cit->rest;
+ }
+ // print the overall numeric coefficient, if present:
+ if (!overall_coeff.is_zero()) {
+ if (overall_coeff.info(info_flags::positive)) os << '+';
+ os << overall_coeff;
+ }
+ if (precedence<=upper_precedence) os << ")";
+}
+
+void add::printraw(ostream & os) const
+{
+ debugmsg("add printraw",LOGLEVEL_PRINT);
+
+ os << "+(";
+ for (epvector::const_iterator it=seq.begin(); it!=seq.end(); ++it) {
+ os << "(";
+ (*it).rest.bp->printraw(os);
+ os << ",";
+ (*it).coeff.bp->printraw(os);
+ os << "),";
+ }
+ os << ",hash=" << hashvalue << ",flags=" << flags;
+ os << ")";
+}
+
+void add::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) const
+{
+ debugmsg("add print csrc", LOGLEVEL_PRINT);
+ if (precedence <= upper_precedence)
+ os << "(";
+
+ // Print arguments, separated by "+"
+ epvector::const_iterator it = seq.begin();
+ epvector::const_iterator itend = seq.end();
+ while (it != itend) {
+
+ // If the coefficient is -1, it is replaced by a single minus sign
+ if (it->coeff.compare(numONE()) == 0) {
+ it->rest.bp->printcsrc(os, type, precedence);
+ } else if (it->coeff.compare(numMINUSONE()) == 0) {
+ os << "-";
+ it->rest.bp->printcsrc(os, type, precedence);
+ } else if (ex_to_numeric(it->coeff).numer().compare(numONE()) == 0) {
+ it->rest.bp->printcsrc(os, type, precedence);
+ os << "/";
+ ex_to_numeric(it->coeff).denom().printcsrc(os, type, precedence);
+ } else if (ex_to_numeric(it->coeff).numer().compare(numMINUSONE()) == 0) {
+ os << "-";
+ it->rest.bp->printcsrc(os, type, precedence);
+ os << "/";
+ ex_to_numeric(it->coeff).denom().printcsrc(os, type, precedence);
+ } else {
+ it->coeff.bp->printcsrc(os, type, precedence);
+ os << "*";
+ it->rest.bp->printcsrc(os, type, precedence);
+ }
+
+ // Separator is "+", except if the following expression would have a leading minus sign
+ it++;
+ if (it != itend && !(it->coeff.compare(numZERO()) < 0 || (it->coeff.compare(numONE()) == 0 && is_ex_exactly_of_type(it->rest, numeric) && it->rest.compare(numZERO()) < 0)))
+ os << "+";
+ }
+
+ if (!overall_coeff.is_equal(exZERO())) {
+ if (overall_coeff.info(info_flags::positive)) os << '+';
+ overall_coeff.bp->printcsrc(os,type,precedence);
+ }
+
+ if (precedence <= upper_precedence)
+ os << ")";
+}
+
bool add::info(unsigned inf) const
{
// TODO: optimize
- if (inf==info_flags::polynomial || inf==info_flags::integer_polynomial || inf==info_flags::rational_polynomial || inf==info_flags::rational_function) {
+ if (inf==info_flags::polynomial ||
+ inf==info_flags::integer_polynomial ||
+ inf==info_flags::cinteger_polynomial ||
+ inf==info_flags::rational_polynomial ||
+ inf==info_flags::crational_polynomial ||
+ inf==info_flags::rational_function) {
for (epvector::const_iterator it=seq.begin(); it!=seq.end(); ++it) {
if (!(recombine_pair_to_ex(*it).info(inf)))
return false;
}
- return true;
+ return overall_coeff.info(inf);
} else {
return expairseq::info(inf);
}
// functions overriding virtual functions from bases classes
public:
basic * duplicate() const;
- void printraw(ostream & os) const;
void print(ostream & os, unsigned upper_precedence=0) const;
+ void printraw(ostream & os) const;
void printcsrc(ostream & os, unsigned type, unsigned upper_precedence=0) const;
bool info(unsigned inf) const;
int degree(symbol const & s) const;
// public
+/** Output to stream formatted to be useful as ginsh input. */
+void basic::print(ostream & os, unsigned upper_precedence) const
+{
+ debugmsg("basic print",LOGLEVEL_PRINT);
+ os << "[basic object]";
+}
+
+/** Output to stream in ugly raw format, so brave developers can have a look
+ * at the underlying structure. */
+void basic::printraw(ostream & os) const
+{
+ debugmsg("basic printraw",LOGLEVEL_PRINT);
+ os << "[basic object]";
+}
+
+/** Output to stream formatted in tree- (indented-) form, so developers can
+ * have a look at the underlying structure. */
+void basic::printtree(ostream & os, unsigned indent) const
+{
+ debugmsg("basic printtree",LOGLEVEL_PRINT);
+ os << string(indent,' ') << "type=" << typeid(*this).name()
+ << ", hash=" << hashvalue << " (0x" << hex << hashvalue << dec << ")"
+ << ", flags=" << flags
+ << ", nops=" << nops() << endl;
+ for (int i=0; i<nops(); ++i) {
+ op(i).printtree(os,indent+delta_indent);
+ }
+}
+
+/** Output to stream formatted as C-source.
+ *
+ * @param os a stream for output
+ * @param type variable type (one of the csrc_types)
+ * @param upper_precedence operator precedence of caller
+ * @see ex::printcsrc */
+void basic::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) const
+{
+ debugmsg("basic print csrc", LOGLEVEL_PRINT);
+}
+
+/** Little wrapper arount print to be called within a debugger. */
+void basic::dbgprint(void) const
+{
+ print(cerr);
+ cerr << endl;
+}
+
+/** Little wrapper arount printtree to be called within a debugger. */
+void basic::dbgprinttree(void) const
+{
+ printtree(cerr,0);
+}
+
basic * basic::duplicate() const
{
debugmsg("basic duplicate",LOGLEVEL_DUPLICATE);
// new virtual functions which can be overridden by derived classes
public: // only const functions please (may break reference counting)
virtual basic * duplicate() const;
+ virtual void print(ostream & os,unsigned upper_precedence=0) const;
virtual void printraw(ostream & os) const;
virtual void printtree(ostream & os, unsigned indent) const;
- virtual void print(ostream & os,unsigned upper_precedence=0) const;
virtual void printcsrc(ostream & os, unsigned type, unsigned upper_precedence=0) const;
virtual void dbgprint(void) const;
virtual void dbgprinttree(void) const;
return new constant(*this);
}
+void constant::print(ostream & os, unsigned upper_precedence) const
+{
+ debugmsg("constant print",LOGLEVEL_PRINT);
+ os << name;
+}
+
+void constant::printraw(ostream & os) const
+{
+ debugmsg("constant printraw",LOGLEVEL_PRINT);
+ os << "constant(" << name << ")";
+}
+
+void constant::printtree(ostream & os, unsigned indent) const
+{
+ debugmsg("constant printtree",LOGLEVEL_PRINT);
+ os << string(indent,' ') << name
+ << ", type=" << typeid(*this).name()
+ << ", hash=" << hashvalue << " (0x" << hex << hashvalue << dec << ")"
+ << ", flags=" << flags << endl;
+}
+
+void constant::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) const
+{
+ debugmsg("constant print csrc",LOGLEVEL_PRINT);
+ os << name;
+}
+
ex constant::evalf(int level) const
{
if (ef!=0) {
// functions overriding virtual functions from bases classes
public:
basic * duplicate() const;
- void printraw(ostream & os) const;
void print(ostream & os, unsigned upper_precedence=0) const;
+ void printraw(ostream & os) const;
+ void printtree(ostream & os, unsigned indent) const;
void printcsrc(ostream & os, unsigned type, unsigned upper_precedence=0) const;
ex evalf(int level=0) const;
ex diff(symbol const & s) const;
{
GINAC_ASSERT(bp!=0);
- if ( nth==0 ) {
+ if (nth==0) {
return *this;
}
ex ndiff = bp->diff(s);
- while ( nth>1 ) {
+ while (nth>1) {
ndiff = ndiff.diff(s);
--nth;
}
// public
+/** Swap the contents of two expressions. */
void ex::swap(ex & other)
{
debugmsg("ex swap",LOGLEVEL_MEMBER_FUNCTION);
other.bp=tmpbp;
}
+/** Output formatted to be useful as ginsh input. */
+void ex::print(ostream & os, unsigned upper_precedence) const
+{
+ debugmsg("ex print",LOGLEVEL_PRINT);
+ GINAC_ASSERT(bp!=0);
+ bp->print(os,upper_precedence);
+}
+
+void ex::printraw(ostream & os) const
+{
+ debugmsg("ex printraw",LOGLEVEL_PRINT);
+ GINAC_ASSERT(bp!=0);
+ os << "ex(";
+ bp->printraw(os);
+ os << ")";
+}
+
+void ex::printtree(ostream & os, unsigned indent) const
+{
+ debugmsg("ex printtree",LOGLEVEL_PRINT);
+ GINAC_ASSERT(bp!=0);
+ // os << "refcount=" << bp->refcount << " ";
+ bp->printtree(os,indent);
+}
+
+/** Print expression as a C++ statement. The output looks like
+ * "<type> <var_name> = <expression>;". The "type" parameter has an effect
+ * on how number literals are printed.
+ *
+ * @param os output stream
+ * @param type variable type (one of the csrc_types)
+ * @param var_name variable name to be printed */
+void ex::printcsrc(ostream & os, unsigned type, const char *var_name) const
+{
+ debugmsg("ex print csrc", LOGLEVEL_PRINT);
+ GINAC_ASSERT(bp!=0);
+ switch (type) {
+ case csrc_types::ctype_float:
+ os << "float ";
+ break;
+ case csrc_types::ctype_double:
+ os << "double ";
+ break;
+ case csrc_types::ctype_cl_N:
+ os << "cl_N ";
+ break;
+ }
+ os << var_name << " = ";
+ bp->printcsrc(os, type, 0);
+ os << ";\n";
+}
+
+/** Little wrapper arount print to be called within a debugger. */
+void ex::dbgprint(void) const
+{
+ debugmsg("ex dbgprint",LOGLEVEL_PRINT);
+ GINAC_ASSERT(bp!=0);
+ bp->dbgprint();
+}
+
+/** Little wrapper arount printtree to be called within a debugger. */
+void ex::dbgprinttree(void) const
+{
+ debugmsg("ex dbgprinttree",LOGLEVEL_PRINT);
+ GINAC_ASSERT(bp!=0);
+ bp->dbgprinttree();
+}
+
bool ex::info(unsigned inf) const
{
if (inf == info_flags::normal_form) {
void ex::construct_from_basic(basic const & other)
{
- if ( (other.flags & status_flags::evaluated)==0 ) {
+ if ((other.flags & status_flags::evaluated)==0) {
// cf. copy constructor
ex const & tmpex = other.eval(1); // evaluate only one (top) level
bp = tmpex.bp;
#include <algorithm>
#include <string>
#include <stdexcept>
+#include <cmath>
#include "expairseq.h"
#include "lst.h"
return new expairseq(*this);
}
+void expairseq::print(ostream & os, unsigned upper_precedence) const
+{
+ debugmsg("expairseq print",LOGLEVEL_PRINT);
+ os << "[[";
+ printseq(os,',',precedence,upper_precedence);
+ os << "]]";
+}
+
+void expairseq::printraw(ostream & os) const
+{
+ debugmsg("expairseq printraw",LOGLEVEL_PRINT);
+
+ os << "expairseq(";
+ for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
+ os << "(";
+ (*cit).rest.printraw(os);
+ os << ",";
+ (*cit).coeff.printraw(os);
+ os << "),";
+ }
+ os << ")";
+}
+
+void expairseq::printtree(ostream & os, unsigned indent) const
+{
+ debugmsg("expairseq printtree",LOGLEVEL_PRINT);
+
+ os << string(indent,' ') << "type=" << typeid(*this).name()
+ << ", hash=" << hashvalue << " (0x" << hex << hashvalue << dec << ")"
+ << ", flags=" << flags
+ << ", nops=" << nops() << endl;
+ for (unsigned i=0; i<seq.size(); ++i) {
+ seq[i].rest.printtree(os,indent+delta_indent);
+ seq[i].coeff.printtree(os,indent+delta_indent);
+ if (i!=seq.size()-1) {
+ os << string(indent+delta_indent,' ') << "-----" << endl;
+ }
+ }
+ if (!overall_coeff.is_equal(default_overall_coeff())) {
+ os << string(indent+delta_indent,' ') << "-----" << endl;
+ os << string(indent+delta_indent,' ') << "overall_coeff" << endl;
+ overall_coeff.printtree(os,indent+delta_indent);
+ }
+ os << string(indent+delta_indent,' ') << "=====" << endl;
+#ifdef EXPAIRSEQ_USE_HASHTAB
+ os << string(indent+delta_indent,' ')
+ << "hashtab size " << hashtabsize << endl;
+ if (hashtabsize==0) return;
+#define MAXCOUNT 5
+ unsigned count[MAXCOUNT+1];
+ for (int i=0; i<MAXCOUNT+1; ++i) count[i]=0;
+ unsigned this_bin_fill;
+ unsigned cum_fill_sq=0;
+ unsigned cum_fill=0;
+ for (unsigned i=0; i<hashtabsize; ++i) {
+ this_bin_fill=0;
+ if (hashtab[i].size()>0) {
+ os << string(indent+delta_indent,' ')
+ << "bin " << i << " with entries ";
+ for (epplist::const_iterator it=hashtab[i].begin();
+ it!=hashtab[i].end(); ++it) {
+ os << *it-seq.begin() << " ";
+ this_bin_fill++;
+ }
+ os << endl;
+ cum_fill += this_bin_fill;
+ cum_fill_sq += this_bin_fill*this_bin_fill;
+ }
+ if (this_bin_fill<MAXCOUNT) {
+ ++count[this_bin_fill];
+ } else {
+ ++count[MAXCOUNT];
+ }
+ }
+ unsigned fact=1;
+ double cum_prob=0;
+ double lambda=(1.0*seq.size())/hashtabsize;
+ for (int k=0; k<MAXCOUNT; ++k) {
+ if (k>0) fact *= k;
+ double prob=pow(lambda,k)/fact*exp(-lambda);
+ cum_prob += prob;
+ os << string(indent+delta_indent,' ') << "bins with " << k << " entries: "
+ << int(1000.0*count[k]/hashtabsize)/10.0 << "% (expected: "
+ << int(prob*1000)/10.0 << ")" << endl;
+ }
+ os << string(indent+delta_indent,' ') << "bins with more entries: "
+ << int(1000.0*count[MAXCOUNT]/hashtabsize)/10.0 << "% (expected: "
+ << int((1-cum_prob)*1000)/10.0 << ")" << endl;
+
+ os << string(indent+delta_indent,' ') << "variance: "
+ << 1.0/hashtabsize*cum_fill_sq-(1.0/hashtabsize*cum_fill)*(1.0/hashtabsize*cum_fill)
+ << endl;
+ os << string(indent+delta_indent,' ') << "average fill: "
+ << (1.0*cum_fill)/hashtabsize
+ << " (should be equal to " << (1.0*seq.size())/hashtabsize << ")" << endl;
+#endif // def EXPAIRSEQ_USE_HASHTAB
+}
+
bool expairseq::info(unsigned inf) const
{
return basic::info(inf);
return expairseq(vp,oc);
}
+void expairseq::printpair(ostream & os, expair const & p, unsigned upper_precedence) const
+{
+ os << "[[";
+ p.rest.bp->print(os,precedence);
+ os << ",";
+ p.coeff.bp->print(os,precedence);
+ os << "]]";
+}
+
+void expairseq::printseq(ostream & os, char delim, unsigned this_precedence,
+ unsigned upper_precedence) const
+{
+ if (this_precedence<=upper_precedence) os << "(";
+ epvector::const_iterator it,it_last;
+ it_last=seq.end();
+ --it_last;
+ for (it=seq.begin(); it!=it_last; ++it) {
+ printpair(os,*it,this_precedence);
+ os << delim;
+ }
+ printpair(os,*it,this_precedence);
+ if (!overall_coeff.is_equal(default_overall_coeff())) {
+ os << delim << overall_coeff;
+ }
+ if (this_precedence<=upper_precedence) os << ")";
+}
+
expair expairseq::split_ex_to_pair(ex const & e) const
{
return expair(e,exONE());
ex const & evaled_ex=(*cit).rest.eval(level);
if (!are_ex_trivially_equal((*cit).rest,evaled_ex)) {
- // something changed, copy seq, eval and return it
+ // something changed, copy seq, eval and return it
epvector *s=new epvector;
s->reserve(seq.size());
// functions overriding virtual functions from bases classes
public:
basic * duplicate() const;
+ void print(ostream & os, unsigned upper_precedence=0) const;
void printraw(ostream & os) const;
void printtree(ostream & os, unsigned indent) const;
- void print(ostream & os, unsigned upper_precedence=0) const;
bool info(unsigned inf) const;
int nops() const;
ex op(int const i) const;
return new fail(*this);
}
+void fail::print(ostream & os, unsigned upper_precedence) const
+{
+ debugmsg("fail print",LOGLEVEL_PRINT);
+ os << "FAIL";
+}
+
+void fail::printraw(ostream & os) const
+{
+ debugmsg("fail printraw",LOGLEVEL_PRINT);
+ os << "FAIL";
+}
+
// protected
int fail::compare_same_type(basic const & other) const
// functions overriding virtual functions from bases classes
public:
basic * duplicate() const;
- void printraw(ostream & os) const;
void print(ostream & os, unsigned upper_precedence=0) const;
+ void printraw(ostream & os) const;
protected:
int compare_same_type(basic const & other) const;
unsigned return_type(void) const { return return_types::noncommutative_composite; };
real,
rational,
integer,
+ crational,
+ cinteger,
positive,
negative,
nonnegative,
// answered by classes numeric, symbol, add, mul, power
polynomial,
integer_polynomial,
+ cinteger_polynomial,
rational_polynomial,
+ crational_polynomial,
rational_function,
// answered by class ex
/** Force inclusion of functions from initcns_gamma and inifcns_zeta
* for static lib (so ginsh will see them). */
unsigned force_include_gamma = function_index_gamma;
-unsigned force_include_zeta = function_index_zeta;
+unsigned force_include_zeta1 = function_index_zeta1;
#ifndef NO_GINAC_NAMESPACE
} // namespace GiNaC
/** Trilogarithm. */
DECLARE_FUNCTION_1P(Li3)
+// overloading at work: we cannot use the macros
/** Riemann's Zeta-function. */
-DECLARE_FUNCTION_1P(zeta)
-//DECLARE_FUNCTION_2P(zeta)
+extern const unsigned function_index_zeta1;
+inline function zeta(ex const & p1) {
+ return function(function_index_zeta1, p1);
+}
+/** Derivatives of Riemann's Zeta-function. */
+extern const unsigned function_index_zeta2;
+inline function zeta(ex const & p1, ex const & p2) {
+ return function(function_index_zeta2, p1, p2);
+}
/** Gamma-function. */
DECLARE_FUNCTION_1P(gamma)
/** Beta-function. */
DECLARE_FUNCTION_2P(beta)
+// overloading at work: we cannot use the macros
/** Psi-function (aka polygamma-function). */
-// overloading @ work: we cannot use the macros
extern const unsigned function_index_psi1;
inline function psi(ex const & p1) {
return function(function_index_psi1, p1);
}
+/** Derivatives of Psi-function (aka polygamma-functions). */
extern const unsigned function_index_psi2;
inline function psi(ex const & p1, ex const & p2) {
return function(function_index_psi2, p1, p2);
}
-//DECLARE_FUNCTION_1P(psi)
-//DECLARE_FUNCTION_2P(psi)
/** Factorial function. */
DECLARE_FUNCTION_1P(factorial)
inline bool is_order_function(ex const & e)
{
- return is_ex_the_function(e, Order);
+ return is_ex_the_function(e, Order);
}
#ifndef NO_GINAC_NAMESPACE
// Gamma-function
//////////
+static ex gamma_evalf(ex const & x)
+{
+ BEGIN_TYPECHECK
+ TYPECHECK(x,numeric)
+ END_TYPECHECK(gamma(x))
+
+ return gamma(ex_to_numeric(x));
+}
+
/** Evaluation of gamma(x). Knows about integer arguments, half-integer
* arguments and that's it. Somebody ought to provide some good numerical
* evaluation some day...
return gamma(x).hold();
}
-static ex gamma_evalf(ex const & x)
-{
- BEGIN_TYPECHECK
- TYPECHECK(x,numeric)
- END_TYPECHECK(gamma(x))
-
- return gamma(ex_to_numeric(x));
-}
-
static ex gamma_diff(ex const & x, unsigned diff_param)
{
GINAC_ASSERT(diff_param==0);
// Beta-function
//////////
+static ex beta_evalf(ex const & x, ex const & y)
+{
+ BEGIN_TYPECHECK
+ TYPECHECK(x,numeric)
+ TYPECHECK(y,numeric)
+ END_TYPECHECK(beta(x,y))
+
+ return gamma(ex_to_numeric(x))*gamma(ex_to_numeric(y))
+ / gamma(ex_to_numeric(x+y));
+}
+
static ex beta_eval(ex const & x, ex const & y)
{
if (x.info(info_flags::numeric) && y.info(info_flags::numeric)) {
return beta(x,y).hold();
}
-static ex beta_evalf(ex const & x, ex const & y)
-{
- BEGIN_TYPECHECK
- TYPECHECK(x,numeric)
- TYPECHECK(y,numeric)
- END_TYPECHECK(beta(x,y))
-
- return gamma(ex_to_numeric(x))*gamma(ex_to_numeric(y))
- / gamma(ex_to_numeric(x+y));
-}
-
static ex beta_diff(ex const & x, ex const & y, unsigned diff_param)
{
GINAC_ASSERT(diff_param<2);
// Psi-function (aka polygamma-function)
//////////
+static ex psi1_evalf(ex const & x)
+{
+ BEGIN_TYPECHECK
+ TYPECHECK(x,numeric)
+ END_TYPECHECK(psi(x))
+
+ return psi(ex_to_numeric(x));
+}
+
/** Evaluation of polygamma-function psi(x).
* Somebody ought to provide some good numerical evaluation some day... */
static ex psi1_eval(ex const & x)
return psi(x).hold();
}
-static ex psi1_evalf(ex const & x)
-{
- BEGIN_TYPECHECK
- TYPECHECK(x,numeric)
- END_TYPECHECK(psi(x))
-
- return psi(ex_to_numeric(x));
-}
-
static ex psi1_diff(ex const & x, unsigned diff_param)
{
GINAC_ASSERT(diff_param==0);
// Psi-functions (aka polygamma-functions) psi(0,x)==psi(x)
//////////
+static ex psi2_evalf(ex const & n, ex const & x)
+{
+ BEGIN_TYPECHECK
+ TYPECHECK(n,numeric)
+ TYPECHECK(x,numeric)
+ END_TYPECHECK(psi(n,x))
+
+ return psi(ex_to_numeric(n), ex_to_numeric(x));
+}
+
/** Evaluation of polygamma-function psi(n,x).
* Somebody ought to provide some good numerical evaluation some day... */
static ex psi2_eval(ex const & n, ex const & x)
return psi(n, x).hold();
}
-static ex psi2_evalf(ex const & n, ex const & x)
-{
- BEGIN_TYPECHECK
- TYPECHECK(n,numeric)
- TYPECHECK(x,numeric)
- END_TYPECHECK(psi(n,x))
-
- return psi(ex_to_numeric(n), ex_to_numeric(x));
-}
-
static ex psi2_diff(ex const & n, ex const & x, unsigned diff_param)
{
GINAC_ASSERT(diff_param<2);
#include "constant.h"
#include "numeric.h"
#include "power.h"
+#include "relational.h"
+#include "symbol.h"
+#include "utils.h"
#ifndef NO_GINAC_NAMESPACE
namespace GiNaC {
return x.op(0);
// exp(float)
- if (x.info(info_flags::numeric) && !x.info(info_flags::rational))
+ if (x.info(info_flags::numeric) && !x.info(info_flags::crational))
return exp_evalf(x);
return exp(x).hold();
-}
+}
static ex exp_diff(ex const & x, unsigned diff_param)
{
if (x.is_equal(exZERO()))
throw(std::domain_error("log_eval(): log(0)"));
// log(float)
- if (!x.info(info_flags::rational))
+ if (!x.info(info_flags::crational))
return log_evalf(x);
}
+ // log(exp(t)) -> t (for real-valued t):
+ if (is_ex_the_function(x, exp)) {
+ ex t=x.op(0);
+ if (t.info(info_flags::real))
+ return t;
+ }
return log(x).hold();
-}
+}
static ex log_diff(ex const & x, unsigned diff_param)
{
static ex sin_eval(ex const & x)
{
- // sin(n*Pi) -> 0
ex xOverPi=x/Pi;
- if (xOverPi.info(info_flags::integer))
- return exZERO();
-
- // sin((2n+1)*Pi/2) -> {+|-}1
- ex xOverPiMinusHalf=xOverPi-exHALF();
- if (xOverPiMinusHalf.info(info_flags::even))
- return exONE();
- else if (xOverPiMinusHalf.info(info_flags::odd))
- return exMINUSONE();
+ if (xOverPi.info(info_flags::numeric)) {
+ // sin(n*Pi) -> 0
+ if (xOverPi.info(info_flags::integer))
+ return exZERO();
+
+ // sin((2n+1)*Pi/2) -> {+|-}1
+ ex xOverPiMinusHalf=xOverPi-exHALF();
+ if (xOverPiMinusHalf.info(info_flags::even))
+ return exONE();
+ else if (xOverPiMinusHalf.info(info_flags::odd))
+ return exMINUSONE();
+ }
if (is_ex_exactly_of_type(x, function)) {
ex t=x.op(0);
}
// sin(float) -> float
- if (x.info(info_flags::numeric) && !x.info(info_flags::rational))
+ if (x.info(info_flags::numeric) && !x.info(info_flags::crational))
return sin_evalf(x);
return sin(x).hold();
static ex cos_eval(ex const & x)
{
- // cos(n*Pi) -> {+|-}1
ex xOverPi=x/Pi;
- if (xOverPi.info(info_flags::even))
- return exONE();
- else if (xOverPi.info(info_flags::odd))
- return exMINUSONE();
-
- // cos((2n+1)*Pi/2) -> 0
- ex xOverPiMinusHalf=xOverPi-exHALF();
- if (xOverPiMinusHalf.info(info_flags::integer))
- return exZERO();
+ if (xOverPi.info(info_flags::numeric)) {
+ // cos(n*Pi) -> {+|-}1
+ if (xOverPi.info(info_flags::even))
+ return exONE();
+ else if (xOverPi.info(info_flags::odd))
+ return exMINUSONE();
+
+ // cos((2n+1)*Pi/2) -> 0
+ ex xOverPiMinusHalf=xOverPi-exHALF();
+ if (xOverPiMinusHalf.info(info_flags::integer))
+ return exZERO();
+ }
if (is_ex_exactly_of_type(x, function)) {
ex t=x.op(0);
}
// cos(float) -> float
- if (x.info(info_flags::numeric) && !x.info(info_flags::rational))
+ if (x.info(info_flags::numeric) && !x.info(info_flags::crational))
return cos_evalf(x);
return cos(x).hold();
}
// tan(float) -> float
- if (x.info(info_flags::numeric) && !x.info(info_flags::rational)) {
+ if (x.info(info_flags::numeric) && !x.info(info_flags::crational)) {
return tan_evalf(x);
}
return (1+power(tan(x),exTWO()));
}
-REGISTER_FUNCTION(tan, tan_eval, tan_evalf, tan_diff, NULL);
+static ex tan_series(ex const & x, symbol const & s, ex const & point, int order)
+{
+ // method:
+ // Taylor series where there is no pole falls back to tan_diff.
+ // On a pole simply expand sin(x)/cos(x).
+ ex xpoint = x.subs(s==point);
+ if (!(2*xpoint/Pi).info(info_flags::odd))
+ throw do_taylor();
+ // if we got here we have to care for a simple pole
+ return (sin(x)/cos(x)).series(s, point, order+2);
+}
+
+REGISTER_FUNCTION(tan, tan_eval, tan_evalf, tan_diff, tan_series);
//////////
// inverse sine (arc sine)
if (x.is_equal(exMINUSONE()))
return numeric(-1,2)*Pi;
// asin(float) -> float
- if (!x.info(info_flags::rational))
+ if (!x.info(info_flags::crational))
return asin_evalf(x);
}
if (x.is_equal(exMINUSONE()))
return Pi;
// acos(float) -> float
- if (!x.info(info_flags::rational))
+ if (!x.info(info_flags::crational))
return acos_evalf(x);
}
if (x.is_equal(exZERO()))
return exZERO();
// atan(float) -> float
- if (!x.info(info_flags::rational))
+ if (!x.info(info_flags::crational))
return atan_evalf(x);
}
static ex atan2_eval(ex const & y, ex const & x)
{
- if (y.info(info_flags::numeric) && !y.info(info_flags::rational) &&
- x.info(info_flags::numeric) && !x.info(info_flags::rational)) {
+ if (y.info(info_flags::numeric) && !y.info(info_flags::crational) &&
+ x.info(info_flags::numeric) && !x.info(info_flags::crational)) {
return atan2_evalf(y,x);
}
if (x.is_zero())
return exZERO();
// sinh(float) -> float
- if (!x.info(info_flags::rational))
+ if (!x.info(info_flags::crational))
return sinh_evalf(x);
}
+ ex xOverPiI=x/Pi/I;
+ if (xOverPiI.info(info_flags::numeric)) {
+ // sinh(n*Pi*I) -> 0
+ if (xOverPiI.info(info_flags::integer))
+ return exZERO();
+
+ // sinh((2n+1)*Pi*I/2) -> {+|-}I
+ ex xOverPiIMinusHalf=xOverPiI-exHALF();
+ if (xOverPiIMinusHalf.info(info_flags::even))
+ return I;
+ else if (xOverPiIMinusHalf.info(info_flags::odd))
+ return -I;
+ }
+
if (is_ex_exactly_of_type(x, function)) {
ex t=x.op(0);
// sinh(asinh(x)) -> x
if (x.is_zero())
return exONE();
// cosh(float) -> float
- if (!x.info(info_flags::rational))
+ if (!x.info(info_flags::crational))
return cosh_evalf(x);
}
+ ex xOverPiI=x/Pi/I;
+ if (xOverPiI.info(info_flags::numeric)) {
+ // cosh(n*Pi*I) -> {+|-}1
+ if (xOverPiI.info(info_flags::even))
+ return exONE();
+ else if (xOverPiI.info(info_flags::odd))
+ return exMINUSONE();
+
+ // cosh((2n+1)*Pi*I/2) -> 0
+ ex xOverPiIMinusHalf=xOverPiI-exHALF();
+ if (xOverPiIMinusHalf.info(info_flags::integer))
+ return exZERO();
+ }
+
if (is_ex_exactly_of_type(x, function)) {
ex t=x.op(0);
// cosh(acosh(x)) -> x
if (x.is_zero())
return exZERO();
// tanh(float) -> float
- if (!x.info(info_flags::rational))
+ if (!x.info(info_flags::crational))
return tanh_evalf(x);
}
return exONE()-power(tanh(x),exTWO());
}
-REGISTER_FUNCTION(tanh, tanh_eval, tanh_evalf, tanh_diff, NULL);
+static ex tanh_series(ex const & x, symbol const & s, ex const & point, int order)
+{
+ // method:
+ // Taylor series where there is no pole falls back to tanh_diff.
+ // On a pole simply expand sinh(x)/cosh(x).
+ ex xpoint = x.subs(s==point);
+ if (!(2*I*xpoint/Pi).info(info_flags::odd))
+ throw do_taylor();
+ // if we got here we have to care for a simple pole
+ return (sinh(x)/cosh(x)).series(s, point, order+2);
+}
+
+REGISTER_FUNCTION(tanh, tanh_eval, tanh_evalf, tanh_diff, tanh_series);
//////////
// inverse hyperbolic sine (trigonometric function)
if (x.is_zero())
return exZERO();
// asinh(float) -> float
- if (!x.info(info_flags::rational))
+ if (!x.info(info_flags::crational))
return asinh_evalf(x);
}
if (x.is_equal(exMINUSONE()))
return Pi*I;
// acosh(float) -> float
- if (!x.info(info_flags::rational))
+ if (!x.info(info_flags::crational))
return acosh_evalf(x);
}
if (x.is_equal(exONE()) || x.is_equal(exONE()))
throw (std::domain_error("atanh_eval(): infinity"));
// atanh(float) -> float
- if (!x.info(info_flags::rational))
+ if (!x.info(info_flags::crational))
return atanh_evalf(x);
}
// Riemann's Zeta-function
//////////
-static ex zeta_eval(ex const & x)
+static ex zeta1_evalf(ex const & x)
+{
+ BEGIN_TYPECHECK
+ TYPECHECK(x,numeric)
+ END_TYPECHECK(zeta(x))
+
+ return zeta(ex_to_numeric(x));
+}
+
+static ex zeta1_eval(ex const & x)
{
if (x.info(info_flags::numeric)) {
numeric y = ex_to_numeric(x);
return zeta(x).hold();
}
-static ex zeta_evalf(ex const & x)
+static ex zeta1_diff(ex const & x, unsigned diff_param)
{
- BEGIN_TYPECHECK
- TYPECHECK(x,numeric)
- END_TYPECHECK(zeta(x))
+ GINAC_ASSERT(diff_param==0);
- return zeta(ex_to_numeric(x));
+ return zeta(exONE(), x);
+}
+
+const unsigned function_index_zeta1 = function::register_new("zeta", zeta1_eval, zeta1_evalf, zeta1_diff, NULL);
+
+//////////
+// Derivatives of Riemann's Zeta-function zeta(0,x)==zeta(x)
+//////////
+
+static ex zeta2_eval(ex const & n, ex const & x)
+{
+ if (n.info(info_flags::numeric)) {
+ // zeta(0,x) -> zeta(x)
+ if (n.is_zero())
+ return zeta(x);
+ }
+
+ return zeta(n, x).hold();
+}
+
+static ex zeta2_diff(ex const & n, ex const & x, unsigned diff_param)
+{
+ GINAC_ASSERT(diff_param<2);
+
+ if (diff_param==0) {
+ // d/dn zeta(n,x)
+ throw(std::logic_error("cannot diff zeta(n,x) with respect to n"));
+ }
+ // d/dx psi(n,x)
+ return zeta(n+1,x);
}
-REGISTER_FUNCTION(zeta, zeta_eval, zeta_evalf, NULL, NULL);
+const unsigned function_index_zeta2 = function::register_new("zeta", zeta2_eval, NULL, zeta2_diff, NULL);
#ifndef NO_GINAC_NAMESPACE
} // namespace GiNaC
return new matrix(*this);
}
+void matrix::print(ostream & os, unsigned upper_precedence) const
+{
+ debugmsg("matrix print",LOGLEVEL_PRINT);
+ os << "[[ ";
+ for (int r=0; r<row-1; ++r) {
+ os << "[[";
+ for (int c=0; c<col-1; ++c) {
+ os << m[r*col+c] << ",";
+ }
+ os << m[col*(r+1)-1] << "]], ";
+ }
+ os << "[[";
+ for (int c=0; c<col-1; ++c) {
+ os << m[(row-1)*col+c] << ",";
+ }
+ os << m[row*col-1] << "]] ]]";
+}
+
+void matrix::printraw(ostream & os) const
+{
+ debugmsg("matrix printraw",LOGLEVEL_PRINT);
+ os << "matrix(" << row << "," << col <<",";
+ for (int r=0; r<row-1; ++r) {
+ os << "(";
+ for (int c=0; c<col-1; ++c) {
+ os << m[r*col+c] << ",";
+ }
+ os << m[col*(r-1)-1] << "),";
+ }
+ os << "(";
+ for (int c=0; c<col-1; ++c) {
+ os << m[(row-1)*col+c] << ",";
+ }
+ os << m[row*col-1] << "))";
+}
+
/** nops is defined to be rows x columns. */
int matrix::nops() const
{
// functions overriding virtual functions from bases classes
public:
basic * duplicate() const;
- void printraw(ostream & os) const;
void print(ostream & os, unsigned upper_precedence=0) const;
+ void printraw(ostream & os) const;
int nops() const;
ex & let_op(int const i);
ex expand(unsigned options=0) const;
return new mul(*this);
}
+void mul::print(ostream & os, unsigned upper_precedence) const
+{
+ debugmsg("mul print",LOGLEVEL_PRINT);
+ if (precedence<=upper_precedence) os << "(";
+ bool first=true;
+ // first print the overall numeric coefficient:
+ if (ex_to_numeric(overall_coeff).csgn()==-1) os << '-';
+ if (!overall_coeff.is_equal(exONE()) &&
+ !overall_coeff.is_equal(exMINUSONE())) {
+ if (ex_to_numeric(overall_coeff).csgn()==-1)
+ (numMINUSONE()*overall_coeff).print(os, precedence);
+ else
+ overall_coeff.print(os, precedence);
+ os << '*';
+ }
+ // then proceed with the remaining factors:
+ for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
+ if (!first) {
+ os << '*';
+ } else {
+ first=false;
+ }
+ recombine_pair_to_ex(*cit).print(os,precedence);
+ }
+ if (precedence<=upper_precedence) os << ")";
+}
+
+void mul::printraw(ostream & os) const
+{
+ debugmsg("mul printraw",LOGLEVEL_PRINT);
+
+ os << "*(";
+ for (epvector::const_iterator it=seq.begin(); it!=seq.end(); ++it) {
+ os << "(";
+ (*it).rest.bp->printraw(os);
+ os << ",";
+ (*it).coeff.bp->printraw(os);
+ os << "),";
+ }
+ os << ",hash=" << hashvalue << ",flags=" << flags;
+ os << ")";
+}
+
+void mul::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) const
+{
+ debugmsg("mul print csrc", LOGLEVEL_PRINT);
+ if (precedence <= upper_precedence)
+ os << "(";
+
+ if (!overall_coeff.is_equal(exONE())) {
+ overall_coeff.bp->printcsrc(os,type,precedence);
+ os << "*";
+ }
+
+ // Print arguments, separated by "*" or "/"
+ epvector::const_iterator it = seq.begin();
+ epvector::const_iterator itend = seq.end();
+ while (it != itend) {
+
+ // If the first argument is a negative integer power, it gets printed as "1.0/<expr>"
+ if (it == seq.begin() && ex_to_numeric(it->coeff).is_integer() && it->coeff.compare(numZERO()) < 0) {
+ if (type == csrc_types::ctype_cl_N)
+ os << "recip(";
+ else
+ os << "1.0/";
+ }
+
+ // If the exponent is 1 or -1, it is left out
+ if (it->coeff.compare(exONE()) == 0 || it->coeff.compare(numMINUSONE()) == 0)
+ it->rest.bp->printcsrc(os, type, precedence);
+ else
+ // outer parens around ex needed for broken gcc-2.95 parser:
+ (ex(power(it->rest, abs(ex_to_numeric(it->coeff))))).bp->printcsrc(os, type, upper_precedence);
+
+ // Separator is "/" for negative integer powers, "*" otherwise
+ it++;
+ if (it != itend) {
+ if (ex_to_numeric(it->coeff).is_integer() && it->coeff.compare(numZERO()) < 0)
+ os << "/";
+ else
+ os << "*";
+ }
+ }
+ if (precedence <= upper_precedence)
+ os << ")";
+}
+
bool mul::info(unsigned inf) const
{
// TODO: optimize
- if (inf==info_flags::polynomial || inf==info_flags::integer_polynomial || inf==info_flags::rational_polynomial || inf==info_flags::rational_function) {
+ if (inf==info_flags::polynomial ||
+ inf==info_flags::integer_polynomial ||
+ inf==info_flags::cinteger_polynomial ||
+ inf==info_flags::rational_polynomial ||
+ inf==info_flags::crational_polynomial ||
+ inf==info_flags::rational_function) {
for (epvector::const_iterator it=seq.begin(); it!=seq.end(); ++it) {
if (!(recombine_pair_to_ex(*it).info(inf)))
return false;
}
- return true;
+ return overall_coeff.info(inf);
} else {
return expairseq::info(inf);
}
all_commutative=0;
}
if ((rt==return_types::noncommutative)&&(!all_commutative)) {
- // another nc element found, compare type_infos
+ // another nc element found, compare type_infos
if ((*cit_noncommutative_element).rest.return_type_tinfo()!=(*cit).rest.return_type_tinfo()) {
- // diffent types -> mul is ncc
- return return_types::noncommutative_composite;
+ // diffent types -> mul is ncc
+ return return_types::noncommutative_composite;
}
}
}
// functions overriding virtual functions from bases classes
public:
basic * duplicate() const;
- void printraw(ostream & os) const;
void print(ostream & os, unsigned upper_precedence) const;
+ void printraw(ostream & os) const;
void printcsrc(ostream & os, unsigned type, unsigned upper_precedence) const;
bool info(unsigned inf) const;
int degree(symbol const & s) const;
return new ncmul(*this);
}
+void ncmul::print(ostream & os, unsigned upper_precedence) const
+{
+ debugmsg("ncmul print",LOGLEVEL_PRINT);
+ printseq(os,'(','%',')',precedence,upper_precedence);
+}
+
+void ncmul::printraw(ostream & os) const
+{
+ debugmsg("ncmul printraw",LOGLEVEL_PRINT);
+
+ os << "%(";
+ for (exvector::const_iterator it=seq.begin(); it!=seq.end(); ++it) {
+ (*it).bp->printraw(os);
+ os << ",";
+ }
+ os << ",hash=" << hashvalue << ",flags=" << flags;
+ os << ")";
+}
+
+void ncmul::printcsrc(ostream & os, unsigned upper_precedence) const
+{
+ debugmsg("ncmul print csrc",LOGLEVEL_PRINT);
+ exvector::const_iterator it;
+ exvector::const_iterator itend = seq.end()-1;
+ os << "ncmul(";
+ for (it=seq.begin(); it!=itend; ++it) {
+ (*it).bp->printcsrc(os,precedence);
+ os << ",";
+ }
+ (*it).bp->printcsrc(os,precedence);
+ os << ")";
+}
+
bool ncmul::info(unsigned inf) const
{
throw(std::logic_error("which flags have to be implemented in ncmul::info()?"));
// functions overriding virtual functions from bases classes
public:
basic * duplicate() const;
- void printraw(ostream & os) const;
void print(ostream & os, unsigned upper_precedence) const;
+ void printraw(ostream & os) const;
void printcsrc(ostream & os, unsigned upper_precedence) const;
bool info(unsigned inf) const;
int degree(symbol const & s) const;
if (is_real())
if (is_rational())
return *this;
- else
- return replace_with_symbol(*this, sym_lst, repl_lst);
+ else
+ return replace_with_symbol(*this, sym_lst, repl_lst);
else { // complex
numeric re = real(), im = imag();
- ex re_ex = re.is_rational() ? re : replace_with_symbol(re, sym_lst, repl_lst);
- ex im_ex = im.is_rational() ? im : replace_with_symbol(im, sym_lst, repl_lst);
- return re_ex + im_ex * replace_with_symbol(I, sym_lst, repl_lst);
- }
+ ex re_ex = re.is_rational() ? re : replace_with_symbol(re, sym_lst, repl_lst);
+ ex im_ex = im.is_rational() ? im : replace_with_symbol(im, sym_lst, repl_lst);
+ return re_ex + im_ex * replace_with_symbol(I, sym_lst, repl_lst);
+ }
}
/*
* Helper function for fraction cancellation (returns cancelled fraction n/d)
*/
-
static ex frac_cancel(const ex &n, const ex &d)
{
ex num = n;
}
}
+void numeric::printtree(ostream & os, unsigned indent) const
+{
+ debugmsg("numeric printtree", LOGLEVEL_PRINT);
+ os << string(indent,' ') << *value
+ << " (numeric): "
+ << "hash=" << hashvalue << " (0x" << hex << hashvalue << dec << ")"
+ << ", flags=" << flags << endl;
+}
+
+void numeric::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) const
+{
+ debugmsg("numeric print csrc", LOGLEVEL_PRINT);
+ ios::fmtflags oldflags = os.flags();
+ os.setf(ios::scientific);
+ if (is_rational() && !is_integer()) {
+ if (compare(numZERO()) > 0) {
+ os << "(";
+ if (type == csrc_types::ctype_cl_N)
+ os << "cl_F(\"" << numer().evalf() << "\")";
+ else
+ os << numer().to_double();
+ } else {
+ os << "-(";
+ if (type == csrc_types::ctype_cl_N)
+ os << "cl_F(\"" << -numer().evalf() << "\")";
+ else
+ os << -numer().to_double();
+ }
+ os << "/";
+ if (type == csrc_types::ctype_cl_N)
+ os << "cl_F(\"" << denom().evalf() << "\")";
+ else
+ os << denom().to_double();
+ os << ")";
+ } else {
+ if (type == csrc_types::ctype_cl_N)
+ os << "cl_F(\"" << evalf() << "\")";
+ else
+ os << to_double();
+ }
+ os.flags(oldflags);
+}
+
bool numeric::info(unsigned inf) const
{
switch (inf) {
case info_flags::rational:
case info_flags::rational_polynomial:
return is_rational();
+ case info_flags::crational:
+ case info_flags::crational_polynomial:
+ return is_crational();
case info_flags::integer:
case info_flags::integer_polynomial:
return is_integer();
+ case info_flags::cinteger:
+ case info_flags::cinteger_polynomial:
+ return is_cinteger();
case info_flags::positive:
return is_positive();
case info_flags::negative:
* currently set.
*
* @param level ignored, but needed for overriding basic::evalf.
- * @return an ex-handle to a numeric. */
+ * @return an ex-handle to a numeric. */
ex numeric::evalf(int level) const
{
// level can safely be discarded for numeric objects.
* @exception overflow_error (division by zero) */
numeric numeric::div(numeric const & other) const
{
- if (zerop(*other.value))
+ if (::zerop(*other.value))
throw (std::overflow_error("division by zero"));
return numeric((*value)/(*other.value));
}
numeric numeric::power(numeric const & other) const
{
static const numeric * numONEp=&numONE();
- if (&other==numONEp) {
+ if (&other==numONEp)
return *this;
- }
- if (zerop(*value) && other.is_real() && minusp(realpart(*other.value)))
+ if (::zerop(*value) && other.is_real() && ::minusp(realpart(*other.value)))
throw (std::overflow_error("division by zero"));
- return numeric(expt(*value,*other.value));
+ return numeric(::expt(*value,*other.value));
}
/** Inverse of a number. */
numeric numeric::inverse(void) const
{
- return numeric(recip(*value)); // -> CLN
+ return numeric(::recip(*value)); // -> CLN
}
numeric const & numeric::add_dyn(numeric const & other) const
numeric const & numeric::div_dyn(numeric const & other) const
{
- if (zerop(*other.value))
+ if (::zerop(*other.value))
throw (std::overflow_error("division by zero"));
return static_cast<numeric const &>((new numeric((*value)/(*other.value)))->
setflag(status_flags::dynallocated));
numeric const & numeric::power_dyn(numeric const & other) const
{
static const numeric * numONEp=&numONE();
- if (&other==numONEp) {
+ if (&other==numONEp)
return *this;
- }
- // The ifs are only a workaround for a bug in CLN. It gets stuck otherwise:
- if ( !other.is_integer() &&
- other.is_rational() &&
- (*this).is_nonneg_integer() ) {
- if ( !zerop(*value) ) {
- return static_cast<numeric const &>((new numeric(exp(*other.value * log(*value))))->
- setflag(status_flags::dynallocated));
- } else {
- if ( !zerop(*other.value) ) { // 0^(n/m)
- return static_cast<numeric const &>((new numeric(0))->
- setflag(status_flags::dynallocated));
- } else { // raise FPE (0^0 requested)
- return static_cast<numeric const &>((new numeric(1/(*other.value)))->
- setflag(status_flags::dynallocated));
- }
- }
- } else { // default -> CLN
- return static_cast<numeric const &>((new numeric(expt(*value,*other.value)))->
- setflag(status_flags::dynallocated));
- }
+ if (::zerop(*value) && other.is_real() && ::minusp(realpart(*other.value)))
+ throw (std::overflow_error("division by zero"));
+ return static_cast<numeric const &>((new numeric(::expt(*value,*other.value)))->
+ setflag(status_flags::dynallocated));
}
numeric const & numeric::operator=(int i)
{
if (is_zero())
return 0;
- if (!zerop(realpart(*value))) {
- if (plusp(realpart(*value)))
+ if (!::zerop(realpart(*value))) {
+ if (::plusp(realpart(*value)))
return 1;
else
return -1;
} else {
- if (plusp(imagpart(*value)))
+ if (::plusp(imagpart(*value)))
return 1;
else
return -1;
// Comparing two real numbers?
if (is_real() && other.is_real())
// Yes, just compare them
- return cl_compare(The(cl_R)(*value), The(cl_R)(*other.value));
+ return ::cl_compare(The(cl_R)(*value), The(cl_R)(*other.value));
else {
// No, first compare real parts
- cl_signean real_cmp = cl_compare(realpart(*value), realpart(*other.value));
+ cl_signean real_cmp = ::cl_compare(realpart(*value), realpart(*other.value));
if (real_cmp)
return real_cmp;
- return cl_compare(imagpart(*value), imagpart(*other.value));
+ return ::cl_compare(imagpart(*value), imagpart(*other.value));
}
}
/** True if object is zero. */
bool numeric::is_zero(void) const
{
- return zerop(*value); // -> CLN
+ return ::zerop(*value); // -> CLN
}
/** True if object is not complex and greater than zero. */
bool numeric::is_positive(void) const
{
- if (is_real()) {
- return plusp(The(cl_R)(*value)); // -> CLN
- }
+ if (is_real())
+ return ::plusp(The(cl_R)(*value)); // -> CLN
return false;
}
/** True if object is not complex and less than zero. */
bool numeric::is_negative(void) const
{
- if (is_real()) {
- return minusp(The(cl_R)(*value)); // -> CLN
- }
+ if (is_real())
+ return ::minusp(The(cl_R)(*value)); // -> CLN
return false;
}
/** True if object is a non-complex integer. */
bool numeric::is_integer(void) const
{
- return instanceof(*value, cl_I_ring); // -> CLN
+ return ::instanceof(*value, cl_I_ring); // -> CLN
}
/** True if object is an exact integer greater than zero. */
bool numeric::is_pos_integer(void) const
{
- return (is_integer() &&
- plusp(The(cl_I)(*value))); // -> CLN
+ return (is_integer() && ::plusp(The(cl_I)(*value))); // -> CLN
}
/** True if object is an exact integer greater or equal zero. */
bool numeric::is_nonneg_integer(void) const
{
- return (is_integer() &&
- !minusp(The(cl_I)(*value))); // -> CLN
+ return (is_integer() && !::minusp(The(cl_I)(*value))); // -> CLN
}
/** True if object is an exact even integer. */
bool numeric::is_even(void) const
{
- return (is_integer() &&
- evenp(The(cl_I)(*value))); // -> CLN
+ return (is_integer() && ::evenp(The(cl_I)(*value))); // -> CLN
}
/** True if object is an exact odd integer. */
bool numeric::is_odd(void) const
{
- return (is_integer() &&
- oddp(The(cl_I)(*value))); // -> CLN
+ return (is_integer() && ::oddp(The(cl_I)(*value))); // -> CLN
}
/** Probabilistic primality test.
* @return true if object is exact integer and prime. */
bool numeric::is_prime(void) const
{
- return (is_integer() &&
- isprobprime(The(cl_I)(*value))); // -> CLN
+ return (is_integer() && ::isprobprime(The(cl_I)(*value))); // -> CLN
}
/** True if object is an exact rational number, may even be complex
* (denominator may be unity). */
bool numeric::is_rational(void) const
{
- return instanceof(*value, cl_RA_ring);
+ return ::instanceof(*value, cl_RA_ring); // -> CLN
}
/** True if object is a real integer, rational or float (but not complex). */
bool numeric::is_real(void) const
{
- return instanceof(*value, cl_R_ring); // -> CLN
+ return ::instanceof(*value, cl_R_ring); // -> CLN
}
bool numeric::operator==(numeric const & other) const
* of the form a+b*I, where a and b are integers. */
bool numeric::is_cinteger(void) const
{
- if (instanceof(*value, cl_I_ring))
+ if (::instanceof(*value, cl_I_ring))
return true;
else if (!is_real()) { // complex case, handle n+m*I
- if (instanceof(realpart(*value), cl_I_ring) &&
- instanceof(imagpart(*value), cl_I_ring))
+ if (::instanceof(realpart(*value), cl_I_ring) &&
+ ::instanceof(imagpart(*value), cl_I_ring))
return true;
}
return false;
* (denominator may be unity). */
bool numeric::is_crational(void) const
{
- if (instanceof(*value, cl_RA_ring))
+ if (::instanceof(*value, cl_RA_ring))
return true;
else if (!is_real()) { // complex case, handle Q(i):
- if (instanceof(realpart(*value), cl_RA_ring) &&
- instanceof(imagpart(*value), cl_RA_ring))
+ if (::instanceof(realpart(*value), cl_RA_ring) &&
+ ::instanceof(imagpart(*value), cl_RA_ring))
return true;
}
return false;
* @exception invalid_argument (complex inequality) */
bool numeric::operator<(numeric const & other) const
{
- if ( is_real() && other.is_real() ) {
+ if (is_real() && other.is_real())
return (bool)(The(cl_R)(*value) < The(cl_R)(*other.value)); // -> CLN
- }
throw (std::invalid_argument("numeric::operator<(): complex inequality"));
return false; // make compiler shut up
}
* @exception invalid_argument (complex inequality) */
bool numeric::operator<=(numeric const & other) const
{
- if ( is_real() && other.is_real() ) {
+ if (is_real() && other.is_real())
return (bool)(The(cl_R)(*value) <= The(cl_R)(*other.value)); // -> CLN
- }
throw (std::invalid_argument("numeric::operator<=(): complex inequality"));
return false; // make compiler shut up
}
* @exception invalid_argument (complex inequality) */
bool numeric::operator>(numeric const & other) const
{
- if ( is_real() && other.is_real() ) {
+ if (is_real() && other.is_real())
return (bool)(The(cl_R)(*value) > The(cl_R)(*other.value)); // -> CLN
- }
throw (std::invalid_argument("numeric::operator>(): complex inequality"));
return false; // make compiler shut up
}
* @exception invalid_argument (complex inequality) */
bool numeric::operator>=(numeric const & other) const
{
- if ( is_real() && other.is_real() ) {
+ if (is_real() && other.is_real())
return (bool)(The(cl_R)(*value) >= The(cl_R)(*other.value)); // -> CLN
- }
throw (std::invalid_argument("numeric::operator>=(): complex inequality"));
return false; // make compiler shut up
}
int numeric::to_int(void) const
{
GINAC_ASSERT(is_integer());
- return cl_I_to_int(The(cl_I)(*value));
+ return ::cl_I_to_int(The(cl_I)(*value)); // -> CLN
}
/** Converts numeric types to machine's double. You should check with is_real()
double numeric::to_double(void) const
{
GINAC_ASSERT(is_real());
- return cl_double_approx(realpart(*value));
+ return ::cl_double_approx(realpart(*value)); // -> CLN
}
/** Real part of a number. */
numeric numeric::real(void) const
{
- return numeric(realpart(*value)); // -> CLN
+ return numeric(::realpart(*value)); // -> CLN
}
/** Imaginary part of a number. */
numeric numeric::imag(void) const
{
- return numeric(imagpart(*value)); // -> CLN
+ return numeric(::imagpart(*value)); // -> CLN
}
#ifndef SANE_LINKER
return numeric(*this);
}
#ifdef SANE_LINKER
- else if (instanceof(*value, cl_RA_ring)) {
- return numeric(numerator(The(cl_RA)(*value)));
+ else if (::instanceof(*value, cl_RA_ring)) {
+ return numeric(::numerator(The(cl_RA)(*value)));
}
else if (!is_real()) { // complex case, handle Q(i):
- cl_R r = realpart(*value);
- cl_R i = imagpart(*value);
- if (instanceof(r, cl_I_ring) && instanceof(i, cl_I_ring))
+ cl_R r = ::realpart(*value);
+ cl_R i = ::imagpart(*value);
+ if (::instanceof(r, cl_I_ring) && ::instanceof(i, cl_I_ring))
return numeric(*this);
- if (instanceof(r, cl_I_ring) && instanceof(i, cl_RA_ring))
- return numeric(complex(r*denominator(The(cl_RA)(i)), numerator(The(cl_RA)(i))));
- if (instanceof(r, cl_RA_ring) && instanceof(i, cl_I_ring))
- return numeric(complex(numerator(The(cl_RA)(r)), i*denominator(The(cl_RA)(r))));
- if (instanceof(r, cl_RA_ring) && instanceof(i, cl_RA_ring)) {
- cl_I s = lcm(denominator(The(cl_RA)(r)), denominator(The(cl_RA)(i)));
- return numeric(complex(numerator(The(cl_RA)(r))*(exquo(s,denominator(The(cl_RA)(r)))),
- numerator(The(cl_RA)(i))*(exquo(s,denominator(The(cl_RA)(i))))));
+ if (::instanceof(r, cl_I_ring) && ::instanceof(i, cl_RA_ring))
+ return numeric(complex(r*::denominator(The(cl_RA)(i)), ::numerator(The(cl_RA)(i))));
+ if (::instanceof(r, cl_RA_ring) && ::instanceof(i, cl_I_ring))
+ return numeric(complex(::numerator(The(cl_RA)(r)), i*::denominator(The(cl_RA)(r))));
+ if (::instanceof(r, cl_RA_ring) && ::instanceof(i, cl_RA_ring)) {
+ cl_I s = lcm(::denominator(The(cl_RA)(r)), ::denominator(The(cl_RA)(i)));
+ return numeric(complex(::numerator(The(cl_RA)(r))*(exquo(s,::denominator(The(cl_RA)(r)))),
+ ::numerator(The(cl_RA)(i))*(exquo(s,::denominator(The(cl_RA)(i))))));
}
}
#else
}
#ifdef SANE_LINKER
if (instanceof(*value, cl_RA_ring)) {
- return numeric(denominator(The(cl_RA)(*value)));
+ return numeric(::denominator(The(cl_RA)(*value)));
}
if (!is_real()) { // complex case, handle Q(i):
cl_R r = realpart(*value);
cl_R i = imagpart(*value);
- if (instanceof(r, cl_I_ring) && instanceof(i, cl_I_ring))
+ if (::instanceof(r, cl_I_ring) && ::instanceof(i, cl_I_ring))
return numONE();
- if (instanceof(r, cl_I_ring) && instanceof(i, cl_RA_ring))
- return numeric(denominator(The(cl_RA)(i)));
- if (instanceof(r, cl_RA_ring) && instanceof(i, cl_I_ring))
- return numeric(denominator(The(cl_RA)(r)));
- if (instanceof(r, cl_RA_ring) && instanceof(i, cl_RA_ring))
- return numeric(lcm(denominator(The(cl_RA)(r)), denominator(The(cl_RA)(i))));
+ if (::instanceof(r, cl_I_ring) && ::instanceof(i, cl_RA_ring))
+ return numeric(::denominator(The(cl_RA)(i)));
+ if (::instanceof(r, cl_RA_ring) && ::instanceof(i, cl_I_ring))
+ return numeric(::denominator(The(cl_RA)(r)));
+ if (::instanceof(r, cl_RA_ring) && ::instanceof(i, cl_RA_ring))
+ return numeric(lcm(::denominator(The(cl_RA)(r)), ::denominator(The(cl_RA)(i))));
}
#else
if (instanceof(*value, cl_RA_ring)) {
* in two's complement if it is an integer, 0 otherwise. */
int numeric::int_length(void) const
{
- if (is_integer()) {
- return integer_length(The(cl_I)(*value)); // -> CLN
- } else {
+ if (is_integer())
+ return ::integer_length(The(cl_I)(*value)); // -> CLN
+ else
return 0;
- }
}
* integer arguments. */
numeric zeta(numeric const & x)
{
- if (x.is_integer())
- return ::cl_zeta(x.to_int()); // -> CLN
- else
- clog << "zeta(): Does anybody know good way to calculate this numerically?" << endl;
+ // A dirty hack to allow for things like zeta(3.0), since CLN currently
+ // only knows about integer arguments and zeta(3).evalf() automatically
+ // cascades down to zeta(3.0).evalf(). The trick is to rely on 3.0-3
+ // being an exact zero for CLN, which can be tested and then we can just
+ // pass the number casted to an int:
+ if (x.is_real()) {
+ int aux = (int)(::cl_double_approx(realpart(*x.value)));
+ if (zerop(*x.value-aux))
+ return ::cl_zeta(aux); // -> CLN
+ }
+ clog << "zeta(" << x
+ << "): Does anybody know good way to calculate this numerically?"
+ << endl;
return numeric(0);
}
* This is only a stub! */
numeric gamma(numeric const & x)
{
- clog << "gamma(): Does anybody know good way to calculate this numerically?" << endl;
+ clog << "gamma(" << x
+ << "): Does anybody know good way to calculate this numerically?"
+ << endl;
return numeric(0);
}
* This is only a stub! */
numeric psi(numeric const & x)
{
- clog << "psi(): Does anybody know good way to calculate this numerically?" << endl;
+ clog << "psi(" << x
+ << "): Does anybody know good way to calculate this numerically?"
+ << endl;
return numeric(0);
}
* This is only a stub! */
numeric psi(numeric const & n, numeric const & x)
{
- clog << "psi(): Does anybody know good way to calculate this numerically?" << endl;
+ clog << "psi(" << n << "," << x
+ << "): Does anybody know good way to calculate this numerically?"
+ << endl;
return numeric(0);
}
* @exception range_error (argument must be integer >= 0) */
numeric factorial(numeric const & nn)
{
- if ( !nn.is_nonneg_integer() ) {
+ if (!nn.is_nonneg_integer())
throw (std::range_error("numeric::factorial(): argument must be integer >= 0"));
- }
-
return numeric(::factorial(nn.to_int())); // -> CLN
}
return numZERO();
} else {
return numMINUSONE().power(k)*binomial(k-n-numONE(),k);
- }
+ }
}
// should really be gamma(n+1)/(gamma(r+1)/gamma(n-r+1) or a suitable limit
* integer, 0 otherwise. */
numeric mod(numeric const & a, numeric const & b)
{
- if (a.is_integer() && b.is_integer()) {
+ if (a.is_integer() && b.is_integer())
return ::mod(The(cl_I)(*a.value), The(cl_I)(*b.value)); // -> CLN
- }
- else {
+ else
return numZERO(); // Throw?
- }
}
/** Modulus (in symmetric representation).
if (a.is_integer() && b.is_integer()) {
cl_I b2 = The(cl_I)(ceiling1(The(cl_I)(*b.value) / 2)) - 1;
return ::mod(The(cl_I)(*a.value) + b2, The(cl_I)(*b.value)) - b2;
- } else {
+ } else
return numZERO(); // Throw?
- }
}
/** Numeric integer remainder.
* @return remainder of a/b if both are integer, 0 otherwise. */
numeric irem(numeric const & a, numeric const & b)
{
- if (a.is_integer() && b.is_integer()) {
+ if (a.is_integer() && b.is_integer())
return ::rem(The(cl_I)(*a.value), The(cl_I)(*b.value)); // -> CLN
- }
- else {
+ else
return numZERO(); // Throw?
- }
}
/** Numeric integer remainder.
* @return truncated quotient of a/b if both are integer, 0 otherwise. */
numeric iquo(numeric const & a, numeric const & b)
{
- if (a.is_integer() && b.is_integer()) {
+ if (a.is_integer() && b.is_integer())
return truncate1(The(cl_I)(*a.value), The(cl_I)(*b.value)); // -> CLN
- } else {
+ else
return numZERO(); // Throw?
- }
}
/** Numeric integer quotient.
/** Integer numeric square root. */
numeric isqrt(numeric const & x)
{
- if (x.is_integer()) {
- cl_I root;
- ::isqrt(The(cl_I)(*x.value), &root); // -> CLN
- return root;
- } else
- return numZERO(); // Throw?
+ if (x.is_integer()) {
+ cl_I root;
+ ::isqrt(The(cl_I)(*x.value), &root); // -> CLN
+ return root;
+ } else
+ return numZERO(); // Throw?
}
/** Greatest Common Divisor.
numeric gcd(numeric const & a, numeric const & b)
{
if (a.is_integer() && b.is_integer())
- return ::gcd(The(cl_I)(*a.value), The(cl_I)(*b.value)); // -> CLN
+ return ::gcd(The(cl_I)(*a.value), The(cl_I)(*b.value)); // -> CLN
else
return numONE();
}
numeric lcm(numeric const & a, numeric const & b)
{
if (a.is_integer() && b.is_integer())
- return ::lcm(The(cl_I)(*a.value), The(cl_I)(*b.value)); // -> CLN
+ return ::lcm(The(cl_I)(*a.value), The(cl_I)(*b.value)); // -> CLN
else
return *a.value * *b.value;
}
friend numeric asinh(numeric const & x);
friend numeric acosh(numeric const & x);
friend numeric atanh(numeric const & x);
+ friend numeric zeta(numeric const & x);
friend numeric bernoulli(numeric const & n);
friend numeric abs(numeric const & x);
friend numeric mod(numeric const & a, numeric const & b);
// functions overriding virtual functions from bases classes
public:
basic * duplicate() const;
+ void print(ostream & os, unsigned precedence=0) const;
void printraw(ostream & os) const;
void printtree(ostream & os, unsigned indent) const;
- void print(ostream & os, unsigned precedence=0) const;
void printcsrc(ostream & os, unsigned type, unsigned precedence=0) const;
bool info(unsigned inf) const;
ex evalf(int level=0) const;
inline bool is_real(numeric const & x)
{ return x.is_real(); }
+inline bool is_cinteger(numeric const & x)
+{ return x.is_cinteger(); }
+
+inline bool is_crational(numeric const & x)
+{ return x.is_crational(); }
+
inline numeric real(numeric const & x)
{ return x.real(); }
// utility functions
inline const numeric &ex_to_numeric(const ex &e)
{
- return static_cast<const numeric &>(*e.bp);
+ return static_cast<const numeric &>(*e.bp);
}
#ifndef NO_GINAC_NAMESPACE
ex operator/(ex const & lh, ex const & rh);
ex operator%(ex const & lh, ex const & rh); // non-commutative multiplication
-/*
-
-// binary arithmetic operators ex with numeric
-ex operator+(ex const & lh, numeric const & rh);
-ex operator-(ex const & lh, numeric const & rh);
-ex operator*(ex const & lh, numeric const & rh);
-ex operator/(ex const & lh, numeric const & rh);
-ex operator%(ex const & lh, numeric const & rh); // non-commutative multiplication
-
-// binary arithmetic operators numeric with ex
-ex operator+(numeric const & lh, ex const & rh);
-ex operator-(numeric const & lh, ex const & rh);
-ex operator*(numeric const & lh, ex const & rh);
-ex operator/(numeric const & lh, ex const & rh);
-ex operator%(numeric const & lh, ex const & rh); // non-commutative multiplication
-
-*/
-
// binary arithmetic operators numeric with numeric
numeric operator+(numeric const & lh, numeric const & rh);
numeric operator-(numeric const & lh, numeric const & rh);
ex const & operator/=(ex & lh, ex const & rh);
ex const & operator%=(ex & lh, ex const & rh); // non-commutative multiplication
-/*
-
-// binary arithmetic assignment operators with numeric
-ex const & operator+=(ex & lh, numeric const & rh);
-ex const & operator-=(ex & lh, numeric const & rh);
-ex const & operator*=(ex & lh, numeric const & rh);
-ex const & operator/=(ex & lh, numeric const & rh);
-ex const & operator%=(ex & lh, numeric const & rh); // non-commutative multiplication
-
-*/
-
// binary arithmetic assignment operators with numeric
numeric const & operator+=(numeric & lh, numeric const & rh);
numeric const & operator-=(numeric & lh, numeric const & rh);
relational operator>(ex const & lh, ex const & rh);
relational operator>=(ex const & lh, ex const & rh);
-/*
-
-// binary relational operators ex with numeric
-relational operator==(ex const & lh, numeric const & rh);
-relational operator!=(ex const & lh, numeric const & rh);
-relational operator<(ex const & lh, numeric const & rh);
-relational operator<=(ex const & lh, numeric const & rh);
-relational operator>(ex const & lh, numeric const & rh);
-relational operator>=(ex const & lh, numeric const & rh);
-
-// binary relational operators numeric with ex
-relational operator==(numeric const & lh, ex const & rh);
-relational operator!=(numeric const & lh, ex const & rh);
-relational operator<(numeric const & lh, ex const & rh);
-relational operator<=(numeric const & lh, ex const & rh);
-relational operator>(numeric const & lh, ex const & rh);
-relational operator>=(numeric const & lh, ex const & rh);
-
-*/
-
// input/output stream operators
ostream & operator<<(ostream & os, ex const & e);
istream & operator>>(istream & is, ex & e);
return new power(*this);
}
+void power::print(ostream & os, unsigned upper_precedence) const
+{
+ debugmsg("power print",LOGLEVEL_PRINT);
+ if (precedence<=upper_precedence) os << "(";
+ basis.print(os,precedence);
+ os << "^";
+ exponent.print(os,precedence);
+ if (precedence<=upper_precedence) os << ")";
+}
+
+void power::printraw(ostream & os) const
+{
+ debugmsg("power printraw",LOGLEVEL_PRINT);
+
+ os << "power(";
+ basis.printraw(os);
+ os << ",";
+ exponent.printraw(os);
+ os << ",hash=" << hashvalue << ",flags=" << flags << ")";
+}
+
+void power::printtree(ostream & os, unsigned indent) const
+{
+ debugmsg("power printtree",LOGLEVEL_PRINT);
+
+ os << string(indent,' ') << "power: "
+ << "hash=" << hashvalue << " (0x" << hex << hashvalue << dec << ")"
+ << ", flags=" << flags << endl;
+ basis.printtree(os,indent+delta_indent);
+ exponent.printtree(os,indent+delta_indent);
+}
+
+static void print_sym_pow(ostream & os, unsigned type, const symbol &x, int exp)
+{
+ // Optimal output of integer powers of symbols to aid compiler CSE
+ if (exp == 1) {
+ x.printcsrc(os, type, 0);
+ } else if (exp == 2) {
+ x.printcsrc(os, type, 0);
+ os << "*";
+ x.printcsrc(os, type, 0);
+ } else if (exp & 1) {
+ x.printcsrc(os, 0);
+ os << "*";
+ print_sym_pow(os, type, x, exp-1);
+ } else {
+ os << "(";
+ print_sym_pow(os, type, x, exp >> 1);
+ os << ")*(";
+ print_sym_pow(os, type, x, exp >> 1);
+ os << ")";
+ }
+}
+
+void power::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) const
+{
+ debugmsg("power print csrc", LOGLEVEL_PRINT);
+
+ // Integer powers of symbols are printed in a special, optimized way
+ if (exponent.info(info_flags::integer) &&
+ (is_ex_exactly_of_type(basis, symbol) ||
+ is_ex_exactly_of_type(basis, constant))) {
+ int exp = ex_to_numeric(exponent).to_int();
+ if (exp > 0)
+ os << "(";
+ else {
+ exp = -exp;
+ if (type == csrc_types::ctype_cl_N)
+ os << "recip(";
+ else
+ os << "1.0/(";
+ }
+ print_sym_pow(os, type, static_cast<const symbol &>(*basis.bp), exp);
+ os << ")";
+
+ // <expr>^-1 is printed as "1.0/<expr>" or with the recip() function of CLN
+ } else if (exponent.compare(numMINUSONE()) == 0) {
+ if (type == csrc_types::ctype_cl_N)
+ os << "recip(";
+ else
+ os << "1.0/(";
+ basis.bp->printcsrc(os, type, 0);
+ os << ")";
+
+ // Otherwise, use the pow() or expt() (CLN) functions
+ } else {
+ if (type == csrc_types::ctype_cl_N)
+ os << "expt(";
+ else
+ os << "pow(";
+ basis.bp->printcsrc(os, type, 0);
+ os << ",";
+ exponent.bp->printcsrc(os, type, 0);
+ os << ")";
+ }
+}
+
bool power::info(unsigned inf) const
{
- if (inf==info_flags::polynomial || inf==info_flags::integer_polynomial || inf==info_flags::rational_polynomial) {
+ if (inf==info_flags::polynomial ||
+ inf==info_flags::integer_polynomial ||
+ inf==info_flags::cinteger_polynomial ||
+ inf==info_flags::rational_polynomial ||
+ inf==info_flags::crational_polynomial) {
return exponent.info(info_flags::nonnegint);
} else if (inf==info_flags::rational_function) {
return exponent.info(info_flags::integer);
if (basis_is_numerical && exponent_is_numerical) {
// ^(c1,c2) -> c1^c2 (c1, c2 numeric(),
// except if c1,c2 are rational, but c1^c2 is not)
- bool basis_is_rational = num_basis->is_rational();
- bool exponent_is_rational = num_exponent->is_rational();
+ bool basis_is_crational = num_basis->is_crational();
+ bool exponent_is_crational = num_exponent->is_crational();
numeric res = (*num_basis).power(*num_exponent);
- if ((!basis_is_rational || !exponent_is_rational)
- || res.is_rational()) {
+ if ((!basis_is_crational || !exponent_is_crational)
+ || res.is_crational()) {
return res;
}
GINAC_ASSERT(!num_exponent->is_integer()); // has been handled by now
// ^(c1,n/m) -> *(c1^q,c1^(n/m-q)), 0<(n/m-h)<1, q integer
- if (basis_is_rational && exponent_is_rational
+ if (basis_is_crational && exponent_is_crational
&& num_exponent->is_real()
&& !num_exponent->is_integer()) {
numeric r, q, n, m;
// functions overriding virtual functions from bases classes
public:
basic * duplicate() const;
+ void print(ostream & os, unsigned upper_precedence=0) const;
void printraw(ostream & os) const;
void printtree(ostream & os, unsigned indent) const;
- void print(ostream & os, unsigned upper_precedence=0) const;
void printcsrc(ostream & os, unsigned type, unsigned upper_precedence=0) const;
bool info(unsigned inf) const;
int nops() const;
+++ /dev/null
-/** @file print.cpp
- *
- * The methods .print() are responsible for the nice default-output of
- * objects. All related helper-functions go in here as well. */
-
-/*
- * GiNaC Copyright (C) 1999 Johannes Gutenberg University Mainz, Germany
- *
- * This program is free software; you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation; either version 2 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program; if not, write to the Free Software
- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
- */
-
-#include <iostream>
-
-#include "basic.h"
-#include "ex.h"
-#include "add.h"
-#include "constant.h"
-#include "expairseq.h"
-#include "fail.h"
-#include "indexed.h"
-#include "inifcns.h"
-#include "matrix.h"
-#include "mul.h"
-#include "ncmul.h"
-#include "numeric.h"
-#include "power.h"
-#include "relational.h"
-#include "series.h"
-#include "symbol.h"
-#include "debugmsg.h"
-
-#ifndef NO_GINAC_NAMESPACE
-namespace GiNaC {
-#endif // ndef NO_GINAC_NAMESPACE
-
-void ex::print(ostream & os, unsigned upper_precedence) const
-{
- debugmsg("ex print",LOGLEVEL_PRINT);
- GINAC_ASSERT(bp!=0);
- bp->print(os,upper_precedence);
-}
-
-void ex::dbgprint(void) const
-{
- debugmsg("ex dbgprint",LOGLEVEL_PRINT);
- GINAC_ASSERT(bp!=0);
- bp->dbgprint();
-}
-
-void basic::print(ostream & os, unsigned upper_precedence) const
-{
- debugmsg("basic print",LOGLEVEL_PRINT);
- os << "[basic object]";
-}
-
-void basic::dbgprint(void) const
-{
- print(cerr);
- cerr << endl;
-}
-
-void symbol::print(ostream & os, unsigned upper_precedence) const
-{
- debugmsg("symbol print",LOGLEVEL_PRINT);
- os << name;
-}
-
-void constant::print(ostream & os, unsigned upper_precedence) const
-{
- debugmsg("constant print",LOGLEVEL_PRINT);
- os << name;
-}
-
-void power::print(ostream & os, unsigned upper_precedence) const
-{
- debugmsg("power print",LOGLEVEL_PRINT);
- if (precedence<=upper_precedence) os << "(";
- basis.print(os,precedence);
- os << "^";
- exponent.print(os,precedence);
- if (precedence<=upper_precedence) os << ")";
-}
-
-void fail::print(ostream & os, unsigned upper_precedence) const
-{
- debugmsg("fail print",LOGLEVEL_PRINT);
- os << "FAIL";
-}
-
-void expairseq::printpair(ostream & os, expair const & p, unsigned upper_precedence) const
-{
- os << "[[";
- p.rest.bp->print(os,precedence);
- os << ",";
- p.coeff.bp->print(os,precedence);
- os << "]]";
-}
-
-void expairseq::printseq(ostream & os, char delim, unsigned this_precedence,
- unsigned upper_precedence) const
-{
- if (this_precedence<=upper_precedence) os << "(";
- epvector::const_iterator it,it_last;
- it_last=seq.end();
- --it_last;
- for (it=seq.begin(); it!=it_last; ++it) {
- printpair(os,*it,this_precedence);
- os << delim;
- }
- printpair(os,*it,this_precedence);
- if (!overall_coeff.is_equal(default_overall_coeff())) {
- os << delim << overall_coeff;
- }
- if (this_precedence<=upper_precedence) os << ")";
-}
-
-void expairseq::print(ostream & os, unsigned upper_precedence) const
-{
- debugmsg("expairseq print",LOGLEVEL_PRINT);
- os << "[[";
- printseq(os,',',precedence,upper_precedence);
- os << "]]";
-}
-
-void add::print(ostream & os, unsigned upper_precedence) const
-{
- debugmsg("add print",LOGLEVEL_PRINT);
- if (precedence<=upper_precedence) os << "(";
- numeric coeff;
- bool first=true;
- for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
- coeff = ex_to_numeric(cit->coeff);
- if (!first) {
- if (coeff.csgn()==-1) os << '-'; else os << '+';
- } else {
- if (coeff.csgn()==-1) os << '-';
- first=false;
- }
- if (!coeff.is_equal(numONE()) &&
- !coeff.is_equal(numMINUSONE())) {
- if (coeff.csgn()==-1)
- (numMINUSONE()*coeff).print(os, precedence);
- else
- coeff.print(os, precedence);
- os << '*';
- }
- os << cit->rest;
- }
- // print the overall numeric coefficient, if present:
- if (!overall_coeff.is_zero()) {
- if (overall_coeff > 0) os << '+';
- os << overall_coeff;
- }
- if (precedence<=upper_precedence) os << ")";
-}
-
-void mul::print(ostream & os, unsigned upper_precedence) const
-{
- debugmsg("mul print",LOGLEVEL_PRINT);
- if (precedence<=upper_precedence) os << "(";
- bool first=true;
- // first print the overall numeric coefficient:
- if (ex_to_numeric(overall_coeff).csgn()==-1) os << '-';
- if (!overall_coeff.is_equal(exONE()) &&
- !overall_coeff.is_equal(exMINUSONE())) {
- if (ex_to_numeric(overall_coeff).csgn()==-1)
- (numMINUSONE()*overall_coeff).print(os, precedence);
- else
- overall_coeff.print(os, precedence);
- os << '*';
- }
- // then proceed with the remaining factors:
- for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
- if (!first) {
- os << '*';
- } else {
- first=false;
- }
- recombine_pair_to_ex(*cit).print(os,precedence);
- }
- if (precedence<=upper_precedence) os << ")";
-}
-
-void ncmul::print(ostream & os, unsigned upper_precedence) const
-{
- debugmsg("ncmul print",LOGLEVEL_PRINT);
- printseq(os,'(','%',')',precedence,upper_precedence);
-}
-
-/*void function::print(ostream & os, unsigned upper_precedence) const
- *{
- * debugmsg("function print",LOGLEVEL_PRINT);
- * os << name;
- * printseq(os,'(',',',')',exprseq::precedence,function::precedence);
- *}*/
-
-void series::print(ostream &os, unsigned upper_precedence) const
-{
- debugmsg("symbol print", LOGLEVEL_PRINT);
- convert_to_poly().print(os, upper_precedence);
-}
-
-void relational::print(ostream & os, unsigned upper_precedence) const
-{
- debugmsg("relational print",LOGLEVEL_PRINT);
- if (precedence<=upper_precedence) os << "(";
- lh.print(os,precedence);
- switch (o) {
- case equal:
- os << "==";
- break;
- case not_equal:
- os << "!=";
- break;
- case less:
- os << "<";
- break;
- case less_or_equal:
- os << "<=";
- break;
- case greater:
- os << ">";
- break;
- case greater_or_equal:
- os << ">=";
- break;
- default:
- os << "(INVALID RELATIONAL OPERATOR)";
- }
- rh.print(os,precedence);
- if (precedence<=upper_precedence) os << ")";
-}
-
-void matrix::print(ostream & os, unsigned upper_precedence) const
-{
- debugmsg("matrix print",LOGLEVEL_PRINT);
- os << "[[ ";
- for (int r=0; r<row-1; ++r) {
- os << "[[";
- for (int c=0; c<col-1; ++c) {
- os << m[r*col+c] << ",";
- }
- os << m[col*(r+1)-1] << "]], ";
- }
- os << "[[";
- for (int c=0; c<col-1; ++c) {
- os << m[(row-1)*col+c] << ",";
- }
- os << m[row*col-1] << "]] ]]";
-}
-
-#ifndef NO_GINAC_NAMESPACE
-} // namespace GiNaC
-#endif // ndef NO_GINAC_NAMESPACE
+++ /dev/null
-/** @file printcsrc.cpp
- *
- * The methods .printcsrc() are responsible for C-source output of
- * objects. All related helper-functions go in here as well. */
-
-/*
- * GiNaC Copyright (C) 1999 Johannes Gutenberg University Mainz, Germany
- *
- * This program is free software; you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation; either version 2 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program; if not, write to the Free Software
- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
- */
-
-#include <iostream>
-
-#include "basic.h"
-#include "ex.h"
-#include "add.h"
-#include "constant.h"
-#include "expairseq.h"
-#include "indexed.h"
-#include "inifcns.h"
-#include "mul.h"
-#include "ncmul.h"
-#include "numeric.h"
-#include "power.h"
-#include "relational.h"
-#include "series.h"
-#include "symbol.h"
-#include "debugmsg.h"
-
-#ifndef NO_GINAC_NAMESPACE
-namespace GiNaC {
-#endif // ndef NO_GINAC_NAMESPACE
-
-/** Print expression as a C++ statement. The output looks like
- * "<type> <var_name> = <expression>;". The "type" parameter has an effect
- * on how number literals are printed.
- *
- * @param os output stream
- * @param type variable type (one of the csrc_types)
- * @param var_name variable name to be printed */
-void ex::printcsrc(ostream & os, unsigned type, const char *var_name) const
-{
- debugmsg("ex print csrc", LOGLEVEL_PRINT);
- GINAC_ASSERT(bp!=0);
- switch (type) {
- case csrc_types::ctype_float:
- os << "float ";
- break;
- case csrc_types::ctype_double:
- os << "double ";
- break;
- case csrc_types::ctype_cl_N:
- os << "cl_N ";
- break;
- }
- os << var_name << " = ";
- bp->printcsrc(os, type, 0);
- os << ";\n";
-}
-
-void basic::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) const
-{
- debugmsg("basic print csrc", LOGLEVEL_PRINT);
-}
-
-void numeric::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) const
-{
- debugmsg("numeric print csrc", LOGLEVEL_PRINT);
- ios::fmtflags oldflags = os.flags();
- os.setf(ios::scientific);
- if (is_rational() && !is_integer()) {
- if (compare(numZERO()) > 0) {
- os << "(";
- if (type == csrc_types::ctype_cl_N)
- os << "cl_F(\"" << numer().evalf() << "\")";
- else
- os << numer().to_double();
- } else {
- os << "-(";
- if (type == csrc_types::ctype_cl_N)
- os << "cl_F(\"" << -numer().evalf() << "\")";
- else
- os << -numer().to_double();
- }
- os << "/";
- if (type == csrc_types::ctype_cl_N)
- os << "cl_F(\"" << denom().evalf() << "\")";
- else
- os << denom().to_double();
- os << ")";
- } else {
- if (type == csrc_types::ctype_cl_N)
- os << "cl_F(\"" << evalf() << "\")";
- else
- os << to_double();
- }
- os.flags(oldflags);
-}
-
-void symbol::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) const
-{
- debugmsg("symbol print csrc", LOGLEVEL_PRINT);
- os << name;
-}
-
-void constant::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) const
-{
- debugmsg("constant print csrc",LOGLEVEL_PRINT);
- os << name;
-}
-
-static void print_sym_pow(ostream & os, unsigned type, const symbol &x, int exp)
-{
- // Optimal output of integer powers of symbols to aid compiler CSE
- if (exp == 1) {
- x.printcsrc(os, type, 0);
- } else if (exp == 2) {
- x.printcsrc(os, type, 0);
- os << "*";
- x.printcsrc(os, type, 0);
- } else if (exp & 1) {
- x.printcsrc(os, 0);
- os << "*";
- print_sym_pow(os, type, x, exp-1);
- } else {
- os << "(";
- print_sym_pow(os, type, x, exp >> 1);
- os << ")*(";
- print_sym_pow(os, type, x, exp >> 1);
- os << ")";
- }
-}
-
-void power::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) const
-{
- debugmsg("power print csrc", LOGLEVEL_PRINT);
-
- // Integer powers of symbols are printed in a special, optimized way
- if (exponent.info(info_flags::integer) &&
- (is_ex_exactly_of_type(basis, symbol) ||
- is_ex_exactly_of_type(basis, constant))) {
- int exp = ex_to_numeric(exponent).to_int();
- if (exp > 0)
- os << "(";
- else {
- exp = -exp;
- if (type == csrc_types::ctype_cl_N)
- os << "recip(";
- else
- os << "1.0/(";
- }
- print_sym_pow(os, type, static_cast<const symbol &>(*basis.bp), exp);
- os << ")";
-
- // <expr>^-1 is printed as "1.0/<expr>" or with the recip() function of CLN
- } else if (exponent.compare(numMINUSONE()) == 0) {
- if (type == csrc_types::ctype_cl_N)
- os << "recip(";
- else
- os << "1.0/(";
- basis.bp->printcsrc(os, type, 0);
- os << ")";
-
- // Otherwise, use the pow() or expt() (CLN) functions
- } else {
- if (type == csrc_types::ctype_cl_N)
- os << "expt(";
- else
- os << "pow(";
- basis.bp->printcsrc(os, type, 0);
- os << ",";
- exponent.bp->printcsrc(os, type, 0);
- os << ")";
- }
-}
-
-void add::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) const
-{
- debugmsg("add print csrc", LOGLEVEL_PRINT);
- if (precedence <= upper_precedence)
- os << "(";
-
- // Print arguments, separated by "+"
- epvector::const_iterator it = seq.begin();
- epvector::const_iterator itend = seq.end();
- while (it != itend) {
-
- // If the coefficient is -1, it is replaced by a single minus sign
- if (it->coeff.compare(numONE()) == 0) {
- it->rest.bp->printcsrc(os, type, precedence);
- } else if (it->coeff.compare(numMINUSONE()) == 0) {
- os << "-";
- it->rest.bp->printcsrc(os, type, precedence);
- } else if (ex_to_numeric(it->coeff).numer().compare(numONE()) == 0) {
- it->rest.bp->printcsrc(os, type, precedence);
- os << "/";
- ex_to_numeric(it->coeff).denom().printcsrc(os, type, precedence);
- } else if (ex_to_numeric(it->coeff).numer().compare(numMINUSONE()) == 0) {
- os << "-";
- it->rest.bp->printcsrc(os, type, precedence);
- os << "/";
- ex_to_numeric(it->coeff).denom().printcsrc(os, type, precedence);
- } else {
- it->coeff.bp->printcsrc(os, type, precedence);
- os << "*";
- it->rest.bp->printcsrc(os, type, precedence);
- }
-
- // Separator is "+", except it the following expression would have a leading minus sign
- it++;
- if (it != itend && !(it->coeff.compare(numZERO()) < 0 || (it->coeff.compare(numONE()) == 0 && is_ex_exactly_of_type(it->rest, numeric) && it->rest.compare(numZERO()) < 0)))
- os << "+";
- }
-
- if (!overall_coeff.is_equal(exZERO())) {
- if (overall_coeff > 0) os << '+';
- overall_coeff.bp->printcsrc(os,type,precedence);
- }
-
- if (precedence <= upper_precedence)
- os << ")";
-}
-
-void mul::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) const
-{
- debugmsg("mul print csrc", LOGLEVEL_PRINT);
- if (precedence <= upper_precedence)
- os << "(";
-
- if (!overall_coeff.is_equal(exONE())) {
- overall_coeff.bp->printcsrc(os,type,precedence);
- os << "*";
- }
-
- // Print arguments, separated by "*" or "/"
- epvector::const_iterator it = seq.begin();
- epvector::const_iterator itend = seq.end();
- while (it != itend) {
-
- // If the first argument is a negative integer power, it gets printed as "1.0/<expr>"
- if (it == seq.begin() && ex_to_numeric(it->coeff).is_integer() && it->coeff.compare(numZERO()) < 0) {
- if (type == csrc_types::ctype_cl_N)
- os << "recip(";
- else
- os << "1.0/";
- }
-
- // If the exponent is 1 or -1, it is left out
- if (it->coeff.compare(exONE()) == 0 || it->coeff.compare(numMINUSONE()) == 0)
- it->rest.bp->printcsrc(os, type, precedence);
- else
- // outer parens around ex needed for broken gcc-2.95 parser:
- (ex(power(it->rest, abs(ex_to_numeric(it->coeff))))).bp->printcsrc(os, type, upper_precedence);
-
- // Separator is "/" for negative integer powers, "*" otherwise
- it++;
- if (it != itend) {
- if (ex_to_numeric(it->coeff).is_integer() && it->coeff.compare(numZERO()) < 0)
- os << "/";
- else
- os << "*";
- }
- }
- if (precedence <= upper_precedence)
- os << ")";
-}
-
-void ncmul::printcsrc(ostream & os, unsigned upper_precedence) const
-{
- debugmsg("ncmul print csrc",LOGLEVEL_PRINT);
- exvector::const_iterator it;
- exvector::const_iterator itend = seq.end()-1;
- os << "ncmul(";
- for (it=seq.begin(); it!=itend; ++it) {
- (*it).bp->printcsrc(os,precedence);
- os << ",";
- }
- (*it).bp->printcsrc(os,precedence);
- os << ")";
-}
-
-void relational::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) const
-{
- debugmsg("relational print csrc", LOGLEVEL_PRINT);
- if (precedence<=upper_precedence)
- os << "(";
-
- // Print left-hand expression
- lh.bp->printcsrc(os, type, precedence);
-
- // Print relational operator
- switch (o) {
- case equal:
- os << "==";
- break;
- case not_equal:
- os << "!=";
- break;
- case less:
- os << "<";
- break;
- case less_or_equal:
- os << "<=";
- break;
- case greater:
- os << ">";
- break;
- case greater_or_equal:
- os << ">=";
- break;
- default:
- os << "(INVALID RELATIONAL OPERATOR)";
- break;
- }
-
- // Print right-hand operator
- rh.bp->printcsrc(os, type, precedence);
-
- if (precedence <= upper_precedence)
- os << ")";
-}
-
-#ifndef NO_GINAC_NAMESPACE
-} // namespace GiNaC
-#endif // ndef NO_GINAC_NAMESPACE
+++ /dev/null
-/** @file printraw.cpp
- *
- * print in ugly raw format, so brave developers can have a look at the
- * underlying structure. */
-
-/*
- * GiNaC Copyright (C) 1999 Johannes Gutenberg University Mainz, Germany
- *
- * This program is free software; you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation; either version 2 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program; if not, write to the Free Software
- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
- */
-
-/* We are cheating here, because we don't want to include the underlying
- * bignum package's headers again, so in this file we omit the definition of
- * void numeric::printraw(ostream & os) const; */
-
-#include <iostream>
-
-#include "basic.h"
-#include "ex.h"
-#include "add.h"
-#include "constant.h"
-#include "expairseq.h"
-#include "fail.h"
-#include "indexed.h"
-#include "inifcns.h"
-#include "matrix.h"
-#include "mul.h"
-#include "ncmul.h"
-#include "numeric.h"
-#include "power.h"
-#include "relational.h"
-#include "series.h"
-#include "symbol.h"
-#include "debugmsg.h"
-
-#ifndef NO_GINAC_NAMESPACE
-namespace GiNaC {
-#endif // ndef NO_GINAC_NAMESPACE
-
-void ex::printraw(ostream & os) const
-{
- debugmsg("ex printraw",LOGLEVEL_PRINT);
- GINAC_ASSERT(bp!=0);
- os << "ex(";
- bp->printraw(os);
- os << ")";
-}
-
-void basic::printraw(ostream & os) const
-{
- debugmsg("basic printraw",LOGLEVEL_PRINT);
- os << "[basic object]";
-}
-
-void symbol::printraw(ostream & os) const
-{
- debugmsg("symbol printraw",LOGLEVEL_PRINT);
- os << "symbol(" << "name=" << name << ",serial=" << serial
- << ",hash=" << hashvalue << ",flags=" << flags << ")";
-}
-
-void constant::printraw(ostream & os) const
-{
- debugmsg("constant printraw",LOGLEVEL_PRINT);
- os << "constant(" << name << ")";
-}
-
-void power::printraw(ostream & os) const
-{
- debugmsg("power printraw",LOGLEVEL_PRINT);
-
- os << "power(";
- basis.printraw(os);
- os << ",";
- exponent.printraw(os);
- os << ",hash=" << hashvalue << ",flags=" << flags << ")";
-}
-
-void fail::printraw(ostream & os) const
-{
- debugmsg("fail printraw",LOGLEVEL_PRINT);
- os << "FAIL";
-}
-
-void expairseq::printraw(ostream & os) const
-{
- debugmsg("expairseq printraw",LOGLEVEL_PRINT);
-
- os << "expairseq(";
- for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
- os << "(";
- (*cit).rest.printraw(os);
- os << ",";
- (*cit).coeff.printraw(os);
- os << "),";
- }
- os << ")";
-}
-
-void add::printraw(ostream & os) const
-{
- debugmsg("add printraw",LOGLEVEL_PRINT);
-
- os << "+(";
- for (epvector::const_iterator it=seq.begin(); it!=seq.end(); ++it) {
- os << "(";
- (*it).rest.bp->printraw(os);
- os << ",";
- (*it).coeff.bp->printraw(os);
- os << "),";
- }
- os << ",hash=" << hashvalue << ",flags=" << flags;
- os << ")";
-}
-
-void mul::printraw(ostream & os) const
-{
- debugmsg("mul printraw",LOGLEVEL_PRINT);
-
- os << "*(";
- for (epvector::const_iterator it=seq.begin(); it!=seq.end(); ++it) {
- os << "(";
- (*it).rest.bp->printraw(os);
- os << ",";
- (*it).coeff.bp->printraw(os);
- os << "),";
- }
- os << ",hash=" << hashvalue << ",flags=" << flags;
- os << ")";
-}
-
-void ncmul::printraw(ostream & os) const
-{
- debugmsg("ncmul printraw",LOGLEVEL_PRINT);
-
- os << "%(";
- for (exvector::const_iterator it=seq.begin(); it!=seq.end(); ++it) {
- (*it).bp->printraw(os);
- os << ",";
- }
- os << ",hash=" << hashvalue << ",flags=" << flags;
- os << ")";
-}
-
-/*void function::printraw(ostream & os) const
- *{
- * debugmsg("function printraw",LOGLEVEL_PRINT);
- *
- * os << "function." << name << "(";
- * for (exvector::const_iterator it=seq.begin(); it!=seq.end(); ++it) {
- * (*it).bp->print(os);
- * os << ",";
- * }
- * os << ")";
- *}*/
-
-void series::printraw(ostream &os) const
-{
- debugmsg("symbol printraw", LOGLEVEL_PRINT);
- os << "series(" << var << ";" << point << ";";
- for (epvector::const_iterator i=seq.begin(); i!=seq.end(); i++) {
- os << "(" << (*i).rest << "," << (*i).coeff << "),";
- }
- os << ")";
-}
-
-void relational::printraw(ostream & os) const
-{
- debugmsg("relational printraw",LOGLEVEL_PRINT);
- os << "RELATIONAL(";
- lh.printraw(os);
- os << ",";
- rh.printraw(os);
- os << ",";
- switch (o) {
- case equal:
- os << "==";
- break;
- case not_equal:
- os << "!=";
- break;
- case less:
- os << "<";
- break;
- case less_or_equal:
- os << "<=";
- break;
- case greater:
- os << ">";
- break;
- case greater_or_equal:
- os << ">=";
- break;
- default:
- os << "(INVALID RELATIONAL OPERATOR)";
- }
- os << ")";
-}
-
-void matrix::printraw(ostream & os) const
-{
- debugmsg("matrix printraw",LOGLEVEL_PRINT);
- os << "matrix(" << row << "," << col <<",";
- for (int r=0; r<row-1; ++r) {
- os << "(";
- for (int c=0; c<col-1; ++c) {
- os << m[r*col+c] << ",";
- }
- os << m[col*(r-1)-1] << "),";
- }
- os << "(";
- for (int c=0; c<col-1; ++c) {
- os << m[(row-1)*col+c] << ",";
- }
- os << m[row*col-1] << "))";
-}
-
-#ifndef NO_GINAC_NAMESPACE
-} // namespace GiNaC
-#endif // ndef NO_GINAC_NAMESPACE
+++ /dev/null
-/** @file printtree.cpp
- *
- * print in tree- (indented-) form, so developers can have a look at the
- * underlying structure. */
-
-/*
- * GiNaC Copyright (C) 1999 Johannes Gutenberg University Mainz, Germany
- *
- * This program is free software; you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation; either version 2 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program; if not, write to the Free Software
- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
- */
-
-#include <iostream>
-#include <math.h>
-
-#include "basic.h"
-#include "ex.h"
-#include "add.h"
-#include "constant.h"
-#include "expairseq.h"
-#include "indexed.h"
-#include "inifcns.h"
-#include "mul.h"
-#include "ncmul.h"
-#include "numeric.h"
-#include "power.h"
-#include "relational.h"
-#include "series.h"
-#include "symbol.h"
-#include "debugmsg.h"
-
-#ifndef NO_GINAC_NAMESPACE
-namespace GiNaC {
-#endif // ndef NO_GINAC_NAMESPACE
-
-void ex::printtree(ostream & os, unsigned indent) const
-{
- debugmsg("ex printtree",LOGLEVEL_PRINT);
- GINAC_ASSERT(bp!=0);
- // os << "refcount=" << bp->refcount << " ";
- bp->printtree(os,indent);
-}
-
-void ex::dbgprinttree(void) const
-{
- debugmsg("ex dbgprinttree",LOGLEVEL_PRINT);
- GINAC_ASSERT(bp!=0);
- bp->dbgprinttree();
-}
-
-void basic::printtree(ostream & os, unsigned indent) const
-{
- debugmsg("basic printtree",LOGLEVEL_PRINT);
- os << string(indent,' ') << "type=" << typeid(*this).name()
- << ", hash=" << hashvalue << " (0x" << hex << hashvalue << dec << ")"
- << ", flags=" << flags
- << ", nops=" << nops() << endl;
- for (int i=0; i<nops(); ++i) {
- op(i).printtree(os,indent+delta_indent);
- }
-}
-
-void basic::dbgprinttree(void) const
-{
- printtree(cerr,0);
-}
-
-void numeric::printtree(ostream & os, unsigned indent) const
-{
- debugmsg("numeric printtree", LOGLEVEL_PRINT);
- // We are cheating here, because we don't want to include the underlying
- // bignum package's headers again, so we use ostream::operator<<(numeric):
- os << string(indent,' ');
- (*this).print(os);
- os << " (numeric): "
- << "hash=" << hashvalue << " (0x" << hex << hashvalue << dec << ")"
- << ", flags=" << flags << endl;
-}
-
-void symbol::printtree(ostream & os, unsigned indent) const
-{
- debugmsg("symbol printtree",LOGLEVEL_PRINT);
- os << string(indent,' ') << name << " (symbol): "
- << "serial=" << serial
- << ", hash=" << hashvalue << " (0x" << hex << hashvalue << dec << ")"
- << ", flags=" << flags << endl;
-}
-
-void power::printtree(ostream & os, unsigned indent) const
-{
- debugmsg("power printtree",LOGLEVEL_PRINT);
-
- os << string(indent,' ') << "power: "
- << "hash=" << hashvalue << " (0x" << hex << hashvalue << dec << ")"
- << ", flags=" << flags << endl;
- basis.printtree(os,indent+delta_indent);
- exponent.printtree(os,indent+delta_indent);
-}
-
-void expairseq::printtree(ostream & os, unsigned indent) const
-{
- debugmsg("expairseq printtree",LOGLEVEL_PRINT);
-
- os << string(indent,' ') << "type=" << typeid(*this).name()
- << ", hash=" << hashvalue << " (0x" << hex << hashvalue << dec << ")"
- << ", flags=" << flags
- << ", nops=" << nops() << endl;
- for (unsigned i=0; i<seq.size(); ++i) {
- seq[i].rest.printtree(os,indent+delta_indent);
- seq[i].coeff.printtree(os,indent+delta_indent);
- if (i!=seq.size()-1) {
- os << string(indent+delta_indent,' ') << "-----" << endl;
- }
- }
- if (!overall_coeff.is_equal(default_overall_coeff())) {
- os << string(indent+delta_indent,' ') << "-----" << endl;
- os << string(indent+delta_indent,' ') << "overall_coeff" << endl;
- overall_coeff.printtree(os,indent+delta_indent);
- }
- os << string(indent+delta_indent,' ') << "=====" << endl;
-#ifdef EXPAIRSEQ_USE_HASHTAB
- os << string(indent+delta_indent,' ')
- << "hashtab size " << hashtabsize << endl;
- if (hashtabsize==0) return;
-#define MAXCOUNT 5
- unsigned count[MAXCOUNT+1];
- for (int i=0; i<MAXCOUNT+1; ++i) count[i]=0;
- unsigned this_bin_fill;
- unsigned cum_fill_sq=0;
- unsigned cum_fill=0;
- for (unsigned i=0; i<hashtabsize; ++i) {
- this_bin_fill=0;
- if (hashtab[i].size()>0) {
- os << string(indent+delta_indent,' ')
- << "bin " << i << " with entries ";
- for (epplist::const_iterator it=hashtab[i].begin();
- it!=hashtab[i].end(); ++it) {
- os << *it-seq.begin() << " ";
- this_bin_fill++;
- }
- os << endl;
- cum_fill += this_bin_fill;
- cum_fill_sq += this_bin_fill*this_bin_fill;
- }
- if (this_bin_fill<MAXCOUNT) {
- ++count[this_bin_fill];
- } else {
- ++count[MAXCOUNT];
- }
- }
- unsigned fact=1;
- double cum_prob=0;
- double lambda=(1.0*seq.size())/hashtabsize;
- for (int k=0; k<MAXCOUNT; ++k) {
- if (k>0) fact *= k;
- double prob=pow(lambda,k)/fact*exp(-lambda);
- cum_prob += prob;
- os << string(indent+delta_indent,' ') << "bins with " << k << " entries: "
- << int(1000.0*count[k]/hashtabsize)/10.0 << "% (expected: "
- << int(prob*1000)/10.0 << ")" << endl;
- }
- os << string(indent+delta_indent,' ') << "bins with more entries: "
- << int(1000.0*count[MAXCOUNT]/hashtabsize)/10.0 << "% (expected: "
- << int((1-cum_prob)*1000)/10.0 << ")" << endl;
-
- os << string(indent+delta_indent,' ') << "variance: "
- << 1.0/hashtabsize*cum_fill_sq-(1.0/hashtabsize*cum_fill)*(1.0/hashtabsize*cum_fill)
- << endl;
- os << string(indent+delta_indent,' ') << "average fill: "
- << (1.0*cum_fill)/hashtabsize
- << " (should be equal to " << (1.0*seq.size())/hashtabsize << ")" << endl;
-#endif // def EXPAIRSEQ_USE_HASHTAB
-}
-
-#ifndef NO_GINAC_NAMESPACE
-} // namespace GiNaC
-#endif // ndef NO_GINAC_NAMESPACE
return new relational(*this);
}
+void relational::print(ostream & os, unsigned upper_precedence) const
+{
+ debugmsg("relational print",LOGLEVEL_PRINT);
+ if (precedence<=upper_precedence) os << "(";
+ lh.print(os,precedence);
+ switch (o) {
+ case equal:
+ os << "==";
+ break;
+ case not_equal:
+ os << "!=";
+ break;
+ case less:
+ os << "<";
+ break;
+ case less_or_equal:
+ os << "<=";
+ break;
+ case greater:
+ os << ">";
+ break;
+ case greater_or_equal:
+ os << ">=";
+ break;
+ default:
+ os << "(INVALID RELATIONAL OPERATOR)";
+ }
+ rh.print(os,precedence);
+ if (precedence<=upper_precedence) os << ")";
+}
+
+void relational::printraw(ostream & os) const
+{
+ debugmsg("relational printraw",LOGLEVEL_PRINT);
+ os << "RELATIONAL(";
+ lh.printraw(os);
+ os << ",";
+ rh.printraw(os);
+ os << ",";
+ switch (o) {
+ case equal:
+ os << "==";
+ break;
+ case not_equal:
+ os << "!=";
+ break;
+ case less:
+ os << "<";
+ break;
+ case less_or_equal:
+ os << "<=";
+ break;
+ case greater:
+ os << ">";
+ break;
+ case greater_or_equal:
+ os << ">=";
+ break;
+ default:
+ os << "(INVALID RELATIONAL OPERATOR)";
+ }
+ os << ")";
+}
+
+void relational::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) const
+{
+ debugmsg("relational print csrc", LOGLEVEL_PRINT);
+ if (precedence<=upper_precedence)
+ os << "(";
+
+ // Print left-hand expression
+ lh.bp->printcsrc(os, type, precedence);
+
+ // Print relational operator
+ switch (o) {
+ case equal:
+ os << "==";
+ break;
+ case not_equal:
+ os << "!=";
+ break;
+ case less:
+ os << "<";
+ break;
+ case less_or_equal:
+ os << "<=";
+ break;
+ case greater:
+ os << ">";
+ break;
+ case greater_or_equal:
+ os << ">=";
+ break;
+ default:
+ os << "(INVALID RELATIONAL OPERATOR)";
+ break;
+ }
+
+ // Print right-hand operator
+ rh.bp->printcsrc(os, type, precedence);
+
+ if (precedence <= upper_precedence)
+ os << ")";
+}
+
bool relational::info(unsigned inf) const
{
switch (inf) {
// functions overriding virtual functions from bases classes
public:
basic * duplicate() const;
- void printraw(ostream & os) const;
void print(ostream & os, unsigned upper_precedence=0) const;
+ void printraw(ostream & os) const;
void printcsrc(ostream & os, unsigned type, unsigned upper_precedence=0) const;
bool info(unsigned inf) const;
int nops() const;
return new series(*this);
}
+void series::print(ostream &os, unsigned upper_precedence) const
+{
+ debugmsg("symbol print", LOGLEVEL_PRINT);
+ convert_to_poly().print(os, upper_precedence);
+}
+
+void series::printraw(ostream &os) const
+{
+ debugmsg("symbol printraw", LOGLEVEL_PRINT);
+ os << "series(" << var << ";" << point << ";";
+ for (epvector::const_iterator i=seq.begin(); i!=seq.end(); i++) {
+ os << "(" << (*i).rest << "," << (*i).coeff << "),";
+ }
+ os << ")";
+}
+
// Highest degree of variable
int series::degree(symbol const &s) const
{
// functions overriding virtual functions from base classes
public:
basic *duplicate() const;
- void printraw(ostream &os) const;
void print(ostream &os, unsigned upper_precedence=0) const;
+ void printraw(ostream &os) const;
int degree(symbol const &s) const;
int ldegree(symbol const &s) const;
ex coeff(symbol const &s, int const n=1) const;
return new symbol(*this);
}
+void symbol::print(ostream & os, unsigned upper_precedence) const
+{
+ debugmsg("symbol print",LOGLEVEL_PRINT);
+ os << name;
+}
+
+void symbol::printraw(ostream & os) const
+{
+ debugmsg("symbol printraw",LOGLEVEL_PRINT);
+ os << "symbol(" << "name=" << name << ",serial=" << serial
+ << ",hash=" << hashvalue << ",flags=" << flags << ")";
+}
+
+void symbol::printtree(ostream & os, unsigned indent) const
+{
+ debugmsg("symbol printtree",LOGLEVEL_PRINT);
+ os << string(indent,' ') << name << " (symbol): "
+ << "serial=" << serial
+ << ", hash=" << hashvalue << " (0x" << hex << hashvalue << dec << ")"
+ << ", flags=" << flags << endl;
+}
+
+void symbol::printcsrc(ostream & os, unsigned type, unsigned upper_precedence) const
+{
+ debugmsg("symbol print csrc", LOGLEVEL_PRINT);
+ os << name;
+}
+
bool symbol::info(unsigned inf) const
{
if (inf==info_flags::symbol) return true;
- if (inf==info_flags::polynomial || inf==info_flags::integer_polynomial || inf==info_flags::rational_polynomial || inf==info_flags::rational_function) {
+ if (inf==info_flags::polynomial ||
+ inf==info_flags::integer_polynomial ||
+ inf==info_flags::cinteger_polynomial ||
+ inf==info_flags::rational_polynomial ||
+ inf==info_flags::crational_polynomial ||
+ inf==info_flags::rational_function) {
return true;
} else {
return basic::info(inf);
// functions overriding virtual functions from base classes
public:
basic * duplicate() const;
+ void print(ostream & os, unsigned upper_precedence=0) const;
void printraw(ostream & os) const;
void printtree(ostream & os, unsigned indent) const;
- void print(ostream & os, unsigned upper_precedence=0) const;
void printcsrc(ostream & os, unsigned type, unsigned upper_precedence=0) const;
bool info(unsigned inf) const;
ex expand(unsigned options=0) const;
@for file in $(DISTFILES); do \
d=$(srcdir); \
if test -d $$d/$$file; then \
- cp -pr $$d/$$file $(distdir)/$$file; \
+ cp -pr $$/$$file $(distdir)/$$file; \
else \
test -f $(distdir)/$$file \
|| ln $$d/$$file $(distdir)/$$file 2> /dev/null \