* Ctors of class numeric are not explicit any more. All built-in callers for
authorRichard Kreckel <Richard.Kreckel@uni-mainz.de>
Tue, 26 Jun 2001 21:19:34 +0000 (21:19 +0000)
committerRichard Kreckel <Richard.Kreckel@uni-mainz.de>
Tue, 26 Jun 2001 21:19:34 +0000 (21:19 +0000)
  pseudofunctions are now templated and default to ex arguments which relaxes
  the need for explicit ctors.

13 files changed:
NEWS
doc/tutorial/ginac.texi
ginac/expairseq.cpp
ginac/function.pl
ginac/inifcns.cpp
ginac/inifcns.h
ginac/inifcns_gamma.cpp
ginac/inifcns_trans.cpp
ginac/mul.cpp
ginac/ncmul.cpp
ginac/numeric.h
ginac/power.cpp
ginac/power.h

diff --git a/NEWS b/NEWS
index de5777b..bc8ab41 100644 (file)
--- a/NEWS
+++ b/NEWS
@@ -1,6 +1,9 @@
 This file records noteworthy changes.
 
 0.9.1 (<date>)
+* Ctors of class numeric are not explicit any more.  All built-in callers for
+  pseudofunctions are now templated and default to ex arguments which relaxes
+  the need for explicit ctors.
 * New functions/methods:
    - find()
    - remove_first(), remove_last(), sort() and unique() for lists
index 6ed6902..5934347 100644 (file)
@@ -856,10 +856,10 @@ using namespace GiNaC;
 
 int main()
 @{
-    numeric two(2);                       // exact integer 2
+    numeric two = 2;                      // exact integer 2
     numeric r(2,3);                       // exact fraction 2/3
     numeric e(2.71828);                   // floating point number
-    numeric p("3.1415926535897932385");   // floating point number
+    numeric p = "3.14159265358979323846"; // constructor from string
     // Trott's constant in scientific notation:
     numeric trott("1.0841015122311136151E-2");
     
@@ -867,18 +867,6 @@ int main()
 @}
 @end example
 
