- enforced GiNaC coding standards :-)
[ginac.git] / check / series_expansion.cpp
1 /** @file series_expansion.cpp
2  *
3  *  Series expansion test (Laurent and Taylor series).
4  *
5  *  GiNaC Copyright (C) 1999 Johannes Gutenberg University Mainz, Germany
6  *
7  *  This program is free software; you can redistribute it and/or modify
8  *  it under the terms of the GNU General Public License as published by
9  *  the Free Software Foundation; either version 2 of the License, or
10  *  (at your option) any later version.
11  *
12  *  This program is distributed in the hope that it will be useful,
13  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
14  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15  *  GNU General Public License for more details.
16  *
17  *  You should have received a copy of the GNU General Public License
18  *  along with this program; if not, write to the Free Software
19  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
20  */
21
22 #include <ginac/ginac.h>
23
24 static symbol x("x");
25
26 static unsigned check_series(const ex &e, const ex &point, const ex &d)
27 {
28         ex es = e.series(x, point, 8);
29         ex ep = static_cast<series *>(es.bp)->convert_to_poly();
30         if ((ep - d).compare(exZERO()) != 0) {
31                 clog << "series expansion of " << e << " at " << point
32              << " erroneously returned " << ep << " (instead of " << d
33              << ")" << endl;
34                 (ep-d).printtree(clog);
35                 return 1;
36         }
37         return 0;
38 }
39
40 // Series expansion
41 static unsigned series1(void)
42 {
43         unsigned result = 0;
44         ex e, d;
45
46         e = sin(x);
47         d = x - pow(x, 3) / 6 + pow(x, 5) / 120 - pow(x, 7) / 5040 + Order(pow(x, 8));
48         result += check_series(e, exZERO(), d);
49
50         e = cos(x);
51         d = 1 - pow(x, 2) / 2 + pow(x, 4) / 24 - pow(x, 6) / 720 + Order(pow(x, 8));
52         result += check_series(e, exZERO(), d);
53
54         e = exp(x);
55         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));
56         result += check_series(e, exZERO(), d);
57
58         e = pow(1 - x, -1);
59         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));
60         result += check_series(e, exZERO(), d);
61
62         e = x + pow(x, -1);
63         d = x + pow(x, -1);
64         result += check_series(e, exZERO(), d);
65
66         e = x + pow(x, -1);
67         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));
68         result += check_series(e, exONE(), d);
69
70         e = pow(x + pow(x, 3), -1);
71         d = pow(x, -1) - x + pow(x, 3) - pow(x, 5) + Order(pow(x, 7));
72         result += check_series(e, exZERO(), d);
73
74         e = pow(pow(x, 2) + pow(x, 4), -1);
75         d = pow(x, -2) - 1 + pow(x, 2) - pow(x, 4) + Order(pow(x, 6));
76         result += check_series(e, exZERO(), d);
77
78         e = pow(sin(x), -2);
79         d = pow(x, -2) + numeric(1,3) + pow(x, 2) / 15 + pow(x, 4) * 2/189 + Order(pow(x, 5));
80         result += check_series(e, exZERO(), d);
81
82         e = sin(x) / cos(x);
83         d = x + pow(x, 3) / 3 + pow(x, 5) * 2/15 + pow(x, 7) * 17/315 + Order(pow(x, 8));
84         result += check_series(e, exZERO(), d);
85
86         e = cos(x) / sin(x);
87         d = pow(x, -1) - x / 3 - pow(x, 3) / 45 - pow(x, 5) * 2/945 + Order(pow(x, 6));
88         result += check_series(e, exZERO(), d);
89
90         e = pow(numeric(2), x);
91         ex t = log(ex(2)) * x;
92         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));
93         result += check_series(e, exZERO(), d.expand());
94
95         e = pow(Pi, x);
96         t = log(Pi) * x;
97         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));
98         result += check_series(e, exZERO(), d.expand());
99
100         return result;
101 }
102
103 // Series addition
104 static unsigned series2(void)
105 {
106         unsigned result = 0;
107         ex e, d;
108
109         e = pow(sin(x), -1).series(x, exZERO(), 8) + pow(sin(-x), -1).series(x, exZERO(), 12);
110         d = Order(pow(x, 6));
111         result += check_series(e, exZERO(), d);
112
113         return result;
114 }
115
116 // Series multiplication
117 static unsigned series3(void)
118 {
119         unsigned result = 0;
120         ex e, d;
121
122         e = sin(x).series(x, exZERO(), 8) * pow(sin(x), -1).series(x, exZERO(), 12);
123         d = 1 + Order(pow(x, 7));
124         result += check_series(e, exZERO(), d);
125
126         return result;
127 }
128
129 unsigned series_expansion(void)
130 {
131         unsigned result = 0;
132
133         cout << "checking series expansion..." << flush;
134         clog << "---------series expansion:" << endl;
135
136         result += series1();
137         result += series2();
138         result += series3();
139
140         if (!result) {
141                 cout << " passed ";
142                 clog << "(no output)" << endl;
143         } else {
144                 cout << " failed ";
145         }
146         return result;
147 }