check/exam_sqrfree.cpp

   1 /** @file exam_sqrfree.cpp
   2  *
   3  */
   4
   5 /*
   6  *  GiNaC Copyright (C) 1999-2022 Johannes Gutenberg University Mainz, Germany
   7  *
   8  *  This program is free software; you can redistribute it and/or modify
   9  *  it under the terms of the GNU General Public License as published by
  10  *  the Free Software Foundation; either version 2 of the License, or
  11  *  (at your option) any later version.
  12  *
  13  *  This program is distributed in the hope that it will be useful,
  14  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  15  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  16  *  GNU General Public License for more details.
  17  *
  18  *  You should have received a copy of the GNU General Public License
  19  *  along with this program; if not, write to the Free Software
  20  *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  21  */
  22
  23 #include "ginac.h"
  24 using namespace GiNaC;
  25
  26 #include <iostream>
  27 using namespace std;
  28
  29 static unsigned exam_sqrfree1()
  30 {
  31         unsigned result = 0;
  32         symbol x("x");
  33         ex e1, e2;
  34
  35         e1 = (1+x)*pow((2+x),2)*pow((3+x),3)*pow((4+x),4);
  36         e2 = sqrfree(expand(e1), lst{x});
  37         if (e1 != e2) {
  38                 clog << "sqrfree(expand(" << e1 << ")) erroneously returned "
  39                      << e2 << endl;
  40                 ++result;
  41         }
  42
  43         return result;
  44 }
  45
  46 static unsigned exam_sqrfree2()
  47 {
  48         unsigned result = 0;
  49         symbol x("x"), y("y");
  50         ex e1, e2;
  51
  52         e1 = (x+y)*pow((x+2*y),2)*pow((x+3*y),3)*pow((x+4*y),4);
  53         e2 = sqrfree(expand(e1));
  54         if (e1 != e2) {
  55                 clog << "sqrfree(expand(" << e1 << ")) erroneously returned "
  56                      << e2 << endl;
  57                 ++result;
  58         }
  59         e2 = sqrfree(expand(e1), lst{x});
  60         if (e1 != e2) {
  61                 clog << "sqrfree(expand(" << e1 << "),[x]) erroneously returned "
  62                      << e2 << endl;
  63                 ++result;
  64         }
  65         e2 = sqrfree(expand(e1), lst{y});
  66         if (e1 != e2) {
  67                 clog << "sqrfree(expand(" << e1 << "),[y]) erroneously returned "
  68                      << e2 << endl;
  69                 ++result;
  70         }
  71         e2 = sqrfree(expand(e1), lst{x,y});
  72         if (e1 != e2) {
  73                 clog << "sqrfree(expand(" << e1 << "),[x,y]) erroneously returned "
  74                      << e2 << endl;
  75                 ++result;
  76         }
  77
  78         return result;
  79 }
  80
  81 static unsigned exam_sqrfree3()
  82 {
  83         unsigned result = 0;
  84         symbol x("x"), y("y"), z("z");
  85         ex e1, e2;
  86
  87         e1 = (x+y)*pow(x, 2)*(-z-1);
  88         e2 = sqrfree(expand(e1));
  89         if (!expand(e1 - e2).is_zero()) {
  90                 clog << "sqrfree(expand(" << e1 << ")) erroneously returned "
  91                      << e2 << endl;
  92                 ++result;
  93         }
  94
  95         e1 = (x+y)*pow(x, 3)*(-z-1);
  96         e2 = sqrfree(expand(e1));
  97         if (!expand(e1 - e2).is_zero()) {
  98                 clog << "sqrfree(expand(" << e1 << ")) erroneously returned "
  99                      << e2 << endl;
 100                 ++result;
 101         }
 102
 103         return result;
 104 }
 105
 106 // Bug in sqrfree_yun (fixed 2016-02-02).
 107 static unsigned exam_hidden_zero1()
 108 {
 109         unsigned result = 0;
 110         symbol x("x");
 111         ex e;
 112
 113         e = (x-1)*(x+1) - x*x + 1;  // an unexpanded 0...
 114         try {
 115                 ex f = sqrfree(e);
 116                 if (!f.is_zero()) {
 117                         clog << "sqrfree(" << e << ") returns " << f << " instead of 0\n";
 118                         ++result;
 119                 }
 120         } catch (const exception &err) {
 121                 clog << "sqrfree(" << e << ") throws " << err.what() << endl;
 122                 ++result;
 123         }
 124
 125         e = pow(x-1,3) - expand(pow(x-1,3));  // ...still after differentiating...
 126         try {
 127                 ex f = sqrfree(e);
 128                 if (!f.is_zero()) {
 129                         clog << "sqrfree(" << e << ") returns " << f << " instead of 0\n";
 130                         ++result;
 131                 }
 132         } catch (const exception &err) {
 133                 clog << "sqrfree(" << e << ") throws " << err.what() << endl;
 134                 ++result;
 135         }
 136
 137         e = pow(x-1,4) - expand(pow(x-1,4));  // ...and after differentiating twice.
 138         try {
 139                 ex f = sqrfree(e);
 140                 if (!f.is_zero()) {
 141                         clog << "sqrfree(" << e << ") returns " << f << " instead of 0\n";
 142                         ++result;
 143                 }
 144         } catch (const exception &err) {
 145                 clog << "sqrfree(" << e << ") throws " << err.what() << endl;
 146                 ++result;
 147         }
 148
 149         return result;
 150 }
 151
 152 static unsigned exam_hidden_zero2()
 153 {
 154         unsigned result = 0;
 155         symbol x("x"), y("y");
 156         ex e1, e2;
 157
 158         e1 = (1 + 3*x + 3*pow(x,2) + pow(x,3) - pow(1+x,3)) * y;
 159         e2 = sqrfree(e1);
 160         if (!e2.is_zero()) {
 161                 clog << "sqrfree(" << e1 << ") erroneously returned "
 162                      << e2 << endl;
 163                 ++result;
 164         }
 165
 166         e1 = (pow(x,2)-2*x*y+pow(y,2)-pow(x-y,2)) * x;
 167         e2 = sqrfree(e1);
 168         if (!e2.is_zero()) {
 169                 clog << "sqrfree(" << e1 << ") erroneously returned "
 170                      << e2 << endl;
 171                 ++result;
 172         }
 173
 174         e1 = (pow(x,2)-2*x*y+pow(y,2)-pow(x-y,2)) * (x+y);
 175         e2 = sqrfree(e1);
 176         if (!e2.is_zero()) {
 177                 clog << "sqrfree(" << e1 << ") erroneously returned "
 178                      << e2 << endl;
 179                 ++result;
 180         }
 181
 182         return result;
 183 }
 184
 185 unsigned exam_sqrfree()
 186 {
 187         unsigned result = 0;
 188
 189         cout << "examining square-free factorization" << flush;
 190
 191         result += exam_sqrfree1();  cout << '.' << flush;
 192         result += exam_sqrfree2();  cout << '.' << flush;
 193         result += exam_sqrfree3();  cout << '.' << flush;
 194         result += exam_hidden_zero1();  cout << '.' << flush;
 195         result += exam_hidden_zero2();  cout << '.' << flush;
 196
 197         return result;
 198 }
 199
 200 unsigned exam_sqrfree_parfrac()
 201 {
 202         symbol x("x");
 203         // (expression, number of terms after partial fraction decomposition)
 204         vector<pair<ex, unsigned>> exams = {
 205                 {ex("(x - 1) / (x^2*(x^2 + 2))", lst{x}), 3},
 206                 {ex("(1 - x^10) / x", lst{x}), 2},
 207                 {ex("(2*x^3 + x + 3) / ((x^2 + 1)^2)", lst{x}), 2},
 208                 {ex("1 / (x * (x+1)^2 * (x+2)^3)", lst{x}), 6},
 209                 {ex("(x*x + 3*x - 1) / (x^2*(x^2 + 2)^3)", lst{x}), 5},
 210                 {ex("(1 - x^10) / (x + 2)", lst{x}), 11},
 211                 {ex("(1 - x + 3*x^2) / (x^3 * (2+x)^2)", lst{x}), 5},
 212                 {ex("(1 - x) / (x^4 * (x - 2)^3)", lst{x}), 6},
 213                 {ex("(1 - 2*x + x^9) / (x^5 * (1 - x + x^2)^6)", lst{x}), 11}
 214         };
 215         unsigned result = 0;
 216
 217         cout << "\n"
 218              << "examining square-free partial fraction decomposition" << flush;
 219         for (auto e: exams) {
 220                 ex e1 = e.first;
 221                 ex e2 = sqrfree_parfrac(e1, x);
 222                 if (e2.nops() != e.second ||
 223                     !is_a<add>(e2) ||
 224                     !normal(e1-e2).is_zero()) {
 225                         clog << "sqrfree_parfrac(" << e1 << ", " << x << ") erroneously returned "
 226                              << e2 << endl;
 227                         ++result;
 228                 }
 229                 cout << '.' << flush;
 230         }
 231
 232         return result;
 233 }
 234
 235 int main(int argc, char** argv)
 236 {
 237         unsigned result = 0;
 238
 239         result += exam_sqrfree();
 240         result += exam_sqrfree_parfrac();
 241
 242         return result;
 243 }