Another restructuring: moved include/*.h -> include/GiNaC/*.h in order to
[ginac.git] / check / series_expansion.cpp
1 // check/series_expansion.cpp
2
3 /* Series expansion test (Laurent and Taylor series). */
4
5 #include <GiNaC/ginac.h>
6
7 static symbol x("x");
8
9 static unsigned check_series(const ex &e, const ex &point, const ex &d)
10 {
11         ex es = e.series(x, point, 8);
12         ex ep = static_cast<series *>(es.bp)->convert_to_poly();
13         if ((ep - d).compare(exZERO()) != 0) {
14                 clog << "series expansion of " << e << " at " << point
15              << " erroneously returned " << ep << " (instead of " << d
16              << ")" << endl;
17                 (ep-d).printtree(clog);
18                 return 1;
19         }
20         return 0;
21 }
22
23 // Series expansion
24 static unsigned series1(void)
25 {
26         unsigned result = 0;
27         ex e, d;
28
29         e = sin(x);
30         d = x - pow(x, 3) / 6 + pow(x, 5) / 120 - pow(x, 7) / 5040 + Order(pow(x, 8));
31         result += check_series(e, exZERO(), d);
32
33         e = cos(x);
34         d = 1 - pow(x, 2) / 2 + pow(x, 4) / 24 - pow(x, 6) / 720 + Order(pow(x, 8));
35         result += check_series(e, exZERO(), d);
36
37         e = exp(x);
38         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));
39         result += check_series(e, exZERO(), d);
40
41         e = pow(1 - x, -1);
42         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));
43         result += check_series(e, exZERO(), d);
44
45         e = x + pow(x, -1);
46         d = x + pow(x, -1);
47         result += check_series(e, exZERO(), d);
48
49         e = x + pow(x, -1);
50         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));
51         result += check_series(e, exONE(), d);
52
53         e = pow(x + pow(x, 3), -1);
54         d = pow(x, -1) - x + pow(x, 3) - pow(x, 5) + Order(pow(x, 7));
55         result += check_series(e, exZERO(), d);
56
57         e = pow(pow(x, 2) + pow(x, 4), -1);
58         d = pow(x, -2) - 1 + pow(x, 2) - pow(x, 4) + Order(pow(x, 6));
59         result += check_series(e, exZERO(), d);
60
61         e = pow(sin(x), -2);
62         d = pow(x, -2) + numeric(1,3) + pow(x, 2) / 15 + pow(x, 4) * 2/189 + Order(pow(x, 5));
63         result += check_series(e, exZERO(), d);
64
65         e = sin(x) / cos(x);
66         d = x + pow(x, 3) / 3 + pow(x, 5) * 2/15 + pow(x, 7) * 17/315 + Order(pow(x, 8));
67         result += check_series(e, exZERO(), d);
68
69         e = cos(x) / sin(x);
70         d = pow(x, -1) - x / 3 - pow(x, 3) / 45 - pow(x, 5) * 2/945 + Order(pow(x, 6));
71         result += check_series(e, exZERO(), d);
72
73         e = pow(numeric(2), x);
74         ex t = log(ex(2)) * x;
75         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));
76         result += check_series(e, exZERO(), d.expand());
77
78         e = pow(Pi, x);
79         t = log(Pi) * x;
80         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));
81         result += check_series(e, exZERO(), d.expand());
82
83         return result;
84 }
85
86 // Series addition
87 static unsigned series2(void)
88 {
89         unsigned result = 0;
90         ex e, d;
91
92         e = pow(sin(x), -1).series(x, exZERO(), 8) + pow(sin(-x), -1).series(x, exZERO(), 12);
93         d = Order(pow(x, 6));
94         result += check_series(e, exZERO(), d);
95
96         return result;
97 }
98
99 // Series multiplication
100 static unsigned series3(void)
101 {
102         unsigned result = 0;
103         ex e, d;
104
105         e = sin(x).series(x, exZERO(), 8) * pow(sin(x), -1).series(x, exZERO(), 12);
106         d = 1 + Order(pow(x, 7));
107         result += check_series(e, exZERO(), d);
108
109         return result;
110 }
111
112 unsigned series_expansion(void)
113 {
114         unsigned result = 0;
115
116         cout << "checking series expansion..." << flush;
117         clog << "---------series expansion:" << endl;
118
119         result += series1();
120         result += series2();
121         result += series3();
122
123         if (!result) {
124                 cout << " passed ";
125                 clog << "(no output)" << endl;
126         } else {
127                 cout << " failed ";
128         }
129         return result;
130 }