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