-Note that all those constructors are @emph{explicit} which means you are
-not allowed to write @code{numeric two=2;}.  This is because the basic
-objects to be handled by GiNaC are the expressions @code{ex} and we want
-to keep things simple and wish objects like @code{pow(x,2)} to be
-handled the same way as @code{pow(x,a)}, which means that we need to
-allow a general @code{ex} as base and exponent.  Therefore there is an
-implicit constructor from C-integers directly to expressions handling
-numerics at work in most of our examples.  This design really becomes
-convenient when one declares own functions having more than one
-parameter but it forbids using implicit constructors because that would
-lead to compile-time ambiguities.
-
 It may be tempting to construct numbers writing @code{numeric r(3/2)}.
 This would, however, call C's built-in operator @code{/} for integers
 first and result in a numeric holding a plain integer 1.  @strong{Never
@@ -1196,11 +1184,12 @@ There are quite a number of useful functions hard-wired into GiNaC.  For
 instance, all trigonometric and hyperbolic functions are implemented
 (@xref{Built-in Functions}, for a complete list).
 
-These functions are all objects of class @code{function}.  They accept
-one or more expressions as arguments and return one expression.  If the
-arguments are not numerical, the evaluation of the function may be
-halted, as it does in the next example, showing how a function returns
-itself twice and finally an expression that may be really useful:
+These functions (better called @emph{pseudofunctions}) are all objects
+of class @code{function}.  They accept one or more expressions as
+arguments and return one expression.  If the arguments are not
+numerical, the evaluation of the function may be halted, as it does in
+the next example, showing how a function returns itself twice and
+finally an expression that may be really useful:
 
 @cindex Gamma function
 @cindex @code{subs()}
@@ -1223,6 +1212,18 @@ Besides evaluation most of these functions allow differentiation, series
 expansion and so on.  Read the next chapter in order to learn more about
 this.
 
+It must be noted that these pseudofunctions are created by inline
+functions, where the argument list is templated.  This means that
+whenever you call @code{GiNaC::sin(1)} it is equivalent to
+@code{sin(ex(1))} and will therefore not result in a floating point
+numeber.  Unless of course the function prototype is explicitly
+overridden -- which is the case for arguments of type @code{numeric}
+(not wrapped inside an @code{ex}).  Hence, in order to obtain a floating
+point number of class @code{numeric} you should call
+@code{sin(numeric(1))}.  This is almost the same as calling
+@code{sin(1).evalf()} except that the latter will return a numeric
+wrapped inside an @code{ex}.
+
 
 @node Relations, Matrices, Mathematical functions, Basic Concepts
 @c    node-name, next, previous, up
index f49fefc..87b19e1 100644 (file)
@@ -1082,9 +1082,9 @@ void expairseq::combine_same_terms_sorted_seq(void)
                // possible from then on the sequence has changed and must be compacted
                bool must_copy = false;
                while (itin2!=last) {
-                       if ((*itin1).rest.compare(itin2->rest)==0) {
-                               (*itin1).coeff = ex_to<numeric>(itin1->coeff).
-                                                add_dyn(ex_to<numeric>(itin2->coeff));
+                       if (itin1->rest.compare(itin2->rest)==0) {
+                               itin1->coeff = ex_to<numeric>(itin1->coeff).
+                                              add_dyn(ex_to<numeric>(itin2->coeff));
                                if (expair_needs_further_processing(itin1))
                                        needs_further_processing = true;
                                must_copy = true;
index 565b0f9..de7433e 100755 (executable)
@@ -15,65 +15,68 @@ sub generate_seq {
 }
 
 sub generate_from_to {
-       my ($template,$seq_template1,$seq_template2,$from,$to)=@_;
+       my ($template,$seq_template1,$seq_template2,$seq_template3,$from,$to)=@_;
        my ($res,$N,$SEQ);
 
        $res='';
        for ($N=$from; $N<=$to; $N++) {
                $SEQ1=generate_seq($seq_template1,$N);
                $SEQ2=generate_seq($seq_template2,$N);
+               $SEQ3=generate_seq($seq_template3,$N);
                $res .= eval('"' . $template . '"');
                $SEQ1=''; # to avoid main::SEQ1 used only once warning
                $SEQ2=''; # same as above
+               $SEQ3=''; # same as above
        }
        return $res;
 }
 
 sub generate {
-       my ($template,$seq_template1,$seq_template2)=@_;
-       return generate_from_to($template,$seq_template1,$seq_template2,1,$maxargs);
+       my ($template,$seq_template1,$seq_template2,$seq_template3)=@_;
+       return generate_from_to($template,$seq_template1,$seq_template2,$seq_template3,1,$maxargs);
 }
 
-$declare_function_macro = generate_from_to(
-       <<'END_OF_DECLARE_FUNCTION_MACRO','const GiNaC::ex & p${N}','p${N}',1,$maxargs);
+$declare_function_macro = generate(
+       <<'END_OF_DECLARE_FUNCTION_MACRO','typename T${N}','const T${N} & p${N}','GiNaC::ex(p${N})');
 #define DECLARE_FUNCTION_${N}P(NAME) \\
 extern const unsigned function_index_##NAME; \\
-inline const GiNaC::function NAME(${SEQ1}) { \\
-       return GiNaC::function(function_index_##NAME, ${SEQ2}); \\
+template<${SEQ1}> \\
+inline const GiNaC::function NAME(${SEQ2}) { \\
+       return GiNaC::function(function_index_##NAME, ${SEQ3}); \\
 }
 
 END_OF_DECLARE_FUNCTION_MACRO
 
 $typedef_eval_funcp=generate(
 'typedef ex (* eval_funcp_${N})(${SEQ1});'."\n",
-'const ex &','');
+'const ex &','','');
 
 $typedef_evalf_funcp=generate(
 'typedef ex (* evalf_funcp_${N})(${SEQ1});'."\n",
-'const ex &','');
+'const ex &','','');
 
 $typedef_derivative_funcp=generate(
 'typedef ex (* derivative_funcp_${N})(${SEQ1}, unsigned);'."\n",
-'const ex &','');
+'const ex &','','');
 
 $typedef_series_funcp=generate(
 'typedef ex (* series_funcp_${N})(${SEQ1}, const relational &, int, unsigned);'."\n",
-'const ex &','');
+'const ex &','','');
 
-$eval_func_interface=generate('    function_options & eval_func(eval_funcp_${N} e);'."\n",'','');
+$eval_func_interface=generate('    function_options & eval_func(eval_funcp_${N} e);'."\n",'','','');
 
-$evalf_func_interface=generate('    function_options & evalf_func(evalf_funcp_${N} ef);'."\n",'','');
+$evalf_func_interface=generate('    function_options & evalf_func(evalf_funcp_${N} ef);'."\n",'','','');
 
-$derivative_func_interface=generate('    function_options & derivative_func(derivative_funcp_${N} d);'."\n",'','');
+$derivative_func_interface=generate('    function_options & derivative_func(derivative_funcp_${N} d);'."\n",'','','');
 
-$series_func_interface=generate('    function_options & series_func(series_funcp_${N} s);'."\n",'','');
+$series_func_interface=generate('    function_options & series_func(series_funcp_${N} s);'."\n",'','','');
 
 $constructors_interface=generate(
 '    function(unsigned ser, ${SEQ1});'."\n",
-'const ex & param${N}','');
+'const ex & param${N}','','');
 
 $constructors_implementation=generate(
-       <<'END_OF_CONSTRUCTORS_IMPLEMENTATION','const ex & param${N}','param${N}');
+       <<'END_OF_CONSTRUCTORS_IMPLEMENTATION','const ex & param${N}','param${N}','');
 function::function(unsigned ser, ${SEQ1})
        : exprseq(${SEQ2}), serial(ser)
 {
@@ -83,26 +86,26 @@ function::function(unsigned ser, ${SEQ1})
 END_OF_CONSTRUCTORS_IMPLEMENTATION
 
 $eval_switch_statement=generate(
-       <<'END_OF_EVAL_SWITCH_STATEMENT','seq[${N}-1]','');
+       <<'END_OF_EVAL_SWITCH_STATEMENT','seq[${N}-1]','','');
        case ${N}:
                eval_result = ((eval_funcp_${N})(registered_functions()[serial].eval_f))(${SEQ1});
                break;
 END_OF_EVAL_SWITCH_STATEMENT
 
 $evalf_switch_statement=generate(
-       <<'END_OF_EVALF_SWITCH_STATEMENT','eseq[${N}-1]','');
+       <<'END_OF_EVALF_SWITCH_STATEMENT','eseq[${N}-1]','','');
        case ${N}:
                return ((evalf_funcp_${N})(registered_functions()[serial].evalf_f))(${SEQ1});
 END_OF_EVALF_SWITCH_STATEMENT
 
 $diff_switch_statement=generate(
-       <<'END_OF_DIFF_SWITCH_STATEMENT','seq[${N}-1]','');
+       <<'END_OF_DIFF_SWITCH_STATEMENT','seq[${N}-1]','','');
        case ${N}:
                return ((derivative_funcp_${N})(registered_functions()[serial].derivative_f))(${SEQ1},diff_param);
 END_OF_DIFF_SWITCH_STATEMENT
 
 $series_switch_statement=generate(
-       <<'END_OF_SERIES_SWITCH_STATEMENT','seq[${N}-1]','');
+       <<'END_OF_SERIES_SWITCH_STATEMENT','seq[${N}-1]','','');
        case ${N}:
                try {
                        res = ((series_funcp_${N})(registered_functions()[serial].series_f))(${SEQ1},r,order,options);
@@ -113,7 +116,7 @@ $series_switch_statement=generate(
 END_OF_SERIES_SWITCH_STATEMENT
 
 $eval_func_implementation=generate(
-       <<'END_OF_EVAL_FUNC_IMPLEMENTATION','','');
+       <<'END_OF_EVAL_FUNC_IMPLEMENTATION','','','');
 function_options & function_options::eval_func(eval_funcp_${N} e)
 {
        test_and_set_nparams(${N});
@@ -123,7 +126,7 @@ function_options & function_options::eval_func(eval_funcp_${N} e)
 END_OF_EVAL_FUNC_IMPLEMENTATION
 
 $evalf_func_implementation=generate(
-       <<'END_OF_EVALF_FUNC_IMPLEMENTATION','','');
+       <<'END_OF_EVALF_FUNC_IMPLEMENTATION','','','');
 function_options & function_options::evalf_func(evalf_funcp_${N} ef)
 {
        test_and_set_nparams(${N});
@@ -133,7 +136,7 @@ function_options & function_options::evalf_func(evalf_funcp_${N} ef)
 END_OF_EVALF_FUNC_IMPLEMENTATION
 
 $derivative_func_implementation=generate(
-       <<'END_OF_DERIVATIVE_FUNC_IMPLEMENTATION','','');
+       <<'END_OF_DERIVATIVE_FUNC_IMPLEMENTATION','','','');
 function_options & function_options::derivative_func(derivative_funcp_${N} d)
 {
        test_and_set_nparams(${N});
@@ -143,7 +146,7 @@ function_options & function_options::derivative_func(derivative_funcp_${N} d)
 END_OF_DERIVATIVE_FUNC_IMPLEMENTATION
 
 $series_func_implementation=generate(
-       <<'END_OF_SERIES_FUNC_IMPLEMENTATION','','');
+       <<'END_OF_SERIES_FUNC_IMPLEMENTATION','','','');
 function_options & function_options::series_func(series_funcp_${N} s)
 {
        test_and_set_nparams(${N});
@@ -984,14 +987,14 @@ unsigned function::register_new(function_options const & opt)
                // usually executed before main(), so the exception could not
                // caught anyhow
                std::cerr << "WARNING: function name " << opt.name
-                                 << " already in use!" << std::endl;
+                         << " already in use!" << std::endl;
        }
        registered_functions().push_back(opt);
        if (opt.use_remember) {
                remember_table::remember_tables().
                        push_back(remember_table(opt.remember_size,
-                                                                        opt.remember_assoc_size,
-                                                                        opt.remember_strategy));
+                                                opt.remember_assoc_size,
+                                                opt.remember_strategy));
        } else {
                remember_table::remember_tables().push_back(remember_table());
        }
