GEOS  3.8.0dev
ttmath.h
Go to the documentation of this file.
1 /*
2  * This file is a part of TTMath Bignum Library
3  * and is distributed under the 3-Clause BSD Licence.
4  * Author: Tomasz Sowa <t.sowa@ttmath.org>
5  */
6 
7 /*
8  * Copyright (c) 2006-2017, Tomasz Sowa
9  * All rights reserved.
10  *
11  * Redistribution and use in source and binary forms, with or without
12  * modification, are permitted provided that the following conditions are met:
13  *
14  * * Redistributions of source code must retain the above copyright notice,
15  * this list of conditions and the following disclaimer.
16  *
17  * * Redistributions in binary form must reproduce the above copyright
18  * notice, this list of conditions and the following disclaimer in the
19  * documentation and/or other materials provided with the distribution.
20  *
21  * * Neither the name Tomasz Sowa nor the names of contributors to this
22  * project may be used to endorse or promote products derived
23  * from this software without specific prior written permission.
24  *
25  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
26  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
28  * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
29  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
30  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
31  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
32  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
33  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
34  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
35  * THE POSSIBILITY OF SUCH DAMAGE.
36  */
37 
38 
39 
40 #ifndef headerfilettmathmathtt
41 #define headerfilettmathmathtt
42 
48 #ifdef _MSC_VER
49 //warning C4127: conditional expression is constant
50 #pragma warning( disable: 4127 )
51 //warning C4702: unreachable code
52 #pragma warning( disable: 4702 )
53 //warning C4800: forcing value to bool 'true' or 'false' (performance warning)
54 #pragma warning( disable: 4800 )
55 #endif
56 
57 
58 #include "ttmathbig.h"
59 #include "ttmathobjects.h"
60 
61 
62 namespace ttmath
63 {
64  /*
65  *
66  * functions defined here are used only with Big<> types
67  *
68  *
69  */
70 
71 
72  /*
73  *
74  * functions for rounding
75  *
76  *
77  */
78 
79 
89  template<class ValueType>
90  ValueType SkipFraction(const ValueType & x)
91  {
92  ValueType result( x );
93  result.SkipFraction();
94 
95  return result;
96  }
97 
98 
108  template<class ValueType>
109  ValueType Round(const ValueType & x, ErrorCode * err = 0)
110  {
111  if( x.IsNan() )
112  {
113  if( err )
114  *err = err_improper_argument;
115 
116  return x; // NaN
117  }
118 
119  ValueType result( x );
120  uint c = result.Round();
121 
122  if( err )
123  *err = c ? err_overflow : err_ok;
124 
125  return result;
126  }
127 
128 
129 
141  template<class ValueType>
142  ValueType Ceil(const ValueType & x, ErrorCode * err = 0)
143  {
144  if( x.IsNan() )
145  {
146  if( err )
147  *err = err_improper_argument;
148 
149  return x; // NaN
150  }
151 
152  ValueType result(x);
153  uint c = 0;
154 
155  result.SkipFraction();
156 
157  if( result != x )
158  {
159  // x is with fraction
160  // if x is negative we don't have to do anything
161  if( !x.IsSign() )
162  {
163  ValueType one;
164  one.SetOne();
165 
166  c += result.Add(one);
167  }
168  }
169 
170  if( err )
171  *err = c ? err_overflow : err_ok;
172 
173  return result;
174  }
175 
176 
188  template<class ValueType>
189  ValueType Floor(const ValueType & x, ErrorCode * err = 0)
190  {
191  if( x.IsNan() )
192  {
193  if( err )
194  *err = err_improper_argument;
195 
196  return x; // NaN
197  }
198 
199  ValueType result(x);
200  uint c = 0;
201 
202  result.SkipFraction();
203 
204  if( result != x )
205  {
206  // x is with fraction
207  // if x is positive we don't have to do anything
208  if( x.IsSign() )
209  {
210  ValueType one;
211  one.SetOne();
212 
213  c += result.Sub(one);
214  }
215  }
216 
217  if( err )
218  *err = c ? err_overflow : err_ok;
219 
220  return result;
221  }
222 
223 
224 
225  /*
226  *
227  * logarithms and the exponent
228  *
229  *
230  */
231 
232 
236  template<class ValueType>
237  ValueType Ln(const ValueType & x, ErrorCode * err = 0)
238  {
239  if( x.IsNan() )
240  {
241  if( err )
242  *err = err_improper_argument;
243 
244  return x; // NaN
245  }
246 
247  ValueType result;
248  uint state = result.Ln(x);
249 
250  if( err )
251  {
252  switch( state )
253  {
254  case 0:
255  *err = err_ok;
256  break;
257  case 1:
258  *err = err_overflow;
259  break;
260  case 2:
261  *err = err_improper_argument;
262  break;
263  default:
264  *err = err_internal_error;
265  break;
266  }
267  }
268 
269 
270  return result;
271  }
272 
273 
277  template<class ValueType>
278  ValueType Log(const ValueType & x, const ValueType & base, ErrorCode * err = 0)
279  {
280  if( x.IsNan() )
281  {
282  if( err ) *err = err_improper_argument;
283  return x;
284  }
285 
286  if( base.IsNan() )
287  {
288  if( err ) *err = err_improper_argument;
289  return base;
290  }
291 
292  ValueType result;
293  uint state = result.Log(x, base);
294 
295  if( err )
296  {
297  switch( state )
298  {
299  case 0:
300  *err = err_ok;
301  break;
302  case 1:
303  *err = err_overflow;
304  break;
305  case 2:
306  case 3:
307  *err = err_improper_argument;
308  break;
309  default:
310  *err = err_internal_error;
311  break;
312  }
313  }
314 
315  return result;
316  }
317 
318 
322  template<class ValueType>
323  ValueType Exp(const ValueType & x, ErrorCode * err = 0)
324  {
325  if( x.IsNan() )
326  {
327  if( err )
328  *err = err_improper_argument;
329 
330  return x; // NaN
331  }
332 
333  ValueType result;
334  uint c = result.Exp(x);
335 
336  if( err )
337  *err = c ? err_overflow : err_ok;
338 
339  return result;
340  }
341 
342 
350  /*
351  this namespace consists of auxiliary functions
352  (something like 'private' in a class)
353 
354  this is excluded from doxygen documentation
355  (option EXCLUDE_SYMBOLS in doxygen.cfg)
356  */
357  namespace auxiliaryfunctions
358  {
359 
364  template<class ValueType>
365  uint PrepareSin(ValueType & x, bool & change_sign)
366  {
367  ValueType temp;
368 
369  change_sign = false;
370 
371  if( x.IsSign() )
372  {
373  // we're using the formula 'sin(-x) = -sin(x)'
374  change_sign = !change_sign;
375  x.ChangeSign();
376  }
377 
378  // we're reducing the period 2*PI
379  // (for big values there'll always be zero)
380  temp.Set2Pi();
381 
382  if( x.Mod(temp) )
383  return 1;
384 
385 
386  // we're setting 'x' as being in the range of <0, 0.5PI>
387 
388  temp.SetPi();
389 
390  if( x > temp )
391  {
392  // x is in (pi, 2*pi>
393  x.Sub( temp );
394  change_sign = !change_sign;
395  }
396 
397  temp.Set05Pi();
398 
399  if( x > temp )
400  {
401  // x is in (0.5pi, pi>
402  x.Sub( temp );
403  x = temp - x;
404  }
405 
406  return 0;
407  }
408 
409 
429  template<class ValueType>
430  ValueType Sin0pi05(const ValueType & x)
431  {
432  ValueType result;
433  ValueType numerator, denominator;
434  ValueType d_numerator, d_denominator;
435  ValueType one, temp, old_result;
436 
437  // temp = pi/4
438  temp.Set05Pi();
439  temp.exponent.SubOne();
440 
441  one.SetOne();
442 
443  if( x < temp )
444  {
445  // we're using the Taylor series with a=0
446  result = x;
447  numerator = x;
448  denominator = one;
449 
450  // d_numerator = x^2
451  d_numerator = x;
452  d_numerator.Mul(x);
453 
454  d_denominator = 2;
455  }
456  else
457  {
458  // we're using the Taylor series with a=PI/2
459  result = one;
460  numerator = one;
461  denominator = one;
462 
463  // d_numerator = (x-pi/2)^2
464  ValueType pi05;
465  pi05.Set05Pi();
466 
467  temp = x;
468  temp.Sub( pi05 );
469  d_numerator = temp;
470  d_numerator.Mul( temp );
471 
472  d_denominator = one;
473  }
474 
475  uint c = 0;
476  bool addition = false;
477 
478  old_result = result;
479  for(uint i=1 ; i<=TTMATH_ARITHMETIC_MAX_LOOP ; ++i)
480  {
481  // we're starting from a second part of the formula
482  c += numerator. Mul( d_numerator );
483  c += denominator. Mul( d_denominator );
484  c += d_denominator.Add( one );
485  c += denominator. Mul( d_denominator );
486  c += d_denominator.Add( one );
487  temp = numerator;
488  c += temp.Div(denominator);
489 
490  if( c )
491  // Sin is from <-1,1> and cannot make an overflow
492  // but the carry can be from the Taylor series
493  // (then we only break our calculations)
494  break;
495 
496  if( addition )
497  result.Add( temp );
498  else
499  result.Sub( temp );
500 
501 
502  addition = !addition;
503 
504  // we're testing whether the result has changed after adding
505  // the next part of the Taylor formula, if not we end the loop
506  // (it means 'x' is zero or 'x' is PI/2 or this part of the formula
507  // is too small)
508  if( result == old_result )
509  break;
510 
511  old_result = result;
512  }
513 
514  return result;
515  }
516 
517  } // namespace auxiliaryfunctions
518 
519 
520 
524  template<class ValueType>
525  ValueType Sin(ValueType x, ErrorCode * err = 0)
526  {
527  using namespace auxiliaryfunctions;
528 
529  ValueType one, result;
530  bool change_sign;
531 
532  if( x.IsNan() )
533  {
534  if( err )
535  *err = err_improper_argument;
536 
537  return x;
538  }
539 
540  if( err )
541  *err = err_ok;
542 
543  if( PrepareSin( x, change_sign ) )
544  {
545  // x is too big, we cannnot reduce the 2*PI period
546  // prior to version 0.8.5 the result was zero
547 
548  // result has NaN flag set by default
549 
550  if( err )
551  *err = err_overflow; // maybe another error code? err_improper_argument?
552 
553  return result; // NaN is set by default
554  }
555 
556  result = Sin0pi05( x );
557 
558  one.SetOne();
559 
560  // after calculations there can be small distortions in the result
561  if( result > one )
562  result = one;
563  else
564  if( result.IsSign() )
565  // we've calculated the sin from <0, pi/2> and the result
566  // should be positive
567  result.SetZero();
568 
569  if( change_sign )
570  result.ChangeSign();
571 
572  return result;
573  }
574 
575 
580  template<class ValueType>
581  ValueType Cos(ValueType x, ErrorCode * err = 0)
582  {
583  if( x.IsNan() )
584  {
585  if( err )
586  *err = err_improper_argument;
587 
588  return x; // NaN
589  }
590 
591  ValueType pi05;
592  pi05.Set05Pi();
593 
594  uint c = x.Add( pi05 );
595 
596  if( c )
597  {
598  if( err )
599  *err = err_overflow;
600 
601  return ValueType(); // result is undefined (NaN is set by default)
602  }
603 
604  return Sin(x, err);
605  }
606 
607 
618  template<class ValueType>
619  ValueType Tan(const ValueType & x, ErrorCode * err = 0)
620  {
621  ValueType result = Cos(x, err);
622 
623  if( err && *err != err_ok )
624  return result;
625 
626  if( result.IsZero() )
627  {
628  if( err )
629  *err = err_improper_argument;
630 
631  result.SetNan();
632 
633  return result;
634  }
635 
636  return Sin(x, err) / result;
637  }
638 
639 
646  template<class ValueType>
647  ValueType Tg(const ValueType & x, ErrorCode * err = 0)
648  {
649  return Tan(x, err);
650  }
651 
652 
660  template<class ValueType>
661  ValueType Cot(const ValueType & x, ErrorCode * err = 0)
662  {
663  ValueType result = Sin(x, err);
664 
665  if( err && *err != err_ok )
666  return result;
667 
668  if( result.IsZero() )
669  {
670  if( err )
671  *err = err_improper_argument;
672 
673  result.SetNan();
674 
675  return result;
676  }
677 
678  return Cos(x, err) / result;
679  }
680 
681 
688  template<class ValueType>
689  ValueType Ctg(const ValueType & x, ErrorCode * err = 0)
690  {
691  return Cot(x, err);
692  }
693 
694 
695  /*
696  *
697  * inverse trigonometric functions
698  *
699  *
700  */
701 
702  namespace auxiliaryfunctions
703  {
704 
714  template<class ValueType>
715  ValueType ASin_0(const ValueType & x)
716  {
717  ValueType nominator, denominator, nominator_add, nominator_x, denominator_add, denominator_x;
718  ValueType two, result(x), x2(x);
719  ValueType nominator_temp, denominator_temp, old_result = result;
720  uint c = 0;
721 
722  x2.Mul(x);
723  two = 2;
724 
725  nominator.SetOne();
726  denominator = two;
727  nominator_add = nominator;
728  denominator_add = denominator;
729  nominator_x = x;
730  denominator_x = 3;
731 
732  for(uint i=1 ; i<=TTMATH_ARITHMETIC_MAX_LOOP ; ++i)
733  {
734  c += nominator_x.Mul(x2);
735  nominator_temp = nominator_x;
736  c += nominator_temp.Mul(nominator);
737  denominator_temp = denominator;
738  c += denominator_temp.Mul(denominator_x);
739  c += nominator_temp.Div(denominator_temp);
740 
741  // if there is a carry somewhere we only break the calculating
742  // the result should be ok -- it's from <-pi/2, pi/2>
743  if( c )
744  break;
745 
746  result.Add(nominator_temp);
747 
748  if( result == old_result )
749  // there's no sense to calculate more
750  break;
751 
752  old_result = result;
753 
754 
755  c += nominator_add.Add(two);
756  c += denominator_add.Add(two);
757  c += nominator.Mul(nominator_add);
758  c += denominator.Mul(denominator_add);
759  c += denominator_x.Add(two);
760  }
761 
762  return result;
763  }
764 
765 
766 
778  template<class ValueType>
779  ValueType ASin_1(const ValueType & x)
780  {
781  ValueType nominator, denominator, nominator_add, nominator_x, nominator_x_add, denominator_add, denominator_x;
782  ValueType denominator2;
783  ValueType one, two, result;
784  ValueType nominator_temp, denominator_temp, old_result;
785  uint c = 0;
786 
787  two = 2;
788 
789  one.SetOne();
790  nominator = one;
791  result = one;
792  old_result = result;
793  denominator = two;
794  nominator_add = nominator;
795  denominator_add = denominator;
796  nominator_x = one;
797  nominator_x.Sub(x);
798  nominator_x_add = nominator_x;
799  denominator_x = 3;
800  denominator2 = two;
801 
802 
803  for(uint i=1 ; i<=TTMATH_ARITHMETIC_MAX_LOOP ; ++i)
804  {
805  nominator_temp = nominator_x;
806  c += nominator_temp.Mul(nominator);
807  denominator_temp = denominator;
808  c += denominator_temp.Mul(denominator_x);
809  c += denominator_temp.Mul(denominator2);
810  c += nominator_temp.Div(denominator_temp);
811 
812  // if there is a carry somewhere we only break the calculating
813  // the result should be ok -- it's from <-pi/2, pi/2>
814  if( c )
815  break;
816 
817  result.Add(nominator_temp);
818 
819  if( result == old_result )
820  // there's no sense to calculate more
821  break;
822 
823  old_result = result;
824 
825  c += nominator_x.Mul(nominator_x_add);
826  c += nominator_add.Add(two);
827  c += denominator_add.Add(two);
828  c += nominator.Mul(nominator_add);
829  c += denominator.Mul(denominator_add);
830  c += denominator_x.Add(two);
831  c += denominator2.Mul(two);
832  }
833 
834 
835  nominator_x_add.exponent.AddOne(); // *2
836  one.exponent.SubOne(); // =0.5
837  nominator_x_add.Pow(one); // =sqrt(nominator_x_add)
838  result.Mul(nominator_x_add);
839 
840  one.Set05Pi();
841  one.Sub(result);
842 
843  return one;
844  }
845 
846 
847  } // namespace auxiliaryfunctions
848 
849 
854  template<class ValueType>
855  ValueType ASin(ValueType x, ErrorCode * err = 0)
856  {
857  using namespace auxiliaryfunctions;
858 
859  ValueType result, one;
860  one.SetOne();
861  bool change_sign = false;
862 
863  if( x.IsNan() )
864  {
865  if( err )
866  *err = err_improper_argument;
867 
868  return x;
869  }
870 
871  if( x.GreaterWithoutSignThan(one) )
872  {
873  if( err )
874  *err = err_improper_argument;
875 
876  return result; // NaN is set by default
877  }
878 
879  if( x.IsSign() )
880  {
881  change_sign = true;
882  x.Abs();
883  }
884 
885  one.exponent.SubOne(); // =0.5
886 
887  // asin(-x) = -asin(x)
888  if( x.GreaterWithoutSignThan(one) )
889  result = ASin_1(x);
890  else
891  result = ASin_0(x);
892 
893  if( change_sign )
894  result.ChangeSign();
895 
896  if( err )
897  *err = err_ok;
898 
899  return result;
900  }
901 
902 
909  template<class ValueType>
910  ValueType ACos(const ValueType & x, ErrorCode * err = 0)
911  {
912  ValueType temp;
913 
914  temp.Set05Pi();
915  temp.Sub(ASin(x, err));
916 
917  return temp;
918  }
919 
920 
921 
922  namespace auxiliaryfunctions
923  {
924 
934  template<class ValueType>
935  ValueType ATan0(const ValueType & x)
936  {
937  ValueType nominator, denominator, nominator_add, denominator_add, temp;
938  ValueType result, old_result;
939  bool adding = false;
940  uint c = 0;
941 
942  result = x;
943  old_result = result;
944  nominator = x;
945  nominator_add = x;
946  nominator_add.Mul(x);
947 
948  denominator.SetOne();
949  denominator_add = 2;
950 
951  for(uint i=1 ; i<=TTMATH_ARITHMETIC_MAX_LOOP ; ++i)
952  {
953  c += nominator.Mul(nominator_add);
954  c += denominator.Add(denominator_add);
955 
956  temp = nominator;
957  c += temp.Div(denominator);
958 
959  if( c )
960  // the result should be ok
961  break;
962 
963  if( adding )
964  result.Add(temp);
965  else
966  result.Sub(temp);
967 
968  if( result == old_result )
969  // there's no sense to calculate more
970  break;
971 
972  old_result = result;
973  adding = !adding;
974  }
975 
976  return result;
977  }
978 
979 
985  template<class ValueType>
986  ValueType ATan01(const ValueType & x)
987  {
988  ValueType half;
989  half.Set05();
990 
991  /*
992  it would be better if we chose about sqrt(2)-1=0.41... instead of 0.5 here
993 
994  because as you can see below:
995  when x = sqrt(2)-1
996  abs(x) = abs( (x-1)/(1+x) )
997  so when we're calculating values around x
998  then they will be better converged to each other
999 
1000  for example if we have x=0.4999 then during calculating ATan0(0.4999)
1001  we have to make about 141 iterations but when we have x=0.5
1002  then during calculating ATan0( (x-1)/(1+x) ) we have to make
1003  only about 89 iterations (both for Big<3,9>)
1004 
1005  in the future this 0.5 can be changed
1006  */
1007  if( x.SmallerWithoutSignThan(half) )
1008  return ATan0(x);
1009 
1010 
1011  /*
1012  x>=0.5 and x<=1
1013  (x can be even smaller than 0.5)
1014 
1015  y = atac(x)
1016  x = tan(y)
1017 
1018  tan(y-b) = (tan(y)-tab(b)) / (1+tan(y)*tan(b))
1019  y-b = atan( (tan(y)-tab(b)) / (1+tan(y)*tan(b)) )
1020  y = b + atan( (x-tab(b)) / (1+x*tan(b)) )
1021 
1022  let b = pi/4
1023  tan(b) = tan(pi/4) = 1
1024  y = pi/4 + atan( (x-1)/(1+x) )
1025 
1026  so
1027  atac(x) = pi/4 + atan( (x-1)/(1+x) )
1028  when x->1 (x converges to 1) the (x-1)/(1+x) -> 0
1029  and we can use ATan0() function here
1030  */
1031 
1032  ValueType n(x),d(x),one,result;
1033 
1034  one.SetOne();
1035  n.Sub(one);
1036  d.Add(one);
1037  n.Div(d);
1038 
1039  result = ATan0(n);
1040 
1041  n.Set05Pi();
1042  n.exponent.SubOne(); // =pi/4
1043  result.Add(n);
1044 
1045  return result;
1046  }
1047 
1048 
1056  template<class ValueType>
1057  ValueType ATanGreaterThanPlusOne(const ValueType & x)
1058  {
1059  ValueType temp, atan;
1060 
1061  temp.SetOne();
1062 
1063  if( temp.Div(x) )
1064  {
1065  // if there was a carry here that means x is very big
1066  // and atan(1/x) fast converged to 0
1067  atan.SetZero();
1068  }
1069  else
1070  atan = ATan01(temp);
1071 
1072  temp.Set05Pi();
1073  temp.Sub(atan);
1074 
1075  return temp;
1076  }
1077 
1078  } // namespace auxiliaryfunctions
1079 
1080 
1084  template<class ValueType>
1085  ValueType ATan(ValueType x)
1086  {
1087  using namespace auxiliaryfunctions;
1088 
1089  ValueType one, result;
1090  one.SetOne();
1091  bool change_sign = false;
1092 
1093  if( x.IsNan() )
1094  return x;
1095 
1096  // if x is negative we're using the formula:
1097  // atan(-x) = -atan(x)
1098  if( x.IsSign() )
1099  {
1100  change_sign = true;
1101  x.Abs();
1102  }
1103 
1104  if( x.GreaterWithoutSignThan(one) )
1105  result = ATanGreaterThanPlusOne(x);
1106  else
1107  result = ATan01(x);
1108 
1109  if( change_sign )
1110  result.ChangeSign();
1111 
1112  return result;
1113  }
1114 
1115 
1122  template<class ValueType>
1123  ValueType ATg(const ValueType & x)
1124  {
1125  return ATan(x);
1126  }
1127 
1128 
1135  template<class ValueType>
1136  ValueType ACot(const ValueType & x)
1137  {
1138  ValueType result;
1139 
1140  result.Set05Pi();
1141  result.Sub(ATan(x));
1142 
1143  return result;
1144  }
1145 
1146 
1153  template<class ValueType>
1154  ValueType ACtg(const ValueType & x)
1155  {
1156  return ACot(x);
1157  }
1158 
1159 
1160  /*
1161  *
1162  * hyperbolic functions
1163  *
1164  *
1165  */
1166 
1167 
1173  template<class ValueType>
1174  ValueType Sinh(const ValueType & x, ErrorCode * err = 0)
1175  {
1176  if( x.IsNan() )
1177  {
1178  if( err )
1179  *err = err_improper_argument;
1180 
1181  return x; // NaN
1182  }
1183 
1184  ValueType ex, emx;
1185  uint c = 0;
1186 
1187  c += ex.Exp(x);
1188  c += emx.Exp(-x);
1189 
1190  c += ex.Sub(emx);
1191  c += ex.exponent.SubOne();
1192 
1193  if( err )
1194  *err = c ? err_overflow : err_ok;
1195 
1196  return ex;
1197  }
1198 
1199 
1205  template<class ValueType>
1206  ValueType Cosh(const ValueType & x, ErrorCode * err = 0)
1207  {
1208  if( x.IsNan() )
1209  {
1210  if( err )
1211  *err = err_improper_argument;
1212 
1213  return x; // NaN
1214  }
1215 
1216  ValueType ex, emx;
1217  uint c = 0;
1218 
1219  c += ex.Exp(x);
1220  c += emx.Exp(-x);
1221 
1222  c += ex.Add(emx);
1223  c += ex.exponent.SubOne();
1224 
1225  if( err )
1226  *err = c ? err_overflow : err_ok;
1227 
1228  return ex;
1229  }
1230 
1231 
1237  template<class ValueType>
1238  ValueType Tanh(const ValueType & x, ErrorCode * err = 0)
1239  {
1240  if( x.IsNan() )
1241  {
1242  if( err )
1243  *err = err_improper_argument;
1244 
1245  return x; // NaN
1246  }
1247 
1248  ValueType ex, emx, nominator, denominator;
1249  uint c = 0;
1250 
1251  c += ex.Exp(x);
1252  c += emx.Exp(-x);
1253 
1254  nominator = ex;
1255  c += nominator.Sub(emx);
1256  denominator = ex;
1257  c += denominator.Add(emx);
1258 
1259  c += nominator.Div(denominator);
1260 
1261  if( err )
1262  *err = c ? err_overflow : err_ok;
1263 
1264  return nominator;
1265  }
1266 
1267 
1274  template<class ValueType>
1275  ValueType Tgh(const ValueType & x, ErrorCode * err = 0)
1276  {
1277  return Tanh(x, err);
1278  }
1279 
1285  template<class ValueType>
1286  ValueType Coth(const ValueType & x, ErrorCode * err = 0)
1287  {
1288  if( x.IsNan() )
1289  {
1290  if( err )
1291  *err = err_improper_argument;
1292 
1293  return x; // NaN
1294  }
1295 
1296  if( x.IsZero() )
1297  {
1298  if( err )
1299  *err = err_improper_argument;
1300 
1301  return ValueType(); // NaN is set by default
1302  }
1303 
1304  ValueType ex, emx, nominator, denominator;
1305  uint c = 0;
1306 
1307  c += ex.Exp(x);
1308  c += emx.Exp(-x);
1309 
1310  nominator = ex;
1311  c += nominator.Add(emx);
1312  denominator = ex;
1313  c += denominator.Sub(emx);
1314 
1315  c += nominator.Div(denominator);
1316 
1317  if( err )
1318  *err = c ? err_overflow : err_ok;
1319 
1320  return nominator;
1321  }
1322 
1323 
1330  template<class ValueType>
1331  ValueType Ctgh(const ValueType & x, ErrorCode * err = 0)
1332  {
1333  return Coth(x, err);
1334  }
1335 
1336 
1337  /*
1338  *
1339  * inverse hyperbolic functions
1340  *
1341  *
1342  */
1343 
1344 
1350  template<class ValueType>
1351  ValueType ASinh(const ValueType & x, ErrorCode * err = 0)
1352  {
1353  if( x.IsNan() )
1354  {
1355  if( err )
1356  *err = err_improper_argument;
1357 
1358  return x; // NaN
1359  }
1360 
1361  ValueType xx(x), one, result;
1362  uint c = 0;
1363  one.SetOne();
1364 
1365  c += xx.Mul(x);
1366  c += xx.Add(one);
1367  one.exponent.SubOne(); // one=0.5
1368  // xx is >= 1
1369  c += xx.PowFrac(one); // xx=sqrt(xx)
1370  c += xx.Add(x);
1371  c += result.Ln(xx); // xx > 0
1372 
1373  // here can only be a carry
1374  if( err )
1375  *err = c ? err_overflow : err_ok;
1376 
1377  return result;
1378  }
1379 
1380 
1386  template<class ValueType>
1387  ValueType ACosh(const ValueType & x, ErrorCode * err = 0)
1388  {
1389  if( x.IsNan() )
1390  {
1391  if( err )
1392  *err = err_improper_argument;
1393 
1394  return x; // NaN
1395  }
1396 
1397  ValueType xx(x), one, result;
1398  uint c = 0;
1399  one.SetOne();
1400 
1401  if( x < one )
1402  {
1403  if( err )
1404  *err = err_improper_argument;
1405 
1406  return result; // NaN is set by default
1407  }
1408 
1409  c += xx.Mul(x);
1410  c += xx.Sub(one);
1411  // xx is >= 0
1412  // we can't call a PowFrac when the 'x' is zero
1413  // if x is 0 the sqrt(0) is 0
1414  if( !xx.IsZero() )
1415  {
1416  one.exponent.SubOne(); // one=0.5
1417  c += xx.PowFrac(one); // xx=sqrt(xx)
1418  }
1419  c += xx.Add(x);
1420  c += result.Ln(xx); // xx >= 1
1421 
1422  // here can only be a carry
1423  if( err )
1424  *err = c ? err_overflow : err_ok;
1425 
1426  return result;
1427  }
1428 
1429 
1435  template<class ValueType>
1436  ValueType ATanh(const ValueType & x, ErrorCode * err = 0)
1437  {
1438  if( x.IsNan() )
1439  {
1440  if( err )
1441  *err = err_improper_argument;
1442 
1443  return x; // NaN
1444  }
1445 
1446  ValueType nominator(x), denominator, one, result;
1447  uint c = 0;
1448  one.SetOne();
1449 
1450  if( !x.SmallerWithoutSignThan(one) )
1451  {
1452  if( err )
1453  *err = err_improper_argument;
1454 
1455  return result; // NaN is set by default
1456  }
1457 
1458  c += nominator.Add(one);
1459  denominator = one;
1460  c += denominator.Sub(x);
1461  c += nominator.Div(denominator);
1462  c += result.Ln(nominator);
1463  c += result.exponent.SubOne();
1464 
1465  // here can only be a carry
1466  if( err )
1467  *err = c ? err_overflow : err_ok;
1468 
1469  return result;
1470  }
1471 
1472 
1476  template<class ValueType>
1477  ValueType ATgh(const ValueType & x, ErrorCode * err = 0)
1478  {
1479  return ATanh(x, err);
1480  }
1481 
1482 
1488  template<class ValueType>
1489  ValueType ACoth(const ValueType & x, ErrorCode * err = 0)
1490  {
1491  if( x.IsNan() )
1492  {
1493  if( err )
1494  *err = err_improper_argument;
1495 
1496  return x; // NaN
1497  }
1498 
1499  ValueType nominator(x), denominator(x), one, result;
1500  uint c = 0;
1501  one.SetOne();
1502 
1503  if( !x.GreaterWithoutSignThan(one) )
1504  {
1505  if( err )
1506  *err = err_improper_argument;
1507 
1508  return result; // NaN is set by default
1509  }
1510 
1511  c += nominator.Add(one);
1512  c += denominator.Sub(one);
1513  c += nominator.Div(denominator);
1514  c += result.Ln(nominator);
1515  c += result.exponent.SubOne();
1516 
1517  // here can only be a carry
1518  if( err )
1519  *err = c ? err_overflow : err_ok;
1520 
1521  return result;
1522  }
1523 
1524 
1528  template<class ValueType>
1529  ValueType ACtgh(const ValueType & x, ErrorCode * err = 0)
1530  {
1531  return ACoth(x, err);
1532  }
1533 
1534 
1535 
1536 
1537 
1538  /*
1539  *
1540  * functions for converting between degrees, radians and gradians
1541  *
1542  *
1543  */
1544 
1545 
1551  template<class ValueType>
1552  ValueType DegToRad(const ValueType & x, ErrorCode * err = 0)
1553  {
1554  ValueType result, temp;
1555  uint c = 0;
1556 
1557  if( x.IsNan() )
1558  {
1559  if( err )
1560  *err = err_improper_argument;
1561 
1562  return x;
1563  }
1564 
1565  result = x;
1566 
1567  // it is better to make division first and then multiplication
1568  // the result is more accurate especially when x is: 90,180,270 or 360
1569  temp = 180;
1570  c += result.Div(temp);
1571 
1572  temp.SetPi();
1573  c += result.Mul(temp);
1574 
1575  if( err )
1576  *err = c ? err_overflow : err_ok;
1577 
1578  return result;
1579  }
1580 
1581 
1587  template<class ValueType>
1588  ValueType RadToDeg(const ValueType & x, ErrorCode * err = 0)
1589  {
1590  ValueType result, delimiter;
1591  uint c = 0;
1592 
1593  if( x.IsNan() )
1594  {
1595  if( err )
1596  *err = err_improper_argument;
1597 
1598  return x;
1599  }
1600 
1601  result = 180;
1602  c += result.Mul(x);
1603 
1604  delimiter.SetPi();
1605  c += result.Div(delimiter);
1606 
1607  if( err )
1608  *err = c ? err_overflow : err_ok;
1609 
1610  return result;
1611  }
1612 
1613 
1632  template<class ValueType>
1633  ValueType DegToDeg( const ValueType & d, const ValueType & m, const ValueType & s,
1634  ErrorCode * err = 0)
1635  {
1636  ValueType delimiter, multipler;
1637  uint c = 0;
1638 
1639  if( d.IsNan() || m.IsNan() || s.IsNan() || m.IsSign() || s.IsSign() )
1640  {
1641  if( err )
1642  *err = err_improper_argument;
1643 
1644  delimiter.SetZeroNan(); // not needed, only to get rid of GCC warning about an uninitialized variable
1645 
1646  return delimiter;
1647  }
1648 
1649  multipler = 60;
1650  delimiter = 3600;
1651 
1652  c += multipler.Mul(m);
1653  c += multipler.Add(s);
1654  c += multipler.Div(delimiter);
1655 
1656  if( d.IsSign() )
1657  multipler.ChangeSign();
1658 
1659  c += multipler.Add(d);
1660 
1661  if( err )
1662  *err = c ? err_overflow : err_ok;
1663 
1664  return multipler;
1665  }
1666 
1667 
1671  template<class ValueType>
1672  ValueType DegToRad( const ValueType & d, const ValueType & m, const ValueType & s,
1673  ErrorCode * err = 0)
1674  {
1675  ValueType temp_deg = DegToDeg(d,m,s,err);
1676 
1677  if( err && *err!=err_ok )
1678  return temp_deg;
1679 
1680  return DegToRad(temp_deg, err);
1681  }
1682 
1683 
1689  template<class ValueType>
1690  ValueType GradToRad(const ValueType & x, ErrorCode * err = 0)
1691  {
1692  ValueType result, temp;
1693  uint c = 0;
1694 
1695  if( x.IsNan() )
1696  {
1697  if( err )
1698  *err = err_improper_argument;
1699 
1700  return x;
1701  }
1702 
1703  result = x;
1704 
1705  // it is better to make division first and then multiplication
1706  // the result is more accurate especially when x is: 100,200,300 or 400
1707  temp = 200;
1708  c += result.Div(temp);
1709 
1710  temp.SetPi();
1711  c += result.Mul(temp);
1712 
1713  if( err )
1714  *err = c ? err_overflow : err_ok;
1715 
1716  return result;
1717  }
1718 
1719 
1725  template<class ValueType>
1726  ValueType RadToGrad(const ValueType & x, ErrorCode * err = 0)
1727  {
1728  ValueType result, delimiter;
1729  uint c = 0;
1730 
1731  if( x.IsNan() )
1732  {
1733  if( err )
1734  *err = err_improper_argument;
1735 
1736  return x;
1737  }
1738 
1739  result = 200;
1740  c += result.Mul(x);
1741 
1742  delimiter.SetPi();
1743  c += result.Div(delimiter);
1744 
1745  if( err )
1746  *err = c ? err_overflow : err_ok;
1747 
1748  return result;
1749  }
1750 
1751 
1757  template<class ValueType>
1758  ValueType DegToGrad(const ValueType & x, ErrorCode * err = 0)
1759  {
1760  ValueType result, temp;
1761  uint c = 0;
1762 
1763  if( x.IsNan() )
1764  {
1765  if( err )
1766  *err = err_improper_argument;
1767 
1768  return x;
1769  }
1770 
1771  result = x;
1772 
1773  temp = 200;
1774  c += result.Mul(temp);
1775 
1776  temp = 180;
1777  c += result.Div(temp);
1778 
1779  if( err )
1780  *err = c ? err_overflow : err_ok;
1781 
1782  return result;
1783  }
1784 
1785 
1789  template<class ValueType>
1790  ValueType DegToGrad( const ValueType & d, const ValueType & m, const ValueType & s,
1791  ErrorCode * err = 0)
1792  {
1793  ValueType temp_deg = DegToDeg(d,m,s,err);
1794 
1795  if( err && *err!=err_ok )
1796  return temp_deg;
1797 
1798  return DegToGrad(temp_deg, err);
1799  }
1800 
1801 
1807  template<class ValueType>
1808  ValueType GradToDeg(const ValueType & x, ErrorCode * err = 0)
1809  {
1810  ValueType result, temp;
1811  uint c = 0;
1812 
1813  if( x.IsNan() )
1814  {
1815  if( err )
1816  *err = err_improper_argument;
1817 
1818  return x;
1819  }
1820 
1821  result = x;
1822 
1823  temp = 180;
1824  c += result.Mul(temp);
1825 
1826  temp = 200;
1827  c += result.Div(temp);
1828 
1829  if( err )
1830  *err = c ? err_overflow : err_ok;
1831 
1832  return result;
1833  }
1834 
1835 
1836 
1837 
1838  /*
1839  *
1840  * another functions
1841  *
1842  *
1843  */
1844 
1845 
1851  template<class ValueType>
1852  ValueType Sqrt(ValueType x, ErrorCode * err = 0)
1853  {
1854  if( x.IsNan() || x.IsSign() )
1855  {
1856  if( err )
1857  *err = err_improper_argument;
1858 
1859  x.SetNan();
1860 
1861  return x;
1862  }
1863 
1864  uint c = x.Sqrt();
1865 
1866  if( err )
1867  *err = c ? err_overflow : err_ok;
1868 
1869  return x;
1870  }
1871 
1872 
1873 
1874  namespace auxiliaryfunctions
1875  {
1876 
1877  template<class ValueType>
1878  bool RootCheckIndexSign(ValueType & x, const ValueType & index, ErrorCode * err)
1879  {
1880  if( index.IsSign() )
1881  {
1882  // index cannot be negative
1883  if( err )
1884  *err = err_improper_argument;
1885 
1886  x.SetNan();
1887 
1888  return true;
1889  }
1890 
1891  return false;
1892  }
1893 
1894 
1895  template<class ValueType>
1896  bool RootCheckIndexZero(ValueType & x, const ValueType & index, ErrorCode * err)
1897  {
1898  if( index.IsZero() )
1899  {
1900  if( x.IsZero() )
1901  {
1902  // there isn't root(0;0) - we assume it's not defined
1903  if( err )
1904  *err = err_improper_argument;
1905 
1906  x.SetNan();
1907 
1908  return true;
1909  }
1910 
1911  // root(x;0) is 1 (if x!=0)
1912  x.SetOne();
1913 
1914  if( err )
1915  *err = err_ok;
1916 
1917  return true;
1918  }
1919 
1920  return false;
1921  }
1922 
1923 
1924  template<class ValueType>
1925  bool RootCheckIndexOne(const ValueType & index, ErrorCode * err)
1926  {
1927  ValueType one;
1928  one.SetOne();
1929 
1930  if( index == one )
1931  {
1932  //root(x;1) is x
1933  // we do it because if we used the PowFrac function
1934  // we would lose the precision
1935  if( err )
1936  *err = err_ok;
1937 
1938  return true;
1939  }
1940 
1941  return false;
1942  }
1943 
1944 
1945  template<class ValueType>
1946  bool RootCheckIndexTwo(ValueType & x, const ValueType & index, ErrorCode * err)
1947  {
1948  if( index == 2 )
1949  {
1950  x = Sqrt(x, err);
1951 
1952  return true;
1953  }
1954 
1955  return false;
1956  }
1957 
1958 
1959  template<class ValueType>
1960  bool RootCheckIndexFrac(ValueType & x, const ValueType & index, ErrorCode * err)
1961  {
1962  if( !index.IsInteger() )
1963  {
1964  // index must be integer
1965  if( err )
1966  *err = err_improper_argument;
1967 
1968  x.SetNan();
1969 
1970  return true;
1971  }
1972 
1973  return false;
1974  }
1975 
1976 
1977  template<class ValueType>
1978  bool RootCheckXZero(ValueType & x, ErrorCode * err)
1979  {
1980  if( x.IsZero() )
1981  {
1982  // root(0;index) is zero (if index!=0)
1983  // RootCheckIndexZero() must be called beforehand
1984  x.SetZero();
1985 
1986  if( err )
1987  *err = err_ok;
1988 
1989  return true;
1990  }
1991 
1992  return false;
1993  }
1994 
1995 
1996  template<class ValueType>
1997  bool RootCheckIndex(ValueType & x, const ValueType & index, ErrorCode * err, bool * change_sign)
1998  {
1999  *change_sign = false;
2000 
2001  if( index.Mod2() )
2002  {
2003  // index is odd (1,3,5...)
2004  if( x.IsSign() )
2005  {
2006  *change_sign = true;
2007  x.Abs();
2008  }
2009  }
2010  else
2011  {
2012  // index is even
2013  // x cannot be negative
2014  if( x.IsSign() )
2015  {
2016  if( err )
2017  *err = err_improper_argument;
2018 
2019  x.SetNan();
2020 
2021  return true;
2022  }
2023  }
2024 
2025  return false;
2026  }
2027 
2028 
2029  template<class ValueType>
2030  uint RootCorrectInteger(ValueType & old_x, ValueType & x, const ValueType & index)
2031  {
2032  if( !old_x.IsInteger() || x.IsInteger() || !index.exponent.IsSign() )
2033  return 0;
2034 
2035  // old_x is integer,
2036  // x is not integer,
2037  // index is relatively small (index.exponent<0 or index.exponent<=0)
2038  // (because we're using a special powering algorithm Big::PowUInt())
2039 
2040  uint c = 0;
2041 
2042  ValueType temp(x);
2043  c += temp.Round();
2044 
2045  ValueType temp_round(temp);
2046  c += temp.PowUInt(index);
2047 
2048  if( temp == old_x )
2049  x = temp_round;
2050 
2051  return (c==0)? 0 : 1;
2052  }
2053 
2054 
2055 
2056  } // namespace auxiliaryfunctions
2057 
2058 
2059 
2075  template<class ValueType>
2076  ValueType Root(ValueType x, const ValueType & index, ErrorCode * err = 0)
2077  {
2078  using namespace auxiliaryfunctions;
2079 
2080  if( x.IsNan() || index.IsNan() )
2081  {
2082  if( err )
2083  *err = err_improper_argument;
2084 
2085  x.SetNan();
2086 
2087  return x;
2088  }
2089 
2090  if( RootCheckIndexSign(x, index, err) ) return x;
2091  if( RootCheckIndexZero(x, index, err) ) return x;
2092  if( RootCheckIndexOne ( index, err) ) return x;
2093  if( RootCheckIndexTwo (x, index, err) ) return x;
2094  if( RootCheckIndexFrac(x, index, err) ) return x;
2095  if( RootCheckXZero (x, err) ) return x;
2096 
2097  // index integer and index!=0
2098  // x!=0
2099 
2100  ValueType old_x(x);
2101  bool change_sign;
2102 
2103  if( RootCheckIndex(x, index, err, &change_sign ) ) return x;
2104 
2105  ValueType temp;
2106  uint c = 0;
2107 
2108  // we're using the formula: root(x ; n) = exp( ln(x) / n )
2109  c += temp.Ln(x);
2110  c += temp.Div(index);
2111  c += x.Exp(temp);
2112 
2113  if( change_sign )
2114  {
2115  // x is different from zero
2116  x.SetSign();
2117  }
2118 
2119  c += RootCorrectInteger(old_x, x, index);
2120 
2121  if( err )
2122  *err = c ? err_overflow : err_ok;
2123 
2124  return x;
2125  }
2126 
2127 
2128 
2136  template<class ValueType>
2137  ValueType Abs(const ValueType & x)
2138  {
2139  ValueType result( x );
2140  result.Abs();
2141 
2142  return result;
2143  }
2144 
2145 
2154  template<class ValueType>
2155  ValueType Sgn(ValueType x)
2156  {
2157  x.Sgn();
2158 
2159  return x;
2160  }
2161 
2162 
2172  template<class ValueType>
2173  ValueType Mod(ValueType a, const ValueType & b, ErrorCode * err = 0)
2174  {
2175  if( a.IsNan() || b.IsNan() )
2176  {
2177  if( err )
2178  *err = err_improper_argument;
2179 
2180  a.SetNan();
2181 
2182  return a;
2183  }
2184 
2185  uint c = a.Mod(b);
2186 
2187  if( err )
2188  *err = c ? err_overflow : err_ok;
2189 
2190  return a;
2191  }
2192 
2193 
2194 
2195  namespace auxiliaryfunctions
2196  {
2197 
2211  template<class ValueType>
2212  void SetFactorialSequence(std::vector<ValueType> & fact, uint more = 20)
2213  {
2214  if( more == 0 )
2215  more = 1;
2216 
2217  uint start = static_cast<uint>(fact.size());
2218  fact.resize(fact.size() + more);
2219 
2220  if( start == 0 )
2221  {
2222  fact[0] = 1;
2223  ++start;
2224  }
2225 
2226  for(uint i=start ; i<fact.size() ; ++i)
2227  {
2228  fact[i] = fact[i-1];
2229  fact[i].MulInt(i);
2230  }
2231  }
2232 
2233 
2246  template<class ValueType>
2247  ValueType SetBernoulliNumbersSum(CGamma<ValueType> & cgamma, const ValueType & n_, uint m,
2248  const volatile StopCalculating * stop = 0)
2249  {
2250  ValueType k_, temp, temp2, temp3, sum;
2251 
2252  sum.SetZero();
2253 
2254  for(uint k=0 ; k<m ; ++k) // k<m means k<=m-1
2255  {
2256  if( stop && (k & 15)==0 ) // means: k % 16 == 0
2257  if( stop->WasStopSignal() )
2258  return ValueType(); // NaN
2259 
2260  if( k>1 && (k & 1) == 1 ) // for that k the Bernoulli number is zero
2261  continue;
2262 
2263  k_ = k;
2264 
2265  temp = n_; // n_ is equal 2
2266  temp.Pow(k_);
2267  // temp = 2^k
2268 
2269  temp2 = cgamma.fact[m];
2270  temp3 = cgamma.fact[k];
2271  temp3.Mul(cgamma.fact[m-k]);
2272  temp2.Div(temp3);
2273  // temp2 = (m k) = m! / ( k! * (m-k)! )
2274 
2275  temp.Mul(temp2);
2276  temp.Mul(cgamma.bern[k]);
2277 
2278  sum.Add(temp);
2279  // sum += 2^k * (m k) * B(k)
2280 
2281  if( sum.IsNan() )
2282  break;
2283  }
2284 
2285  return sum;
2286  }
2287 
2288 
2298  template<class ValueType>
2299  bool SetBernoulliNumbersMore(CGamma<ValueType> & cgamma, uint start, const volatile StopCalculating * stop = 0)
2300  {
2301  ValueType denominator, temp, temp2, temp3, m_, sum, sum2, n_, k_;
2302 
2303  const uint n = 2;
2304  n_ = n;
2305 
2306  // start is >= 2
2307  for(uint m=start ; m<cgamma.bern.size() ; ++m)
2308  {
2309  if( (m & 1) == 1 )
2310  {
2311  cgamma.bern[m].SetZero();
2312  }
2313  else
2314  {
2315  m_ = m;
2316 
2317  temp = n_; // n_ = 2
2318  temp.Pow(m_);
2319  // temp = 2^m
2320 
2321  denominator.SetOne();
2322  denominator.Sub(temp);
2323  if( denominator.exponent.AddOne() ) // it means: denominator.MulInt(2)
2324  denominator.SetNan();
2325 
2326  // denominator = 2 * (1 - 2^m)
2327 
2328  cgamma.bern[m] = SetBernoulliNumbersSum(cgamma, n_, m, stop);
2329 
2330  if( stop && stop->WasStopSignal() )
2331  {
2332  cgamma.bern.resize(m); // valid numbers are in [0, m-1]
2333  return false;
2334  }
2335 
2336  cgamma.bern[m].Div(denominator);
2337  }
2338  }
2339 
2340  return true;
2341  }
2342 
2343 
2359  template<class ValueType>
2360  bool SetBernoulliNumbers(CGamma<ValueType> & cgamma, uint more = 20, const volatile StopCalculating * stop = 0)
2361  {
2362  if( more == 0 )
2363  more = 1;
2364 
2365  uint start = static_cast<uint>(cgamma.bern.size());
2366  cgamma.bern.resize(cgamma.bern.size() + more);
2367 
2368  if( start == 0 )
2369  {
2370  cgamma.bern[0].SetOne();
2371  ++start;
2372  }
2373 
2374  if( cgamma.bern.size() == 1 )
2375  return true;
2376 
2377  if( start == 1 )
2378  {
2379  cgamma.bern[1].Set05();
2380  cgamma.bern[1].ChangeSign();
2381  ++start;
2382  }
2383 
2384  // we should have sufficient factorials in cgamma.fact
2385  if( cgamma.fact.size() < cgamma.bern.size() )
2386  SetFactorialSequence(cgamma.fact, static_cast<uint>(cgamma.bern.size() - cgamma.fact.size()));
2387 
2388 
2389  return SetBernoulliNumbersMore(cgamma, start, stop);
2390  }
2391 
2392 
2403  template<class ValueType>
2404  ValueType GammaFactorialHighSum(const ValueType & n, CGamma<ValueType> & cgamma, ErrorCode & err,
2405  const volatile StopCalculating * stop)
2406  {
2407  ValueType temp, temp2, denominator, sum, oldsum;
2408 
2409  sum.SetZero();
2410 
2411  for(uint m=2 ; m<TTMATH_ARITHMETIC_MAX_LOOP ; m+=2)
2412  {
2413  if( stop && (m & 3)==0 ) // (m & 3)==0 means: (m % 4)==0
2414  if( stop->WasStopSignal() )
2415  {
2416  err = err_interrupt;
2417  return ValueType(); // NaN
2418  }
2419 
2420  temp = (m-1);
2421  denominator = n;
2422  denominator.Pow(temp);
2423  // denominator = n ^ (m-1)
2424 
2425  temp = m;
2426  temp2 = temp;
2427  temp.Mul(temp2);
2428  temp.Sub(temp2);
2429  // temp = m^2 - m
2430 
2431  denominator.Mul(temp);
2432  // denominator = (m^2 - m) * n ^ (m-1)
2433 
2434  if( m >= cgamma.bern.size() )
2435  {
2436  if( !SetBernoulliNumbers(cgamma, m - cgamma.bern.size() + 1 + 3, stop) ) // 3 more than needed
2437  {
2438  // there was the stop signal
2439  err = err_interrupt;
2440  return ValueType(); // NaN
2441  }
2442  }
2443 
2444  temp = cgamma.bern[m];
2445  temp.Div(denominator);
2446 
2447  oldsum = sum;
2448  sum.Add(temp);
2449 
2450  if( sum.IsNan() || oldsum==sum )
2451  break;
2452  }
2453 
2454  return sum;
2455  }
2456 
2457 
2468  template<class ValueType>
2469  ValueType GammaFactorialHigh(const ValueType & n, CGamma<ValueType> & cgamma, ErrorCode & err,
2470  const volatile StopCalculating * stop)
2471  {
2472  ValueType temp, temp2, temp3, denominator, sum;
2473 
2474  temp.Set2Pi();
2475  temp.Mul(n);
2476  temp2 = Sqrt(temp);
2477  // temp2 = sqrt(2*pi*n)
2478 
2479  temp = n;
2480  temp3.SetE();
2481  temp.Div(temp3);
2482  temp.Pow(n);
2483  // temp = (n/e)^n
2484 
2485  sum = GammaFactorialHighSum(n, cgamma, err, stop);
2486  temp3.Exp(sum);
2487  // temp3 = exp(sum)
2488 
2489  temp.Mul(temp2);
2490  temp.Mul(temp3);
2491 
2492  return temp;
2493  }
2494 
2495 
2501  template<class ValueType>
2502  ValueType GammaPlusHigh(ValueType n, CGamma<ValueType> & cgamma, ErrorCode & err, const volatile StopCalculating * stop)
2503  {
2504  ValueType one;
2505 
2506  one.SetOne();
2507  n.Sub(one);
2508 
2509  return GammaFactorialHigh(n, cgamma, err, stop);
2510  }
2511 
2512 
2521  template<class ValueType>
2523  {
2524  TTMATH_ASSERT( n > 0 )
2525 
2526  if( n - 1 < static_cast<uint>(cgamma.fact.size()) )
2527  return cgamma.fact[n - 1];
2528 
2529  ValueType res;
2530  uint start = 2;
2531 
2532  if( cgamma.fact.size() < 2 )
2533  {
2534  res.SetOne();
2535  }
2536  else
2537  {
2538  start = static_cast<uint>(cgamma.fact.size());
2539  res = cgamma.fact[start-1];
2540  }
2541 
2542  for(uint i=start ; i<n ; ++i)
2543  res.MulInt(i);
2544 
2545  return res;
2546  }
2547 
2548 
2554  template<class ValueType>
2555  ValueType GammaPlusLowInteger(const ValueType & n, CGamma<ValueType> & cgamma)
2556  {
2557  sint n_;
2558 
2559  n.ToInt(n_);
2560 
2561  return GammaPlusLowIntegerInt(n_, cgamma);
2562  }
2563 
2564 
2577  template<class ValueType>
2578  ValueType GammaPlusLow(ValueType n, CGamma<ValueType> & cgamma, ErrorCode & err, const volatile StopCalculating * stop)
2579  {
2580  ValueType one, denominator, temp, boundary;
2581 
2582  if( n.IsInteger() )
2583  return GammaPlusLowInteger(n, cgamma);
2584 
2585  one.SetOne();
2586  denominator = n;
2587  boundary = TTMATH_GAMMA_BOUNDARY;
2588 
2589  while( n < boundary )
2590  {
2591  n.Add(one);
2592  denominator.Mul(n);
2593  }
2594 
2595  n.Add(one);
2596 
2597  // now n is sufficient big
2598  temp = GammaPlusHigh(n, cgamma, err, stop);
2599  temp.Div(denominator);
2600 
2601  return temp;
2602  }
2603 
2604 
2608  template<class ValueType>
2609  ValueType GammaPlus(const ValueType & n, CGamma<ValueType> & cgamma, ErrorCode & err, const volatile StopCalculating * stop)
2610  {
2611  if( n > TTMATH_GAMMA_BOUNDARY )
2612  return GammaPlusHigh(n, cgamma, err, stop);
2613 
2614  return GammaPlusLow(n, cgamma, err, stop);
2615  }
2616 
2617 
2627  template<class ValueType>
2628  ValueType GammaMinus(const ValueType & n, CGamma<ValueType> & cgamma, ErrorCode & err, const volatile StopCalculating * stop)
2629  {
2630  ValueType pi, denominator, temp, temp2;
2631 
2632  if( n.IsInteger() )
2633  {
2634  // gamma function is not defined when n is negative and integer
2635  err = err_improper_argument;
2636  return temp; // NaN
2637  }
2638 
2639  pi.SetPi();
2640 
2641  temp = pi;
2642  temp.Mul(n);
2643  temp2 = Sin(temp);
2644  // temp2 = sin(pi * n)
2645 
2646  temp.SetOne();
2647  temp.Sub(n);
2648  temp = GammaPlus(temp, cgamma, err, stop);
2649  // temp = gamma(1 - n)
2650 
2651  temp.Mul(temp2);
2652  pi.Div(temp);
2653 
2654  return pi;
2655  }
2656 
2657  } // namespace auxiliaryfunctions
2658 
2659 
2660 
2678  template<class ValueType>
2679  ValueType Gamma(const ValueType & n, CGamma<ValueType> & cgamma, ErrorCode * err = 0,
2680  const volatile StopCalculating * stop = 0)
2681  {
2682  using namespace auxiliaryfunctions;
2683 
2684  ValueType result;
2685  ErrorCode err_tmp;
2686 
2687  if( n.IsNan() )
2688  {
2689  if( err )
2690  *err = err_improper_argument;
2691 
2692  return n;
2693  }
2694 
2695  if( cgamma.history.Get(n, result, err_tmp) )
2696  {
2697  if( err )
2698  *err = err_tmp;
2699 
2700  return result;
2701  }
2702 
2703  err_tmp = err_ok;
2704 
2705  if( n.IsSign() )
2706  {
2707  result = GammaMinus(n, cgamma, err_tmp, stop);
2708  }
2709  else
2710  if( n.IsZero() )
2711  {
2712  err_tmp = err_improper_argument;
2713  result.SetNan();
2714  }
2715  else
2716  {
2717  result = GammaPlus(n, cgamma, err_tmp, stop);
2718  }
2719 
2720  if( result.IsNan() && err_tmp==err_ok )
2721  err_tmp = err_overflow;
2722 
2723  if( err )
2724  *err = err_tmp;
2725 
2726  if( stop && !stop->WasStopSignal() )
2727  cgamma.history.Add(n, result, err_tmp);
2728 
2729  return result;
2730  }
2731 
2732 
2738  template<class ValueType>
2739  ValueType Gamma(const ValueType & n, ErrorCode * err = 0)
2740  {
2741  // warning: this static object is not thread safe
2742  static CGamma<ValueType> cgamma;
2743 
2744  return Gamma(n, cgamma, err);
2745  }
2746 
2747 
2748 
2749  namespace auxiliaryfunctions
2750  {
2751 
2758  template<class ValueType>
2759  ValueType Factorial2(ValueType x,
2760  CGamma<ValueType> * cgamma = 0,
2761  ErrorCode * err = 0,
2762  const volatile StopCalculating * stop = 0)
2763  {
2764  ValueType result, one;
2765 
2766  if( x.IsNan() || x.IsSign() || !x.IsInteger() )
2767  {
2768  if( err )
2769  *err = err_improper_argument;
2770 
2771  x.SetNan();
2772 
2773  return x;
2774  }
2775 
2776  one.SetOne();
2777  x.Add(one);
2778 
2779  if( cgamma )
2780  return Gamma(x, *cgamma, err, stop);
2781 
2782  return Gamma(x, err);
2783  }
2784 
2785  } // namespace auxiliaryfunctions
2786 
2787 
2788 
2808  template<class ValueType>
2809  ValueType Factorial(const ValueType & x, CGamma<ValueType> & cgamma, ErrorCode * err = 0,
2810  const volatile StopCalculating * stop = 0)
2811  {
2812  return auxiliaryfunctions::Factorial2(x, &cgamma, err, stop);
2813  }
2814 
2815 
2823  template<class ValueType>
2824  ValueType Factorial(const ValueType & x, ErrorCode * err = 0)
2825  {
2826  return auxiliaryfunctions::Factorial2(x, (CGamma<ValueType>*)0, err, 0);
2827  }
2828 
2829 
2839  template<class ValueType>
2841  {
2842  ValueType x = TTMATH_GAMMA_BOUNDARY + 1;
2843 
2844  // history.Remove(x) removes only one object
2845  // we must be sure that there are not others objects with the key 'x'
2846  while( history.Remove(x) )
2847  {
2848  }
2849 
2850  // the simplest way to initialize is to call the Gamma function with (TTMATH_GAMMA_BOUNDARY + 1)
2851  // when x is larger then fewer coefficients we need
2852  Gamma(x, *this);
2853  }
2854 
2855 
2856 
2857 } // namespace
2858 
2859 
2864 #include "ttmathparser.h"
2865 
2866 // Dec is not finished yet
2867 //#include "ttmathdec.h"
2868 
2869 
2870 
2871 #ifdef _MSC_VER
2872 //warning C4127: conditional expression is constant
2873 #pragma warning( default: 4127 )
2874 //warning C4702: unreachable code
2875 #pragma warning( default: 4702 )
2876 //warning C4800: forcing value to bool 'true' or 'false' (performance warning)
2877 #pragma warning( default: 4800 )
2878 #endif
2879 
2880 #endif
ValueType Sgn(ValueType x)
Definition: ttmath.h:2155
ValueType ACot(const ValueType &x)
Definition: ttmath.h:1136
ValueType GammaMinus(const ValueType &n, CGamma< ValueType > &cgamma, ErrorCode &err, const volatile StopCalculating *stop)
Definition: ttmath.h:2628
ValueType ACoth(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:1489
bool SetBernoulliNumbers(CGamma< ValueType > &cgamma, uint more=20, const volatile StopCalculating *stop=0)
Definition: ttmath.h:2360
ValueType Tan(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:619
ValueType Log(const ValueType &x, const ValueType &base, ErrorCode *err=0)
Definition: ttmath.h:278
void InitAll()
Definition: ttmath.h:2840
ValueType Sin(ValueType x, ErrorCode *err=0)
Definition: ttmath.h:525
ValueType ATg(const ValueType &x)
Definition: ttmath.h:1123
#define TTMATH_ARITHMETIC_MAX_LOOP
Definition: ttmathtypes.h:322
ValueType ATan01(const ValueType &x)
Definition: ttmath.h:986
Definition: ttmathtypes.h:554
A mathematical parser.
ValueType GammaPlusHigh(ValueType n, CGamma< ValueType > &cgamma, ErrorCode &err, const volatile StopCalculating *stop)
Definition: ttmath.h:2502
ValueType Exp(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:323
ValueType ACosh(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:1387
ValueType RadToDeg(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:1588
ValueType ASin(ValueType x, ErrorCode *err=0)
Definition: ttmath.h:855
ValueType ASinh(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:1351
ValueType SkipFraction(const ValueType &x)
Definition: ttmath.h:90
ValueType Tg(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:647
ValueType Cos(ValueType x, ErrorCode *err=0)
Definition: ttmath.h:581
ValueType ATanh(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:1436
History< ValueType > history
Definition: ttmathobjects.h:782
Definition: ttmathobjects.h:742
ValueType Ln(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:237
ValueType RadToGrad(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:1726
ValueType Abs(const ValueType &x)
Definition: ttmath.h:2137
std::vector< ValueType > fact
Definition: ttmathobjects.h:755
ValueType Sin0pi05(const ValueType &x)
Definition: ttmath.h:430
ValueType ACtgh(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:1529
ValueType GradToDeg(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:1808
ErrorCode
Definition: ttmathtypes.h:382
ValueType ATan0(const ValueType &x)
Definition: ttmath.h:935
ValueType Tgh(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:1275
a namespace for the TTMath library
Definition: ttmath.h:62
ValueType GammaPlusLowInteger(const ValueType &n, CGamma< ValueType > &cgamma)
Definition: ttmath.h:2555
ValueType Root(ValueType x, const ValueType &index, ErrorCode *err=0)
Definition: ttmath.h:2076
ValueType DegToRad(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:1552
ValueType GradToRad(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:1690
ValueType GammaPlusLow(ValueType n, CGamma< ValueType > &cgamma, ErrorCode &err, const volatile StopCalculating *stop)
Definition: ttmath.h:2578
ValueType Mod(ValueType a, const ValueType &b, ErrorCode *err=0)
Definition: ttmath.h:2173
ValueType ATanGreaterThanPlusOne(const ValueType &x)
Definition: ttmath.h:1057
ValueType Ctg(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:689
ValueType ATan(ValueType x)
Definition: ttmath.h:1085
ValueType Coth(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:1286
ValueType Factorial(const ValueType &x, CGamma< ValueType > &cgamma, ErrorCode *err=0, const volatile StopCalculating *stop=0)
Definition: ttmath.h:2809
ValueType GammaPlusLowIntegerInt(uint n, CGamma< ValueType > &cgamma)
Definition: ttmath.h:2522
ValueType Floor(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:189
ValueType Tanh(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:1238
ValueType Sinh(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:1174
ValueType Cot(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:661
ValueType ASin_1(const ValueType &x)
Definition: ttmath.h:779
ValueType DegToGrad(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:1758
ValueType Gamma(const ValueType &n, CGamma< ValueType > &cgamma, ErrorCode *err=0, const volatile StopCalculating *stop=0)
Definition: ttmath.h:2679
Mathematic functions.
ValueType Cosh(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:1206
#define TTMATH_GAMMA_BOUNDARY
Definition: ttmathtypes.h:350
ValueType Factorial2(ValueType x, CGamma< ValueType > *cgamma=0, ErrorCode *err=0, const volatile StopCalculating *stop=0)
Definition: ttmath.h:2759
ValueType Round(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:109
ValueType ACtg(const ValueType &x)
Definition: ttmath.h:1154
ValueType DegToDeg(const ValueType &d, const ValueType &m, const ValueType &s, ErrorCode *err=0)
Definition: ttmath.h:1633
ValueType ATgh(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:1477
A Class for representing floating point numbers.
void SetFactorialSequence(std::vector< ValueType > &fact, uint more=20)
Definition: ttmath.h:2212
ValueType Ctgh(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:1331
unsigned int uint
Definition: ttmathtypes.h:186
ValueType GammaFactorialHigh(const ValueType &n, CGamma< ValueType > &cgamma, ErrorCode &err, const volatile StopCalculating *stop)
Definition: ttmath.h:2469
ValueType GammaPlus(const ValueType &n, CGamma< ValueType > &cgamma, ErrorCode &err, const volatile StopCalculating *stop)
Definition: ttmath.h:2609
ValueType Ceil(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:142
ValueType ASin_0(const ValueType &x)
Definition: ttmath.h:715
ValueType GammaFactorialHighSum(const ValueType &n, CGamma< ValueType > &cgamma, ErrorCode &err, const volatile StopCalculating *stop)
Definition: ttmath.h:2404
ValueType Sqrt(ValueType x, ErrorCode *err=0)
Definition: ttmath.h:1852
ValueType ACos(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:910
std::vector< ValueType > bern
Definition: ttmathobjects.h:773
uint PrepareSin(ValueType &x, bool &change_sign)
Definition: ttmath.h:365
bool SetBernoulliNumbersMore(CGamma< ValueType > &cgamma, uint start, const volatile StopCalculating *stop=0)
Definition: ttmath.h:2299
ValueType SetBernoulliNumbersSum(CGamma< ValueType > &cgamma, const ValueType &n_, uint m, const volatile StopCalculating *stop=0)
Definition: ttmath.h:2247