[CALC] Improve multi-precision support, and powers/roots. CORE-8486
[reactos.git] / base / applications / calc / fun_ieee.c
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 }