All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 (new) BSD licence.
4  * Author: Tomasz Sowa <t.sowa@ttmath.org>
5  */
6 
7 /*
8  * Copyright (c) 2006-2012, 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 
87  template<class ValueType>
88  ValueType SkipFraction(const ValueType & x)
89  {
90  ValueType result( x );
91  result.SkipFraction();
92 
93  return result;
94  }
95 
96 
104  template<class ValueType>
105  ValueType Round(const ValueType & x, ErrorCode * err = 0)
106  {
107  if( x.IsNan() )
108  {
109  if( err )
110  *err = err_improper_argument;
111 
112  return x; // NaN
113  }
114 
115  ValueType result( x );
116  uint c = result.Round();
117 
118  if( err )
119  *err = c ? err_overflow : err_ok;
120 
121  return result;
122  }
123 
124 
125 
137  template<class ValueType>
138  ValueType Ceil(const ValueType & x, ErrorCode * err = 0)
139  {
140  if( x.IsNan() )
141  {
142  if( err )
143  *err = err_improper_argument;
144 
145  return x; // NaN
146  }
147 
148  ValueType result(x);
149  uint c = 0;
150 
151  result.SkipFraction();
152 
153  if( result != x )
154  {
155  // x is with fraction
156  // if x is negative we don't have to do anything
157  if( !x.IsSign() )
158  {
159  ValueType one;
160  one.SetOne();
161 
162  c += result.Add(one);
163  }
164  }
165 
166  if( err )
167  *err = c ? err_overflow : err_ok;
168 
169  return result;
170  }
171 
172 
184  template<class ValueType>
185  ValueType Floor(const ValueType & x, ErrorCode * err = 0)
186  {
187  if( x.IsNan() )
188  {
189  if( err )
190  *err = err_improper_argument;
191 
192  return x; // NaN
193  }
194 
195  ValueType result(x);
196  uint c = 0;
197 
198  result.SkipFraction();
199 
200  if( result != x )
201  {
202  // x is with fraction
203  // if x is positive we don't have to do anything
204  if( x.IsSign() )
205  {
206  ValueType one;
207  one.SetOne();
208 
209  c += result.Sub(one);
210  }
211  }
212 
213  if( err )
214  *err = c ? err_overflow : err_ok;
215 
216  return result;
217  }
218 
219 
220 
221  /*
222  *
223  * logarithms and the exponent
224  *
225  *
226  */
227 
228 
232  template<class ValueType>
233  ValueType Ln(const ValueType & x, ErrorCode * err = 0)
234  {
235  if( x.IsNan() )
236  {
237  if( err )
238  *err = err_improper_argument;
239 
240  return x; // NaN
241  }
242 
243  ValueType result;
244  uint state = result.Ln(x);
245 
246  if( err )
247  {
248  switch( state )
249  {
250  case 0:
251  *err = err_ok;
252  break;
253  case 1:
254  *err = err_overflow;
255  break;
256  case 2:
257  *err = err_improper_argument;
258  break;
259  default:
260  *err = err_internal_error;
261  break;
262  }
263  }
264 
265 
266  return result;
267  }
268 
269 
273  template<class ValueType>
274  ValueType Log(const ValueType & x, const ValueType & base, ErrorCode * err = 0)
275  {
276  if( x.IsNan() )
277  {
278  if( err ) *err = err_improper_argument;
279  return x;
280  }
281 
282  if( base.IsNan() )
283  {
284  if( err ) *err = err_improper_argument;
285  return base;
286  }
287 
288  ValueType result;
289  uint state = result.Log(x, base);
290 
291  if( err )
292  {
293  switch( state )
294  {
295  case 0:
296  *err = err_ok;
297  break;
298  case 1:
299  *err = err_overflow;
300  break;
301  case 2:
302  case 3:
303  *err = err_improper_argument;
304  break;
305  default:
306  *err = err_internal_error;
307  break;
308  }
309  }
310 
311  return result;
312  }
313 
314 
318  template<class ValueType>
319  ValueType Exp(const ValueType & x, ErrorCode * err = 0)
320  {
321  if( x.IsNan() )
322  {
323  if( err )
324  *err = err_improper_argument;
325 
326  return x; // NaN
327  }
328 
329  ValueType result;
330  uint c = result.Exp(x);
331 
332  if( err )
333  *err = c ? err_overflow : err_ok;
334 
335  return result;
336  }
337 
338 
346  /*
347  this namespace consists of auxiliary functions
348  (something like 'private' in a class)
349  */
350  namespace auxiliaryfunctions
351  {
352 
357  template<class ValueType>
358  uint PrepareSin(ValueType & x, bool & change_sign)
359  {
360  ValueType temp;
361 
362  change_sign = false;
363 
364  if( x.IsSign() )
365  {
366  // we're using the formula 'sin(-x) = -sin(x)'
367  change_sign = !change_sign;
368  x.ChangeSign();
369  }
370 
371  // we're reducing the period 2*PI
372  // (for big values there'll always be zero)
373  temp.Set2Pi();
374 
375  if( x.Mod(temp) )
376  return 1;
377 
378 
379  // we're setting 'x' as being in the range of <0, 0.5PI>
380 
381  temp.SetPi();
382 
383  if( x > temp )
384  {
385  // x is in (pi, 2*pi>
386  x.Sub( temp );
387  change_sign = !change_sign;
388  }
389 
390  temp.Set05Pi();
391 
392  if( x > temp )
393  {
394  // x is in (0.5pi, pi>
395  x.Sub( temp );
396  x = temp - x;
397  }
398 
399  return 0;
400  }
401 
402 
422  template<class ValueType>
423  ValueType Sin0pi05(const ValueType & x)
424  {
425  ValueType result;
426  ValueType numerator, denominator;
427  ValueType d_numerator, d_denominator;
428  ValueType one, temp, old_result;
429 
430  // temp = pi/4
431  temp.Set05Pi();
432  temp.exponent.SubOne();
433 
434  one.SetOne();
435 
436  if( x < temp )
437  {
438  // we're using the Taylor series with a=0
439  result = x;
440  numerator = x;
441  denominator = one;
442 
443  // d_numerator = x^2
444  d_numerator = x;
445  d_numerator.Mul(x);
446 
447  d_denominator = 2;
448  }
449  else
450  {
451  // we're using the Taylor series with a=PI/2
452  result = one;
453  numerator = one;
454  denominator = one;
455 
456  // d_numerator = (x-pi/2)^2
457  ValueType pi05;
458  pi05.Set05Pi();
459 
460  temp = x;
461  temp.Sub( pi05 );
462  d_numerator = temp;
463  d_numerator.Mul( temp );
464 
465  d_denominator = one;
466  }
467 
468  uint c = 0;
469  bool addition = false;
470 
471  old_result = result;
472  for(uint i=1 ; i<=TTMATH_ARITHMETIC_MAX_LOOP ; ++i)
473  {
474  // we're starting from a second part of the formula
475  c += numerator. Mul( d_numerator );
476  c += denominator. Mul( d_denominator );
477  c += d_denominator.Add( one );
478  c += denominator. Mul( d_denominator );
479  c += d_denominator.Add( one );
480  temp = numerator;
481  c += temp.Div(denominator);
482 
483  if( c )
484  // Sin is from <-1,1> and cannot make an overflow
485  // but the carry can be from the Taylor series
486  // (then we only break our calculations)
487  break;
488 
489  if( addition )
490  result.Add( temp );
491  else
492  result.Sub( temp );
493 
494 
495  addition = !addition;
496 
497  // we're testing whether the result has changed after adding
498  // the next part of the Taylor formula, if not we end the loop
499  // (it means 'x' is zero or 'x' is PI/2 or this part of the formula
500  // is too small)
501  if( result == old_result )
502  break;
503 
504  old_result = result;
505  }
506 
507  return result;
508  }
509 
510  } // namespace auxiliaryfunctions
511 
512 
513 
517  template<class ValueType>
518  ValueType Sin(ValueType x, ErrorCode * err = 0)
519  {
520  using namespace auxiliaryfunctions;
521 
522  ValueType one, result;
523  bool change_sign;
524 
525  if( x.IsNan() )
526  {
527  if( err )
528  *err = err_improper_argument;
529 
530  return x;
531  }
532 
533  if( err )
534  *err = err_ok;
535 
536  if( PrepareSin( x, change_sign ) )
537  {
538  // x is too big, we cannnot reduce the 2*PI period
539  // prior to version 0.8.5 the result was zero
540 
541  // result has NaN flag set by default
542 
543  if( err )
544  *err = err_overflow; // maybe another error code? err_improper_argument?
545 
546  return result; // NaN is set by default
547  }
548 
549  result = Sin0pi05( x );
550 
551  one.SetOne();
552 
553  // after calculations there can be small distortions in the result
554  if( result > one )
555  result = one;
556  else
557  if( result.IsSign() )
558  // we've calculated the sin from <0, pi/2> and the result
559  // should be positive
560  result.SetZero();
561 
562  if( change_sign )
563  result.ChangeSign();
564 
565  return result;
566  }
567 
568 
573  template<class ValueType>
574  ValueType Cos(ValueType x, ErrorCode * err = 0)
575  {
576  if( x.IsNan() )
577  {
578  if( err )
579  *err = err_improper_argument;
580 
581  return x; // NaN
582  }
583 
584  ValueType pi05;
585  pi05.Set05Pi();
586 
587  uint c = x.Add( pi05 );
588 
589  if( c )
590  {
591  if( err )
592  *err = err_overflow;
593 
594  return ValueType(); // result is undefined (NaN is set by default)
595  }
596 
597  return Sin(x, err);
598  }
599 
600 
611  template<class ValueType>
612  ValueType Tan(const ValueType & x, ErrorCode * err = 0)
613  {
614  ValueType result = Cos(x, err);
615 
616  if( err && *err != err_ok )
617  return result;
618 
619  if( result.IsZero() )
620  {
621  if( err )
622  *err = err_improper_argument;
623 
624  result.SetNan();
625 
626  return result;
627  }
628 
629  return Sin(x, err) / result;
630  }
631 
632 
639  template<class ValueType>
640  ValueType Tg(const ValueType & x, ErrorCode * err = 0)
641  {
642  return Tan(x, err);
643  }
644 
645 
653  template<class ValueType>
654  ValueType Cot(const ValueType & x, ErrorCode * err = 0)
655  {
656  ValueType result = Sin(x, err);
657 
658  if( err && *err != err_ok )
659  return result;
660 
661  if( result.IsZero() )
662  {
663  if( err )
664  *err = err_improper_argument;
665 
666  result.SetNan();
667 
668  return result;
669  }
670 
671  return Cos(x, err) / result;
672  }
673 
674 
681  template<class ValueType>
682  ValueType Ctg(const ValueType & x, ErrorCode * err = 0)
683  {
684  return Cot(x, err);
685  }
686 
687 
688  /*
689  *
690  * inverse trigonometric functions
691  *
692  *
693  */
694 
695  namespace auxiliaryfunctions
696  {
697 
707  template<class ValueType>
708  ValueType ASin_0(const ValueType & x)
709  {
710  ValueType nominator, denominator, nominator_add, nominator_x, denominator_add, denominator_x;
711  ValueType two, result(x), x2(x);
712  ValueType nominator_temp, denominator_temp, old_result = result;
713  uint c = 0;
714 
715  x2.Mul(x);
716  two = 2;
717 
718  nominator.SetOne();
719  denominator = two;
720  nominator_add = nominator;
721  denominator_add = denominator;
722  nominator_x = x;
723  denominator_x = 3;
724 
725  for(uint i=1 ; i<=TTMATH_ARITHMETIC_MAX_LOOP ; ++i)
726  {
727  c += nominator_x.Mul(x2);
728  nominator_temp = nominator_x;
729  c += nominator_temp.Mul(nominator);
730  denominator_temp = denominator;
731  c += denominator_temp.Mul(denominator_x);
732  c += nominator_temp.Div(denominator_temp);
733 
734  // if there is a carry somewhere we only break the calculating
735  // the result should be ok -- it's from <-pi/2, pi/2>
736  if( c )
737  break;
738 
739  result.Add(nominator_temp);
740 
741  if( result == old_result )
742  // there's no sense to calculate more
743  break;
744 
745  old_result = result;
746 
747 
748  c += nominator_add.Add(two);
749  c += denominator_add.Add(two);
750  c += nominator.Mul(nominator_add);
751  c += denominator.Mul(denominator_add);
752  c += denominator_x.Add(two);
753  }
754 
755  return result;
756  }
757 
758 
759 
771  template<class ValueType>
772  ValueType ASin_1(const ValueType & x)
773  {
774  ValueType nominator, denominator, nominator_add, nominator_x, nominator_x_add, denominator_add, denominator_x;
775  ValueType denominator2;
776  ValueType one, two, result;
777  ValueType nominator_temp, denominator_temp, old_result;
778  uint c = 0;
779 
780  two = 2;
781 
782  one.SetOne();
783  nominator = one;
784  result = one;
785  old_result = result;
786  denominator = two;
787  nominator_add = nominator;
788  denominator_add = denominator;
789  nominator_x = one;
790  nominator_x.Sub(x);
791  nominator_x_add = nominator_x;
792  denominator_x = 3;
793  denominator2 = two;
794 
795 
796  for(uint i=1 ; i<=TTMATH_ARITHMETIC_MAX_LOOP ; ++i)
797  {
798  nominator_temp = nominator_x;
799  c += nominator_temp.Mul(nominator);
800  denominator_temp = denominator;
801  c += denominator_temp.Mul(denominator_x);
802  c += denominator_temp.Mul(denominator2);
803  c += nominator_temp.Div(denominator_temp);
804 
805  // if there is a carry somewhere we only break the calculating
806  // the result should be ok -- it's from <-pi/2, pi/2>
807  if( c )
808  break;
809 
810  result.Add(nominator_temp);
811 
812  if( result == old_result )
813  // there's no sense to calculate more
814  break;
815 
816  old_result = result;
817 
818  c += nominator_x.Mul(nominator_x_add);
819  c += nominator_add.Add(two);
820  c += denominator_add.Add(two);
821  c += nominator.Mul(nominator_add);
822  c += denominator.Mul(denominator_add);
823  c += denominator_x.Add(two);
824  c += denominator2.Mul(two);
825  }
826 
827 
828  nominator_x_add.exponent.AddOne(); // *2
829  one.exponent.SubOne(); // =0.5
830  nominator_x_add.Pow(one); // =sqrt(nominator_x_add)
831  result.Mul(nominator_x_add);
832 
833  one.Set05Pi();
834  one.Sub(result);
835 
836  return one;
837  }
838 
839 
840  } // namespace auxiliaryfunctions
841 
842 
847  template<class ValueType>
848  ValueType ASin(ValueType x, ErrorCode * err = 0)
849  {
850  using namespace auxiliaryfunctions;
851 
852  ValueType result, one;
853  one.SetOne();
854  bool change_sign = false;
855 
856  if( x.IsNan() )
857  {
858  if( err )
859  *err = err_improper_argument;
860 
861  return x;
862  }
863 
864  if( x.GreaterWithoutSignThan(one) )
865  {
866  if( err )
867  *err = err_improper_argument;
868 
869  return result; // NaN is set by default
870  }
871 
872  if( x.IsSign() )
873  {
874  change_sign = true;
875  x.Abs();
876  }
877 
878  one.exponent.SubOne(); // =0.5
879 
880  // asin(-x) = -asin(x)
881  if( x.GreaterWithoutSignThan(one) )
882  result = ASin_1(x);
883  else
884  result = ASin_0(x);
885 
886  if( change_sign )
887  result.ChangeSign();
888 
889  if( err )
890  *err = err_ok;
891 
892  return result;
893  }
894 
895 
902  template<class ValueType>
903  ValueType ACos(const ValueType & x, ErrorCode * err = 0)
904  {
905  ValueType temp;
906 
907  temp.Set05Pi();
908  temp.Sub(ASin(x, err));
909 
910  return temp;
911  }
912 
913 
914 
915  namespace auxiliaryfunctions
916  {
917 
927  template<class ValueType>
928  ValueType ATan0(const ValueType & x)
929  {
930  ValueType nominator, denominator, nominator_add, denominator_add, temp;
931  ValueType result, old_result;
932  bool adding = false;
933  uint c = 0;
934 
935  result = x;
936  old_result = result;
937  nominator = x;
938  nominator_add = x;
939  nominator_add.Mul(x);
940 
941  denominator.SetOne();
942  denominator_add = 2;
943 
944  for(uint i=1 ; i<=TTMATH_ARITHMETIC_MAX_LOOP ; ++i)
945  {
946  c += nominator.Mul(nominator_add);
947  c += denominator.Add(denominator_add);
948 
949  temp = nominator;
950  c += temp.Div(denominator);
951 
952  if( c )
953  // the result should be ok
954  break;
955 
956  if( adding )
957  result.Add(temp);
958  else
959  result.Sub(temp);
960 
961  if( result == old_result )
962  // there's no sense to calculate more
963  break;
964 
965  old_result = result;
966  adding = !adding;
967  }
968 
969  return result;
970  }
971 
972 
978  template<class ValueType>
979  ValueType ATan01(const ValueType & x)
980  {
981  ValueType half;
982  half.Set05();
983 
984  /*
985  it would be better if we chose about sqrt(2)-1=0.41... instead of 0.5 here
986 
987  because as you can see below:
988  when x = sqrt(2)-1
989  abs(x) = abs( (x-1)/(1+x) )
990  so when we're calculating values around x
991  then they will be better converged to each other
992 
993  for example if we have x=0.4999 then during calculating ATan0(0.4999)
994  we have to make about 141 iterations but when we have x=0.5
995  then during calculating ATan0( (x-1)/(1+x) ) we have to make
996  only about 89 iterations (both for Big<3,9>)
997 
998  in the future this 0.5 can be changed
999  */
1000  if( x.SmallerWithoutSignThan(half) )
1001  return ATan0(x);
1002 
1003 
1004  /*
1005  x>=0.5 and x<=1
1006  (x can be even smaller than 0.5)
1007 
1008  y = atac(x)
1009  x = tan(y)
1010 
1011  tan(y-b) = (tan(y)-tab(b)) / (1+tan(y)*tan(b))
1012  y-b = atan( (tan(y)-tab(b)) / (1+tan(y)*tan(b)) )
1013  y = b + atan( (x-tab(b)) / (1+x*tan(b)) )
1014 
1015  let b = pi/4
1016  tan(b) = tan(pi/4) = 1
1017  y = pi/4 + atan( (x-1)/(1+x) )
1018 
1019  so
1020  atac(x) = pi/4 + atan( (x-1)/(1+x) )
1021  when x->1 (x converges to 1) the (x-1)/(1+x) -> 0
1022  and we can use ATan0() function here
1023  */
1024 
1025  ValueType n(x),d(x),one,result;
1026 
1027  one.SetOne();
1028  n.Sub(one);
1029  d.Add(one);
1030  n.Div(d);
1031 
1032  result = ATan0(n);
1033 
1034  n.Set05Pi();
1035  n.exponent.SubOne(); // =pi/4
1036  result.Add(n);
1037 
1038  return result;
1039  }
1040 
1041 
1049  template<class ValueType>
1050  ValueType ATanGreaterThanPlusOne(const ValueType & x)
1051  {
1052  ValueType temp, atan;
1053 
1054  temp.SetOne();
1055 
1056  if( temp.Div(x) )
1057  {
1058  // if there was a carry here that means x is very big
1059  // and atan(1/x) fast converged to 0
1060  atan.SetZero();
1061  }
1062  else
1063  atan = ATan01(temp);
1064 
1065  temp.Set05Pi();
1066  temp.Sub(atan);
1067 
1068  return temp;
1069  }
1070 
1071  } // namespace auxiliaryfunctions
1072 
1073 
1077  template<class ValueType>
1078  ValueType ATan(ValueType x)
1079  {
1080  using namespace auxiliaryfunctions;
1081 
1082  ValueType one, result;
1083  one.SetOne();
1084  bool change_sign = false;
1085 
1086  if( x.IsNan() )
1087  return x;
1088 
1089  // if x is negative we're using the formula:
1090  // atan(-x) = -atan(x)
1091  if( x.IsSign() )
1092  {
1093  change_sign = true;
1094  x.Abs();
1095  }
1096 
1097  if( x.GreaterWithoutSignThan(one) )
1098  result = ATanGreaterThanPlusOne(x);
1099  else
1100  result = ATan01(x);
1101 
1102  if( change_sign )
1103  result.ChangeSign();
1104 
1105  return result;
1106  }
1107 
1108 
1115  template<class ValueType>
1116  ValueType ATg(const ValueType & x)
1117  {
1118  return ATan(x);
1119  }
1120 
1121 
1128  template<class ValueType>
1129  ValueType ACot(const ValueType & x)
1130  {
1131  ValueType result;
1132 
1133  result.Set05Pi();
1134  result.Sub(ATan(x));
1135 
1136  return result;
1137  }
1138 
1139 
1146  template<class ValueType>
1147  ValueType ACtg(const ValueType & x)
1148  {
1149  return ACot(x);
1150  }
1151 
1152 
1153  /*
1154  *
1155  * hyperbolic functions
1156  *
1157  *
1158  */
1159 
1160 
1166  template<class ValueType>
1167  ValueType Sinh(const ValueType & x, ErrorCode * err = 0)
1168  {
1169  if( x.IsNan() )
1170  {
1171  if( err )
1172  *err = err_improper_argument;
1173 
1174  return x; // NaN
1175  }
1176 
1177  ValueType ex, emx;
1178  uint c = 0;
1179 
1180  c += ex.Exp(x);
1181  c += emx.Exp(-x);
1182 
1183  c += ex.Sub(emx);
1184  c += ex.exponent.SubOne();
1185 
1186  if( err )
1187  *err = c ? err_overflow : err_ok;
1188 
1189  return ex;
1190  }
1191 
1192 
1198  template<class ValueType>
1199  ValueType Cosh(const ValueType & x, ErrorCode * err = 0)
1200  {
1201  if( x.IsNan() )
1202  {
1203  if( err )
1204  *err = err_improper_argument;
1205 
1206  return x; // NaN
1207  }
1208 
1209  ValueType ex, emx;
1210  uint c = 0;
1211 
1212  c += ex.Exp(x);
1213  c += emx.Exp(-x);
1214 
1215  c += ex.Add(emx);
1216  c += ex.exponent.SubOne();
1217 
1218  if( err )
1219  *err = c ? err_overflow : err_ok;
1220 
1221  return ex;
1222  }
1223 
1224 
1230  template<class ValueType>
1231  ValueType Tanh(const ValueType & x, ErrorCode * err = 0)
1232  {
1233  if( x.IsNan() )
1234  {
1235  if( err )
1236  *err = err_improper_argument;
1237 
1238  return x; // NaN
1239  }
1240 
1241  ValueType ex, emx, nominator, denominator;
1242  uint c = 0;
1243 
1244  c += ex.Exp(x);
1245  c += emx.Exp(-x);
1246 
1247  nominator = ex;
1248  c += nominator.Sub(emx);
1249  denominator = ex;
1250  c += denominator.Add(emx);
1251 
1252  c += nominator.Div(denominator);
1253 
1254  if( err )
1255  *err = c ? err_overflow : err_ok;
1256 
1257  return nominator;
1258  }
1259 
1260 
1267  template<class ValueType>
1268  ValueType Tgh(const ValueType & x, ErrorCode * err = 0)
1269  {
1270  return Tanh(x, err);
1271  }
1272 
1278  template<class ValueType>
1279  ValueType Coth(const ValueType & x, ErrorCode * err = 0)
1280  {
1281  if( x.IsNan() )
1282  {
1283  if( err )
1284  *err = err_improper_argument;
1285 
1286  return x; // NaN
1287  }
1288 
1289  if( x.IsZero() )
1290  {
1291  if( err )
1292  *err = err_improper_argument;
1293 
1294  return ValueType(); // NaN is set by default
1295  }
1296 
1297  ValueType ex, emx, nominator, denominator;
1298  uint c = 0;
1299 
1300  c += ex.Exp(x);
1301  c += emx.Exp(-x);
1302 
1303  nominator = ex;
1304  c += nominator.Add(emx);
1305  denominator = ex;
1306  c += denominator.Sub(emx);
1307 
1308  c += nominator.Div(denominator);
1309 
1310  if( err )
1311  *err = c ? err_overflow : err_ok;
1312 
1313  return nominator;
1314  }
1315 
1316 
1323  template<class ValueType>
1324  ValueType Ctgh(const ValueType & x, ErrorCode * err = 0)
1325  {
1326  return Coth(x, err);
1327  }
1328 
1329 
1330  /*
1331  *
1332  * inverse hyperbolic functions
1333  *
1334  *
1335  */
1336 
1337 
1343  template<class ValueType>
1344  ValueType ASinh(const ValueType & x, ErrorCode * err = 0)
1345  {
1346  if( x.IsNan() )
1347  {
1348  if( err )
1349  *err = err_improper_argument;
1350 
1351  return x; // NaN
1352  }
1353 
1354  ValueType xx(x), one, result;
1355  uint c = 0;
1356  one.SetOne();
1357 
1358  c += xx.Mul(x);
1359  c += xx.Add(one);
1360  one.exponent.SubOne(); // one=0.5
1361  // xx is >= 1
1362  c += xx.PowFrac(one); // xx=sqrt(xx)
1363  c += xx.Add(x);
1364  c += result.Ln(xx); // xx > 0
1365 
1366  // here can only be a carry
1367  if( err )
1368  *err = c ? err_overflow : err_ok;
1369 
1370  return result;
1371  }
1372 
1373 
1379  template<class ValueType>
1380  ValueType ACosh(const ValueType & x, ErrorCode * err = 0)
1381  {
1382  if( x.IsNan() )
1383  {
1384  if( err )
1385  *err = err_improper_argument;
1386 
1387  return x; // NaN
1388  }
1389 
1390  ValueType xx(x), one, result;
1391  uint c = 0;
1392  one.SetOne();
1393 
1394  if( x < one )
1395  {
1396  if( err )
1397  *err = err_improper_argument;
1398 
1399  return result; // NaN is set by default
1400  }
1401 
1402  c += xx.Mul(x);
1403  c += xx.Sub(one);
1404  // xx is >= 0
1405  // we can't call a PowFrac when the 'x' is zero
1406  // if x is 0 the sqrt(0) is 0
1407  if( !xx.IsZero() )
1408  {
1409  one.exponent.SubOne(); // one=0.5
1410  c += xx.PowFrac(one); // xx=sqrt(xx)
1411  }
1412  c += xx.Add(x);
1413  c += result.Ln(xx); // xx >= 1
1414 
1415  // here can only be a carry
1416  if( err )
1417  *err = c ? err_overflow : err_ok;
1418 
1419  return result;
1420  }
1421 
1422 
1428  template<class ValueType>
1429  ValueType ATanh(const ValueType & x, ErrorCode * err = 0)
1430  {
1431  if( x.IsNan() )
1432  {
1433  if( err )
1434  *err = err_improper_argument;
1435 
1436  return x; // NaN
1437  }
1438 
1439  ValueType nominator(x), denominator, one, result;
1440  uint c = 0;
1441  one.SetOne();
1442 
1443  if( !x.SmallerWithoutSignThan(one) )
1444  {
1445  if( err )
1446  *err = err_improper_argument;
1447 
1448  return result; // NaN is set by default
1449  }
1450 
1451  c += nominator.Add(one);
1452  denominator = one;
1453  c += denominator.Sub(x);
1454  c += nominator.Div(denominator);
1455  c += result.Ln(nominator);
1456  c += result.exponent.SubOne();
1457 
1458  // here can only be a carry
1459  if( err )
1460  *err = c ? err_overflow : err_ok;
1461 
1462  return result;
1463  }
1464 
1465 
1469  template<class ValueType>
1470  ValueType ATgh(const ValueType & x, ErrorCode * err = 0)
1471  {
1472  return ATanh(x, err);
1473  }
1474 
1475 
1481  template<class ValueType>
1482  ValueType ACoth(const ValueType & x, ErrorCode * err = 0)
1483  {
1484  if( x.IsNan() )
1485  {
1486  if( err )
1487  *err = err_improper_argument;
1488 
1489  return x; // NaN
1490  }
1491 
1492  ValueType nominator(x), denominator(x), one, result;
1493  uint c = 0;
1494  one.SetOne();
1495 
1496  if( !x.GreaterWithoutSignThan(one) )
1497  {
1498  if( err )
1499  *err = err_improper_argument;
1500 
1501  return result; // NaN is set by default
1502  }
1503 
1504  c += nominator.Add(one);
1505  c += denominator.Sub(one);
1506  c += nominator.Div(denominator);
1507  c += result.Ln(nominator);
1508  c += result.exponent.SubOne();
1509 
1510  // here can only be a carry
1511  if( err )
1512  *err = c ? err_overflow : err_ok;
1513 
1514  return result;
1515  }
1516 
1517 
1521  template<class ValueType>
1522  ValueType ACtgh(const ValueType & x, ErrorCode * err = 0)
1523  {
1524  return ACoth(x, err);
1525  }
1526 
1527 
1528 
1529 
1530 
1531  /*
1532  *
1533  * functions for converting between degrees, radians and gradians
1534  *
1535  *
1536  */
1537 
1538 
1544  template<class ValueType>
1545  ValueType DegToRad(const ValueType & x, ErrorCode * err = 0)
1546  {
1547  ValueType result, temp;
1548  uint c = 0;
1549 
1550  if( x.IsNan() )
1551  {
1552  if( err )
1553  *err = err_improper_argument;
1554 
1555  return x;
1556  }
1557 
1558  result = x;
1559 
1560  // it is better to make division first and then multiplication
1561  // the result is more accurate especially when x is: 90,180,270 or 360
1562  temp = 180;
1563  c += result.Div(temp);
1564 
1565  temp.SetPi();
1566  c += result.Mul(temp);
1567 
1568  if( err )
1569  *err = c ? err_overflow : err_ok;
1570 
1571  return result;
1572  }
1573 
1574 
1580  template<class ValueType>
1581  ValueType RadToDeg(const ValueType & x, ErrorCode * err = 0)
1582  {
1583  ValueType result, delimiter;
1584  uint c = 0;
1585 
1586  if( x.IsNan() )
1587  {
1588  if( err )
1589  *err = err_improper_argument;
1590 
1591  return x;
1592  }
1593 
1594  result = 180;
1595  c += result.Mul(x);
1596 
1597  delimiter.SetPi();
1598  c += result.Div(delimiter);
1599 
1600  if( err )
1601  *err = c ? err_overflow : err_ok;
1602 
1603  return result;
1604  }
1605 
1606 
1624  template<class ValueType>
1625  ValueType DegToDeg( const ValueType & d, const ValueType & m, const ValueType & s,
1626  ErrorCode * err = 0)
1627  {
1628  ValueType delimiter, multipler;
1629  uint c = 0;
1630 
1631  if( d.IsNan() || m.IsNan() || s.IsNan() || m.IsSign() || s.IsSign() )
1632  {
1633  if( err )
1634  *err = err_improper_argument;
1635 
1636  delimiter.SetZeroNan(); // not needed, only to get rid of GCC warning about an uninitialized variable
1637 
1638  return delimiter;
1639  }
1640 
1641  multipler = 60;
1642  delimiter = 3600;
1643 
1644  c += multipler.Mul(m);
1645  c += multipler.Add(s);
1646  c += multipler.Div(delimiter);
1647 
1648  if( d.IsSign() )
1649  multipler.ChangeSign();
1650 
1651  c += multipler.Add(d);
1652 
1653  if( err )
1654  *err = c ? err_overflow : err_ok;
1655 
1656  return multipler;
1657  }
1658 
1659 
1663  template<class ValueType>
1664  ValueType DegToRad( const ValueType & d, const ValueType & m, const ValueType & s,
1665  ErrorCode * err = 0)
1666  {
1667  ValueType temp_deg = DegToDeg(d,m,s,err);
1668 
1669  if( err && *err!=err_ok )
1670  return temp_deg;
1671 
1672  return DegToRad(temp_deg, err);
1673  }
1674 
1675 
1681  template<class ValueType>
1682  ValueType GradToRad(const ValueType & x, ErrorCode * err = 0)
1683  {
1684  ValueType result, temp;
1685  uint c = 0;
1686 
1687  if( x.IsNan() )
1688  {
1689  if( err )
1690  *err = err_improper_argument;
1691 
1692  return x;
1693  }
1694 
1695  result = x;
1696 
1697  // it is better to make division first and then multiplication
1698  // the result is more accurate especially when x is: 100,200,300 or 400
1699  temp = 200;
1700  c += result.Div(temp);
1701 
1702  temp.SetPi();
1703  c += result.Mul(temp);
1704 
1705  if( err )
1706  *err = c ? err_overflow : err_ok;
1707 
1708  return result;
1709  }
1710 
1711 
1717  template<class ValueType>
1718  ValueType RadToGrad(const ValueType & x, ErrorCode * err = 0)
1719  {
1720  ValueType result, delimiter;
1721  uint c = 0;
1722 
1723  if( x.IsNan() )
1724  {
1725  if( err )
1726  *err = err_improper_argument;
1727 
1728  return x;
1729  }
1730 
1731  result = 200;
1732  c += result.Mul(x);
1733 
1734  delimiter.SetPi();
1735  c += result.Div(delimiter);
1736 
1737  if( err )
1738  *err = c ? err_overflow : err_ok;
1739 
1740  return result;
1741  }
1742 
1743 
1749  template<class ValueType>
1750  ValueType DegToGrad(const ValueType & x, ErrorCode * err = 0)
1751  {
1752  ValueType result, temp;
1753  uint c = 0;
1754 
1755  if( x.IsNan() )
1756  {
1757  if( err )
1758  *err = err_improper_argument;
1759 
1760  return x;
1761  }
1762 
1763  result = x;
1764 
1765  temp = 200;
1766  c += result.Mul(temp);
1767 
1768  temp = 180;
1769  c += result.Div(temp);
1770 
1771  if( err )
1772  *err = c ? err_overflow : err_ok;
1773 
1774  return result;
1775  }
1776 
1777 
1781  template<class ValueType>
1782  ValueType DegToGrad( const ValueType & d, const ValueType & m, const ValueType & s,
1783  ErrorCode * err = 0)
1784  {
1785  ValueType temp_deg = DegToDeg(d,m,s,err);
1786 
1787  if( err && *err!=err_ok )
1788  return temp_deg;
1789 
1790  return DegToGrad(temp_deg, err);
1791  }
1792 
1793 
1799  template<class ValueType>
1800  ValueType GradToDeg(const ValueType & x, ErrorCode * err = 0)
1801  {
1802  ValueType result, temp;
1803  uint c = 0;
1804 
1805  if( x.IsNan() )
1806  {
1807  if( err )
1808  *err = err_improper_argument;
1809 
1810  return x;
1811  }
1812 
1813  result = x;
1814 
1815  temp = 180;
1816  c += result.Mul(temp);
1817 
1818  temp = 200;
1819  c += result.Div(temp);
1820 
1821  if( err )
1822  *err = c ? err_overflow : err_ok;
1823 
1824  return result;
1825  }
1826 
1827 
1828 
1829 
1830  /*
1831  *
1832  * another functions
1833  *
1834  *
1835  */
1836 
1837 
1843  template<class ValueType>
1844  ValueType Sqrt(ValueType x, ErrorCode * err = 0)
1845  {
1846  if( x.IsNan() || x.IsSign() )
1847  {
1848  if( err )
1849  *err = err_improper_argument;
1850 
1851  x.SetNan();
1852 
1853  return x;
1854  }
1855 
1856  uint c = x.Sqrt();
1857 
1858  if( err )
1859  *err = c ? err_overflow : err_ok;
1860 
1861  return x;
1862  }
1863 
1864 
1865 
1866  namespace auxiliaryfunctions
1867  {
1868 
1869  template<class ValueType>
1870  bool RootCheckIndexSign(ValueType & x, const ValueType & index, ErrorCode * err)
1871  {
1872  if( index.IsSign() )
1873  {
1874  // index cannot be negative
1875  if( err )
1876  *err = err_improper_argument;
1877 
1878  x.SetNan();
1879 
1880  return true;
1881  }
1882 
1883  return false;
1884  }
1885 
1886 
1887  template<class ValueType>
1888  bool RootCheckIndexZero(ValueType & x, const ValueType & index, ErrorCode * err)
1889  {
1890  if( index.IsZero() )
1891  {
1892  if( x.IsZero() )
1893  {
1894  // there isn't root(0;0) - we assume it's not defined
1895  if( err )
1896  *err = err_improper_argument;
1897 
1898  x.SetNan();
1899 
1900  return true;
1901  }
1902 
1903  // root(x;0) is 1 (if x!=0)
1904  x.SetOne();
1905 
1906  if( err )
1907  *err = err_ok;
1908 
1909  return true;
1910  }
1911 
1912  return false;
1913  }
1914 
1915 
1916  template<class ValueType>
1917  bool RootCheckIndexOne(const ValueType & index, ErrorCode * err)
1918  {
1919  ValueType one;
1920  one.SetOne();
1921 
1922  if( index == one )
1923  {
1924  //root(x;1) is x
1925  // we do it because if we used the PowFrac function
1926  // we would lose the precision
1927  if( err )
1928  *err = err_ok;
1929 
1930  return true;
1931  }
1932 
1933  return false;
1934  }
1935 
1936 
1937  template<class ValueType>
1938  bool RootCheckIndexTwo(ValueType & x, const ValueType & index, ErrorCode * err)
1939  {
1940  if( index == 2 )
1941  {
1942  x = Sqrt(x, err);
1943 
1944  return true;
1945  }
1946 
1947  return false;
1948  }
1949 
1950 
1951  template<class ValueType>
1952  bool RootCheckIndexFrac(ValueType & x, const ValueType & index, ErrorCode * err)
1953  {
1954  if( !index.IsInteger() )
1955  {
1956  // index must be integer
1957  if( err )
1958  *err = err_improper_argument;
1959 
1960  x.SetNan();
1961 
1962  return true;
1963  }
1964 
1965  return false;
1966  }
1967 
1968 
1969  template<class ValueType>
1970  bool RootCheckXZero(ValueType & x, ErrorCode * err)
1971  {
1972  if( x.IsZero() )
1973  {
1974  // root(0;index) is zero (if index!=0)
1975  // RootCheckIndexZero() must be called beforehand
1976  x.SetZero();
1977 
1978  if( err )
1979  *err = err_ok;
1980 
1981  return true;
1982  }
1983 
1984  return false;
1985  }
1986 
1987 
1988  template<class ValueType>
1989  bool RootCheckIndex(ValueType & x, const ValueType & index, ErrorCode * err, bool * change_sign)
1990  {
1991  *change_sign = false;
1992 
1993  if( index.Mod2() )
1994  {
1995  // index is odd (1,3,5...)
1996  if( x.IsSign() )
1997  {
1998  *change_sign = true;
1999  x.Abs();
2000  }
2001  }
2002  else
2003  {
2004  // index is even
2005  // x cannot be negative
2006  if( x.IsSign() )
2007  {
2008  if( err )
2009  *err = err_improper_argument;
2010 
2011  x.SetNan();
2012 
2013  return true;
2014  }
2015  }
2016 
2017  return false;
2018  }
2019 
2020 
2021  template<class ValueType>
2022  uint RootCorrectInteger(ValueType & old_x, ValueType & x, const ValueType & index)
2023  {
2024  if( !old_x.IsInteger() || x.IsInteger() || !index.exponent.IsSign() )
2025  return 0;
2026 
2027  // old_x is integer,
2028  // x is not integer,
2029  // index is relatively small (index.exponent<0 or index.exponent<=0)
2030  // (because we're using a special powering algorithm Big::PowUInt())
2031 
2032  uint c = 0;
2033 
2034  ValueType temp(x);
2035  c += temp.Round();
2036 
2037  ValueType temp_round(temp);
2038  c += temp.PowUInt(index);
2039 
2040  if( temp == old_x )
2041  x = temp_round;
2042 
2043  return (c==0)? 0 : 1;
2044  }
2045 
2046 
2047 
2048  } // namespace auxiliaryfunctions
2049 
2050 
2051 
2066  template<class ValueType>
2067  ValueType Root(ValueType x, const ValueType & index, ErrorCode * err = 0)
2068  {
2069  using namespace auxiliaryfunctions;
2070 
2071  if( x.IsNan() || index.IsNan() )
2072  {
2073  if( err )
2074  *err = err_improper_argument;
2075 
2076  x.SetNan();
2077 
2078  return x;
2079  }
2080 
2081  if( RootCheckIndexSign(x, index, err) ) return x;
2082  if( RootCheckIndexZero(x, index, err) ) return x;
2083  if( RootCheckIndexOne ( index, err) ) return x;
2084  if( RootCheckIndexTwo (x, index, err) ) return x;
2085  if( RootCheckIndexFrac(x, index, err) ) return x;
2086  if( RootCheckXZero (x, err) ) return x;
2087 
2088  // index integer and index!=0
2089  // x!=0
2090 
2091  ValueType old_x(x);
2092  bool change_sign;
2093 
2094  if( RootCheckIndex(x, index, err, &change_sign ) ) return x;
2095 
2096  ValueType temp;
2097  uint c = 0;
2098 
2099  // we're using the formula: root(x ; n) = exp( ln(x) / n )
2100  c += temp.Ln(x);
2101  c += temp.Div(index);
2102  c += x.Exp(temp);
2103 
2104  if( change_sign )
2105  {
2106  // x is different from zero
2107  x.SetSign();
2108  }
2109 
2110  c += RootCorrectInteger(old_x, x, index);
2111 
2112  if( err )
2113  *err = c ? err_overflow : err_ok;
2114 
2115  return x;
2116  }
2117 
2118 
2119 
2125  template<class ValueType>
2126  ValueType Abs(const ValueType & x)
2127  {
2128  ValueType result( x );
2129  result.Abs();
2130 
2131  return result;
2132  }
2133 
2134 
2141  template<class ValueType>
2142  ValueType Sgn(ValueType x)
2143  {
2144  x.Sgn();
2145 
2146  return x;
2147  }
2148 
2149 
2159  template<class ValueType>
2160  ValueType Mod(ValueType a, const ValueType & b, ErrorCode * err = 0)
2161  {
2162  if( a.IsNan() || b.IsNan() )
2163  {
2164  if( err )
2165  *err = err_improper_argument;
2166 
2167  a.SetNan();
2168 
2169  return a;
2170  }
2171 
2172  uint c = a.Mod(b);
2173 
2174  if( err )
2175  *err = c ? err_overflow : err_ok;
2176 
2177  return a;
2178  }
2179 
2180 
2181 
2182  namespace auxiliaryfunctions
2183  {
2184 
2197  template<class ValueType>
2198  void SetFactorialSequence(std::vector<ValueType> & fact, uint more = 20)
2199  {
2200  if( more == 0 )
2201  more = 1;
2202 
2203  uint start = static_cast<uint>(fact.size());
2204  fact.resize(fact.size() + more);
2205 
2206  if( start == 0 )
2207  {
2208  fact[0] = 1;
2209  ++start;
2210  }
2211 
2212  for(uint i=start ; i<fact.size() ; ++i)
2213  {
2214  fact[i] = fact[i-1];
2215  fact[i].MulInt(i);
2216  }
2217  }
2218 
2219 
2231  template<class ValueType>
2232  ValueType SetBernoulliNumbersSum(CGamma<ValueType> & cgamma, const ValueType & n_, uint m,
2233  const volatile StopCalculating * stop = 0)
2234  {
2235  ValueType k_, temp, temp2, temp3, sum;
2236 
2237  sum.SetZero();
2238 
2239  for(uint k=0 ; k<m ; ++k) // k<m means k<=m-1
2240  {
2241  if( stop && (k & 15)==0 ) // means: k % 16 == 0
2242  if( stop->WasStopSignal() )
2243  return ValueType(); // NaN
2244 
2245  if( k>1 && (k & 1) == 1 ) // for that k the Bernoulli number is zero
2246  continue;
2247 
2248  k_ = k;
2249 
2250  temp = n_; // n_ is equal 2
2251  temp.Pow(k_);
2252  // temp = 2^k
2253 
2254  temp2 = cgamma.fact[m];
2255  temp3 = cgamma.fact[k];
2256  temp3.Mul(cgamma.fact[m-k]);
2257  temp2.Div(temp3);
2258  // temp2 = (m k) = m! / ( k! * (m-k)! )
2259 
2260  temp.Mul(temp2);
2261  temp.Mul(cgamma.bern[k]);
2262 
2263  sum.Add(temp);
2264  // sum += 2^k * (m k) * B(k)
2265 
2266  if( sum.IsNan() )
2267  break;
2268  }
2269 
2270  return sum;
2271  }
2272 
2273 
2282  template<class ValueType>
2283  bool SetBernoulliNumbersMore(CGamma<ValueType> & cgamma, uint start, const volatile StopCalculating * stop = 0)
2284  {
2285  ValueType denominator, temp, temp2, temp3, m_, sum, sum2, n_, k_;
2286 
2287  const uint n = 2;
2288  n_ = n;
2289 
2290  // start is >= 2
2291  for(uint m=start ; m<cgamma.bern.size() ; ++m)
2292  {
2293  if( (m & 1) == 1 )
2294  {
2295  cgamma.bern[m].SetZero();
2296  }
2297  else
2298  {
2299  m_ = m;
2300 
2301  temp = n_; // n_ = 2
2302  temp.Pow(m_);
2303  // temp = 2^m
2304 
2305  denominator.SetOne();
2306  denominator.Sub(temp);
2307  if( denominator.exponent.AddOne() ) // it means: denominator.MulInt(2)
2308  denominator.SetNan();
2309 
2310  // denominator = 2 * (1 - 2^m)
2311 
2312  cgamma.bern[m] = SetBernoulliNumbersSum(cgamma, n_, m, stop);
2313 
2314  if( stop && stop->WasStopSignal() )
2315  {
2316  cgamma.bern.resize(m); // valid numbers are in [0, m-1]
2317  return false;
2318  }
2319 
2320  cgamma.bern[m].Div(denominator);
2321  }
2322  }
2323 
2324  return true;
2325  }
2326 
2327 
2342  template<class ValueType>
2343  bool SetBernoulliNumbers(CGamma<ValueType> & cgamma, uint more = 20, const volatile StopCalculating * stop = 0)
2344  {
2345  if( more == 0 )
2346  more = 1;
2347 
2348  uint start = static_cast<uint>(cgamma.bern.size());
2349  cgamma.bern.resize(cgamma.bern.size() + more);
2350 
2351  if( start == 0 )
2352  {
2353  cgamma.bern[0].SetOne();
2354  ++start;
2355  }
2356 
2357  if( cgamma.bern.size() == 1 )
2358  return true;
2359 
2360  if( start == 1 )
2361  {
2362  cgamma.bern[1].Set05();
2363  cgamma.bern[1].ChangeSign();
2364  ++start;
2365  }
2366 
2367  // we should have sufficient factorials in cgamma.fact
2368  if( cgamma.fact.size() < cgamma.bern.size() )
2369  SetFactorialSequence(cgamma.fact, static_cast<uint>(cgamma.bern.size() - cgamma.fact.size()));
2370 
2371 
2372  return SetBernoulliNumbersMore(cgamma, start, stop);
2373  }
2374 
2375 
2384  template<class ValueType>
2385  ValueType GammaFactorialHighSum(const ValueType & n, CGamma<ValueType> & cgamma, ErrorCode & err,
2386  const volatile StopCalculating * stop)
2387  {
2388  ValueType temp, temp2, denominator, sum, oldsum;
2389 
2390  sum.SetZero();
2391 
2392  for(uint m=2 ; m<TTMATH_ARITHMETIC_MAX_LOOP ; m+=2)
2393  {
2394  if( stop && (m & 3)==0 ) // (m & 3)==0 means: (m % 4)==0
2395  if( stop->WasStopSignal() )
2396  {
2397  err = err_interrupt;
2398  return ValueType(); // NaN
2399  }
2400 
2401  temp = (m-1);
2402  denominator = n;
2403  denominator.Pow(temp);
2404  // denominator = n ^ (m-1)
2405 
2406  temp = m;
2407  temp2 = temp;
2408  temp.Mul(temp2);
2409  temp.Sub(temp2);
2410  // temp = m^2 - m
2411 
2412  denominator.Mul(temp);
2413  // denominator = (m^2 - m) * n ^ (m-1)
2414 
2415  if( m >= cgamma.bern.size() )
2416  {
2417  if( !SetBernoulliNumbers(cgamma, m - cgamma.bern.size() + 1 + 3, stop) ) // 3 more than needed
2418  {
2419  // there was the stop signal
2420  err = err_interrupt;
2421  return ValueType(); // NaN
2422  }
2423  }
2424 
2425  temp = cgamma.bern[m];
2426  temp.Div(denominator);
2427 
2428  oldsum = sum;
2429  sum.Add(temp);
2430 
2431  if( sum.IsNan() || oldsum==sum )
2432  break;
2433  }
2434 
2435  return sum;
2436  }
2437 
2438 
2447  template<class ValueType>
2448  ValueType GammaFactorialHigh(const ValueType & n, CGamma<ValueType> & cgamma, ErrorCode & err,
2449  const volatile StopCalculating * stop)
2450  {
2451  ValueType temp, temp2, temp3, denominator, sum;
2452 
2453  temp.Set2Pi();
2454  temp.Mul(n);
2455  temp2 = Sqrt(temp);
2456  // temp2 = sqrt(2*pi*n)
2457 
2458  temp = n;
2459  temp3.SetE();
2460  temp.Div(temp3);
2461  temp.Pow(n);
2462  // temp = (n/e)^n
2463 
2464  sum = GammaFactorialHighSum(n, cgamma, err, stop);
2465  temp3.Exp(sum);
2466  // temp3 = exp(sum)
2467 
2468  temp.Mul(temp2);
2469  temp.Mul(temp3);
2470 
2471  return temp;
2472  }
2473 
2474 
2480  template<class ValueType>
2481  ValueType GammaPlusHigh(ValueType n, CGamma<ValueType> & cgamma, ErrorCode & err, const volatile StopCalculating * stop)
2482  {
2483  ValueType one;
2484 
2485  one.SetOne();
2486  n.Sub(one);
2487 
2488  return GammaFactorialHigh(n, cgamma, err, stop);
2489  }
2490 
2491 
2499  template<class ValueType>
2501  {
2502  TTMATH_ASSERT( n > 0 )
2503 
2504  if( n - 1 < static_cast<uint>(cgamma.fact.size()) )
2505  return cgamma.fact[n - 1];
2506 
2507  ValueType res;
2508  uint start = 2;
2509 
2510  if( cgamma.fact.size() < 2 )
2511  {
2512  res.SetOne();
2513  }
2514  else
2515  {
2516  start = static_cast<uint>(cgamma.fact.size());
2517  res = cgamma.fact[start-1];
2518  }
2519 
2520  for(uint i=start ; i<n ; ++i)
2521  res.MulInt(i);
2522 
2523  return res;
2524  }
2525 
2526 
2532  template<class ValueType>
2533  ValueType GammaPlusLowInteger(const ValueType & n, CGamma<ValueType> & cgamma)
2534  {
2535  sint n_;
2536 
2537  n.ToInt(n_);
2538 
2539  return GammaPlusLowIntegerInt(n_, cgamma);
2540  }
2541 
2542 
2554  template<class ValueType>
2555  ValueType GammaPlusLow(ValueType n, CGamma<ValueType> & cgamma, ErrorCode & err, const volatile StopCalculating * stop)
2556  {
2557  ValueType one, denominator, temp, boundary;
2558 
2559  if( n.IsInteger() )
2560  return GammaPlusLowInteger(n, cgamma);
2561 
2562  one.SetOne();
2563  denominator = n;
2564  boundary = TTMATH_GAMMA_BOUNDARY;
2565 
2566  while( n < boundary )
2567  {
2568  n.Add(one);
2569  denominator.Mul(n);
2570  }
2571 
2572  n.Add(one);
2573 
2574  // now n is sufficient big
2575  temp = GammaPlusHigh(n, cgamma, err, stop);
2576  temp.Div(denominator);
2577 
2578  return temp;
2579  }
2580 
2581 
2585  template<class ValueType>
2586  ValueType GammaPlus(const ValueType & n, CGamma<ValueType> & cgamma, ErrorCode & err, const volatile StopCalculating * stop)
2587  {
2588  if( n > TTMATH_GAMMA_BOUNDARY )
2589  return GammaPlusHigh(n, cgamma, err, stop);
2590 
2591  return GammaPlusLow(n, cgamma, err, stop);
2592  }
2593 
2594 
2604  template<class ValueType>
2605  ValueType GammaMinus(const ValueType & n, CGamma<ValueType> & cgamma, ErrorCode & err, const volatile StopCalculating * stop)
2606  {
2607  ValueType pi, denominator, temp, temp2;
2608 
2609  if( n.IsInteger() )
2610  {
2611  // gamma function is not defined when n is negative and integer
2612  err = err_improper_argument;
2613  return temp; // NaN
2614  }
2615 
2616  pi.SetPi();
2617 
2618  temp = pi;
2619  temp.Mul(n);
2620  temp2 = Sin(temp);
2621  // temp2 = sin(pi * n)
2622 
2623  temp.SetOne();
2624  temp.Sub(n);
2625  temp = GammaPlus(temp, cgamma, err, stop);
2626  // temp = gamma(1 - n)
2627 
2628  temp.Mul(temp2);
2629  pi.Div(temp);
2630 
2631  return pi;
2632  }
2633 
2634  } // namespace auxiliaryfunctions
2635 
2636 
2637 
2653  template<class ValueType>
2654  ValueType Gamma(const ValueType & n, CGamma<ValueType> & cgamma, ErrorCode * err = 0,
2655  const volatile StopCalculating * stop = 0)
2656  {
2657  using namespace auxiliaryfunctions;
2658 
2659  ValueType result;
2660  ErrorCode err_tmp;
2661 
2662  if( n.IsNan() )
2663  {
2664  if( err )
2665  *err = err_improper_argument;
2666 
2667  return n;
2668  }
2669 
2670  if( cgamma.history.Get(n, result, err_tmp) )
2671  {
2672  if( err )
2673  *err = err_tmp;
2674 
2675  return result;
2676  }
2677 
2678  err_tmp = err_ok;
2679 
2680  if( n.IsSign() )
2681  {
2682  result = GammaMinus(n, cgamma, err_tmp, stop);
2683  }
2684  else
2685  if( n.IsZero() )
2686  {
2687  err_tmp = err_improper_argument;
2688  result.SetNan();
2689  }
2690  else
2691  {
2692  result = GammaPlus(n, cgamma, err_tmp, stop);
2693  }
2694 
2695  if( result.IsNan() && err_tmp==err_ok )
2696  err_tmp = err_overflow;
2697 
2698  if( err )
2699  *err = err_tmp;
2700 
2701  if( stop && !stop->WasStopSignal() )
2702  cgamma.history.Add(n, result, err_tmp);
2703 
2704  return result;
2705  }
2706 
2707 
2713  template<class ValueType>
2714  ValueType Gamma(const ValueType & n, ErrorCode * err = 0)
2715  {
2716  // warning: this static object is not thread safe
2717  static CGamma<ValueType> cgamma;
2718 
2719  return Gamma(n, cgamma, err);
2720  }
2721 
2722 
2723 
2724  namespace auxiliaryfunctions
2725  {
2726 
2733  template<class ValueType>
2734  ValueType Factorial2(ValueType x,
2735  CGamma<ValueType> * cgamma = 0,
2736  ErrorCode * err = 0,
2737  const volatile StopCalculating * stop = 0)
2738  {
2739  ValueType result, one;
2740 
2741  if( x.IsNan() || x.IsSign() || !x.IsInteger() )
2742  {
2743  if( err )
2744  *err = err_improper_argument;
2745 
2746  x.SetNan();
2747 
2748  return x;
2749  }
2750 
2751  one.SetOne();
2752  x.Add(one);
2753 
2754  if( cgamma )
2755  return Gamma(x, *cgamma, err, stop);
2756 
2757  return Gamma(x, err);
2758  }
2759 
2760  } // namespace auxiliaryfunctions
2761 
2762 
2763 
2781  template<class ValueType>
2782  ValueType Factorial(const ValueType & x, CGamma<ValueType> & cgamma, ErrorCode * err = 0,
2783  const volatile StopCalculating * stop = 0)
2784  {
2785  return auxiliaryfunctions::Factorial2(x, &cgamma, err, stop);
2786  }
2787 
2788 
2796  template<class ValueType>
2797  ValueType Factorial(const ValueType & x, ErrorCode * err = 0)
2798  {
2799  return auxiliaryfunctions::Factorial2(x, (CGamma<ValueType>*)0, err, 0);
2800  }
2801 
2802 
2812  template<class ValueType>
2814  {
2815  ValueType x = TTMATH_GAMMA_BOUNDARY + 1;
2816 
2817  // history.Remove(x) removes only one object
2818  // we must be sure that there are not others objects with the key 'x'
2819  while( history.Remove(x) )
2820  {
2821  }
2822 
2823  // the simplest way to initialize is to call the Gamma function with (TTMATH_GAMMA_BOUNDARY + 1)
2824  // when x is larger then fewer coefficients we need
2825  Gamma(x, *this);
2826  }
2827 
2828 
2829 
2830 } // namespace
2831 
2832 
2837 #include "ttmathparser.h"
2838 
2839 // Dec is not finished yet
2840 //#include "ttmathdec.h"
2841 
2842 
2843 
2844 #ifdef _MSC_VER
2845 //warning C4127: conditional expression is constant
2846 #pragma warning( default: 4127 )
2847 //warning C4702: unreachable code
2848 #pragma warning( default: 4702 )
2849 //warning C4800: forcing value to bool 'true' or 'false' (performance warning)
2850 #pragma warning( default: 4800 )
2851 #endif
2852 
2853 #endif
ValueType Sgn(ValueType x)
Definition: ttmath.h:2142
ValueType ACot(const ValueType &x)
Definition: ttmath.h:1129
Definition: ttmathtypes.h:366
ValueType GammaMinus(const ValueType &n, CGamma< ValueType > &cgamma, ErrorCode &err, const volatile StopCalculating *stop)
Definition: ttmath.h:2605
ValueType ACoth(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:1482
bool SetBernoulliNumbers(CGamma< ValueType > &cgamma, uint more=20, const volatile StopCalculating *stop=0)
Definition: ttmath.h:2343
ValueType Tan(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:612
ValueType Log(const ValueType &x, const ValueType &base, ErrorCode *err=0)
Definition: ttmath.h:274
ValueType Sin(ValueType x, ErrorCode *err=0)
Definition: ttmath.h:518
ValueType ATg(const ValueType &x)
Definition: ttmath.h:1116
#define TTMATH_ARITHMETIC_MAX_LOOP
Definition: ttmathtypes.h:289
Definition: ttmathtypes.h:364
ValueType ATan01(const ValueType &x)
Definition: ttmath.h:979
NetworKit::index index
Definition: BloomFilter.h:16
A mathematical parser.
bool RootCheckIndexTwo(ValueType &x, const ValueType &index, ErrorCode *err)
Definition: ttmath.h:1938
ValueType GammaPlusHigh(ValueType n, CGamma< ValueType > &cgamma, ErrorCode &err, const volatile StopCalculating *stop)
Definition: ttmath.h:2481
ValueType Exp(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:319
ValueType ACosh(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:1380
ValueType RadToDeg(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:1581
ValueType ASin(ValueType x, ErrorCode *err=0)
Definition: ttmath.h:848
virtual bool WasStopSignal() const volatile
Definition: ttmathtypes.h:523
std::vector< ValueType > bern
Definition: ttmathobjects.h:770
void InitAll()
Definition: ttmath.h:2813
ValueType ASinh(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:1344
bool RootCheckIndexOne(const ValueType &index, ErrorCode *err)
Definition: ttmath.h:1917
ValueType SkipFraction(const ValueType &x)
Definition: ttmath.h:88
ValueType Tg(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:640
ValueType Cos(ValueType x, ErrorCode *err=0)
Definition: ttmath.h:574
ValueType ATanh(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:1429
bool RootCheckIndexFrac(ValueType &x, const ValueType &index, ErrorCode *err)
Definition: ttmath.h:1952
ValueType Ln(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:233
ValueType RadToGrad(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:1718
ValueType Abs(const ValueType &x)
Definition: ttmath.h:2126
Definition: ttmathtypes.h:359
History< ValueType > history
Definition: ttmathobjects.h:779
ValueType Sin0pi05(const ValueType &x)
Definition: ttmath.h:423
ValueType ACtgh(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:1522
ValueType GradToDeg(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:1800
ErrorCode
Definition: ttmathtypes.h:349
uint RootCorrectInteger(ValueType &old_x, ValueType &x, const ValueType &index)
Definition: ttmath.h:2022
signed int sint
Definition: ttmathtypes.h:164
ValueType ATan0(const ValueType &x)
Definition: ttmath.h:928
ValueType Tgh(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:1268
ValueType GammaPlusLowInteger(const ValueType &n, CGamma< ValueType > &cgamma)
Definition: ttmath.h:2533
ValueType Root(ValueType x, const ValueType &index, ErrorCode *err=0)
Definition: ttmath.h:2067
ValueType DegToRad(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:1545
ValueType GradToRad(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:1682
ValueType GammaPlusLow(ValueType n, CGamma< ValueType > &cgamma, ErrorCode &err, const volatile StopCalculating *stop)
Definition: ttmath.h:2555
ValueType Mod(ValueType a, const ValueType &b, ErrorCode *err=0)
Definition: ttmath.h:2160
ValueType ATanGreaterThanPlusOne(const ValueType &x)
Definition: ttmath.h:1050
ValueType Ctg(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:682
ValueType ATan(ValueType x)
Definition: ttmath.h:1078
Definition: ttmathtypes.h:358
#define TTMATH_ASSERT(expression)
Definition: ttmathtypes.h:660
ValueType Coth(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:1279
ValueType Factorial(const ValueType &x, CGamma< ValueType > &cgamma, ErrorCode *err=0, const volatile StopCalculating *stop=0)
Definition: ttmath.h:2782
ValueType GammaPlusLowIntegerInt(uint n, CGamma< ValueType > &cgamma)
Definition: ttmath.h:2500
ValueType Floor(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:185
ValueType Tanh(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:1231
ValueType Sinh(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:1167
Definition: ttmathobjects.h:739
ValueType Cot(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:654
ValueType ASin_1(const ValueType &x)
Definition: ttmath.h:772
ValueType DegToGrad(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:1750
ValueType Gamma(const ValueType &n, CGamma< ValueType > &cgamma, ErrorCode *err=0, const volatile StopCalculating *stop=0)
Definition: ttmath.h:2654
std::vector< ValueType > fact
Definition: ttmathobjects.h:752
Mathematic functions.
ValueType Cosh(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:1199
#define TTMATH_GAMMA_BOUNDARY
Definition: ttmathtypes.h:317
ValueType Factorial2(ValueType x, CGamma< ValueType > *cgamma=0, ErrorCode *err=0, const volatile StopCalculating *stop=0)
Definition: ttmath.h:2734
ValueType Round(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:105
ValueType ACtg(const ValueType &x)
Definition: ttmath.h:1147
ValueType DegToDeg(const ValueType &d, const ValueType &m, const ValueType &s, ErrorCode *err=0)
Definition: ttmath.h:1625
Definition: ttmathtypes.h:520
ValueType ATgh(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:1470
A Class for representing floating point numbers.
void SetFactorialSequence(std::vector< ValueType > &fact, uint more=20)
Definition: ttmath.h:2198
ValueType Ctgh(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:1324
unsigned int uint
Definition: ttmathtypes.h:163
ValueType GammaFactorialHigh(const ValueType &n, CGamma< ValueType > &cgamma, ErrorCode &err, const volatile StopCalculating *stop)
Definition: ttmath.h:2448
ValueType GammaPlus(const ValueType &n, CGamma< ValueType > &cgamma, ErrorCode &err, const volatile StopCalculating *stop)
Definition: ttmath.h:2586
ValueType Ceil(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:138
ValueType ASin_0(const ValueType &x)
Definition: ttmath.h:708
ValueType GammaFactorialHighSum(const ValueType &n, CGamma< ValueType > &cgamma, ErrorCode &err, const volatile StopCalculating *stop)
Definition: ttmath.h:2385
ValueType Sqrt(ValueType x, ErrorCode *err=0)
Definition: ttmath.h:1844
bool RootCheckXZero(ValueType &x, ErrorCode *err)
Definition: ttmath.h:1970
Definition: ttmathtypes.h:351
ValueType ACos(const ValueType &x, ErrorCode *err=0)
Definition: ttmath.h:903
bool RootCheckIndexZero(ValueType &x, const ValueType &index, ErrorCode *err)
Definition: ttmath.h:1888
uint PrepareSin(ValueType &x, bool &change_sign)
Definition: ttmath.h:358
bool RootCheckIndexSign(ValueType &x, const ValueType &index, ErrorCode *err)
Definition: ttmath.h:1870
bool RootCheckIndex(ValueType &x, const ValueType &index, ErrorCode *err, bool *change_sign)
Definition: ttmath.h:1989
bool SetBernoulliNumbersMore(CGamma< ValueType > &cgamma, uint start, const volatile StopCalculating *stop=0)
Definition: ttmath.h:2283
ValueType SetBernoulliNumbersSum(CGamma< ValueType > &cgamma, const ValueType &n_, uint m, const volatile StopCalculating *stop=0)
Definition: ttmath.h:2232