index 7ac2153..8c93b6f 100644 (file)
@@ -274,7 +274,7 @@ static ex Li2_series(const ex &x, const relational &rel, int order, unsigned opt
                        // method:
                        // construct series manually in a dummy symbol s
                        const symbol s;
-                       ex ser = zeta(2);
+                       ex ser = zeta(_ex2());
                        // manually construct the primitive expansion
                        for (int i=1; i<order; ++i)
                                ser += pow(1-s,i) * (numeric(1,i)*(I*Pi+log(s-1)) - numeric(1,i*i));
index fd0d9e4..8acc2bf 100644 (file)
@@ -91,13 +91,15 @@ DECLARE_FUNCTION_1P(Li3)
 // overloading at work: we cannot use the macros here
 /** Riemann's Zeta-function. */
 extern const unsigned function_index_zeta1;
-inline function zeta(const ex & p1) {
-       return function(function_index_zeta1, p1);
+template<typename T1>
+inline function zeta(const T1 & p1) {
+       return function(function_index_zeta1, ex(p1));
 }
 /** Derivatives of Riemann's Zeta-function. */
 extern const unsigned function_index_zeta2;
-inline function zeta(const ex & p1, const ex & p2) {
-       return function(function_index_zeta2, p1, p2);
+template<typename T1, typename T2>
+inline function zeta(const T1 & p1, const T2 & p2) {
+       return function(function_index_zeta2, ex(p1), ex(p2));
 }
 
 /** Gamma-function. */
@@ -110,13 +112,15 @@ DECLARE_FUNCTION_2P(beta)
 // overloading at work: we cannot use the macros here
 /** Psi-function (aka digamma-function). */
 extern const unsigned function_index_psi1;
-inline function psi(const ex & p1) {
-       return function(function_index_psi1, p1);
+template<typename T1>
+inline function psi(const T1 & p1) {
+       return function(function_index_psi1, ex(p1));
 }
 /** Derivatives of Psi-function (aka polygamma-functions). */
 extern const unsigned function_index_psi2;
-inline function psi(const ex & p1, const ex & p2) {
-       return function(function_index_psi2, p1, p2);
+template<typename T1, typename T2>
+inline function psi(const T1 & p1, const T2 & p2) {
+       return function(function_index_psi2, ex(p1), ex(p2));
 }
        
 /** Factorial function. */
index 5e1b1eb..4377efd 100644 (file)
@@ -101,7 +101,7 @@ static ex lgamma_series(const ex & arg,
        // if we got here we have to care for a simple pole of tgamma(-m):
        numeric m = -ex_to<numeric>(arg_pt);
        ex recur;
-       for (numeric p; p<=m; ++p)
+       for (numeric p = 0; p<=m; ++p)
                recur += log(arg+p);
        return (lgamma(arg+m+_ex1())-recur).series(rel, order, options);
 }
@@ -137,30 +137,27 @@ static ex tgamma_eval(const ex & x)
 {
        if (x.info(info_flags::numeric)) {
                // trap integer arguments:
-               if (x.info(info_flags::integer)) {
+               const numeric two_x = _num2()*ex_to<numeric>(x);
+               if (two_x.is_even()) {
                        // tgamma(n) -> (n-1)! for postitive n
-                       if (x.info(info_flags::posint)) {
+                       if (two_x.is_positive()) {
                                return factorial(ex_to<numeric>(x).sub(_num1()));
                        } else {
                                throw (pole_error("tgamma_eval(): simple pole",1));
                        }
                }
                // trap half integer arguments:
-               if ((x*2).info(info_flags::integer)) {
+               if (two_x.is_integer()) {
                        // trap positive x==(n+1/2)
                        // tgamma(n+1/2) -> Pi^(1/2)*(1*3*..*(2*n-1))/(2^n)
-                       if ((x*_ex2()).info(info_flags::posint)) {
-                               numeric n = ex_to<numeric>(x).sub(_num1_2());
-                               numeric coefficient = doublefactorial(n.mul(_num2()).sub(_num1()));
-                               coefficient = coefficient.div(pow(_num2(),n));
-                               return coefficient * pow(Pi,_ex1_2());
+                       if (two_x.is_positive()) {
+                               const numeric n = ex_to<numeric>(x).sub(_num1_2());
+                               return (doublefactorial(n.mul(_num2()).sub(_num1())).div(pow(_num2(),n))) * pow(Pi,_ex1_2());
                        } else {
                                // trap negative x==(-n+1/2)
                                // tgamma(-n+1/2) -> Pi^(1/2)*(-2)^n/(1*3*..*(2*n-1))
-                               numeric n = abs(ex_to<numeric>(x).sub(_num1_2()));
-                               numeric coefficient = pow(_num_2(), n);
-                               coefficient = coefficient.div(doublefactorial(n.mul(_num2()).sub(_num1())));;
-                               return coefficient*power(Pi,_ex1_2());
+                               const numeric n = abs(ex_to<numeric>(x).sub(_num1_2()));
+                               return (pow(_num_2(), n).div(doublefactorial(n.mul(_num2()).sub(_num1()))))*power(Pi,_ex1_2());
                        }
                }
                //  tgamma_evalf should be called here once it becomes available
@@ -196,7 +193,7 @@ static ex tgamma_series(const ex & arg,
        if (!arg_pt.info(info_flags::integer) || arg_pt.info(info_flags::positive))
                throw do_taylor();  // caught by function::series()
        // if we got here we have to care for a simple pole at -m:
-       numeric m = -ex_to<numeric>(arg_pt);
+       const numeric m = -ex_to<numeric>(arg_pt);
        ex ser_denom = _ex1();
        for (numeric p; p<=m; ++p)
                ser_denom *= arg+p;
@@ -232,8 +229,8 @@ static ex beta_eval(const ex & x, const ex & y)
                // treat all problematic x and y that may not be passed into tgamma,
                // because they would throw there although beta(x,y) is well-defined
                // using the formula beta(x,y) == (-1)^y * beta(1-x-y, y)
-               numeric nx(ex_to<numeric>(x));
-               numeric ny(ex_to<numeric>(y));
+               const numeric nx = ex_to<numeric>(x);
+               const numeric ny = ex_to<numeric>(y);
                if (nx.is_real() && nx.is_integer() &&
                        ny.is_real() && ny.is_integer()) {
                        if (nx.is_negative()) {
@@ -342,13 +339,13 @@ static ex psi1_evalf(const ex & x)
 static ex psi1_eval(const ex & x)
 {
        if (x.info(info_flags::numeric)) {
-               numeric nx = ex_to<numeric>(x);
+               const numeric nx = ex_to<numeric>(x);
                if (nx.is_integer()) {
                        // integer case 
                        if (nx.is_positive()) {
                                // psi(n) -> 1 + 1/2 +...+ 1/(n-1) - Euler
-                               numeric rat(0);
-                               for (numeric i(nx+_num_1()); i.is_positive(); --i)
+                               numeric rat = 0;
+                               for (numeric i(nx+_num_1()); i>0; --i)
                                        rat += i.inverse();
                                return rat-Euler;
                        } else {
@@ -360,8 +357,8 @@ static ex psi1_eval(const ex & x)
                        // half integer case
                        if (nx.is_positive()) {
                                // psi((2m+1)/2) -> 2/(2m+1) + 2/2m +...+ 2/1 - Euler - 2log(2)
-                               numeric rat(0);
-                               for (numeric i((nx+_num_1())*_num2()); i.is_positive(); i-=_num2())
+                               numeric rat = 0;
+                               for (numeric i = (nx+_num_1())*_num2(); i>0; i-=_num2())
                                        rat += _num2()*i.inverse();
                                return rat-Euler-_ex2()*log(_ex2());
                        } else {
@@ -370,8 +367,8 @@ static ex psi1_eval(const ex & x)
                                // to relate psi(-m-1/2) to psi(1/2):
                                //   psi(-m-1/2) == psi(1/2) + r
                                // where r == ((-1/2)^(-1) + ... + (-m-1/2)^(-1))
-                               numeric recur(0);
-                               for (numeric p(nx); p<0; ++p)
+                               numeric recur = 0;
+                               for (numeric p = nx; p<0; ++p)
                                        recur -= pow(p, _num_1());
                                return recur+psi(_ex1_2());
                        }
@@ -407,7 +404,7 @@ static ex psi1_series(const ex & arg,
        if (!arg_pt.info(info_flags::integer) || arg_pt.info(info_flags::positive))
                throw do_taylor();  // caught by function::series()
        // if we got here we have to care for a simple pole at -m:
-       numeric m = -ex_to<numeric>(arg_pt);
+       const numeric m = -ex_to<numeric>(arg_pt);
        ex recur;
        for (numeric p; p<=m; ++p)
                recur += power(arg+p,_ex_1());
@@ -449,8 +446,8 @@ static ex psi2_eval(const ex & n, const ex & x)
                return log(tgamma(x));
        if (n.info(info_flags::numeric) && n.info(info_flags::posint) &&
                x.info(info_flags::numeric)) {
-               numeric nn = ex_to<numeric>(n);
-               numeric nx = ex_to<numeric>(x);
+               const numeric nn = ex_to<numeric>(n);
+               const numeric nx = ex_to<numeric>(x);
                if (nx.is_integer()) {
                        // integer case 
                        if (nx.is_equal(_num1()))
@@ -462,8 +459,8 @@ static ex psi2_eval(const ex & n, const ex & x)
                                // to relate psi(n,m) to psi(n,1):
                                //   psi(n,m) == psi(n,1) + r
                                // where r == (-)^n * n! * (1^(-n-1) + ... + (m-1)^(-n-1))
-                               numeric recur(0);
-                               for (numeric p(1); p<nx; ++p)
+                               numeric recur = 0;
+                               for (numeric p = 1; p<nx; ++p)
                                        recur += pow(p, -nn+_num_1());
                                recur *= factorial(nn)*pow(_num_1(), nn);
                                return recur+psi(n,_ex1());
@@ -478,7 +475,7 @@ static ex psi2_eval(const ex & n, const ex & x)
                                // use psi(n,1/2) == (-)^(n+1) * n! * (2^(n+1)-1) * zeta(n+1)
                                return pow(_num_1(),nn+_num1())*factorial(nn)*(pow(_num2(),nn+_num1()) + _num_1())*zeta(ex(nn+_num1()));
                        if (nx.is_positive()) {
-                               numeric m = nx - _num1_2();
+                               const numeric m = nx - _num1_2();
                                // use the multiplication formula
                                //   psi(n,2*m) == (psi(n,m) + psi(n,m+1/2)) / 2^(n+1)
                                // to revert to positive integer case
@@ -489,8 +486,8 @@ static ex psi2_eval(const ex & n, const ex & x)
                                // to relate psi(n,-m-1/2) to psi(n,1/2):
                                //   psi(n,-m-1/2) == psi(n,1/2) + r
                                // where r == (-)^(n+1) * n! * ((-1/2)^(-n-1) + ... + (-m-1/2)^(-n-1))
-                               numeric recur(0);
-                               for (numeric p(nx); p<0; ++p)
+                               numeric recur = 0;
+                               for (numeric p = nx; p<0; ++p)
                                        recur += pow(p, -nn+_num_1());
                                recur *= factorial(nn)*pow(_num_1(), nn+_num_1());
                                return recur+psi(n,_ex1_2());
@@ -533,7 +530,7 @@ static ex psi2_series(const ex & n,
        if (!arg_pt.info(info_flags::integer) || arg_pt.info(info_flags::positive))
                throw do_taylor();  // caught by function::series()
        // if we got here we have to care for a pole of order n+1 at -m:
-       numeric m = -ex_to<numeric>(arg_pt);
+       const numeric m = -ex_to<numeric>(arg_pt);
        ex recur;
        for (numeric p; p<=m; ++p)
                recur += power(arg+p,-n+_ex_1());
index 89b0454..229d3f7 100644 (file)
@@ -56,7 +56,7 @@ static ex exp_eval(const ex & x)
                return _ex1();
        }
        // exp(n*Pi*I/2) -> {+1|+I|-1|-I}
-       ex TwoExOverPiI=(_ex2()*x)/(Pi*I);
+       const ex TwoExOverPiI=(_ex2()*x)/(Pi*I);
        if (TwoExOverPiI.info(info_flags::integer)) {
                numeric z=mod(ex_to<numeric>(TwoExOverPiI),_num4());
                if (z.is_equal(_num0()))
@@ -199,7 +199,7 @@ static ex log_series(const ex &arg,
                const symbol *s = static_cast<symbol *>(rel.lhs().bp);
                const ex point = rel.rhs();
                const symbol foo;
-               ex replarg = series(log(arg), *s==foo, order).subs(foo==point);
+               const ex replarg = series(log(arg), *s==foo, order).subs(foo==point);
                epvector seq;
                seq.push_back(expair(-I*csgn(arg*I)*Pi, _ex0()));
                seq.push_back(expair(Order(_ex1()), order));
@@ -230,7 +230,7 @@ static ex sin_evalf(const ex & x)
 static ex sin_eval(const ex & x)
 {
        // sin(n/d*Pi) -> { all known non-nested radicals }
-       ex SixtyExOverPi = _ex60()*x/Pi;
+       const ex SixtyExOverPi = _ex60()*x/Pi;
        ex sign = _ex1();
        if (SixtyExOverPi.info(info_flags::integer)) {
                numeric z = mod(ex_to<numeric>(SixtyExOverPi),_num120());
@@ -312,7 +312,7 @@ static ex cos_evalf(const ex & x)
 static ex cos_eval(const ex & x)
 {
        // cos(n/d*Pi) -> { all known non-nested radicals }
-       ex SixtyExOverPi = _ex60()*x/Pi;
+       const ex SixtyExOverPi = _ex60()*x/Pi;
        ex sign = _ex1();
        if (SixtyExOverPi.info(info_flags::integer)) {
                numeric z = mod(ex_to<numeric>(SixtyExOverPi),_num120());
@@ -394,7 +394,7 @@ static ex tan_evalf(const ex & x)
 static ex tan_eval(const ex & x)
 {
        // tan(n/d*Pi) -> { all known non-nested radicals }
-       ex SixtyExOverPi = _ex60()*x/Pi;
+       const ex SixtyExOverPi = _ex60()*x/Pi;
        ex sign = _ex1();
        if (SixtyExOverPi.info(info_flags::integer)) {
                numeric z = mod(ex_to<numeric>(SixtyExOverPi),_num60());
@@ -650,7 +650,7 @@ static ex atan_series(const ex &arg,
                const symbol *s = static_cast<symbol *>(rel.lhs().bp);
                const ex point = rel.rhs();
                const symbol foo;
-               ex replarg = series(atan(arg), *s==foo, order).subs(foo==point);
+               const ex replarg = series(atan(arg), *s==foo, order).subs(foo==point);
                ex Order0correction = replarg.op(0)+csgn(arg)*Pi*_ex_1_2();
                if ((I*arg_pt)<_ex0())
                        Order0correction += log((I*arg_pt+_ex_1())/(I*arg_pt+_ex1()))*I*_ex_1_2();
@@ -1043,7 +1043,7 @@ static ex atanh_series(const ex &arg,
                const symbol *s = static_cast<symbol *>(rel.lhs().bp);
                const ex point = rel.rhs();
                const symbol foo;
-               ex replarg = series(atanh(arg), *s==foo, order).subs(foo==point);
+               const ex replarg = series(atanh(arg), *s==foo, order).subs(foo==point);
                ex Order0correction = replarg.op(0)+csgn(I*arg)*Pi*I*_ex1_2();
                if (arg_pt<_ex0())
                        Order0correction += log((arg_pt+_ex_1())/(arg_pt+_ex1()))*_ex1_2();
index 92e21b3..9bb753a 100644 (file)
@@ -709,11 +709,11 @@ ex mul::expand(unsigned options) const
        }
        if (expanded_seqp)
                delete expanded_seqp;
-
+       
        // Now the only remaining thing to do is to multiply the factors which
        // were not sums into the "last_expanded" sum
        if (is_ex_exactly_of_type(last_expanded, add)) {
-               add const & finaladd = ex_to<add>(last_expanded);
+               const add & finaladd = ex_to<add>(last_expanded);
                exvector distrseq;
                int n = finaladd.nops();
                distrseq.reserve(n);
index c4cfe69..f7d1a41 100644 (file)
@@ -149,7 +149,7 @@ ex ncmul::expand(unsigned options) const
 {
        // First, expand the children
        exvector expanded_seq = expandchildren(options);
-
+       
        // Now, look for all the factors that are sums and remember their
        // position and number of terms. One remark is in order here: we do not
        // take into account the overall_coeff of the add objects. This is
@@ -171,9 +171,9 @@ ex ncmul::expand(unsigned options) const
                        const add & expanded_addref = ex_to<add>(*cit);
                        number_of_add_operands[number_of_adds] = expanded_addref.seq.size();
                        number_of_expanded_terms *= expanded_addref.seq.size();
-                       number_of_adds++;
+                       ++number_of_adds;
                }
-               current_position++;
+               ++current_position;
        }
 
        // If there are no sums, we are done
index 71781f4..251e573 100644 (file)
@@ -72,13 +72,13 @@ class numeric : public basic
        
        // other ctors
 public:
-       explicit numeric(int i);
-       explicit numeric(unsigned int i);
-       explicit numeric(long i);
-       explicit numeric(unsigned long i);
-       explicit numeric(long numer, long denom);
-       explicit numeric(double d);
-       explicit numeric(const char *);
+       numeric(int i);
+       numeric(unsigned int i);
+       numeric(long i);
+       numeric(unsigned long i);
+       numeric(long numer, long denom);
+       numeric(double d);
+       numeric(const char *);
        
        // functions overriding virtual functions from bases classes
 public:
index f911c49..5daa841 100644 (file)
@@ -67,17 +67,7 @@ DEFAULT_DESTROY(power)
 // other ctors
 //////////
 
-power::power(const ex & lh, const ex & rh) : inherited(TINFO_power), basis(lh), exponent(rh)
-{
-       debugmsg("power ctor from ex,ex",LOGLEVEL_CONSTRUCT);
-}
-
-/** Ctor from an ex and a bare numeric.  This is somewhat more efficient than
- *  the normal ctor from two ex whenever it can be used. */
-power::power(const ex & lh, const numeric & rh) : inherited(TINFO_power), basis(lh), exponent(rh)
-{
-       debugmsg("power ctor from ex,numeric",LOGLEVEL_CONSTRUCT);
-}
+// all inlined
 
 //////////
 // archiving
@@ -493,8 +483,8 @@ ex power::evalf(int level) const
 
 ex power::evalm(void) const
 {
-       ex ebasis = basis.evalm();
-       ex eexponent = exponent.evalm();
+       const ex ebasis = basis.evalm();
+       const ex eexponent = exponent.evalm();
        if (is_ex_of_type(ebasis,matrix)) {
                if (is_ex_of_type(eexponent,numeric)) {
                        return (new matrix(ex_to<matrix>(ebasis).pow(eexponent)))->setflag(status_flags::dynallocated);
@@ -658,7 +648,7 @@ ex power::expand_add(const add & a, int n) const
                upper_limit[l] = n;
        }
        
-       while (1) {
+       while (true) {
                exvector term;
                term.reserve(m+1);
                for (l=0; l<m-1; l++) {
@@ -695,23 +685,12 @@ ex power::expand_add(const add & a, int n) const
                
                term.push_back(f);
                
-               /*
-               cout << "begin term" << endl;
-               for (int i=0; i<m-1; i++) {
-                       cout << "k[" << i << "]=" << k[i] << endl;
-                       cout << "k_cum[" << i << "]=" << k_cum[i] << endl;
-                       cout << "upper_limit[" << i << "]=" << upper_limit[i] << endl;
-               }
-               for_each(term.begin(), term.end(), ostream_iterator<ex>(cout, "\n"));
-               cout << "end term" << endl;
-               */
-               
-               // TODO: optimize this
+               // TODO: Can we optimize this?  Alex seemed to think so...
                sum.push_back((new mul(term))->setflag(status_flags::dynallocated));
                
                // increment k[]
                l = m-2;
-               while ((l>=0)&&((++k[l])>upper_limit[l])) {
+               while ((l>=0) && ((++k[l])>upper_limit[l])) {
                        k[l] = 0;    
                        --l;
                }
@@ -730,7 +709,7 @@ ex power::expand_add(const add & a, int n) const
                        upper_limit[i] = n-k_cum[i-1];
        }
        return (new add(sum))->setflag(status_flags::dynallocated |
-                                                                  status_flags::expanded );
+                                      status_flags::expanded );
 }
 
 
@@ -746,8 +725,8 @@ ex power::expand_add_2(const add & a) const
        // power(+(x,...,z;c),2)=power(+(x,...,z;0),2)+2*c*+(x,...,z;0)+c*c
        // first part: ignore overall_coeff and expand other terms
        for (epvector::const_iterator cit0=a.seq.begin(); cit0!=last; ++cit0) {
-               const ex & r = (*cit0).rest;
-               const ex & c = (*cit0).coeff;
+               const ex & r = cit0->rest;
+               const ex & c = cit0->coeff;
                
                GINAC_ASSERT(!is_ex_exactly_of_type(r,add));
                GINAC_ASSERT(!is_ex_exactly_of_type(r,power) ||
@@ -776,8 +755,8 @@ ex power::expand_add_2(const add & a) const
                }
                        
                for (epvector::const_iterator cit1=cit0+1; cit1!=last; ++cit1) {
-                       const ex & r1 = (*cit1).rest;
-                       const ex & c1 = (*cit1).coeff;
+                       const ex & r1 = cit1->rest;
+                       const ex & c1 = cit1->coeff;
                        sum.push_back(a.combine_ex_with_coeff_to_pair((new mul(r,r1))->setflag(status_flags::dynallocated),
                                                                      _num2().mul(ex_to<numeric>(c)).mul_dyn(ex_to<numeric>(c1))));
                }
@@ -824,18 +803,6 @@ ex power::expand_mul(const mul & m, const numeric & n) const
        return (new mul(distrseq,ex_to<numeric>(m.overall_coeff).power_dyn(n)))->setflag(status_flags::dynallocated);
 }
 
-/*
-ex power::expand_noncommutative(const ex & basis, const numeric & exponent,
-                                                               unsigned options) const
-{
-       ex rest_power = ex(power(basis,exponent.add(_num_1()))).
-                       expand(options | expand_options::internal_do_not_expand_power_operands);
-
-       return ex(mul(rest_power,basis),0).
-              expand(options | expand_options::internal_do_not_expand_mul_operands);
-}
-*/
-
 // helper function
 
 ex sqrt(const ex & a)
index 049d55d..952a72a 100644 (file)
@@ -43,8 +43,8 @@ class power : public basic
        
        // other ctors
 public:
-       power(const ex & lh, const ex & rh);
-       power(const ex & lh, const numeric & rh);
+       power(const ex & lh, const ex & rh) : inherited(TINFO_power), basis(lh), exponent(rh) {}
+       template<typename T> power(const ex & lh, const T & rh) : inherited(TINFO_power), basis(lh), exponent(rh) {}
        
        // functions overriding virtual functions from bases classes
 public:
@@ -80,9 +80,6 @@ protected:
        ex expand_add(const add & a, int n) const;
        ex expand_add_2(const add & a) const;
        ex expand_mul(const mul & m, const numeric & n) const;
-       //ex expand_commutative_3(const ex & basis, const numeric & exponent,
-       //                        unsigned options) const;
-       //ex expand_noncommutative(const ex & basis, const numeric & exponent, unsigned options) const;
        
 // member variables
        
@@ -116,6 +113,11 @@ inline ex pow(const ex & b, const ex & e)
 {
        return power(b, e);
 }
+template<typename T1, typename T2>
+inline ex pow(const T1 & b, const T2 & e)
+{
+       return power(ex(b), ex(e));
+}
 
 /** Square root expression.  Returns a power-object with exponent 1/2. */
 ex sqrt(const ex & a);