src/paranoidnumber.cpp

   1 #include "paranoidnumber.h"
   2
   3 #include <sstream>
   4 #include <fenv.h>
   5 #include "log.h"
   6 #include <iostream>
   7
   8 using namespace std;
   9 namespace IPDF
  10 {
  11 int64_t ParanoidNumber::g_count = 0;
  12
  13 ParanoidNumber::ParanoidNumber(const char * str) : m_value(0), m_op(ADD), m_next_term(NULL), m_next_factor(NULL)
  14 {
  15         int dp = 0;
  16         int end = 0;
  17         while (str[dp] != '\0' && str[dp] != '.')
  18         {
  19                 ++dp;
  20                 ++end;
  21         }
  22         while (str[end] != '\0')
  23                 ++end;
  24
  25         ParanoidNumber m(1);
  26         for (int i = dp-1; i >= 0; --i)
  27         {
  28                 ParanoidNumber b(str[i]-'0');
  29                 b*=m;
  30                 //Debug("m is %s", m.Str().c_str());
  31                 //Debug("Add %s", b.Str().c_str());
  32                 this->operator+=(b);
  33                 //Debug("Now at %s", Str().c_str());
  34                 m*=10;
  35         }
  36         ParanoidNumber n(1);
  37         for (int i = dp+1; i < end; ++i)
  38         {
  39                 n/=10;
  40                 ParanoidNumber b(str[i]-'0');
  41                 //Debug("%s * %s", b.Str().c_str(), n.Str().c_str());
  42                 b*=n;
  43                 //Debug("b -> %s", b.Str().c_str());
  44                 //Debug("Add %s", b.Str().c_str());
  45                 this->operator+=(b);
  46                 //Debug("Now at %s", Str().c_str());
  47
  48         }
  49         //Debug("Constructed {%s} from %s (%f)", Str().c_str(), str, ToDouble());
  50 }
  51
  52 ParanoidNumber & ParanoidNumber::operator=(const ParanoidNumber & a)
  53 {
  54         //TODO: Optimise
  55         delete m_next_term;
  56         delete m_next_factor;
  57         m_op = a.m_op;
  58         if (a.m_next_term != NULL)
  59         {
  60                 m_next_term = new ParanoidNumber(*(a.m_next_term));
  61         }
  62         if (a.m_next_factor != NULL)
  63         {
  64                 m_next_factor = new ParanoidNumber(*(a.m_next_factor));
  65         }
  66         return *this;
  67 }
  68
  69 ParanoidNumber & ParanoidNumber::operator+=(const ParanoidNumber & a)
  70 {
  71
  72         if (m_next_factor == NULL && a.Floating())
  73         {
  74                 if (ParanoidOp<digit_t>(m_value, a.m_value, ADD))
  75                 {
  76                         Simplify();
  77                         return *this;
  78                 }
  79         }
  80         ParanoidNumber * nt = m_next_term;
  81         ParanoidNumber * nf = m_next_factor;
  82
  83         ParanoidNumber ca(a);
  84         if (m_next_factor != NULL)
  85         {
  86                 if (m_next_factor->m_op == MULTIPLY)
  87                         ca /= (*m_next_factor);
  88                 else
  89                         ca *= (*m_next_factor);
  90
  91                 if (ca.Floating())
  92                 {
  93                         m_next_factor = NULL;
  94                         m_next_term = NULL;
  95                         operator+=(ca);
  96                         m_next_factor = nf;
  97                         m_next_term = nt;
  98                         Simplify();
  99                         return *this;
 100                 }
 101
 102         }
 103
 104         m_next_term = new ParanoidNumber(a, ADD);
 105         ParanoidNumber * t = m_next_term;
 106         while (t->m_next_term != NULL)
 107                 t = t->m_next_term;
 108         t->m_next_term = nt;
 109         //Debug("Simplify {%s} after add", Str().c_str());
 110         Simplify();
 111         return *this;
 112 }
 113
 114 ParanoidNumber & ParanoidNumber::operator-=(const ParanoidNumber & a)
 115 {
 116         // this = v + t + (a)
 117         // -> v + (a) + t
 118         if (m_next_factor == NULL && a.Floating())
 119         {
 120                 if (ParanoidOp<digit_t>(m_value, a.m_value, ADD))
 121                 {
 122                         Simplify();
 123                         return *this;
 124                 }
 125         }
 126
 127         ParanoidNumber * nt = m_next_term;
 128         ParanoidNumber * nf = m_next_factor;
 129
 130         ParanoidNumber ca(a, SUBTRACT);
 131         if (m_next_factor != NULL)
 132         {
 133                 if (m_next_factor->m_op == MULTIPLY)
 134                         ca /= (*m_next_factor);
 135                 else
 136                         ca *= (*m_next_factor);
 137
 138                 if (ca.Floating())
 139                 {
 140                         m_next_factor = NULL;
 141                         m_next_term = NULL;
 142                         operator-=(ca);
 143                         m_next_factor = nf;
 144                         m_next_term = nt;
 145                         Simplify();
 146                         return *this;
 147                 }
 148
 149         }
 150
 151         m_next_term = new ParanoidNumber(a,SUBTRACT);
 152         ParanoidNumber * t = m_next_term;
 153         while (t->m_next_term != NULL)
 154         {
 155                 t->m_op = SUBTRACT;
 156                 t = t->m_next_term;
 157         }
 158         t->m_op = SUBTRACT;
 159         //Debug("next term {%s}", m_next_term->Str().c_str());
 160         t->m_next_term = nt;
 161         //Debug("Simplify {%s} after sub", Str().c_str());
 162         Simplify();
 163         return *this;
 164 }
 165
 166 ParanoidNumber & ParanoidNumber::operator*=(const ParanoidNumber & a)
 167 {
 168
 169         //if (m_value == 0)
 170         //              return *this;
 171         //Debug("{%s} *= {%s}", Str().c_str(), a.Str().c_str());
 172         // this = (vf + t) * (a)
 173         if (a.Floating() && ParanoidOp<digit_t>(m_value, a.m_value, MULTIPLY))
 174         {
 175                 if (m_next_term != NULL)
 176                         m_next_term->operator*=(a);
 177                 Simplify();
 178                 return *this;
 179         }
 180
 181         ParanoidNumber * t = this;
 182         while (t->m_next_factor != NULL)
 183                 t = t->m_next_factor;
 184         t->m_next_factor = new ParanoidNumber(a, MULTIPLY);
 185
 186         if (m_next_term != NULL)
 187                 m_next_term->operator*=(a);
 188
 189         //Debug("Simplify {%s}", Str().c_str());
 190         Simplify();
 191         //Debug("Simplified to {%s}", Str().c_str());
 192         return *this;
 193 }
 194
 195
 196 ParanoidNumber & ParanoidNumber::operator/=(const ParanoidNumber & a)
 197 {
 198
 199
 200
 201         if (a.Floating() && ParanoidOp<digit_t>(m_value, a.m_value, DIVIDE))
 202         {
 203                 if (m_next_term != NULL)
 204                         m_next_term->operator/=(a);
 205                 Simplify();
 206                 return *this;
 207         }
 208
 209         //Debug("Called %s /= %s", Str().c_str(), a.Str().c_str());
 210         // this = (vf + t) * (a)
 211         ParanoidNumber * t = this;
 212         while (t->m_next_factor != NULL)
 213         {
 214                 t = t->m_next_factor;
 215         }
 216         t->m_next_factor = new ParanoidNumber(a, DIVIDE);
 217
 218         if (m_next_term != NULL)
 219                 m_next_term->operator/=(a);
 220
 221         Simplify();
 222         return *this;
 223 }
 224
 225
 226
 227 void ParanoidNumber::SimplifyTerms()
 228 {
 229
 230         //Debug("Simplify {%s}", Str().c_str());
 231         if (m_next_term == NULL)
 232         {
 233                 //Debug("No terms!");
 234                 return;
 235         }
 236
 237         for (ParanoidNumber * a = this; a != NULL; a = a->m_next_term)
 238         {
 239                 ParanoidNumber * b = a->m_next_term;
 240                 if (a->m_next_factor != NULL)
 241                 {
 242                         continue;
 243                 }
 244
 245                 ParanoidNumber * bprev = a;
 246                 while (b != NULL)
 247                 {
 248                         //Debug("Simplify factors of %s", b->Str().c_str());
 249                         b->SimplifyFactors();
 250                         if (b->m_next_factor != NULL)
 251                         {
 252                                 bprev = b;
 253                                 b = b->m_next_term;
 254                                 continue;
 255                         }
 256                         bool simplify = false;
 257                         simplify = ParanoidOp<digit_t>(a->m_value, b->Head<digit_t>(), ADD);
 258                         if (simplify)
 259                         {
 260                                 bprev->m_next_term = b->m_next_term;
 261                                 b->m_next_term = NULL;
 262                                 delete b;
 263                                 b = bprev;
 264                         }
 265
 266                         bprev = b;
 267                         b = b->m_next_term;
 268                 }
 269         }
 270 }
 271
 272 void ParanoidNumber::SimplifyFactors()
 273 {
 274
 275         //Debug("Simplify {%s}", Str().c_str());
 276         if (m_next_factor == NULL)
 277         {
 278                 //Debug("No factors!");
 279                 return;
 280         }
 281
 282         for (ParanoidNumber * a = this; a != NULL; a = a->m_next_factor)
 283         {
 284                 if ((a->m_op != ADD || a->m_op != SUBTRACT) && a->m_next_term != NULL)
 285                         continue;
 286
 287                 ParanoidNumber * bprev = a;
 288                 ParanoidNumber * b = a->m_next_factor;
 289                 while (b != NULL)
 290                 {
 291                         b->SimplifyTerms();
 292                         if (b->m_next_term != NULL)
 293                         {
 294                                 bprev = b;
 295                                 b = b->m_next_factor;
 296                                 continue;
 297                         }
 298
 299                         Optype op = b->m_op;
 300                         if (a->m_op == DIVIDE)
 301                         {
 302                                 op = (b->m_op == DIVIDE) ? MULTIPLY : DIVIDE;
 303                         }
 304
 305                         if (ParanoidOp<digit_t>(a->m_value, b->m_value, op))
 306                         {
 307
 308                                 bprev->m_next_factor = b->m_next_factor;
 309                                 b->m_next_factor = NULL;
 310                                 delete b;
 311                                 b = bprev;
 312                         }
 313                         bprev = b;
 314                         b = b->m_next_factor;
 315                 }
 316         }
 317 }
 318
 319 void ParanoidNumber::Simplify()
 320 {
 321         SimplifyFactors();
 322         SimplifyTerms();
 323 }
 324
 325 string ParanoidNumber::Str() const
 326 {
 327         string result("");
 328         stringstream s;
 329         s << (double)m_value;
 330
 331         if (m_next_factor != NULL)
 332         {
 333                 result += s.str();
 334                 result += OpChar(m_next_factor->m_op);
 335                 if (m_next_factor->m_next_term != NULL)
 336                         result += "(" + m_next_factor->Str() + ")";
 337                 else
 338                         result += m_next_factor->Str();
 339         }
 340         else
 341         {
 342                 result += s.str();
 343         }
 344
 345         if (m_next_term != NULL)
 346         {
 347                 result += " ";
 348                 result += OpChar(m_next_term->m_op);
 349                 result += m_next_term->Str();
 350         }
 351         return result;
 352 }
 353
 354 template <>
 355 bool TrustingOp<float>(float & a, const float & b, Optype op)
 356 {
 357         feclearexcept(FE_ALL_EXCEPT);
 358         switch (op)
 359         {
 360                 case ADD:
 361                         a += b;
 362                         break;
 363                 case SUBTRACT:
 364                         a -= b;
 365                         break;
 366                 case MULTIPLY:
 367                         a *= b;
 368                         break;
 369                 case DIVIDE:
 370                         a /= b;
 371                         break;
 372         }
 373         return !fetestexcept(FE_ALL_EXCEPT);
 374 }
 375
 376 template <>
 377 bool TrustingOp<double>(double & a, const double & b, Optype op)
 378 {
 379         feclearexcept(FE_ALL_EXCEPT);
 380         switch (op)
 381         {
 382                 case ADD:
 383                         a += b;
 384                         break;
 385                 case SUBTRACT:
 386                         a -= b;
 387                         break;
 388                 case MULTIPLY:
 389                         a *= b;
 390                         break;
 391                 case DIVIDE:
 392                         a /= b;
 393                         break;
 394         }
 395         return !fetestexcept(FE_ALL_EXCEPT);
 396 }
 397
 398 template <>
 399 bool TrustingOp<int8_t>(int8_t & a, const int8_t & b, Optype op)
 400 {
 401         int16_t sa(a);
 402         bool exact = true;
 403         switch (op)
 404         {
 405                 case ADD:
 406                         sa += b;
 407                         exact = (abs(sa) <= 127);
 408                         break;
 409                 case SUBTRACT:
 410                         sa -= b;
 411                         exact = (abs(sa) <= 127);
 412                         break;
 413                 case MULTIPLY:
 414                         sa *= b;
 415                         exact = (abs(sa) <= 127);
 416                         break;
 417                 case DIVIDE:
 418                         exact = (b != 0 && sa > b && sa % b == 0);
 419                         sa /= b;
 420                         break;
 421         }
 422         a = (int8_t)(sa);
 423         return exact;
 424 }
 425
 426 }