1 /* floathmath.c: higher mathematical functions, based on floatnum. */
2 /*
3 Copyright (C) 2007, 2008 Wolf Lammen.
4
5 This program is free software; you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation; either version 2 of the License , or
8 (at your option) any later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License
16 along with this program; see the file COPYING. If not, write to:
17
18 The Free Software Foundation, Inc.
19 59 Temple Place, Suite 330
20 Boston, MA 02111-1307 USA.
21
22
23 You may contact the author by:
24 e-mail: ookami1 <at> gmx <dot> de
25 mail: Wolf Lammen
26 Oertzweg 45
27 22307 Hamburg
28 Germany
29
30 *************************************************************************/
31
32 #include "floathmath.h"
33 #include "floatconst.h"
34 #include "floatcommon.h"
35 #include "floatlog.h"
36 #include "floatexp.h"
37 #include "floattrig.h"
38 #include "floatpower.h"
39 #include "floatipower.h"
40 #include "floatgamma.h"
41 #include "floaterf.h"
42 #include "floatlogic.h"
43
44 static char
_cvtlogic(t_longint * lx,cfloatnum x)45 _cvtlogic(
46 t_longint* lx,
47 cfloatnum x)
48 {
49 if (float_isnan(x))
50 {
51 float_seterror(NoOperand);
52 return 0;
53 }
54 if (_floatnum2logic(lx, x))
55 return 1;
56 float_seterror(OutOfLogicRange);
57 return 0;
58 }
59
60 char
float_lnxplus1(floatnum x,int digits)61 float_lnxplus1(
62 floatnum x,
63 int digits)
64 {
65 if (!chckmathparam(x, digits))
66 return 0;
67 if (float_getsign(x) < 0 && float_getexponent(x) >= 0)
68 return _seterror(x, OutOfDomain);
69 _lnxplus1(x, digits);
70 return 1;
71 }
72
73 char
float_ln(floatnum x,int digits)74 float_ln(floatnum x, int digits)
75 {
76 if (!chckmathparam(x, digits))
77 return 0;
78 if (float_getsign(x) <= 0)
79 return _seterror(x, OutOfDomain);
80 _ln(x, digits);
81 return 1;
82 }
83
84 char
float_artanh(floatnum x,int digits)85 float_artanh(
86 floatnum x,
87 int digits)
88 {
89 if (!chckmathparam(x, digits))
90 return 0;
91 if (float_getexponent(x) >= 0)
92 return _seterror(x, OutOfDomain);
93 _artanh(x, digits);
94 return 1;
95 }
96
97 char
float_artanhxplus1(floatnum x,int digits)98 float_artanhxplus1(
99 floatnum x,
100 int digits)
101 {
102 if (!chckmathparam(x, digits))
103 return 0;
104 if (float_getsign(x) >= 0 || float_abscmp(x, &c2) >= 0)
105 return _seterror(x, OutOfDomain);
106 if (float_cmp(x, &c1Div2) < 0)
107 {
108 float_neg(x);
109 _artanh1minusx(x, digits);
110 }
111 else
112 {
113 float_sub(x, &c1, x, digits+1);
114 _artanh(x, digits);
115 }
116 return 1;
117 }
118
119 char
float_arsinh(floatnum x,int digits)120 float_arsinh(
121 floatnum x,
122 int digits)
123 {
124 if (!chckmathparam(x, digits))
125 return 0;
126 _arsinh(x, digits);
127 return 1;
128 }
129
130 char
float_arcoshxplus1(floatnum x,int digits)131 float_arcoshxplus1(
132 floatnum x,
133 int digits)
134 {
135 if (!chckmathparam(x, digits))
136 return 0;
137 if (float_getsign(x) < 0)
138 return _seterror(x, OutOfDomain);
139 _arcoshxplus1(x, digits);
140 return 1;
141 }
142
143 char
float_arcosh(floatnum x,int digits)144 float_arcosh(
145 floatnum x,
146 int digits)
147 {
148 if (!chckmathparam(x, digits))
149 return 0;
150 float_sub(x, x, &c1, digits+1);
151 return float_arcoshxplus1(x, digits);
152 }
153
154 char
float_lg(floatnum x,int digits)155 float_lg(
156 floatnum x,
157 int digits)
158 {
159 floatstruct tmp;
160 int expx;
161
162 if (!chckmathparam(x, digits))
163 return 0;
164 if (float_getsign(x) <= 0)
165 return _seterror(x, OutOfDomain);
166 float_create(&tmp);
167 expx = float_getexponent(x);
168 float_setexponent(x, 0);
169 _ln(x, digits);
170 float_div(x, x, &cLn10, digits);
171 float_setinteger(&tmp, expx);
172 float_add(x, x, &tmp, digits);
173 float_free(&tmp);
174 return 1;
175 }
176
177 char
float_lb(floatnum x,int digits)178 float_lb(
179 floatnum x,
180 int digits)
181 {
182 if (!chckmathparam(x, digits))
183 return 0;
184 if (float_getsign(x) <= 0)
185 return _seterror(x, OutOfDomain);
186 _ln(x, digits);
187 float_div(x, x, &cLn2, digits);
188 return 1;
189 }
190
191 char
float_exp(floatnum x,int digits)192 float_exp(
193 floatnum x,
194 int digits)
195 {
196 signed char sgn;
197
198 if (!chckmathparam(x, digits))
199 return 0;
200 sgn = float_getsign(x);
201 if (_exp(x, digits))
202 return 1;
203 if (sgn < 0)
204 return _seterror(x, Underflow);
205 return _seterror(x, Overflow);
206 }
207
208 char
float_expminus1(floatnum x,int digits)209 float_expminus1(
210 floatnum x,
211 int digits)
212 {
213 if (!chckmathparam(x, digits))
214 return 0;
215 if (_expminus1(x, digits))
216 return 1;
217 return _seterror(x, Overflow);
218 }
219
220 char
float_cosh(floatnum x,int digits)221 float_cosh(
222 floatnum x,
223 int digits)
224 {
225 int expx;
226
227 if (!chckmathparam(x, digits))
228 return 0;
229 expx = float_getexponent(x);
230 if (2*expx+2 <= -digits || !_coshminus1(x, digits+2*expx))
231 {
232 if (expx > 0)
233 return _seterror(x, Overflow);
234 float_setzero(x);
235 }
236 return float_add(x, x, &c1, digits);
237 }
238
239 char
float_coshminus1(floatnum x,int digits)240 float_coshminus1(
241 floatnum x,
242 int digits)
243 {
244 int expx;
245
246 if (!chckmathparam(x, digits))
247 return 0;
248 expx = float_getexponent(x);
249 if (_coshminus1(x, digits))
250 return 1;
251 if (expx < 0)
252 return _seterror(x, Underflow);
253 return _seterror(x, Overflow);
254 }
255
256 char
float_sinh(floatnum x,int digits)257 float_sinh(
258 floatnum x,
259 int digits)
260 {
261 if (!chckmathparam(x, digits))
262 return 0;
263 if (_sinh(x, digits))
264 return 1;
265 return _seterror(x, Overflow);
266 }
267
268 char
float_tanh(floatnum x,int digits)269 float_tanh(
270 floatnum x,
271 int digits)
272 {
273 signed char sgn;
274
275 if (!chckmathparam(x, digits))
276 return 0;
277 sgn = float_getsign(x);
278 float_abs(x);
279 if (float_cmp(x, &c1Div2) >= 0)
280 _tanhgt0_5(x, digits);
281 else
282 _tanhlt0_5(x, digits);
283 float_setsign(x, sgn);
284 return 1;
285 }
286
287 char
float_tanhminus1(floatnum x,int digits)288 float_tanhminus1(
289 floatnum x,
290 int digits)
291 {
292 if (!chckmathparam(x, digits))
293 return 0;
294 if (float_cmp(x, &c1Div2) >= 0)
295 return _tanhminus1gt0(x, digits)? 1 : _seterror(x, Underflow);
296 if (!float_iszero(x))
297 {
298 if (float_abscmp(x, &c1Div2) <= 0)
299 _tanhlt0_5(x, digits);
300 else
301 {
302 float_setsign(x, 1);
303 _tanhgt0_5(x, digits);
304 float_setsign(x, -1);
305 }
306 }
307 return float_sub(x, x, &c1, digits);
308 }
309
310 char
float_arctan(floatnum x,int digits)311 float_arctan(
312 floatnum x,
313 int digits)
314 {
315 if (!chckmathparam(x, digits))
316 return 0;
317 _arctan(x, digits);
318 return 1;
319 }
320
321 char
float_arcsin(floatnum x,int digits)322 float_arcsin(
323 floatnum x,
324 int digits)
325 {
326 if (!chckmathparam(x, digits))
327 return 0;
328 if (float_abscmp(x, &c1) > 0)
329 return _seterror(x, OutOfDomain);
330 _arcsin(x, digits);
331 return 1;
332 }
333
334 char
float_arccos(floatnum x,int digits)335 float_arccos(
336 floatnum x,
337 int digits)
338 {
339 if (!chckmathparam(x, digits))
340 return 0;
341 if (float_abscmp(x, &c1) > 0)
342 return _seterror(x, OutOfDomain);
343 _arccos(x, digits);
344 return 1;
345 }
346
347 char
float_arccosxplus1(floatnum x,int digits)348 float_arccosxplus1(
349 floatnum x,
350 int digits)
351 {
352 if (!chckmathparam(x, digits))
353 return 0;
354 if (float_getsign(x) > 0 || float_abscmp(x, &c2) > 0)
355 return _seterror(x, OutOfDomain);
356 _arccosxplus1(x, digits);
357 return 1;
358 }
359
360 char
float_cos(floatnum x,int digits)361 float_cos(
362 floatnum x,
363 int digits)
364 {
365 if (!chckmathparam(x, digits))
366 return 0;
367 if (float_getexponent(x) >= DECPRECISION - 1 || !_trigreduce(x, digits))
368 return _seterror(x, EvalUnstable);
369 _cos(x, digits);
370 return 1;
371 }
372
373 char
float_cosminus1(floatnum x,int digits)374 float_cosminus1(
375 floatnum x,
376 int digits)
377 {
378 if (!chckmathparam(x, digits))
379 return 0;
380 if (!_trigreduce(x, digits))
381 return _seterror(x, EvalUnstable);
382 return _cosminus1(x, digits)? 1 : _seterror(x, Underflow);
383 }
384
385 char
float_sin(floatnum x,int digits)386 float_sin(
387 floatnum x,
388 int digits)
389 {
390 if (!chckmathparam(x, digits))
391 return 0;
392 if (float_getexponent(x) >= DECPRECISION - 1 || !_trigreduce(x, digits))
393 return _seterror(x, EvalUnstable);
394 _sin(x, digits);
395 return 1;
396 }
397
398 char
float_tan(floatnum x,int digits)399 float_tan(
400 floatnum x,
401 int digits)
402 {
403 if (!chckmathparam(x, digits))
404 return 0;
405 return float_getexponent(x) >= DECPRECISION - 1
406 || !_trigreduce(x, digits)
407 || !_tan(x, digits)? _seterror(x, EvalUnstable) : 1;
408 }
409
410 char
float_raisei(floatnum power,cfloatnum base,int exponent,int digits)411 float_raisei(
412 floatnum power,
413 cfloatnum base,
414 int exponent,
415 int digits)
416 {
417 if (digits <= 0 || digits > maxdigits)
418 return _seterror(power, InvalidPrecision);
419 if (float_isnan(base))
420 return _seterror(power, NoOperand);
421 if (float_iszero(base))
422 {
423 if (exponent == 0)
424 return _seterror(power, OutOfDomain);
425 if (exponent < 0)
426 return _seterror(power, ZeroDivide);
427 return _setzero(power);
428 }
429 digits += 14;
430 if (digits > maxdigits)
431 digits = maxdigits;
432 float_copy(power, base, digits);
433 if (!_raisei(power, exponent, digits)
434 || !float_isvalidexp(float_getexponent(power)))
435 {
436 if (float_getexponent(base) < 0)
437 return _seterror(power, Underflow);
438 return _seterror(power, Overflow);
439 }
440 return 1;
441 }
442
443 char
float_raise(floatnum power,cfloatnum base,cfloatnum exponent,int digits)444 float_raise(
445 floatnum power,
446 cfloatnum base,
447 cfloatnum exponent,
448 int digits)
449 {
450 signed char sgn;
451
452 if (float_isnan(exponent) || float_isnan(base))
453 return _seterror(power, NoOperand);
454 if (digits <= 0 || digits > MATHPRECISION)
455 return _seterror(power, InvalidPrecision);
456 if (float_iszero(base))
457 {
458 switch(float_getsign(exponent))
459 {
460 case 0:
461 return _seterror(power, OutOfDomain);
462 case -1:
463 return _seterror(power, ZeroDivide);
464 }
465 return _setzero(power);
466 }
467 sgn = float_getsign(base);
468 if (sgn < 0)
469 {
470 if (!float_isinteger(exponent))
471 return _seterror(power, OutOfDomain);
472 if ((float_getdigit(exponent, float_getexponent(exponent)) & 1) == 0)
473 sgn = 1;
474 }
475 float_copy(power, base, digits+1);
476 float_abs(power);
477 if (!_raise(power, exponent, digits))
478 {
479 float_seterror(Overflow);
480 if (float_getexponent(base) * float_getsign(exponent) < 0)
481 float_seterror(Underflow);
482 return _setnan(power);
483 }
484 float_setsign(power, sgn);
485 return 1;
486 }
487
488 char
float_power10(floatnum x,int digits)489 float_power10(
490 floatnum x,
491 int digits)
492 {
493 signed char sign;
494
495 if (!chckmathparam(x, digits))
496 return 0;
497 sign = float_getsign(x);
498 if (_power10(x, digits))
499 return 1;
500 return sign > 0? _seterror(x, Overflow) : _seterror(x, Underflow);
501 }
502
503 char
float_gamma(floatnum x,int digits)504 float_gamma(
505 floatnum x,
506 int digits)
507 {
508 signed char sign;
509 char result;
510
511 if (!chckmathparam(x, digits))
512 return 0;
513 sign = float_getsign(x);
514 if (float_isinteger(x))
515 {
516 if (sign <= 0)
517 return _seterror(x, ZeroDivide);
518 result = _gammaint(x, digits);
519 }
520 else if (float_getlength(x) - float_getexponent(x) == 2
521 && float_getdigit(x, float_getlength(x) - 1) == 5)
522 result = _gamma0_5(x, digits);
523 else
524 result = _gamma(x, digits);
525 if (!result)
526 {
527 if (sign < 0)
528 float_seterror(Underflow);
529 else
530 float_seterror(Overflow);
531 float_setnan(x);
532 }
533 return result;
534 }
535
536 char
float_lngamma(floatnum x,int digits)537 float_lngamma(
538 floatnum x,
539 int digits)
540 {
541 if (!x)
542 return _seterror(x, OutOfDomain);
543 return chckmathparam(x, digits) && _lngamma(x, digits)?
544 1 : _setnan(x);
545 }
546
547 char
float_factorial(floatnum x,int digits)548 float_factorial(
549 floatnum x,
550 int digits)
551 {
552 if (!float_isnan(x))
553 float_add(x, x, &c1, digits);
554 return float_gamma(x, digits);
555 }
556
557 char
float_pochhammer(floatnum x,cfloatnum delta,int digits)558 float_pochhammer(
559 floatnum x,
560 cfloatnum delta,
561 int digits)
562 {
563 if (!chckmathparam(x, digits))
564 return 0;
565 return float_isnan(delta)?
566 _seterror(x, NoOperand)
567 : _pochhammer(x, delta, digits);
568 }
569
570 char
float_erf(floatnum x,int digits)571 float_erf(floatnum x, int digits)
572 {
573 return chckmathparam(x, digits)? _erf(x, digits) : 0;
574 }
575
576 char
float_erfc(floatnum x,int digits)577 float_erfc(floatnum x, int digits)
578 {
579 return chckmathparam(x, digits)? _erfc(x, digits) : 0;
580 }
581
582 char
float_not(floatnum x)583 float_not(
584 floatnum x)
585 {
586 t_longint lx;
587
588 if(!_cvtlogic(&lx, x))
589 return _setnan(x);
590 _not(&lx);
591 _logic2floatnum(x, &lx);
592 return 1;
593 }
594
595 char
float_and(floatnum dest,cfloatnum x,cfloatnum y)596 float_and(
597 floatnum dest,
598 cfloatnum x,
599 cfloatnum y)
600 {
601 t_longint lx, ly;
602
603 if(!_cvtlogic(&lx, x) || !_cvtlogic(&ly, y))
604 return _setnan(dest);
605 _and(&lx, &ly);
606 _logic2floatnum(dest, &lx);
607 return 1;
608 }
609
610 char
float_or(floatnum dest,cfloatnum x,cfloatnum y)611 float_or(
612 floatnum dest,
613 cfloatnum x,
614 cfloatnum y)
615 {
616 t_longint lx, ly;
617
618 if(!_cvtlogic(&lx, x) || !_cvtlogic(&ly, y))
619 return _setnan(dest);
620 _or(&lx, &ly);
621 _logic2floatnum(dest, &lx);
622 return 1;
623 }
624
625 char
float_xor(floatnum dest,cfloatnum x,cfloatnum y)626 float_xor(
627 floatnum dest,
628 cfloatnum x,
629 cfloatnum y)
630 {
631 t_longint lx, ly;
632
633 if(!_cvtlogic(&lx, x) || !_cvtlogic(&ly, y))
634 return _setnan(dest);
635 _xor(&lx, &ly);
636 _logic2floatnum(dest, &lx);
637 return 1;
638 }
639
640 char
_doshift(floatnum dest,cfloatnum x,cfloatnum shift,char right)641 _doshift(
642 floatnum dest,
643 cfloatnum x,
644 cfloatnum shift,
645 char right)
646 {
647 int ishift;
648 t_longint lx;
649
650 if (float_isnan(shift))
651 return _seterror(dest, NoOperand);
652 if (!float_isinteger(shift))
653 return _seterror(dest, OutOfDomain);
654 if(!_cvtlogic(&lx, x))
655 return 0;
656 if (float_iszero(shift))
657 {
658 float_copy(dest, x, EXACT);
659 return 1;
660 }
661 ishift = float_asinteger(shift);
662 if (ishift == 0)
663 ishift = (3*LOGICRANGE) * float_getsign(shift);
664 if (!right)
665 ishift = -ishift;
666 if (ishift > 0)
667 _shr(&lx, ishift);
668 else
669 _shl(&lx, -ishift);
670 _logic2floatnum(dest, &lx);
671 return 1;
672 }
673
674 char
float_shr(floatnum dest,cfloatnum x,cfloatnum shift)675 float_shr(
676 floatnum dest,
677 cfloatnum x,
678 cfloatnum shift)
679 {
680 return _doshift(dest, x, shift, 1);
681 }
682
683 char
float_shl(floatnum dest,cfloatnum x,cfloatnum shift)684 float_shl(
685 floatnum dest,
686 cfloatnum x,
687 cfloatnum shift)
688 {
689 return _doshift(dest, x, shift, 0);
690 }
691