Make .eval() evaluate top-level only.
[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-2015 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 .eval(), g should be zero:
77         if (!g.is_zero()) {
78                 clog << "e = (x + z*x); f = e*y; g = (f - e*y) erroneously returned g == "
79                      << g << endl;
80                 ++result;
81         }
82
83         return result;
84 }
85
86 // The third bug was introduced on May 18th 1999, discovered on May 19 and
87 // fixed that same day.  It worked when x was substituted by 1 but not with
88 // other numbers:
89 static unsigned exam_paranoia3()
90 {
91         unsigned result = 0;
92         symbol x("x"), y("y");
93         ex e, f;
94
95         e = x*y - y;
96         f = e.subs(x == 2);
97
98         if (!f.is_equal(y)) {
99                 clog << "e = x*y - y; f = e.subs(x == 2) erroneously returned "
100                      << f << endl;
101                 ++result;
102         }
103
104         return result;
105 }
106
107 // The fourth bug was also discovered on May 19th 1999 and fixed immediately:
108 static unsigned exam_paranoia4()
109 {
110         unsigned result = 0;
111         symbol x("x");
112         ex e, f, g;
113
114         e = pow(x, 2) + x + 1;
115         f = pow(x, 2) + x + 1;
116         g = e - f;
117
118         if (!g.is_zero()) {
119                 clog << "e = pow(x,2) + x + 1; f = pow(x,2) + x + 1; g = e-f; g erroneously returned "
120                      << g << endl;
121                 ++result;
122         }
123
124         return result;
125 }
126
127 // The fifth oops was discovered on May 20th 1999 and fixed a day later:
128 static unsigned exam_paranoia5()
129 {
130         unsigned result = 0;
131         symbol x("x"), y("y");
132
133         ex e, f;
134         e = pow(x*y + 1, 2);
135         f = pow(x, 2) * pow(y, 2) + 2*x*y + 1;
136
137         if (!(e-f).expand().is_zero()) {
138                 clog << "e = pow(x*y+1,2); f = pow(x,2)*pow(y,2) + 2*x*y + 1; (e-f).expand() erroneously returned "
139                      << (e-f).expand() << endl;
140                 ++result;
141         }
142
143         return result;
144 }
145
146 // This one was discovered on Jun 1st 1999 and fixed the same day:
147 static unsigned exam_paranoia6()
148 {
149         unsigned result = 0;
150         symbol x("x");
151
152         ex e, f;
153         e = pow(x, -5);
154         f = e.denom();
155
156         if (!f.is_equal(pow(x, 5))) {
157                 clog << "e = pow(x, -5); f = e.denom(); f was " << f << " (should be x^5)" << endl;
158                 ++result;
159         }
160         return result;
161 }
162
163 // This one was introduced on June 1st 1999 by some aggressive manual
164 // optimization. Discovered and fixed on June 2nd.
165 static unsigned exam_paranoia7()
166 {
167         unsigned result = 0;
168         symbol x("x"), y("y");
169
170         ex e = y + y*x + 2;
171         ex f = expand(pow(e, 2) - (e*y*(x + 1)));
172
173         if (f.nops() > 3) {
174                 clog << "e=y+y*x+2; f=expand(pow(e,2)-(e*y*(x+1))) has "
175                      << f.nops() << " arguments instead of 3 ( f=="
176                      << f << " )" << endl;
177                 ++result;
178         }
179         return result;
180 }
181
182 // This one was a result of the rewrite of mul::max_coefficient when we
183 // introduced the overall_coefficient field in expairseq objects on Oct 1st
184 // 1999. Fixed on Oct 4th.
185 static unsigned exam_paranoia8()
186 {
187         unsigned result = 0;
188         symbol x("x");
189
190         ex e = -x / (x+1);
191         ex f;
192         
193         try {
194                 f = e.normal();
195                 if (!f.is_equal(e)) {
196                         clog << "normal(-x/(x+1)) returns " << f << " instead of -x/(x+1)\n";
197                         ++result;
198                 }
199         } catch (const exception &err) {
200                 clog << "normal(-x/(x+1) throws " << err.what() << endl;
201                 ++result;
202         }
203         return result;
204 }
205
206 // This one was a result of a modification to frac_cancel() & Co. to avoid
207 // expanding the numerator and denominator when bringing them from Q[X] to
208 // Z[X]. multiply_lcm() forgot to multiply the x-linear term with the LCM of
209 // the coefficient's denominators (2 in this case).  Introduced on Jan 25th
210 // 2000 and fixed on Jan 31th.
211 static unsigned exam_paranoia9()
212 {
213         unsigned result = 0;
214         symbol x("x");
215
216         ex e = (exp(-x)-2*x*exp(-x)+pow(x,2)/2*exp(-x))/exp(-x);
217         ex f = e.normal();
218
219         if (!f.is_equal(1-2*x+pow(x,2)/2)) {
220                 clog << "normal(" << e << ") returns " << f << " instead of 1-2*x+1/2*x^2\n";
221                 ++result;
222         }
223         return result;
224 }
225
226 // I have no idea when this broke.  It has been working long ago, before 0.4.0
227 // and on Feb 13th 2000 I found out that things like 2^(3/2) throw an exception
228 // "power::eval(): pow(0,0) is undefined" instead of simplifying to 2*2^(1/2).
229 // It was fixed that same day.
230 static unsigned exam_paranoia10()
231 {
232         unsigned result = 0;
233         
234         ex b = numeric(2);
235         ex e = numeric(3,2);
236         ex r;
237         
238         try {
239                 r = pow(b, e);
240                 if (!(r-2*sqrt(ex(2))).is_zero()) {
241                         clog << "2^(3/2) erroneously returned " << r << " instead of 2*sqrt(2)" << endl;
242                         ++result;
243                 }
244         } catch (const exception &err) {
245                 clog << "2^(3/2) throws " << err.what() << endl;
246                 ++result;
247         }
248         return result;
249 }
250
251 // After the rewriting of basic::normal() & Co. to return {num, den} lists,
252 // add::normal() forgot to multiply the denominator of the overall_coeff of
253 // its expanded and normalized children with the denominator of the expanded
254 // child (did you get this? Well, never mind...). Fixed on Feb 21th 2000.
255 static unsigned exam_paranoia11()
256 {
257         unsigned result = 0;
258         symbol x("x");
259
260         ex e = ((-5-2*x)-((2-5*x)/(-2+x))*(3+2*x))/(5-4*x);
261         ex f = e.normal();
262         ex d = normal((4+10*x+8*pow(x,2))/(x-2)/(5-4*x));
263
264         if (!(f - d).expand().is_zero()) {
265                 clog << "normal(" << e << ") returns " << f << " instead of " << d << endl;
266                 ++result;
267         }
268         return result;
269 }
270
271 // This one returned 0 because add::normal() incorrectly assumed that if the
272 // common denominator is 1, all the denominators would be 1 (they can in fact
273 // be +/-1). Fixed on Aug 2nd 2000.
274 static unsigned exam_paranoia12()
275 {
276         unsigned result = 0;
277         symbol x("x");
278         
279         ex e = 2-2*(1+x)/(-1-x);
280         ex f = e.normal();
281         ex d = 4;
282         
283         if (!(f - d).expand().is_zero()) {
284                 clog << "normal(" << e << ") returns " << f
285                      << " instead of " << d << endl;
286                 ++result;
287         }
288         return result;
289 }
290
291 // This one caused a division by 0 because heur_gcd() didn't check its
292 // input polynomials against 0. Fixed on Aug 4th 2000.
293 static unsigned exam_paranoia13()
294 {
295         unsigned result = 0;
296         symbol a("a"), b("b"), c("c");
297         
298         ex e = (b*a-c*a)/(4-a);
299         ex d = (c*a-b*a)/(a-4);
300         
301         try {
302                 ex f = e.normal();      
303                 if (!(f - d).expand().is_zero()) {
304                         clog << "normal(" << e << ") returns " << f
305                              << " instead of " << d << endl;
306                         ++result;
307                 }
308         } catch (const exception &err) {
309                 clog << "normal(" << e << ") throws " << err.what() << endl;
310                 ++result;
311         }
312         return result;
313 }
314
315 // A bug introduced on July 19, 2001. quo() and rem() would sometimes call
316 // vector::reserve() with a negative argument. Fixed on Dec 20, 2001.
317 static unsigned exam_paranoia14()
318 {
319         unsigned result = 0;
320         symbol x("x");
321
322         ex q = quo(1, pow(x, 3), x);
323         if (!q.is_zero()) {
324                 clog << "quo(1,x^3,x) erroneously returned " << q << " instead of 0\n";
325                 ++result;
326         }
327
328         return result;
329 }
330
331 // Under certain conditions, power::expand_add_2() could produce non-canonical
332 // numeric expairs. Fixed on Oct 24, 2002.
333 static unsigned exam_paranoia15()
334 {
335         unsigned result = 0;
336
337         ex q = (pow(pow(2, numeric(1, 2))*2+1, 2)).expand();
338         // this used to produce "1+4*sqrt(2)+4*2" which would never evaluate
339         // to "9+4*sqrt(2)"
340
341         if (!(q-9-4*pow(2, numeric(1, 2))).is_zero()) {
342                 clog << "expand((sqrt(2)*2+1)^2) erroneously returned " << q << " instead of 9-4*sqrt(2)\n";
343                 ++result;
344         }
345
346         return result;
347 }
348
349 // Expanding products containing powers of sums could return results that
350 // were not fully expanded. Fixed on Dec 10, 2003.
351 static unsigned exam_paranoia16()
352 {
353         unsigned result = 0;
354         symbol a("a"), b("b"), c("c"), d("d"), e("e");
355         ex e1, e2, e3;
356
357         e1 = pow(1+a*sqrt(b+c), 2);
358         e2 = e1.expand();
359
360         if (e2.has(pow(a, 2)*(b+c))) {
361                 clog << "expand(" << e1 << ") didn't fully expand\n";
362                 ++result;
363         }
364
365         e1 = (d*sqrt(a+b)+a*sqrt(c+d))*(b*sqrt(a+b)+a*sqrt(c+d));
366         e2 = e1.expand();
367
368         if (e2.has(pow(a, 2)*(c+d))) {
369                 clog << "expand(" << e1 << ") didn't fully expand\n";
370                 ++result;
371         }
372
373         e1 = (a+sqrt(b+c))*sqrt(b+c)*(d+sqrt(b+c));
374         e2 = e1.expand();
375
376         if (e2.has(a*(b+c))) {
377                 clog << "expand(" << e1 << ") didn't fully expand\n";
378                 ++result;
379         }
380
381         e1 = pow(sqrt(a+b)+sqrt(c+d), 3);
382         e2 = e1.expand();
383
384         if (e2.has(3*(a+b)*sqrt(c+d)) || e2.has(3*(c+d)*sqrt(a+b))) {
385                 clog << "expand(" << e1 << ") didn't fully expand\n";
386                 ++result;
387         }
388
389         e1 = a*(b+c*(d+e));
390         e2 = e1.expand();
391
392         if (e2.has(c*(d+e))) {
393                 clog << "expand(" << e1 << ") didn't fully expand\n";
394                 ++result;
395         }
396
397         e1 = 2*pow(1+a, 2)/a;
398         e2 = e1.expand();
399
400         if (e2.has(pow(a, 2))) {
401                 clog << "expand(" << e1 << ") didn't fully expand\n";
402                 ++result;
403         }
404
405         e1 = a*(a+b);
406         e2 = pow(pow(e1, -1), -1);
407
408         if (e2.has(a*b)) {
409                 clog << "double reciprocal expanded where it should not\n";
410                 ++result;
411         }
412
413         return result;
414 }
415
416 // Bug in reposition_dummy_indices() could result in correct expression
417 // turned into one with inconsistent indices. Fixed on Aug 29, 2006
418 static unsigned exam_paranoia17()
419 {
420         varidx mu1(symbol("mu1"), 4);
421         varidx mu2(symbol("mu2"), 4);
422         varidx mu3(symbol("mu3"), 4);
423         varidx mu4(symbol("mu4"), 4);
424         varidx mu5(symbol("mu5"), 4);
425         varidx mu6(symbol("mu6"), 4);
426
427         exvector ev2;
428         ev2.push_back(mu3.toggle_variance());
429         ev2.push_back(mu6);
430         ev2.push_back(mu5.toggle_variance());
431         ev2.push_back(mu6.toggle_variance());
432         ev2.push_back(mu5);
433         ev2.push_back(mu3); 
434         // notice: all indices are contracted ...
435
436         ex test_cycl = indexed(symbol("A"), sy_cycl(), ev2);
437         test_cycl = test_cycl.simplify_indexed();
438         // ... so there should be zero free indices in the end.
439         return test_cycl.get_free_indices().size();
440 }
441
442 // Bug in add::eval() could result in numeric terms not being collected into
443 // the overall coefficient. Fixed on Sep 22, 2010
444 static unsigned exam_paranoia18()
445 {
446         ex sqrt2 = sqrt(ex(2));
447         ex e = 1+2*(sqrt2+1)*(sqrt2-1);
448         if ( e.real_part() != 3 ) {
449                 clog << "real_part(1+2*(sqrt(2)+1)*(sqrt(2)-1)) failed to evaluate to 3\n";
450                 return 1;
451         }
452         return 0;
453 }
454
455 // Bug in mul::conjugate when factors are evaluated at branch cuts, reported as
456 // Sage bug #10964.
457 static unsigned exam_paranoia19()
458 {
459         symbol a("a");
460         ex e = conjugate(a*sqrt(ex(-2))*sqrt(ex(-3)));
461         ex c = a*conjugate(sqrt(ex(-2)))*conjugate(sqrt(ex(-3)));
462         if (!subs(e-c, a==42).is_zero()) {
463                 clog << "subs(a*conjugate(sqrt(-2))*conjugate(sqrt(-3))-conjugate(a*sqrt(-2)*sqrt(-3)),a==42) failed to evaluate to 0\n";
464                 return 1;
465         }
466         return 0;
467 }
468
469 // Bugs in is_polynomial (fixed 2011-05-20 and 2014-07-26).
470 static unsigned exam_paranoia20()
471 {
472         unsigned result = 0;
473         symbol x("x"), y("y");
474         ex e1 = sqrt(x*x+1)*sqrt(x+1);
475         if (e1.is_polynomial(x)) {
476                 clog << "sqrt(x*x+1)*sqrt(x+1) is wrongly reported to be a polynomial in x\n";
477                 ++result;
478         }
479         ex e2 = sqrt(Pi)*x;
480         if (!e2.is_polynomial(x)) {
481                 clog << "sqrt(Pi)*x is wrongly reported to be no polynomial in x\n";
482                 ++result;
483         }
484         ex e3 = sqrt(x);
485         if (!e3.is_polynomial(y)) {
486                 clog << "sqrt(x) is wrongly reported to be no polynomial in y\n";
487                 ++result;
488         }
489         ex e4 = (1+y)/(2+x);
490         if (e4.is_polynomial(x)) {
491                 clog << "(1+y)/(2+x) is wrongly reported to be a polynomial in x\n";
492                 ++result;
493         }
494         return result;
495 }
496
497 static unsigned is_polynomial_false_positive()
498 {
499         unsigned result = 0;
500         symbol x("x"), n("n");
501         exvector nonpoly_exprs;
502         nonpoly_exprs.push_back(1/(1-x));
503         nonpoly_exprs.push_back(1/(x+1));
504         nonpoly_exprs.push_back(-1/(x-1));
505         nonpoly_exprs.push_back(1/(1-x*x));
506         nonpoly_exprs.push_back(1/(1-pow(x,n)));
507         nonpoly_exprs.push_back(x-1/(x-1));
508         for (exvector::const_iterator ep = nonpoly_exprs.begin();
509              ep != nonpoly_exprs.end(); ++ep) {
510                 if (ep->is_polynomial(x)) {
511                         clog << "(" << *ep << ").is_polynomial(" << x << ") "
512                                 "erroneously returned true" << endl;
513                         ++result;
514                 }
515         }
516         return result;
517 }
518
519 // Bug in power::expand reported by Isuru Fernando (fixed 2015-05-07).
520 static unsigned exam_paranoia21()
521 {
522         symbol x("x");
523         ex e = pow(x + sqrt(ex(2))*x, 2).expand();
524         if (e.nops() != 2) {
525                 clog << "(x+sqrt(2)*x)^2 was wrongly expanded to " << e << "\n";
526                 return 1;
527         }
528         return 0;
529 }
530
531 // Bug in power::expand (fixed 2015-07-18).
532 static unsigned exam_paranoia22()
533 {
534         symbol x("x"), y("y");
535         ex e = pow(sqrt(1+x)+y*sqrt(1+x), 2).expand();
536         if (e.nops() != 6) {
537                 clog << "(sqrt(1+x)+y*sqrt(1+x))^2 was wrongly expanded to " << e << "\n";
538                 return 1;
539         }
540         return 0;
541 }
542
543 // Bug in expairseq::evalchildren().
544 static unsigned exam_paranoia23()
545 {
546         unsigned result = 0;
547         symbol x("x");
548
549         epvector v1;
550         v1.push_back(expair(1, 1));
551         v1.push_back(expair(2*x, -1));
552         ex e1 = add(v1);  // Should be e==1-2*x,
553         if (!e1.is_equal(1-2*x)) {
554                 clog << "Failure constructing " << e1 << " from add.\n";
555                 ++result;
556         }
557
558         epvector v2;
559         v2.push_back(expair(x, 1));
560         v2.push_back(expair(1,-1));
561         ex e2 = mul(v2);  // Should be e==x;
562         if (!e2.is_equal(x)) {
563                 clog << "Failure constructing " << e2 << " from mul.\n";
564                 ++result;
565         }
566
567         return result;
568 }
569
570 unsigned exam_paranoia()
571 {
572         unsigned result = 0;
573         
574         cout << "examining several historic failures just out of paranoia" << flush;
575         
576         result += exam_paranoia1();  cout << '.' << flush;
577         result += exam_paranoia2();  cout << '.' << flush;
578         result += exam_paranoia3();  cout << '.' << flush;
579         result += exam_paranoia4();  cout << '.' << flush;
580         result += exam_paranoia5();  cout << '.' << flush;
581         result += exam_paranoia6();  cout << '.' << flush;
582         result += exam_paranoia7();  cout << '.' << flush;
583         result += exam_paranoia8();  cout << '.' << flush;
584         result += exam_paranoia9();  cout << '.' << flush;
585         result += exam_paranoia10();  cout << '.' << flush;
586         result += exam_paranoia11();  cout << '.' << flush;
587         result += exam_paranoia12();  cout << '.' << flush;
588         result += exam_paranoia13();  cout << '.' << flush;
589         result += exam_paranoia14();  cout << '.' << flush;
590         result += exam_paranoia15();  cout << '.' << flush;
591         result += exam_paranoia16();  cout << '.' << flush;
592         result += exam_paranoia17();  cout << '.' << flush;
593         result += exam_paranoia18();  cout << '.' << flush;
594         result += exam_paranoia19();  cout << '.' << flush;
595         result += exam_paranoia20();  cout << '.' << flush;
596         result += is_polynomial_false_positive(); cout << '.' << flush;
597         result += exam_paranoia21();  cout << '.' << flush;
598         result += exam_paranoia22();  cout << '.' << flush;
599         result += exam_paranoia23();  cout << '.' << flush;
600         
601         return result;
602 }
603
604 int main(int argc, char** argv)
605 {
606         return exam_paranoia();
607 }