Update copyright statements.
[ginac.git] / check / exam_paranoia.cpp
1 /** @file exam_paranoia.cpp
2  *
3  *  This set of tests checks for some of GiNaC's oopses which showed up during
4  *  development.  Things were evaluated wrongly and so.  Such a sick behaviour
5  *  shouldn't occur any more.  But we are paranoic and we want to exclude these
6  *  these oopses for good, so we run those stupid tests... */
7
8 /*
9  *  GiNaC Copyright (C) 1999-2014 Johannes Gutenberg University Mainz, Germany
10  *
11  *  This program is free software; you can redistribute it and/or modify
12  *  it under the terms of the GNU General Public License as published by
13  *  the Free Software Foundation; either version 2 of the License, or
14  *  (at your option) any later version.
15  *
16  *  This program is distributed in the hope that it will be useful,
17  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
18  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
19  *  GNU General Public License for more details.
20  *
21  *  You should have received a copy of the GNU General Public License
22  *  along with this program; if not, write to the Free Software
23  *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
24  */
25
26 #include "ginac.h"
27 using namespace GiNaC;
28
29 #include <iostream>
30 using namespace std;
31
32 // The very first pair of historic problems had its roots in power.cpp and was
33 // finally resolved on April 27th 1999. (Fixing the first on April 23rd
34 // actually introduced the second.)
35 static unsigned exam_paranoia1()
36 {
37         unsigned result = 0;
38         symbol x("x"), y("y"), z("z");
39         ex e, f, g;
40
41         e = x * y * z;
42         f = y * z;
43         g = e / f;
44
45         // In the first one expand did not do any job at all:
46         if (!g.expand().is_equal(x)) {
47                 clog << "e = x*y*z; f = y*z; expand(e/f) erroneously returned "
48                      << g.expand() << endl;
49                 ++result;
50         }
51
52         // This one somehow used to return 0:
53         e = pow(x + 1, -1);
54         if (!e.expand().is_equal(e)) {
55                 clog << "expand(pow(x + 1, -1)) erroneously returned "
56                      << e.expand() << endl;
57                 ++result;
58         }
59
60         return result;
61 }
62
63 // And here the second oops which showed up until May 17th 1999.  It had to do
64 // with lexicographic canonicalization and thus showed up only if the variables
65 // had the names as given here:
66 static unsigned exam_paranoia2()
67 {
68         unsigned result = 0;
69         symbol x("x"), y("y"), z("z");
70         ex e, f, g;
71
72         e = x + z*x;
73         f = e*y;
74         g = f - e*y;
75
76         // After .expand(), g should be zero:
77         if (!g.expand().is_zero()) {
78                 clog << "e = (x + z*x); f = e*y; expand(f - e*y) erroneously returned "
79                      << g.expand() << endl;
80                 ++result;
81         }
82         // After .eval(), g should be zero:
83         if (!g.eval().is_zero()) {
84                 clog << "e = (x + z*x); f = e*y; eval(f - e*y) erroneously returned "
85                      << g.eval() << endl;
86                 ++result;
87         }
88         // This actually worked already back in April 1999.
89         // But we are *very* paranoic!
90         if (!g.expand().eval().is_zero()) {
91                 clog << "e = (x + z*x); f = e*y; eval(expand(f - e*y)) erroneously returned "
92                      << g.expand().eval() << endl;
93                 ++result;
94         }
95
96         return result;
97 }
98
99 // The third bug was introduced on May 18th 1999, discovered on May 19 and
100 // fixed that same day.  It worked when x was substituted by 1 but not with
101 // other numbers:
102 static unsigned exam_paranoia3()
103 {
104         unsigned result = 0;
105         symbol x("x"), y("y");
106         ex e, f;
107
108         e = x*y - y;
109         f = e.subs(x == 2);
110
111         if (!f.is_equal(y)) {
112                 clog << "e = x*y - y; f = e.subs(x == 2) erroneously returned "
113                      << f << endl;
114                 ++result;
115         }
116         if (!f.eval().is_equal(y)) {
117                 clog << "e = x*y - y; eval(e.subs(x == 2)) erroneously returned "
118                      << f.eval() << endl;
119                 ++result;
120         }
121         if (!f.expand().is_equal(y)) {
122                 clog << "e = x*y - y; expand(e.subs(x == 2)) erroneously returned "
123                      << f.expand() << endl;
124                 ++result;
125         }
126
127         return result;
128 }
129
130 // The fourth bug was also discovered on May 19th 1999 and fixed immediately:
131 static unsigned exam_paranoia4()
132 {
133         unsigned result = 0;
134         symbol x("x");
135         ex e, f, g;
136
137         e = pow(x, 2) + x + 1;
138         f = pow(x, 2) + x + 1;
139         g = e - f;
140
141         if (!g.is_zero()) {
142                 clog << "e = pow(x,2) + x + 1; f = pow(x,2) + x + 1; g = e-f; g erroneously returned "
143                      << g << endl;
144                 ++result;
145         }
146         if (!g.is_zero()) {
147                 clog << "e = pow(x,2) + x + 1; f = pow(x,2) + x + 1; g = e-f; g.eval() erroneously returned "
148                      << g.eval() << endl;
149                 ++result;
150         }
151
152         return result;
153 }
154
155 // The fifth oops was discovered on May 20th 1999 and fixed a day later:
156 static unsigned exam_paranoia5()
157 {
158         unsigned result = 0;
159         symbol x("x"), y("y");
160
161         ex e, f;
162         e = pow(x*y + 1, 2);
163         f = pow(x, 2) * pow(y, 2) + 2*x*y + 1;
164
165         if (!(e-f).expand().is_zero()) {
166                 clog << "e = pow(x*y+1,2); f = pow(x,2)*pow(y,2) + 2*x*y + 1; (e-f).expand() erroneously returned "
167                      << (e-f).expand() << endl;
168                 ++result;
169         }
170
171         return result;
172 }
173
174 // This one was discovered on Jun 1st 1999 and fixed the same day:
175 static unsigned exam_paranoia6()
176 {
177         unsigned result = 0;
178         symbol x("x");
179
180         ex e, f;
181         e = pow(x, -5);
182         f = e.denom();
183
184         if (!f.is_equal(pow(x, 5))) {
185                 clog << "e = pow(x, -5); f = e.denom(); f was " << f << " (should be x^5)" << endl;
186                 ++result;
187         }
188         return result;
189 }
190
191 // This one was introduced on June 1st 1999 by some aggressive manual
192 // optimization. Discovered and fixed on June 2nd.
193 static unsigned exam_paranoia7()
194 {
195         unsigned result = 0;
196         symbol x("x"), y("y");
197
198         ex e = y + y*x + 2;
199         ex f = expand(pow(e, 2) - (e*y*(x + 1)));
200
201         if (f.nops() > 3) {
202                 clog << "e=y+y*x+2; f=expand(pow(e,2)-(e*y*(x+1))) has "
203                      << f.nops() << " arguments instead of 3 ( f=="
204                      << f << " )" << endl;
205                 ++result;
206         }
207         return result;
208 }
209
210 // This one was a result of the rewrite of mul::max_coefficient when we
211 // introduced the overall_coefficient field in expairseq objects on Oct 1st
212 // 1999. Fixed on Oct 4th.
213 static unsigned exam_paranoia8()
214 {
215         unsigned result = 0;
216         symbol x("x");
217
218         ex e = -x / (x+1);
219         ex f;
220         
221         try {
222                 f = e.normal();
223                 if (!f.is_equal(e)) {
224                         clog << "normal(-x/(x+1)) returns " << f << " instead of -x/(x+1)\n";
225                         ++result;
226                 }
227         } catch (const exception &err) {
228                 clog << "normal(-x/(x+1) throws " << err.what() << endl;
229                 ++result;
230         }
231         return result;
232 }
233
234 // This one was a result of a modification to frac_cancel() & Co. to avoid
235 // expanding the numerator and denominator when bringing them from Q[X] to
236 // Z[X]. multiply_lcm() forgot to multiply the x-linear term with the LCM of
237 // the coefficient's denominators (2 in this case).  Introduced on Jan 25th
238 // 2000 and fixed on Jan 31th.
239 static unsigned exam_paranoia9()
240 {
241         unsigned result = 0;
242         symbol x("x");
243
244         ex e = (exp(-x)-2*x*exp(-x)+pow(x,2)/2*exp(-x))/exp(-x);
245         ex f = e.normal();
246
247         if (!f.is_equal(1-2*x+pow(x,2)/2)) {
248                 clog << "normal(" << e << ") returns " << f << " instead of 1-2*x+1/2*x^2\n";
249                 ++result;
250         }
251         return result;
252 }
253
254 // I have no idea when this broke.  It has been working long ago, before 0.4.0
255 // and on Feb 13th 2000 I found out that things like 2^(3/2) throw an exception
256 // "power::eval(): pow(0,0) is undefined" instead of simplifying to 2*2^(1/2).
257 // It was fixed that same day.
258 static unsigned exam_paranoia10()
259 {
260         unsigned result = 0;
261         
262         ex b = numeric(2);
263         ex e = numeric(3,2);
264         ex r;
265         
266         try {
267                 r = pow(b,e).eval();
268                 if (!(r-2*sqrt(ex(2))).is_zero()) {
269                         clog << "2^(3/2) erroneously returned " << r << " instead of 2*sqrt(2)" << endl;
270                         ++result;
271                 }
272         } catch (const exception &err) {
273                 clog << "2^(3/2) throws " << err.what() << endl;
274                 ++result;
275         }
276         return result;
277 }
278
279 // After the rewriting of basic::normal() & Co. to return {num, den} lists,
280 // add::normal() forgot to multiply the denominator of the overall_coeff of
281 // its expanded and normalized children with the denominator of the expanded
282 // child (did you get this? Well, never mind...). Fixed on Feb 21th 2000.
283 static unsigned exam_paranoia11()
284 {
285         unsigned result = 0;
286         symbol x("x");
287
288         ex e = ((-5-2*x)-((2-5*x)/(-2+x))*(3+2*x))/(5-4*x);
289         ex f = e.normal();
290         ex d = normal((4+10*x+8*pow(x,2))/(x-2)/(5-4*x));
291
292         if (!(f - d).expand().is_zero()) {
293                 clog << "normal(" << e << ") returns " << f << " instead of " << d << endl;
294                 ++result;
295         }
296         return result;
297 }
298
299 // This one returned 0 because add::normal() incorrectly assumed that if the
300 // common denominator is 1, all the denominators would be 1 (they can in fact
301 // be +/-1). Fixed on Aug 2nd 2000.
302 static unsigned exam_paranoia12()
303 {
304         unsigned result = 0;
305         symbol x("x");
306         
307         ex e = 2-2*(1+x)/(-1-x);
308         ex f = e.normal();
309         ex d = 4;
310         
311         if (!(f - d).expand().is_zero()) {
312                 clog << "normal(" << e << ") returns " << f
313                      << " instead of " << d << endl;
314                 ++result;
315         }
316         return result;
317 }
318
319 // This one caused a division by 0 because heur_gcd() didn't check its
320 // input polynomials against 0. Fixed on Aug 4th 2000.
321 static unsigned exam_paranoia13()
322 {
323         unsigned result = 0;
324         symbol a("a"), b("b"), c("c");
325         
326         ex e = (b*a-c*a)/(4-a);
327         ex d = (c*a-b*a)/(a-4);
328         
329         try {
330                 ex f = e.normal();      
331                 if (!(f - d).expand().is_zero()) {
332                         clog << "normal(" << e << ") returns " << f
333                              << " instead of " << d << endl;
334                         ++result;
335                 }
336         } catch (const exception &err) {
337                 clog << "normal(" << e << ") throws " << err.what() << endl;
338                 ++result;
339         }
340         return result;
341 }
342
343 // A bug introduced on July 19, 2001. quo() and rem() would sometimes call
344 // vector::reserve() with a negative argument. Fixed on Dec 20, 2001.
345 static unsigned exam_paranoia14()
346 {
347         unsigned result = 0;
348         symbol x("x");
349
350         ex q = quo(1, pow(x, 3), x);
351         if (!q.is_zero()) {
352                 clog << "quo(1,x^3,x) erroneously returned " << q << " instead of 0\n";
353                 ++result;
354         }
355
356         return result;
357 }
358
359 // Under certain conditions, power::expand_add_2() could produce non-canonical
360 // numeric expairs. Fixed on Oct 24, 2002.
361 static unsigned exam_paranoia15()
362 {
363         unsigned result = 0;
364
365         ex q = (pow(pow(2, numeric(1, 2))*2+1, 2)).expand();
366         // this used to produce "1+4*sqrt(2)+4*2" which would never evaluate
367         // to "9+4*sqrt(2)"
368
369         if (!(q-9-4*pow(2, numeric(1, 2))).is_zero()) {
370                 clog << "expand((sqrt(2)*2+1)^2) erroneously returned " << q << " instead of 9-4*sqrt(2)\n";
371                 ++result;
372         }
373
374         return result;
375 }
376
377 // Expanding products containing powers of sums could return results that
378 // were not fully expanded. Fixed on Dec 10, 2003.
379 static unsigned exam_paranoia16()
380 {
381         unsigned result = 0;
382         symbol a("a"), b("b"), c("c"), d("d"), e("e");
383         ex e1, e2, e3;
384
385         e1 = pow(1+a*sqrt(b+c), 2);
386         e2 = e1.expand();
387
388         if (e2.has(pow(a, 2)*(b+c))) {
389                 clog << "expand(" << e1 << ") didn't fully expand\n";
390                 ++result;
391         }
392
393         e1 = (d*sqrt(a+b)+a*sqrt(c+d))*(b*sqrt(a+b)+a*sqrt(c+d));
394         e2 = e1.expand();
395
396         if (e2.has(pow(a, 2)*(c+d))) {
397                 clog << "expand(" << e1 << ") didn't fully expand\n";
398                 ++result;
399         }
400
401         e1 = (a+sqrt(b+c))*sqrt(b+c)*(d+sqrt(b+c));
402         e2 = e1.expand();
403
404         if (e2.has(a*(b+c))) {
405                 clog << "expand(" << e1 << ") didn't fully expand\n";
406                 ++result;
407         }
408
409         e1 = pow(sqrt(a+b)+sqrt(c+d), 3);
410         e2 = e1.expand();
411
412         if (e2.has(3*(a+b)*sqrt(c+d)) || e2.has(3*(c+d)*sqrt(a+b))) {
413                 clog << "expand(" << e1 << ") didn't fully expand\n";
414                 ++result;
415         }
416
417         e1 = a*(b+c*(d+e));
418         e2 = e1.expand();
419
420         if (e2.has(c*(d+e))) {
421                 clog << "expand(" << e1 << ") didn't fully expand\n";
422                 ++result;
423         }
424
425         e1 = 2*pow(1+a, 2)/a;
426         e2 = e1.expand();
427
428         if (e2.has(pow(a, 2))) {
429                 clog << "expand(" << e1 << ") didn't fully expand\n";
430                 ++result;
431         }
432
433         e1 = a*(a+b);
434         e2 = pow(pow(e1, -1), -1);
435
436         if (e2.has(a*b)) {
437                 clog << "double reciprocal expanded where it should not\n";
438                 ++result;
439         }
440
441         return result;
442 }
443
444 // Bug in reposition_dummy_indices() could result in correct expression
445 // turned into one with inconsistent indices. Fixed on Aug 29, 2006
446 static unsigned exam_paranoia17()
447 {
448         varidx mu1(symbol("mu1"), 4);
449         varidx mu2(symbol("mu2"), 4);
450         varidx mu3(symbol("mu3"), 4);
451         varidx mu4(symbol("mu4"), 4);
452         varidx mu5(symbol("mu5"), 4);
453         varidx mu6(symbol("mu6"), 4);
454
455         exvector ev2;
456         ev2.push_back(mu3.toggle_variance());
457         ev2.push_back(mu6);
458         ev2.push_back(mu5.toggle_variance());
459         ev2.push_back(mu6.toggle_variance());
460         ev2.push_back(mu5);
461         ev2.push_back(mu3); 
462         // notice: all indices are contracted ...
463
464         ex test_cycl = indexed(symbol("A"), sy_cycl(), ev2);
465         test_cycl = test_cycl.simplify_indexed();
466         // ... so there should be zero free indices in the end.
467         return test_cycl.get_free_indices().size();
468 }
469
470 // Bug in add::eval() could result in numeric terms not being collected into
471 // the overall coefficient. Fixed on Sep 22, 2010
472 static unsigned exam_paranoia18()
473 {
474         ex sqrt2 = sqrt(ex(2));
475         ex e = 1+2*(sqrt2+1)*(sqrt2-1);
476         if ( e.real_part() != 3 ) {
477                 clog << "real_part(1+2*(sqrt(2)+1)*(sqrt(2)-1)) failed to evaluate to 3\n";
478                 return 1;
479         }
480         return 0;
481 }
482
483 // Bug in mul::conjugate when factors are evaluated at branch cuts, reported as
484 // Sage bug #10964.
485 static unsigned exam_paranoia19()
486 {
487         symbol a("a");
488         ex e = conjugate(a*sqrt(ex(-2))*sqrt(ex(-3)));
489         ex c = a*conjugate(sqrt(ex(-2)))*conjugate(sqrt(ex(-3)));
490         if (!subs(e-c, a==42).is_zero()) {
491                 clog << "subs(a*conjugate(sqrt(-2))*conjugate(sqrt(-3))-conjugate(a*sqrt(-2)*sqrt(-3)),a==42) failed to evaluate to 0\n";
492                 return 1;
493         }
494         return 0;
495 }
496
497 // Bugs in is_polynomial (fixed 2011-05-20 and 2014-07-26).
498 static unsigned exam_paranoia20()
499 {
500         unsigned result = 0;
501         symbol x("x"), y("y");
502         ex e1 = sqrt(x*x+1)*sqrt(x+1);
503         if (e1.is_polynomial(x)) {
504                 clog << "sqrt(x*x+1)*sqrt(x+1) is wrongly reported to be a polynomial in x\n";
505                 ++result;
506         }
507         ex e2 = sqrt(Pi)*x;
508         if (!e2.is_polynomial(x)) {
509                 clog << "sqrt(Pi)*x is wrongly reported to be no polynomial in x\n";
510                 ++result;
511         }
512         ex e3 = sqrt(x);
513         if (!e3.is_polynomial(y)) {
514                 clog << "sqrt(x) is wrongly reported to be no polynomial in y\n";
515                 ++result;
516         }
517         ex e4 = (1+y)/(2+x);
518         if (e4.is_polynomial(x)) {
519                 clog << "(1+y)/(2+x) is wrongly reported to be a polynomial in x\n";
520                 ++result;
521         }
522         return result;
523 }
524
525 static unsigned is_polynomial_false_positive()
526 {
527         unsigned result = 0;
528         symbol x("x"), n("n");
529         exvector nonpoly_exprs;
530         nonpoly_exprs.push_back(1/(1-x));
531         nonpoly_exprs.push_back(1/(x+1));
532         nonpoly_exprs.push_back(-1/(x-1));
533         nonpoly_exprs.push_back(1/(1-x*x));
534         nonpoly_exprs.push_back(1/(1-pow(x,n)));
535         nonpoly_exprs.push_back(x-1/(x-1));
536         for (exvector::const_iterator ep = nonpoly_exprs.begin();
537              ep != nonpoly_exprs.end(); ++ep) {
538                 if (ep->is_polynomial(x)) {
539                         clog << "(" << *ep << ").is_polynomial(" << x << ") "
540                                 "erroneously returned true" << endl;
541                         ++result;
542                 }
543         }
544         return result;
545 }
546
547 unsigned exam_paranoia()
548 {
549         unsigned result = 0;
550         
551         cout << "examining several historic failures just out of paranoia" << flush;
552         
553         result += exam_paranoia1();  cout << '.' << flush;
554         result += exam_paranoia2();  cout << '.' << flush;
555         result += exam_paranoia3();  cout << '.' << flush;
556         result += exam_paranoia4();  cout << '.' << flush;
557         result += exam_paranoia5();  cout << '.' << flush;
558         result += exam_paranoia6();  cout << '.' << flush;
559         result += exam_paranoia7();  cout << '.' << flush;
560         result += exam_paranoia8();  cout << '.' << flush;
561         result += exam_paranoia9();  cout << '.' << flush;
562         result += exam_paranoia10();  cout << '.' << flush;
563         result += exam_paranoia11();  cout << '.' << flush;
564         result += exam_paranoia12();  cout << '.' << flush;
565         result += exam_paranoia13();  cout << '.' << flush;
566         result += exam_paranoia14();  cout << '.' << flush;
567         result += exam_paranoia15();  cout << '.' << flush;
568         result += exam_paranoia16();  cout << '.' << flush;
569         result += exam_paranoia17();  cout << '.' << flush;
570         result += exam_paranoia18();  cout << '.' << flush;
571         result += exam_paranoia19();  cout << '.' << flush;
572         result += exam_paranoia20();  cout << '.' << flush;
573         result += is_polynomial_false_positive(); cout << '.' << flush;
574         
575         return result;
576 }
577
578 int main(int argc, char** argv)
579 {
580         return exam_paranoia();
581 }