[CALC] Improve multi-precision support, and powers/roots. CORE-8486
[reactos.git] / base / applications / calc / fun_mpfr.c
1 /*
2 * ReactOS Calc (Math functions, GMP/MPFR 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 void validate_rad2angle(calc_number_t *c);
24 static void validate_angle2rad(calc_number_t *c);
25
26 void apply_int_mask(calc_number_t *r)
27 {
28 mpz_t a, mask;
29
30 switch (calc.size) {
31 case IDC_RADIO_QWORD:
32 mpz_init_set_str(mask, "FFFFFFFFFFFFFFFF", 16);
33 break;
34 case IDC_RADIO_DWORD:
35 mpz_init_set_str(mask, "00000000FFFFFFFF", 16);
36 break;
37 case IDC_RADIO_WORD:
38 mpz_init_set_str(mask, "000000000000FFFF", 16);
39 break;
40 case IDC_RADIO_BYTE:
41 mpz_init_set_str(mask, "00000000000000FF", 16);
42 break;
43 default:
44 mpz_init_set_si(mask, -1);
45 }
46 mpz_init(a);
47 mpfr_get_z(a, r->mf, MPFR_DEFAULT_RND);
48 mpz_and(a, a, mask);
49 mpfr_set_z(r->mf, a, MPFR_DEFAULT_RND);
50 mpz_clear(a);
51 mpz_clear(mask);
52 }
53
54 static void validate_rad2angle(calc_number_t *r)
55 {
56 mpfr_t mult, divs;
57
58 mpfr_init(mult);
59 mpfr_init(divs);
60 switch (calc.degr) {
61 case IDC_RADIO_DEG:
62 mpfr_set_ui(mult, 180, MPFR_DEFAULT_RND);
63 mpfr_const_pi(divs, MPFR_DEFAULT_RND);
64 break;
65 case IDC_RADIO_RAD:
66 mpfr_set_ui(mult, 1, MPFR_DEFAULT_RND);
67 mpfr_set_ui(divs, 1, MPFR_DEFAULT_RND);
68 break;
69 case IDC_RADIO_GRAD:
70 mpfr_set_ui(mult, 200, MPFR_DEFAULT_RND);
71 mpfr_const_pi(divs, MPFR_DEFAULT_RND);
72 break;
73 }
74 mpfr_mul(r->mf, r->mf, mult, MPFR_DEFAULT_RND);
75 mpfr_div(r->mf, r->mf, divs, MPFR_DEFAULT_RND);
76
77 mpfr_clear(mult);
78 mpfr_clear(divs);
79 }
80
81 static void validate_angle2rad(calc_number_t *r)
82 {
83 mpfr_t mult, divs;
84
85 if (!mpfr_number_p(r->mf)) {
86 calc.is_nan = TRUE;
87 return;
88 }
89 mpfr_init(mult);
90 mpfr_init(divs);
91 switch (calc.degr) {
92 case IDC_RADIO_DEG:
93 mpfr_const_pi(mult, MPFR_DEFAULT_RND);
94 mpfr_set_ui(divs, 180, MPFR_DEFAULT_RND);
95 break;
96 case IDC_RADIO_RAD:
97 mpfr_set_ui(mult, 1, MPFR_DEFAULT_RND);
98 mpfr_set_ui(divs, 1, MPFR_DEFAULT_RND);
99 break;
100 case IDC_RADIO_GRAD:
101 mpfr_const_pi(mult, MPFR_DEFAULT_RND);
102 mpfr_set_ui(divs, 200, MPFR_DEFAULT_RND);
103 break;
104 }
105 mpfr_mul(r->mf, r->mf, mult, MPFR_DEFAULT_RND);
106 mpfr_div(r->mf, r->mf, divs, MPFR_DEFAULT_RND);
107
108 mpfr_clear(mult);
109 mpfr_clear(divs);
110 }
111
112 static void build_rad_const(
113 mpfr_t *mp_pi,
114 mpfr_t *mp_pi_2,
115 mpfr_t *mp_3_pi_2,
116 mpfr_t *mp_2_pi)
117 {
118 mpfr_init(*mp_pi);
119 mpfr_init(*mp_pi_2);
120 mpfr_init(*mp_3_pi_2);
121 mpfr_init(*mp_2_pi);
122 mpfr_const_pi(*mp_pi, MPFR_DEFAULT_RND);
123 mpfr_div_ui(*mp_pi_2, *mp_pi, 2, MPFR_DEFAULT_RND);
124 mpfr_mul_ui(*mp_3_pi_2, *mp_pi, 3, MPFR_DEFAULT_RND);
125 mpfr_div_ui(*mp_3_pi_2, *mp_3_pi_2, 2, MPFR_DEFAULT_RND);
126 mpfr_mul_ui(*mp_2_pi, *mp_pi, 2, MPFR_DEFAULT_RND);
127 }
128
129 void rpn_sin(calc_number_t *c)
130 {
131 mpfr_t mp_pi, mp_pi_2, mp_3_pi_2, mp_2_pi;
132
133 validate_angle2rad(c);
134 build_rad_const(&mp_pi, &mp_pi_2, &mp_3_pi_2, &mp_2_pi);
135
136 if (rpn_is_zero(c) || !mpfr_cmp(c->mf, mp_pi) || !mpfr_cmp(c->mf, mp_2_pi))
137 rpn_zero(c);
138 else
139 if (!mpfr_cmp(c->mf, mp_3_pi_2))
140 mpfr_set_si(c->mf, -1, MPFR_DEFAULT_RND);
141 else
142 if (!mpfr_cmp(c->mf, mp_pi_2))
143 mpfr_set_si(c->mf, 1, MPFR_DEFAULT_RND);
144 else {
145 mpfr_sin(c->mf, c->mf, MPFR_DEFAULT_RND);
146 if (!mpfr_number_p(c->mf)) calc.is_nan = TRUE;
147 }
148 mpfr_clear(mp_pi);
149 mpfr_clear(mp_pi_2);
150 mpfr_clear(mp_3_pi_2);
151 mpfr_clear(mp_2_pi);
152 }
153 void rpn_cos(calc_number_t *c)
154 {
155 mpfr_t mp_pi, mp_pi_2, mp_3_pi_2, mp_2_pi;
156
157 validate_angle2rad(c);
158 build_rad_const(&mp_pi, &mp_pi_2, &mp_3_pi_2, &mp_2_pi);
159
160 if (!mpfr_cmp(c->mf, mp_pi_2) || !mpfr_cmp(c->mf, mp_3_pi_2))
161 rpn_zero(c);
162 else
163 if (!mpfr_cmp(c->mf, mp_pi))
164 mpfr_set_si(c->mf, -1, MPFR_DEFAULT_RND);
165 else
166 if (!mpfr_cmp(c->mf, mp_2_pi))
167 mpfr_set_si(c->mf, 1, MPFR_DEFAULT_RND);
168 else {
169 mpfr_cos(c->mf, c->mf, MPFR_DEFAULT_RND);
170 if (!mpfr_number_p(c->mf)) calc.is_nan = TRUE;
171 }
172 mpfr_clear(mp_pi);
173 mpfr_clear(mp_pi_2);
174 mpfr_clear(mp_3_pi_2);
175 mpfr_clear(mp_2_pi);
176 }
177 void rpn_tan(calc_number_t *c)
178 {
179 mpfr_t mp_pi, mp_pi_2, mp_3_pi_2, mp_2_pi;
180
181 validate_angle2rad(c);
182 build_rad_const(&mp_pi, &mp_pi_2, &mp_3_pi_2, &mp_2_pi);
183
184 if (!mpfr_cmp(c->mf, mp_pi_2) || !mpfr_cmp(c->mf, mp_3_pi_2))
185 calc.is_nan = TRUE;
186 else
187 if (!mpfr_cmp(c->mf, mp_pi) || !mpfr_cmp(c->mf, mp_2_pi))
188 rpn_zero(c);
189 else {
190 mpfr_tan(c->mf, c->mf, MPFR_DEFAULT_RND);
191 if (!mpfr_number_p(c->mf)) calc.is_nan = TRUE;
192 }
193 mpfr_clear(mp_pi);
194 mpfr_clear(mp_pi_2);
195 mpfr_clear(mp_3_pi_2);
196 mpfr_clear(mp_2_pi);
197 }
198
199 void rpn_asin(calc_number_t *c)
200 {
201 mpfr_asin(c->mf, c->mf, MPFR_DEFAULT_RND);
202 validate_rad2angle(c);
203 }
204 void rpn_acos(calc_number_t *c)
205 {
206 mpfr_acos(c->mf, c->mf, MPFR_DEFAULT_RND);
207 validate_rad2angle(c);
208 }
209 void rpn_atan(calc_number_t *c)
210 {
211 mpfr_atan(c->mf, c->mf, MPFR_DEFAULT_RND);
212 validate_rad2angle(c);
213 }
214
215 void rpn_sinh(calc_number_t *c)
216 {
217 mpfr_sinh(c->mf, c->mf, MPFR_DEFAULT_RND);
218 if (!mpfr_number_p(c->mf)) calc.is_nan = TRUE;
219 }
220 void rpn_cosh(calc_number_t *c)
221 {
222 mpfr_cosh(c->mf, c->mf, MPFR_DEFAULT_RND);
223 if (!mpfr_number_p(c->mf)) calc.is_nan = TRUE;
224 }
225 void rpn_tanh(calc_number_t *c)
226 {
227 mpfr_tanh(c->mf, c->mf, MPFR_DEFAULT_RND);
228 if (!mpfr_number_p(c->mf)) calc.is_nan = TRUE;
229 }
230
231 void rpn_asinh(calc_number_t *c)
232 {
233 mpfr_asinh(c->mf, c->mf, MPFR_DEFAULT_RND);
234 if (!mpfr_number_p(c->mf)) calc.is_nan = TRUE;
235 }
236 void rpn_acosh(calc_number_t *c)
237 {
238 mpfr_acosh(c->mf, c->mf, MPFR_DEFAULT_RND);
239 if (!mpfr_number_p(c->mf)) calc.is_nan = TRUE;
240 }
241 void rpn_atanh(calc_number_t *c)
242 {
243 mpfr_atanh(c->mf, c->mf, MPFR_DEFAULT_RND);
244 if (!mpfr_number_p(c->mf)) calc.is_nan = TRUE;
245 }
246
247 void rpn_int(calc_number_t *c)
248 {
249 mpfr_trunc(c->mf, c->mf);
250 }
251
252 void rpn_frac(calc_number_t *c)
253 {
254 mpfr_frac(c->mf, c->mf, MPFR_DEFAULT_RND);
255 }
256
257 void rpn_reci(calc_number_t *c)
258 {
259 if (mpfr_sgn(c->mf) == 0)
260 calc.is_nan = TRUE;
261 else
262 mpfr_ui_div(c->mf, 1, c->mf, MPFR_DEFAULT_RND);
263 }
264
265 void rpn_fact(calc_number_t *c)
266 {
267 if (mpfr_sgn(c->mf) < 0) {
268 calc.is_nan = TRUE;
269 return;
270 }
271
272 mpfr_trunc(c->mf, c->mf);
273 if (mpfr_fits_ulong_p(c->mf, MPFR_DEFAULT_RND) == 0)
274 calc.is_nan = TRUE;
275 else {
276 mpfr_fac_ui(c->mf, mpfr_get_ui(c->mf, MPFR_DEFAULT_RND), MPFR_DEFAULT_RND);
277 if (!mpfr_number_p(c->mf)) calc.is_nan = TRUE;
278 }
279 }
280
281 void rpn_not(calc_number_t *c)
282 {
283 mpz_t a;
284
285 mpz_init(a);
286 mpfr_get_z(a, c->mf, MPFR_DEFAULT_RND);
287 mpz_com(a, a);
288 mpfr_set_z(c->mf, a, MPFR_DEFAULT_RND);
289 mpz_clear(a);
290 }
291
292 void rpn_pi(calc_number_t *c)
293 {
294 mpfr_const_pi(c->mf, MPFR_DEFAULT_RND);
295 }
296
297 void rpn_2pi(calc_number_t *c)
298 {
299 mpfr_const_pi(c->mf, MPFR_DEFAULT_RND);
300 mpfr_mul_ui(c->mf, c->mf, 2, MPFR_DEFAULT_RND);
301 }
302
303 void rpn_sign(calc_number_t *c)
304 {
305 mpfr_mul_si(c->mf, c->mf, -1, MPFR_DEFAULT_RND);
306 }
307
308 void rpn_exp2(calc_number_t *c)
309 {
310 mpfr_sqr(c->mf, c->mf, MPFR_DEFAULT_RND);
311 if (!mpfr_number_p(c->mf)) calc.is_nan = TRUE;
312 }
313
314 void rpn_exp3(calc_number_t *c)
315 {
316 mpfr_pow_ui(c->mf, c->mf, 3, MPFR_DEFAULT_RND);
317 if (!mpfr_number_p(c->mf)) calc.is_nan = TRUE;
318 }
319
320 void rpn_sqrt(calc_number_t *c)
321 {
322 mpfr_sqrt(c->mf, c->mf, MPFR_DEFAULT_RND);
323 if (!mpfr_number_p(c->mf)) calc.is_nan = TRUE;
324 }
325
326 void rpn_cbrt(calc_number_t *c)
327 {
328 mpfr_cbrt(c->mf, c->mf, MPFR_DEFAULT_RND);
329 if (!mpfr_number_p(c->mf)) calc.is_nan = TRUE;
330 }
331
332 void rpn_exp(calc_number_t *c)
333 {
334 mpfr_exp(c->mf, c->mf, MPFR_DEFAULT_RND);
335 if (!mpfr_number_p(c->mf)) calc.is_nan = TRUE;
336 }
337
338 void rpn_exp10(calc_number_t *c)
339 {
340 mpfr_exp10(c->mf, c->mf, MPFR_DEFAULT_RND);
341 if (!mpfr_number_p(c->mf)) calc.is_nan = TRUE;
342 }
343
344 void rpn_ln(calc_number_t *c)
345 {
346 mpfr_log(c->mf, c->mf, MPFR_DEFAULT_RND);
347 if (!mpfr_number_p(c->mf)) calc.is_nan = TRUE;
348 }
349
350 void rpn_log(calc_number_t *c)
351 {
352 mpfr_log10(c->mf, c->mf, MPFR_DEFAULT_RND);
353 if (!mpfr_number_p(c->mf)) calc.is_nan = TRUE;
354 }
355
356 static void stat_sum(mpfr_t sum)
357 {
358 statistic_t *p = calc.stat;
359
360 mpfr_set_ui(sum, 0, MPFR_DEFAULT_RND);
361 while (p != NULL) {
362 mpfr_add(sum, sum, p->num.mf, MPFR_DEFAULT_RND);
363 p = (statistic_t *)(p->next);
364 }
365 }
366
367 static void stat_sum2(mpfr_t sum)
368 {
369 statistic_t *p = calc.stat;
370 mpfr_t sqr;
371
372 mpfr_init(sqr);
373 mpfr_set_ui(sum, 0, MPFR_DEFAULT_RND);
374 while (p != NULL) {
375 mpfr_mul(sqr, p->num.mf, p->num.mf, MPFR_DEFAULT_RND);
376 mpfr_add(sum, sum, sqr, MPFR_DEFAULT_RND);
377 p = (statistic_t *)(p->next);
378 }
379 mpfr_clear(sqr);
380 }
381
382 void rpn_ave(calc_number_t *c)
383 {
384 int n;
385
386 stat_sum(c->mf);
387 n = SendDlgItemMessage(calc.hStatWnd, IDC_LIST_STAT, LB_GETCOUNT, 0, 0);
388
389 if (n)
390 mpfr_div_ui(c->mf, c->mf, n, MPFR_DEFAULT_RND);
391
392 if (calc.base != IDC_RADIO_DEC)
393 mpfr_trunc(c->mf, c->mf);
394 }
395
396 void rpn_ave2(calc_number_t *c)
397 {
398 int n;
399
400 stat_sum2(c->mf);
401 n = SendDlgItemMessage(calc.hStatWnd, IDC_LIST_STAT, LB_GETCOUNT, 0, 0);
402
403 if (n)
404 mpfr_div_ui(c->mf, c->mf, n, MPFR_DEFAULT_RND);
405
406 if (calc.base != IDC_RADIO_DEC)
407 mpfr_trunc(c->mf, c->mf);
408 }
409
410 void rpn_sum(calc_number_t *c)
411 {
412 stat_sum(c->mf);
413
414 if (calc.base != IDC_RADIO_DEC)
415 mpfr_trunc(c->mf, c->mf);
416 }
417
418 void rpn_sum2(calc_number_t *c)
419 {
420 stat_sum2(c->mf);
421
422 if (calc.base != IDC_RADIO_DEC)
423 mpfr_trunc(c->mf, c->mf);
424 }
425
426 static void rpn_s_ex(calc_number_t *c, int pop_type)
427 {
428 mpfr_t dev;
429 mpfr_t num;
430 unsigned long n = 0;
431 statistic_t *p = calc.stat;
432
433 n = SendDlgItemMessage(calc.hStatWnd, IDC_LIST_STAT, LB_GETCOUNT, 0, 0);
434 if (n < 2) {
435 mpfr_set_ui(c->mf, 0, MPFR_DEFAULT_RND);
436 return;
437 }
438
439 stat_sum(c->mf);
440 mpfr_div_ui(c->mf, c->mf, n, MPFR_DEFAULT_RND);
441
442 mpfr_init(dev);
443 mpfr_init(num);
444
445 mpfr_set_ui(dev, 0, MPFR_DEFAULT_RND);
446 p = calc.stat;
447 while (p != NULL) {
448 mpfr_sub(num, p->num.mf, c->mf, MPFR_DEFAULT_RND);
449 mpfr_sqr(num, num, MPFR_DEFAULT_RND);
450 mpfr_add(dev, dev, num, MPFR_DEFAULT_RND);
451 p = (statistic_t *)(p->next);
452 }
453 mpfr_div_ui(c->mf, dev, pop_type ? n-1 : n, MPFR_DEFAULT_RND);
454 mpfr_sqrt(c->mf, c->mf, MPFR_DEFAULT_RND);
455
456 if (calc.base != IDC_RADIO_DEC)
457 mpfr_trunc(c->mf, c->mf);
458
459 mpfr_clear(dev);
460 mpfr_clear(num);
461 }
462
463 void rpn_s(calc_number_t *c)
464 {
465 rpn_s_ex(c, 0);
466 }
467
468 void rpn_s_m1(calc_number_t *c)
469 {
470 rpn_s_ex(c, 1);
471 }
472
473 void rpn_dms2dec(calc_number_t *c)
474 {
475 mpfr_t d, m, s;
476
477 mpfr_init(d);
478 mpfr_init(m);
479 mpfr_init(s);
480
481 mpfr_trunc(d, c->mf);
482 mpfr_frac(m, c->mf, MPFR_DEFAULT_RND);
483 mpfr_mul_ui(m, m, 100, MPFR_DEFAULT_RND);
484
485 mpfr_frac(s, m, MPFR_DEFAULT_RND);
486 mpfr_trunc(m, m);
487 mpfr_mul_ui(s, s, 100, MPFR_DEFAULT_RND);
488 mpfr_ceil(s, s);
489
490 mpfr_div_ui(m, m, 60, MPFR_DEFAULT_RND);
491 mpfr_div_ui(s, s, 3600, MPFR_DEFAULT_RND);
492 mpfr_add(c->mf, d, m, MPFR_DEFAULT_RND);
493 mpfr_add(c->mf, c->mf, s, MPFR_DEFAULT_RND);
494
495 mpfr_clear(d);
496 mpfr_clear(m);
497 mpfr_clear(s);
498 }
499
500 void rpn_dec2dms(calc_number_t *c)
501 {
502 mpfr_t d, m, s;
503
504 mpfr_init(d);
505 mpfr_init(m);
506 mpfr_init(s);
507
508 mpfr_trunc(d, c->mf);
509 mpfr_frac(m, c->mf, MPFR_DEFAULT_RND);
510 mpfr_mul_ui(m, m, 60, MPFR_DEFAULT_RND);
511
512 mpfr_frac(s, m, MPFR_DEFAULT_RND);
513 mpfr_trunc(m, m);
514 mpfr_mul_ui(s, s, 60, MPFR_DEFAULT_RND);
515 mpfr_ceil(s, s);
516
517 mpfr_div_ui(m, m, 100, MPFR_DEFAULT_RND);
518 mpfr_div_ui(s, s, 10000, MPFR_DEFAULT_RND);
519 mpfr_add(c->mf, d, m, MPFR_DEFAULT_RND);
520 mpfr_add(c->mf, c->mf, s, MPFR_DEFAULT_RND);
521
522 mpfr_clear(d);
523 mpfr_clear(m);
524 mpfr_clear(s);
525 }
526
527 void rpn_zero(calc_number_t *c)
528 {
529 mpfr_set_ui(c->mf, 0, MPFR_DEFAULT_RND);
530 }
531
532 void rpn_copy(calc_number_t *dst, calc_number_t *src)
533 {
534 mpfr_set(dst->mf, src->mf, MPFR_DEFAULT_RND);
535 }
536
537 int rpn_is_zero(calc_number_t *c)
538 {
539 return (mpfr_sgn(c->mf) == 0);
540 }
541
542 void rpn_alloc(calc_number_t *c)
543 {
544 mpfr_init(c->mf);
545 }
546
547 void rpn_free(calc_number_t *c)
548 {
549 mpfr_clear(c->mf);
550 }