1 /* 2 * ReactOS Calc (Math functions, IEEE-754 engine) 3 * 4 * Copyright 2007-2017, Carlo Bramini 5 * 6 * This program is free software; you can redistribute it and/or 7 * modify it under the terms of the GNU Lesser General Public 8 * License as published by the Free Software Foundation; either 9 * version 2 of the License, or (at your option) any later version. 10 * 11 * This library is distributed in the hope that it will be useful, 12 * but WITHOUT ANY WARRANTY; without even the implied warranty of 13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 14 * Lesser General Public License for more details. 15 * 16 * You should have received a copy of the GNU Lesser General Public 17 * License along with this library; if not, write to the Free Software 18 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 19 */ 20 21 #include "calc.h" 22 23 static double validate_rad2angle(double a); 24 static double validate_angle2rad(calc_number_t *c); 25 26 void apply_int_mask(calc_number_t *r) 27 { 28 unsigned __int64 mask; 29 30 switch (calc.size) { 31 case IDC_RADIO_QWORD: 32 mask = _UI64_MAX; 33 break; 34 case IDC_RADIO_DWORD: 35 mask = ULONG_MAX; 36 break; 37 case IDC_RADIO_WORD: 38 mask = USHRT_MAX; 39 break; 40 case IDC_RADIO_BYTE: 41 mask = UCHAR_MAX; 42 break; 43 default: 44 mask = (unsigned __int64)-1; 45 } 46 r->i &= mask; 47 } 48 49 double asinh(double x) 50 { 51 return log(x+sqrt(x*x+1)); 52 } 53 54 double acosh(double x) 55 { 56 // must be x>=1, if not return Nan (Not a Number) 57 if(!(x>=1.0)) return sqrt(-1.0); 58 59 // return only the positive result (as sqrt does). 60 return log(x+sqrt(x*x-1.0)); 61 } 62 63 double atanh(double x) 64 { 65 // must be x>-1, x<1, if not return Nan (Not a Number) 66 if(!(x>-1.0 && x<1.0)) return sqrt(-1.0); 67 68 return log((1.0+x)/(1.0-x))/2.0; 69 } 70 71 static double validate_rad2angle(double a) 72 { 73 switch (calc.degr) { 74 case IDC_RADIO_DEG: 75 a = a * (180.0/CALC_PI); 76 break; 77 case IDC_RADIO_RAD: 78 break; 79 case IDC_RADIO_GRAD: 80 a = a * (200.0/CALC_PI); 81 break; 82 } 83 return a; 84 } 85 86 static double validate_angle2rad(calc_number_t *c) 87 { 88 switch (calc.degr) { 89 case IDC_RADIO_DEG: 90 c->f = c->f * (CALC_PI/180.0); 91 break; 92 case IDC_RADIO_RAD: 93 break; 94 case IDC_RADIO_GRAD: 95 c->f = c->f * (CALC_PI/200.0); 96 break; 97 } 98 return c->f; 99 } 100 101 void rpn_sin(calc_number_t *c) 102 { 103 double angle = validate_angle2rad(c); 104 105 if (angle == 0 || angle == CALC_PI) 106 c->f = 0; 107 else 108 if (angle == CALC_3_PI_2) 109 c->f = -1; 110 else 111 if (angle == CALC_2_PI) 112 c->f = 1; 113 else 114 c->f = sin(angle); 115 } 116 void rpn_cos(calc_number_t *c) 117 { 118 double angle = validate_angle2rad(c); 119 120 if (angle == CALC_PI_2 || angle == CALC_3_PI_2) 121 c->f = 0; 122 else 123 if (angle == CALC_PI) 124 c->f = -1; 125 else 126 if (angle == CALC_2_PI) 127 c->f = 1; 128 else 129 c->f = cos(angle); 130 } 131 void rpn_tan(calc_number_t *c) 132 { 133 double angle = validate_angle2rad(c); 134 135 if (angle == CALC_PI_2 || angle == CALC_3_PI_2) 136 calc.is_nan = TRUE; 137 else 138 if (angle == CALC_PI || angle == CALC_2_PI) 139 c->f = 0; 140 else 141 c->f = tan(angle); 142 } 143 144 void rpn_asin(calc_number_t *c) 145 { 146 c->f = validate_rad2angle(asin(c->f)); 147 if (_isnan(c->f)) 148 calc.is_nan = TRUE; 149 } 150 void rpn_acos(calc_number_t *c) 151 { 152 c->f = validate_rad2angle(acos(c->f)); 153 if (_isnan(c->f)) 154 calc.is_nan = TRUE; 155 } 156 void rpn_atan(calc_number_t *c) 157 { 158 c->f = validate_rad2angle(atan(c->f)); 159 if (_isnan(c->f)) 160 calc.is_nan = TRUE; 161 } 162 163 void rpn_sinh(calc_number_t *c) 164 { 165 c->f = sinh(c->f); 166 if (_isnan(c->f)) 167 calc.is_nan = TRUE; 168 } 169 void rpn_cosh(calc_number_t *c) 170 { 171 c->f = cosh(c->f); 172 if (_isnan(c->f)) 173 calc.is_nan = TRUE; 174 } 175 void rpn_tanh(calc_number_t *c) 176 { 177 c->f = tanh(c->f); 178 if (_isnan(c->f)) 179 calc.is_nan = TRUE; 180 } 181 182 void rpn_asinh(calc_number_t *c) 183 { 184 c->f = asinh(c->f); 185 if (_isnan(c->f)) 186 calc.is_nan = TRUE; 187 } 188 void rpn_acosh(calc_number_t *c) 189 { 190 c->f = acosh(c->f); 191 if (_isnan(c->f)) 192 calc.is_nan = TRUE; 193 } 194 void rpn_atanh(calc_number_t *c) 195 { 196 c->f = atanh(c->f); 197 if (_isnan(c->f)) 198 calc.is_nan = TRUE; 199 } 200 201 void rpn_int(calc_number_t *c) 202 { 203 double int_part; 204 205 modf(calc.code.f, &int_part); 206 c->f = int_part; 207 } 208 209 void rpn_frac(calc_number_t *c) 210 { 211 double int_part; 212 213 c->f = modf(calc.code.f, &int_part); 214 } 215 216 void rpn_reci(calc_number_t *c) 217 { 218 if (c->f == 0) 219 calc.is_nan = TRUE; 220 else 221 c->f = 1./c->f; 222 } 223 224 void rpn_fact(calc_number_t *c) 225 { 226 double fact, mult, num; 227 228 if (calc.base == IDC_RADIO_DEC) 229 num = c->f; 230 else 231 num = (double)c->i; 232 if (num > 1000) { 233 calc.is_nan = TRUE; 234 return; 235 } 236 if (num < 0) { 237 calc.is_nan = TRUE; 238 return; 239 } else 240 if (num == 0) 241 fact = 1; 242 else { 243 rpn_int(c); 244 fact = 1; 245 mult = 2; 246 while (mult <= num) { 247 fact *= mult; 248 mult++; 249 } 250 c->f = fact; 251 } 252 if (_finite(fact) == 0) 253 calc.is_nan = TRUE; 254 else 255 if (calc.base == IDC_RADIO_DEC) 256 c->f = fact; 257 else 258 c->i = (__int64)fact; 259 } 260 261 __int64 logic_dbl2int(calc_number_t *a) 262 { 263 double int_part; 264 int width; 265 266 modf(a->f, &int_part); 267 width = (int_part==0) ? 1 : (int)log10(fabs(int_part))+1; 268 if (width > 63) { 269 calc.is_nan = TRUE; 270 return 0; 271 } 272 return (__int64)int_part; 273 } 274 275 double logic_int2dbl(calc_number_t *a) 276 { 277 return (double)a->i; 278 } 279 280 void rpn_not(calc_number_t *c) 281 { 282 if (calc.base == IDC_RADIO_DEC) { 283 calc_number_t n; 284 n.i = logic_dbl2int(c); 285 c->f = (long double)(~n.i); 286 } else 287 c->i = ~c->i; 288 } 289 290 void rpn_pi(calc_number_t *c) 291 { 292 c->f = CALC_PI; 293 } 294 295 void rpn_2pi(calc_number_t *c) 296 { 297 c->f = CALC_PI*2; 298 } 299 300 void rpn_sign(calc_number_t *c) 301 { 302 if (calc.base == IDC_RADIO_DEC) 303 c->f = 0-c->f; 304 else 305 c->i = 0-c->i; 306 } 307 308 void rpn_exp2(calc_number_t *c) 309 { 310 if (calc.base == IDC_RADIO_DEC) { 311 c->f *= c->f; 312 if (_finite(c->f) == 0) 313 calc.is_nan = TRUE; 314 } else 315 c->i *= c->i; 316 } 317 318 void rpn_exp3(calc_number_t *c) 319 { 320 if (calc.base == IDC_RADIO_DEC) { 321 c->f = pow(c->f, 3.); 322 if (_finite(c->f) == 0) 323 calc.is_nan = TRUE; 324 } else 325 c->i *= (c->i*c->i); 326 } 327 328 static __int64 myabs64(__int64 number) 329 { 330 return (number < 0) ? 0-number : number; 331 } 332 333 static unsigned __int64 sqrti(unsigned __int64 number) 334 { 335 /* modified form of Newton's method for approximating roots */ 336 #define NEXT(n, i) (((n) + (i)/(n)) >> 1) 337 unsigned __int64 n, n1; 338 339 #ifdef __GNUC__ 340 if (number == 0xffffffffffffffffULL) 341 #else 342 if (number == 0xffffffffffffffffUI64) 343 #endif 344 return 0xffffffff; 345 346 n = 1; 347 n1 = NEXT(n, number); 348 while (myabs64(n1 - n) > 1) { 349 n = n1; 350 n1 = NEXT(n, number); 351 } 352 while((n1*n1) > number) 353 n1--; 354 return n1; 355 #undef NEXT 356 } 357 358 void rpn_sqrt(calc_number_t *c) 359 { 360 if (calc.base == IDC_RADIO_DEC) { 361 if (c->f < 0) 362 calc.is_nan = TRUE; 363 else 364 c->f = sqrt(c->f); 365 } else { 366 c->i = sqrti(c->i); 367 } 368 } 369 370 static __int64 cbrti(__int64 x) { 371 __int64 s, y, b; 372 373 s = 60; 374 y = 0; 375 while(s >= 0) { 376 y = 2*y; 377 b = (3*y*(y + 1) + 1) << s; 378 s = s - 3; 379 if (x >= b) { 380 x = x - b; 381 y = y + 1; 382 } 383 } 384 return y; 385 } 386 387 void rpn_cbrt(calc_number_t *c) 388 { 389 if (calc.base == IDC_RADIO_DEC) 390 #if defined(__GNUC__) && !defined(__REACTOS__) 391 c->f = cbrt(c->f); 392 #else 393 c->f = pow(c->f,1./3.); 394 #endif 395 else { 396 c->i = cbrti(c->i); 397 } 398 } 399 400 void rpn_exp(calc_number_t *c) 401 { 402 c->f = exp(c->f); 403 if (_finite(c->f) == 0) 404 calc.is_nan = TRUE; 405 } 406 407 void rpn_exp10(calc_number_t *c) 408 { 409 double int_part; 410 411 modf(c->f, &int_part); 412 if (fmod(int_part, 2.) == 0.) 413 calc.is_nan = TRUE; 414 else { 415 c->f = pow(10., c->f); 416 if (_finite(c->f) == 0) 417 calc.is_nan = TRUE; 418 } 419 } 420 421 void rpn_ln(calc_number_t *c) 422 { 423 if (c->f <= 0) 424 calc.is_nan = TRUE; 425 else 426 c->f = log(c->f); 427 } 428 429 void rpn_log(calc_number_t *c) 430 { 431 if (c->f <= 0) 432 calc.is_nan = TRUE; 433 else 434 c->f = log10(c->f); 435 } 436 437 static double stat_sum(void) 438 { 439 double sum = 0; 440 statistic_t *p = calc.stat; 441 442 while (p != NULL) { 443 if (p->base == IDC_RADIO_DEC) 444 sum += p->num.f; 445 else 446 sum += p->num.i; 447 p = (statistic_t *)(p->next); 448 } 449 return sum; 450 } 451 452 static double stat_sum2(void) 453 { 454 double sum = 0; 455 statistic_t *p = calc.stat; 456 457 while (p != NULL) { 458 if (p->base == IDC_RADIO_DEC) 459 sum += p->num.f * p->num.f; 460 else 461 sum += (double)p->num.i * (double)p->num.i; 462 p = (statistic_t *)(p->next); 463 } 464 return sum; 465 } 466 467 void rpn_ave(calc_number_t *c) 468 { 469 double ave = 0; 470 int n; 471 472 ave = stat_sum(); 473 n = SendDlgItemMessage(calc.hStatWnd, IDC_LIST_STAT, LB_GETCOUNT, 0, 0); 474 475 if (n) 476 ave = ave / (double)n; 477 if (calc.base == IDC_RADIO_DEC) 478 c->f = ave; 479 else 480 c->i = (__int64)ave; 481 } 482 483 void rpn_ave2(calc_number_t *c) 484 { 485 double ave = 0; 486 int n; 487 488 ave = stat_sum2(); 489 n = SendDlgItemMessage(calc.hStatWnd, IDC_LIST_STAT, LB_GETCOUNT, 0, 0); 490 491 if (n) 492 ave = ave / (double)n; 493 if (calc.base == IDC_RADIO_DEC) 494 c->f = ave; 495 else 496 c->i = (__int64)ave; 497 } 498 499 void rpn_sum(calc_number_t *c) 500 { 501 double sum = stat_sum(); 502 503 if (calc.base == IDC_RADIO_DEC) 504 c->f = sum; 505 else 506 c->i = (__int64)sum; 507 } 508 509 void rpn_sum2(calc_number_t *c) 510 { 511 double sum = stat_sum2(); 512 513 if (calc.base == IDC_RADIO_DEC) 514 c->f = sum; 515 else 516 c->i = (__int64)sum; 517 } 518 519 static void rpn_s_ex(calc_number_t *c, int pop_type) 520 { 521 double ave = 0; 522 double n = 0; 523 double dev = 0; 524 double num = 0; 525 statistic_t *p = calc.stat; 526 527 ave = stat_sum(); 528 n = (double)SendDlgItemMessage(calc.hStatWnd, IDC_LIST_STAT, LB_GETCOUNT, 0, 0); 529 530 if (n == 0) { 531 c->f = 0; 532 return; 533 } 534 ave = ave / n; 535 536 dev = 0; 537 p = calc.stat; 538 while (p != NULL) { 539 if (p->base == IDC_RADIO_DEC) 540 num = p->num.f; 541 else 542 num = (double)p->num.i; 543 dev += pow(num-ave, 2.); 544 p = (statistic_t *)(p->next); 545 } 546 dev = sqrt(dev/(pop_type ? n-1 : n)); 547 if (calc.base == IDC_RADIO_DEC) 548 c->f = dev; 549 else 550 c->i = (__int64)dev; 551 } 552 553 void rpn_s(calc_number_t *c) 554 { 555 rpn_s_ex(c, 0); 556 } 557 558 void rpn_s_m1(calc_number_t *c) 559 { 560 rpn_s_ex(c, 1); 561 } 562 563 void rpn_dms2dec(calc_number_t *c) 564 { 565 double d, m, s; 566 567 m = modf(c->f, &d) * 100; 568 s = (modf(m, &m) * 100)+.5; 569 modf(s, &s); 570 571 m = m/60; 572 s = s/3600; 573 574 c->f = d + m + s; 575 } 576 577 void rpn_dec2dms(calc_number_t *c) 578 { 579 double d, m, s; 580 581 m = modf(c->f, &d) * 60; 582 s = ceil(modf(m, &m) * 60); 583 c->f = d + m/100. + s/10000.; 584 } 585 586 void rpn_zero(calc_number_t *c) 587 { 588 c->f = 0; 589 } 590 591 void rpn_copy(calc_number_t *dst, calc_number_t *src) 592 { 593 *dst = *src; 594 } 595 596 int rpn_is_zero(calc_number_t *c) 597 { 598 return (c->f == 0); 599 } 600 601 void rpn_alloc(calc_number_t *c) 602 { 603 } 604 605 void rpn_free(calc_number_t *c) 606 { 607 } 608