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         if (m_next_factor == NULL && a.Floating())
  72         {
  73                 if (ParanoidOp<digit_t>(m_value, a.m_value, ADD))
  74                 {
  75                         Simplify();
  76                         return *this;
  77                 }
  78         }
  79         ParanoidNumber * nt = m_next_term;
  80         ParanoidNumber * nf = m_next_factor;
  81
  82         ParanoidNumber ca(a);
  83         if (m_next_factor != NULL)
  84         {
  85                 if (m_next_factor->m_op == MULTIPLY)
  86                         ca /= (*m_next_factor);
  87                 else
  88                         ca *= (*m_next_factor);
  89
  90                 if (ca.Floating())
  91                 {
  92                         m_next_factor = NULL;
  93                         m_next_term = NULL;
  94                         operator+=(ca);
  95                         m_next_factor = nf;
  96                         m_next_term = nt;
  97                         Simplify();
  98                         return *this;
  99                 }
 100
 101         }
 102
 103         m_next_term = new ParanoidNumber(a, ADD);
 104         ParanoidNumber * t = m_next_term;
 105         while (t->m_next_term != NULL)
 106                 t = t->m_next_term;
 107         t->m_next_term = nt;
 108         //Debug("Simplify {%s} after add", Str().c_str());
 109         Simplify();
 110         return *this;
 111 }
 112
 113 ParanoidNumber & ParanoidNumber::operator-=(const ParanoidNumber & a)
 114 {
 115         // this = v + t + (a)
 116         // -> v + (a) + t
 117         if (m_next_factor == NULL && a.Floating())
 118         {
 119                 if (ParanoidOp<digit_t>(m_value, a.m_value, ADD))
 120                 {
 121                         Simplify();
 122                         return *this;
 123                 }
 124         }
 125
 126         ParanoidNumber * nt = m_next_term;
 127         ParanoidNumber * nf = m_next_factor;
 128
 129         ParanoidNumber ca(a, SUBTRACT);
 130         if (m_next_factor != NULL)
 131         {
 132                 if (m_next_factor->m_op == MULTIPLY)
 133                         ca /= (*m_next_factor);
 134                 else
 135                         ca *= (*m_next_factor);
 136
 137                 if (ca.Floating())
 138                 {
 139                         m_next_factor = NULL;
 140                         m_next_term = NULL;
 141                         operator-=(ca);
 142                         m_next_factor = nf;
 143                         m_next_term = nt;
 144                         Simplify();
 145                         return *this;
 146                 }
 147
 148         }
 149
 150         m_next_term = new ParanoidNumber(a,SUBTRACT);
 151         ParanoidNumber * t = m_next_term;
 152         while (t->m_next_term != NULL)
 153         {
 154                 t->m_op = SUBTRACT;
 155                 t = t->m_next_term;
 156         }
 157         t->m_op = SUBTRACT;
 158         //Debug("next term {%s}", m_next_term->Str().c_str());
 159         t->m_next_term = nt;
 160         //Debug("Simplify {%s} after sub", Str().c_str());
 161         Simplify();
 162         return *this;
 163 }
 164
 165 ParanoidNumber & ParanoidNumber::operator*=(const ParanoidNumber & a)
 166 {
 167
 168         //if (m_value == 0)
 169         //              return *this;
 170         //Debug("{%s} *= {%s}", Str().c_str(), a.Str().c_str());
 171         // this = (vf + t) * (a)
 172         if (a.Floating() && ParanoidOp<digit_t>(m_value, a.m_value, MULTIPLY))
 173         {
 174                 if (m_next_term != NULL)
 175                         m_next_term->operator*=(a);
 176                 Simplify();
 177                 return *this;
 178         }
 179
 180         ParanoidNumber * t = this;
 181         while (t->m_next_factor != NULL)
 182                 t = t->m_next_factor;
 183         t->m_next_factor = new ParanoidNumber(a, MULTIPLY);
 184
 185         if (m_next_term != NULL)
 186                 m_next_term->operator*=(a);
 187
 188         //Debug("Simplify {%s}", Str().c_str());
 189         Simplify();
 190         //Debug("Simplified to {%s}", Str().c_str());
 191         return *this;
 192 }
 193
 194
 195 ParanoidNumber & ParanoidNumber::operator/=(const ParanoidNumber & a)
 196 {
 197
 198
 199
 200         if (a.Floating() && ParanoidOp<digit_t>(m_value, a.m_value, DIVIDE))
 201         {
 202                 if (m_next_term != NULL)
 203                         m_next_term->operator/=(a);
 204                 Simplify();
 205                 return *this;
 206         }
 207
 208         //Debug("Called %s /= %s", Str().c_str(), a.Str().c_str());
 209         // this = (vf + t) * (a)
 210         ParanoidNumber * t = this;
 211         while (t->m_next_factor != NULL)
 212         {
 213                 t = t->m_next_factor;
 214         }
 215         t->m_next_factor = new ParanoidNumber(a, DIVIDE);
 216
 217         if (m_next_term != NULL)
 218                 m_next_term->operator/=(a);
 219
 220         Simplify();
 221         return *this;
 222 }
 223
 224
 225
 226 void ParanoidNumber::SimplifyTerms()
 227 {
 228
 229         //Debug("Simplify {%s}", Str().c_str());
 230         if (m_next_term == NULL)
 231         {
 232                 //Debug("No terms!");
 233                 return;
 234         }
 235
 236         for (ParanoidNumber * a = this; a != NULL; a = a->m_next_term)
 237         {
 238                 ParanoidNumber * b = a->m_next_term;
 239                 if (a->m_next_factor != NULL && !a->m_next_factor->Floating())
 240                 {
 241                         continue;
 242                 }
 243
 244                 ParanoidNumber * bprev = a;
 245                 while (b != NULL)
 246                 {
 247                         //Debug("Simplify factors of %s", b->Str().c_str());
 248                         b->SimplifyFactors();
 249                         if (b->m_next_factor != NULL && !b->m_next_factor->Floating())
 250                         {
 251                                 bprev = b;
 252                                 b = b->m_next_term;
 253                                 continue;
 254                         }
 255
 256                         bool simplify = false;
 257                         if (a->m_next_factor != NULL || b->m_next_factor != NULL)
 258                         {
 259                                 digit_t aa(a->Head<digit_t>());
 260                                 digit_t ab = (a->m_next_factor != NULL) ? a->m_next_factor->Head<digit_t>() : 1;
 261                                 digit_t bc(b->Head<digit_t>());
 262                                 digit_t bd = (b->m_next_factor != NULL) ? b->m_next_factor->Head<digit_t>() : 1;
 263                                 Optype aop = (a->m_next_factor != NULL) ? a->m_next_factor->m_op : DIVIDE;
 264                                 Optype cop = (b->m_next_factor != NULL) ? b->m_next_factor->m_op : DIVIDE;
 265                                 simplify = CombineTerms<digit_t>(aa, aop, ab, bc, cop, bd);
 266                                 if (simplify)
 267                                 {
 268                                         a->m_value = aa;
 269                                         if (a->m_next_factor != NULL)
 270                                                 a->m_next_factor->m_value = ab;
 271                                         else if (ab != 1)
 272                                         {
 273                                                 a->m_next_factor = b->m_next_factor;
 274                                                 b->m_next_factor = NULL;
 275                                                 a->m_next_factor->m_value = ab;
 276                                         }
 277                                 }
 278                         }
 279                         else
 280                         {
 281                                 simplify = ParanoidOp<digit_t>(a->m_value, b->Head<digit_t>(), ADD);
 282                         }
 283                         if (simplify)
 284                         {
 285                                 bprev->m_next_term = b->m_next_term;
 286                                 b->m_next_term = NULL;
 287                                 delete b;
 288                                 b = bprev;
 289                         }
 290
 291                         bprev = b;
 292                         b = b->m_next_term;
 293                 }
 294         }
 295 }
 296
 297 void ParanoidNumber::SimplifyFactors()
 298 {
 299
 300         //Debug("Simplify {%s}", Str().c_str());
 301         if (m_next_factor == NULL)
 302         {
 303                 //Debug("No factors!");
 304                 return;
 305         }
 306
 307         for (ParanoidNumber * a = this; a != NULL; a = a->m_next_factor)
 308         {
 309                 if ((a->m_op != ADD || a->m_op != SUBTRACT) && a->m_next_term != NULL)
 310                         continue;
 311
 312                 ParanoidNumber * bprev = a;
 313                 ParanoidNumber * b = a->m_next_factor;
 314                 while (b != NULL)
 315                 {
 316                         b->SimplifyTerms();
 317                         if (b->m_next_term != NULL)
 318                         {
 319                                 bprev = b;
 320                                 b = b->m_next_factor;
 321                                 continue;
 322                         }
 323
 324                         Optype op = b->m_op;
 325                         if (a->m_op == DIVIDE)
 326                         {
 327                                 op = (b->m_op == DIVIDE) ? MULTIPLY : DIVIDE;
 328                         }
 329
 330                         if (ParanoidOp<digit_t>(a->m_value, b->m_value, op))
 331                         {
 332
 333                                 bprev->m_next_factor = b->m_next_factor;
 334                                 b->m_next_factor = NULL;
 335                                 delete b;
 336                                 b = bprev;
 337                         }
 338                         bprev = b;
 339                         b = b->m_next_factor;
 340                 }
 341         }
 342 }
 343
 344 void ParanoidNumber::Simplify()
 345 {
 346         SimplifyFactors();
 347         SimplifyTerms();
 348 }
 349
 350 string ParanoidNumber::Str() const
 351 {
 352         string result("");
 353         stringstream s;
 354         s << (double)m_value;
 355
 356         if (m_next_factor != NULL)
 357         {
 358                 result += s.str();
 359                 result += OpChar(m_next_factor->m_op);
 360                 if (m_next_factor->m_next_term != NULL)
 361                         result += "(" + m_next_factor->Str() + ")";
 362                 else
 363                         result += m_next_factor->Str();
 364         }
 365         else
 366         {
 367                 result += s.str();
 368         }
 369
 370         if (m_next_term != NULL)
 371         {
 372                 result += " ";
 373                 result += OpChar(m_next_term->m_op);
 374                 result += m_next_term->Str();
 375         }
 376         return result;
 377 }
 378
 379 template <>
 380 bool TrustingOp<float>(float & a, const float & b, Optype op)
 381 {
 382         feclearexcept(FE_ALL_EXCEPT);
 383         switch (op)
 384         {
 385                 case ADD:
 386                         a += b;
 387                         break;
 388                 case SUBTRACT:
 389                         a -= b;
 390                         break;
 391                 case MULTIPLY:
 392                         a *= b;
 393                         break;
 394                 case DIVIDE:
 395                         a /= b;
 396                         break;
 397         }
 398         return !fetestexcept(FE_ALL_EXCEPT);
 399 }
 400
 401 template <>
 402 bool TrustingOp<double>(double & a, const double & b, Optype op)
 403 {
 404         feclearexcept(FE_ALL_EXCEPT);
 405         switch (op)
 406         {
 407                 case ADD:
 408                         a += b;
 409                         break;
 410                 case SUBTRACT:
 411                         a -= b;
 412                         break;
 413                 case MULTIPLY:
 414                         a *= b;
 415                         break;
 416                 case DIVIDE:
 417                         a /= b;
 418                         break;
 419         }
 420         return !fetestexcept(FE_ALL_EXCEPT);
 421 }
 422
 423 template <>
 424 bool TrustingOp<int8_t>(int8_t & a, const int8_t & b, Optype op)
 425 {
 426         int16_t sa(a);
 427         bool exact = true;
 428         switch (op)
 429         {
 430                 case ADD:
 431                         sa += b;
 432                         exact = (abs(sa) <= 127);
 433                         break;
 434                 case SUBTRACT:
 435                         sa -= b;
 436                         exact = (abs(sa) <= 127);
 437                         break;
 438                 case MULTIPLY:
 439                         sa *= b;
 440                         exact = (abs(sa) <= 127);
 441                         break;
 442                 case DIVIDE:
 443                         exact = (b != 0 && sa > b && sa % b == 0);
 444                         sa /= b;
 445                         break;
 446         }
 447         a = (int8_t)(sa);
 448         return exact;
 449 }
 450
 451 }