src/tests/realops.cpp

   1 /**
   2  * Test mathematical operations on the Real type and consistency with double
   3  */
   4
   5 #include "main.h"
   6 #include "real.h"
   7 #include "progressbar.h"
   8
   9 using namespace std;
  10 using namespace IPDF;
  11
  12 #define TEST_CASES 1000
  13
  14 static double g_totalerror = 0;
  15
  16 bool NotEqual(double a, double b, double threshold=1e-4)
  17 {
  18         double error = fabs(a-b);
  19         g_totalerror += error;
  20         return (error > threshold);
  21 }
  22
  23 int main(int argc, char ** argv)
  24 {
  25         srand(time(NULL));
  26         DebugRealInfo();
  27
  28         unsigned failures = 0;
  29         Real acumulate(0);
  30         double dacumulate = 0;
  31
  32         for (unsigned i = 0; i < TEST_CASES; ++i)
  33         {
  34                 ProgressBar(i, TEST_CASES, 50);
  35                 //Debug("Test %u of %u", i, TEST_CASES);
  36                 double da = (double)(rand() + 1) / (double)(rand() + 1);
  37                 double db = (double)(rand() + 1) / (double)(rand() + 1);
  38
  39                 if (rand() % 2 == 0)
  40                         da = -da;
  41                 if (rand() % 2 == 0)
  42                         db = -db;
  43
  44                 Real a(da);
  45                 Real b(db);
  46                 Real aa(a);
  47                 Real bb(b);
  48                 unsigned old_failures = failures;
  49                 if (NotEqual(Double(a), Double(aa)))
  50                 {
  51                         failures++;
  52                         Warn("a != Real(a); %f vs %f", Double(a), Double(aa));
  53                 }
  54                 if (NotEqual(Double(b), Double(bb)))
  55                 {
  56                         failures++;
  57                         Warn("b != Real(b); %f vs %f", Double(b), Double(bb));
  58                 }
  59
  60                 if (NotEqual(Double(a), da))
  61                 {
  62                         failures++;
  63                         Warn("a != da; %f vs %f", Double(a), da);
  64                 }
  65                 if (NotEqual(Double(b), db))
  66                 {
  67                         failures++;
  68                         Warn("b != db; %f vs %f", Double(b), db);
  69                 }
  70                 if (NotEqual(Double(a+b), da+db))
  71                 {
  72                         failures++;
  73                         Warn("a + b = %f should be %f", Double(a+b), da+db);
  74                 }
  75                 if (NotEqual(Double(a-b), da-db))
  76                 {
  77                         failures++;
  78                         Warn("a - b = %f should be %f", Double(a-b), da-db);
  79                 }
  80                 if (NotEqual(Double(a*b), da*db))
  81                 {
  82                         failures++;
  83                         Warn("a * b = %f should be %f", Double(a*b), da*db);
  84                 }
  85                 if (NotEqual(Double(a/b), da/db))
  86                 {
  87                         failures++;
  88                         Warn("a / b = %f should be %f", Double(a/b), da/db);
  89                 }
  90
  91                 if (NotEqual(Double(a), da))
  92                 {
  93                         failures++;
  94                         Warn("a has changed after +-*/ from %f to %f", da, Double(a));
  95                 }
  96                 if (NotEqual(Double(b), db))
  97                 {
  98                         failures++;
  99                         Warn("b has changed after +-*/ from %f to %f", db, Double(b));
 100                 }
 101                 Real abeforeop(a);
 102                 if (NotEqual(Double(a+=b), da+=db))
 103                 {
 104                         failures++;
 105                         Warn("a += b = %f should be %f, a before op was %f", Double(a), da, Double(abeforeop));
 106                 }
 107                 abeforeop = a;
 108                 if (NotEqual(Double(a-=b), da-=db))
 109                 {
 110                         failures++;
 111                         Warn("a -= b = %f should be %f, a before op was %f", Double(a), da, Double(abeforeop));
 112                 }
 113                 abeforeop = a;
 114                 if (NotEqual(Double(a*=b), da*=db))
 115                 {
 116                         failures++;
 117                         Warn("a *= b = %f should be %f, a before op was %f", Double(a), da, Double(abeforeop));
 118                 }
 119                 abeforeop = a;
 120                 if (NotEqual(Double(a/=b), da/=db))
 121                 {
 122                         failures++;
 123                         Warn("a /= b = %f should be %f, a before op was %f", Double(a), da, Double(abeforeop));
 124                 }
 125                 if (NotEqual(Double(a*0.0 + 1.0), da*0.0 + 1.0))
 126                 {
 127                         failures++;
 128                         Warn("a * 0 = %f should be %f, a before op was %f", Double(a), da, Double(abeforeop));
 129                 }
 130
 131                 if (NotEqual(Double(a=b), da=db))
 132                 {
 133                         failures++;
 134                         Warn("a = b = %f should be %f, a before op was %f", Double(a), da, Double(abeforeop));
 135                 }
 136
 137                 if (NotEqual(Double(-a), -da))
 138                 {
 139                         failures++;
 140                         Warn("-a = %f should be %f, a before op was %f", Double(-a), -da, Double(abeforeop));
 141                 }
 142
 143                 if (NotEqual(Double(Sqrt(a)), Sqrt(da)))
 144                 {
 145                         failures++;
 146                         Warn("Sqrt(a) = %f should be %f, a before op was %f", Double(Sqrt(a)), Sqrt(da), Double(abeforeop));
 147                 }
 148
 149                 if (NotEqual(Double(a), da))
 150                 {
 151                         failures++;
 152                         Warn("a = %f, should be %f, a before ops was %f", Double(a), da, Double(abeforeop));
 153                 }
 154
 155                 switch (rand() % 4)
 156                 {
 157                         case 0:
 158                                 acumulate += a;
 159                                 dacumulate += da;
 160                                 break;
 161                         case 1:
 162                                 acumulate -= a;
 163                                 dacumulate -= da;
 164                                 break;
 165                         case 2:
 166                                 acumulate *= a;
 167                                 dacumulate *= da;
 168                                 break;
 169                         case 3:
 170                                 acumulate /= a;
 171                                 dacumulate /= da;
 172                                 break;
 173                 }
 174                 if (NotEqual(Double(acumulate), dacumulate))
 175                 {
 176                         Warn("Accumulated result %.30lf vs %.30lf is wrong", Double(acumulate), dacumulate);
 177                         failures++;
 178                 }
 179
 180                 if (failures > old_failures)
 181                 {
 182                         Error("%u failures on case %u da = %f, db = %f, a = %f, b = %f, aa = %f, bb = %f", failures-old_failures, i, da, db, Double(a), Double(b), Double(aa), Double(bb));
 183                         #if REAL == REAL_RATIONAL || REAL == REAL_RATIONAL_ARBINT
 184                                 Debug("\tStrings are a = %s, b = %s, aa = %s, bb = %s", a.Str().c_str(), b.Str().c_str(), aa.Str().c_str(), bb.Str().c_str());
 185                         #endif
 186                 }
 187         }
 188         Debug("Completed %u test cases with total of %u operations, %u failures", TEST_CASES, 18*TEST_CASES, failures);
 189         Debug("Total accumulated difference between Real and Double operations was %f", g_totalerror);
 190         Debug("Real: %.40lf", Double(acumulate));
 191         Debug("Doub: %.40lf", dacumulate);
 192         Debug("Diff: %.40lf", Double(acumulate) - dacumulate);
 193
 194 }