src/paranoidnumber.cpp

   1 #include "paranoidnumber.h"
   2
   3 #include <sstream>
   4 #include <fenv.h>
   5 #include "log.h"
   6 #include <cassert>
   7 #include <iostream>
   8
   9 using namespace std;
  10 namespace IPDF
  11 {
  12 int64_t ParanoidNumber::g_count = 0;
  13
  14
  15 ParanoidNumber::~ParanoidNumber()
  16 {
  17         g_count--;
  18         for (int i = 0; i < NOP; ++i)
  19         {
  20                 for (auto n : m_next[i])
  21                         delete n;
  22         }
  23 }
  24
  25 ParanoidNumber::ParanoidNumber(const char * str) : m_value(0)
  26 {
  27         Construct();
  28         int dp = 0;
  29         int end = 0;
  30         while (str[dp] != '\0' && str[dp] != '.')
  31         {
  32                 ++dp;
  33                 ++end;
  34         }
  35         while (str[end] != '\0')
  36                 ++end;
  37         ParanoidNumber m(1);
  38         for (int i = dp-1; i >= 0; --i)
  39         {
  40                 ParanoidNumber b(str[i]-'0');
  41                 b*=m;
  42                 this->operator+=(b);
  43                 m*=10;
  44         }
  45         ParanoidNumber n(1);
  46         for (int i = dp+1; i < end; ++i)
  47         {
  48                 n/=10;
  49                 ParanoidNumber b(str[i]-'0');
  50                 b*=n;
  51                 this->operator+=(b);
  52         }
  53 }
  54
  55 ParanoidNumber & ParanoidNumber::operator=(const ParanoidNumber & a)
  56 {
  57         m_value = a.m_value;
  58         for (int i = 0; i < NOP; ++i)
  59         {
  60                 for (unsigned j = 0; j < m_next[i].size() && j < a.m_next[i].size(); ++j)
  61                 {
  62                         m_next[i][j]->operator=(*(a.m_next[i][j]));
  63                 }
  64
  65                 for (unsigned j = a.m_next[i].size(); j < m_next[i].size(); ++j)
  66                 {
  67                         delete m_next[i][j];
  68                 }
  69                 m_next[i].resize(a.m_next[i].size());
  70         }
  71         return *this;
  72 }
  73
  74
  75 string ParanoidNumber::Str() const
  76 {
  77         string result("");
  78         stringstream s;
  79         s << (double)m_value;
  80         result += s.str();
  81         for (auto mul : m_next[MULTIPLY])
  82         {
  83                 result += "*";
  84                 if (!mul->Floating())
  85                         result += "(" + mul->Str() + ")";
  86                 else
  87                         result += mul->Str();
  88         }
  89         for (auto div : m_next[DIVIDE])
  90         {
  91                 result += "/";
  92                 if (!div->Floating())
  93                         result += "(" + div->Str() + ")";
  94                 else
  95                         result += div->Str();
  96         }
  97
  98         for (auto add : m_next[ADD])
  99         {
 100                 result += "+";
 101                 if (!add->Floating())
 102                         result += "(" + add->Str() + ")";
 103                 else
 104                         result += add->Str();
 105         }
 106         for (auto sub : m_next[SUBTRACT])
 107         {
 108                 result += "-";
 109                 if (!sub->Floating())
 110                         result += "(" + sub->Str() + ")";
 111                 else
 112                         result += sub->Str();
 113         }
 114
 115
 116         return result;
 117 }
 118
 119 template <>
 120 bool TrustingOp<float>(float & a, const float & b, Optype op)
 121 {
 122         feclearexcept(FE_ALL_EXCEPT);
 123         switch (op)
 124         {
 125                 case ADD:
 126                         a += b;
 127                         break;
 128                 case SUBTRACT:
 129                         a -= b;
 130                         break;
 131                 case MULTIPLY:
 132                         a *= b;
 133                         break;
 134                 case DIVIDE:
 135                         a /= b;
 136                         break;
 137                 case NOP:
 138                         break;
 139         }
 140         return !fetestexcept(FE_ALL_EXCEPT);
 141 }
 142
 143 template <>
 144 bool TrustingOp<double>(double & a, const double & b, Optype op)
 145 {
 146         feclearexcept(FE_ALL_EXCEPT);
 147         switch (op)
 148         {
 149                 case ADD:
 150                         a += b;
 151                         break;
 152                 case SUBTRACT:
 153                         a -= b;
 154                         break;
 155                 case MULTIPLY:
 156                         a *= b;
 157                         break;
 158                 case DIVIDE:
 159                         a /= b;
 160                         break;
 161                 case NOP:
 162                         break;
 163         }
 164         return !fetestexcept(FE_ALL_EXCEPT);
 165 }
 166
 167 template <>
 168 bool TrustingOp<int8_t>(int8_t & a, const int8_t & b, Optype op)
 169 {
 170         int16_t sa(a);
 171         bool exact = true;
 172         switch (op)
 173         {
 174                 case ADD:
 175                         sa += b;
 176                         exact = (abs(sa) <= 127);
 177                         break;
 178                 case SUBTRACT:
 179                         sa -= b;
 180                         exact = (abs(sa) <= 127);
 181                         break;
 182                 case MULTIPLY:
 183                         sa *= b;
 184                         exact = (abs(sa) <= 127);
 185                         break;
 186                 case DIVIDE:
 187                         exact = (b != 0 && sa > b && sa % b == 0);
 188                         sa /= b;
 189                         break;
 190                 case NOP:
 191                         break;
 192         }
 193         a = (int8_t)(sa);
 194         return exact;
 195 }
 196
 197
 198 ParanoidNumber & ParanoidNumber::operator+=(const ParanoidNumber & a)
 199 {
 200         delete Operation(new ParanoidNumber(a), ADD);
 201         return *this;
 202 }
 203
 204
 205 ParanoidNumber & ParanoidNumber::operator-=(const ParanoidNumber & a)
 206 {
 207         delete Operation(new ParanoidNumber(a), SUBTRACT);
 208         return *this;
 209 }
 210
 211 ParanoidNumber & ParanoidNumber::operator*=(const ParanoidNumber & a)
 212 {
 213         delete Operation(new ParanoidNumber(a), MULTIPLY);
 214         return *this;
 215 }
 216
 217
 218 ParanoidNumber & ParanoidNumber::operator/=(const ParanoidNumber & a)
 219 {
 220         delete Operation(new ParanoidNumber(a), DIVIDE);
 221         return *this;
 222 }
 223
 224 // a + b
 225 ParanoidNumber * ParanoidNumber::OperationTerm(ParanoidNumber * b, Optype op, ParanoidNumber ** merge_point, Optype * merge_op)
 226 {
 227
 228         if (Floating() && m_value == 0) // 0 + b = b
 229         {
 230                 m_value = b->m_value;
 231                 if (op == SUBTRACT)
 232                 {
 233                         m_value = -m_value;
 234                         swap(b->m_next[ADD], b->m_next[SUBTRACT]);
 235                 }
 236
 237                 for (int i = 0; i < NOP; ++i)
 238                 {
 239                         m_next[i] = b->m_next[i];
 240                         b->m_next[i].clear();
 241                 }
 242                 return b;
 243         }
 244         if (b->Floating() && b->m_value == 0) // a + 0 = a
 245                 return b;
 246
 247
 248
 249         if (NoFactors() && b->NoFactors())
 250         {
 251                 if (ParanoidOp<digit_t>(m_value, b->m_value, op))
 252                 {
 253                         Optype addop = (op == ADD) ? ADD : SUBTRACT;
 254                         for (auto add : b->m_next[ADD])
 255                         {
 256                                 delete OperationTerm(add, addop);
 257                         }
 258                         Optype subop = (op == ADD) ? SUBTRACT : ADD;
 259                         for (auto sub : b->m_next[SUBTRACT])
 260                                 delete OperationTerm(sub, subop);
 261
 262                         b->m_next[ADD].clear();
 263                         b->m_next[SUBTRACT].clear();
 264                         return b;
 265                 }
 266         }
 267
 268
 269
 270
 271         bool parent = (merge_point == NULL);
 272         ParanoidNumber * merge = this;
 273         Optype mop = op;
 274         assert(mop != NOP); // silence compiler warning
 275         if (parent)
 276         {
 277                 merge_point = &merge;
 278                 merge_op = &mop;
 279         }
 280         else
 281         {
 282                 merge = *merge_point;
 283                 mop = *merge_op;
 284         }
 285
 286         Optype invop = InverseOp(op); // inverse of p
 287         Optype fwd = op;
 288         Optype rev = invop;
 289         if (op == SUBTRACT)
 290         {
 291                 fwd = ADD;
 292                 rev = SUBTRACT;
 293         }
 294
 295         for (auto prev : m_next[invop])
 296         {
 297                 if (prev->OperationTerm(b, rev, merge_point, merge_op) == b)
 298                         return b;
 299
 300         }
 301         for (auto next : m_next[op])
 302         {
 303                 if (next->OperationTerm(b, fwd, merge_point, merge_op) == b)
 304                         return b;
 305         }
 306
 307
 308
 309
 310         if (parent)
 311         {
 312                 merge->m_next[*merge_op].push_back(b);
 313         }
 314         else
 315         {
 316                 if (m_next[op].size() == 0)
 317                 {
 318                         *merge_point = this;
 319                         *merge_op = op;
 320                 }
 321         }
 322         return NULL;
 323 }
 324
 325 ParanoidNumber * ParanoidNumber::OperationFactor(ParanoidNumber * b, Optype op, ParanoidNumber ** merge_point, Optype * merge_op)
 326 {
 327
 328         if (Floating() && m_value == 0)
 329         {
 330                 return b;
 331         }
 332
 333         if (Floating() && m_value == 1 && op == MULTIPLY)
 334         {
 335                 m_value = b->m_value;
 336                 for (int i = 0; i < NOP; ++i)
 337                 {
 338                         for (auto n : m_next[i])
 339                                 delete n;
 340                         m_next[i].clear();
 341                         swap(m_next[i], b->m_next[i]);
 342                 }
 343                 return b;
 344         }
 345         if (b->Floating() && b->m_value == 1)
 346                 return b;
 347
 348         if (NoTerms() && b->NoTerms())
 349         {
 350                 if (ParanoidOp<digit_t>(m_value, b->m_value, op))
 351                 {
 352                         Optype mulop = (op == MULTIPLY) ? MULTIPLY : DIVIDE;
 353                         for (auto mul : b->m_next[MULTIPLY])
 354                         {
 355                                 delete OperationFactor(mul, mulop);
 356                         }
 357                         Optype divop = (op == MULTIPLY) ? DIVIDE : MULTIPLY;
 358                         for (auto div : b->m_next[DIVIDE])
 359                                 delete OperationFactor(div, divop);
 360
 361                         b->m_next[DIVIDE].clear();
 362                         b->m_next[MULTIPLY].clear();
 363                         return b;
 364                 }
 365         }
 366
 367
 368         bool parent = (merge_point == NULL);
 369         ParanoidNumber * merge = this;
 370         Optype mop = op;
 371         if (parent)
 372         {
 373                 merge_point = &merge;
 374                 merge_op = &mop;
 375         }
 376         else
 377         {
 378                 merge = *merge_point;
 379                 mop = *merge_op;
 380         }
 381
 382         Optype invop = InverseOp(op); // inverse of p
 383         Optype fwd = op;
 384         Optype rev = invop;
 385         if (op == DIVIDE)
 386         {
 387                 fwd = MULTIPLY;
 388                 rev = DIVIDE;
 389         }
 390
 391         ParanoidNumber * cpy_b = NULL;
 392
 393         if (m_next[ADD].size() > 0 || m_next[SUBTRACT].size() > 0)
 394         {
 395                 cpy_b = new ParanoidNumber(*b);
 396         }
 397
 398         for (auto prev : m_next[invop])
 399         {
 400                 if (prev->OperationFactor(b, rev, merge_point, merge_op) == b)
 401                 {
 402                         for (auto add : m_next[ADD])
 403                                 delete add->OperationFactor(new ParanoidNumber(*cpy_b), op);
 404                         for (auto sub : m_next[SUBTRACT])
 405                                 delete sub->OperationFactor(new ParanoidNumber(*cpy_b), op);
 406
 407                         delete cpy_b;
 408                         return b;
 409                 }
 410         }
 411         for (auto next : m_next[op])
 412         {
 413                 if (next->OperationFactor(b, fwd, merge_point, merge_op) == b)
 414                 {
 415                         for (auto add : m_next[ADD])
 416                                 delete add->OperationFactor(new ParanoidNumber(*cpy_b), op);
 417                         for (auto sub : m_next[SUBTRACT])
 418                                 delete sub->OperationFactor(new ParanoidNumber(*cpy_b), op);
 419                         delete cpy_b;
 420                         return b;
 421                 }
 422         }
 423
 424         if (parent)
 425         {
 426                 m_next[op].push_back(b);
 427                 for (auto add : m_next[ADD])
 428                         delete add->OperationFactor(new ParanoidNumber(*cpy_b), op);
 429                 for (auto sub : m_next[SUBTRACT])
 430                         delete sub->OperationFactor(new ParanoidNumber(*cpy_b), op);
 431         }
 432         return NULL;
 433 }
 434
 435
 436
 437 /**
 438  * Performs the operation on a with argument b (a += b, a -= b, a *= b, a /= b)
 439  * @returns b if b can safely be deleted
 440  * @returns NULL if b has been merged with a
 441  * append indicates that b should be merged
 442  */
 443 ParanoidNumber * ParanoidNumber::Operation(ParanoidNumber * b, Optype op, ParanoidNumber ** merge_point, Optype * merge_op)
 444 {
 445
 446         if (b == NULL)
 447                 return NULL;
 448
 449
 450         if (op == SUBTRACT || op == ADD)
 451                 return OperationTerm(b, op, merge_point, merge_op);
 452         if (op == MULTIPLY || op == DIVIDE)
 453                 return OperationFactor(b, op, merge_point, merge_op);
 454         return b;
 455 }
 456
 457
 458
 459 string ParanoidNumber::PStr() const
 460 {
 461         stringstream s;
 462         for (int i = 0; i < NOP; ++i)
 463         {
 464                 Optype f = Optype(i);
 465                 s << this;
 466                 for (auto n : m_next[f])
 467                 {
 468                         s << OpChar(f) << n->PStr();
 469                 }
 470         }
 471         return s.str();
 472 }
 473
 474 bool ParanoidNumber::Simplify(Optype op)
 475 {
 476         vector<ParanoidNumber*> next(0);
 477         swap(m_next[op], next);
 478         for (auto n : next)
 479         {
 480                 ParanoidNumber * result = Operation(n, op);
 481                 if (result != NULL)
 482                         delete result;
 483                 else
 484                         m_next[op].push_back(n);
 485         }
 486         return (next.size() > m_next[op].size());
 487 }
 488
 489
 490
 491
 492 }