$declare_function_macro=generate(
<<'END_OF_DECLARE_FUNCTION_MACRO','GiNaC::ex const & p${N}','p${N}');
#define DECLARE_FUNCTION_${N}P(NAME) \\
-extern unsigned function_index_##NAME; \\
+extern const unsigned function_index_##NAME; \\
inline GiNaC::function NAME(${SEQ1}) { \\
return GiNaC::function(function_index_##NAME, ${SEQ2}); \\
}
// end of generated lines
#define REGISTER_FUNCTION(NAME,E,EF,D,S) \\
-unsigned function_index_##NAME=GiNaC::function::register_new(#NAME,E,EF,D,S);
+const unsigned function_index_##NAME=GiNaC::function::register_new(#NAME,E,EF,D,S);
#define BEGIN_TYPECHECK \\
bool automatic_typecheck=true;
DECLARE_FUNCTION_1P(gamma)
/** Psi-function (aka polygamma-function). */
-extern unsigned function_index_psi1;
+extern const unsigned function_index_psi1;
inline function psi(ex const & p1) {
return function(function_index_psi1, p1);
}
-extern unsigned function_index_psi2;
+extern const unsigned function_index_psi2;
inline function psi(ex const & p1, ex const & p2) {
return function(function_index_psi2, p1, p2);
}
static ex psi1_eval(ex const & x)
{
if (x.info(info_flags::numeric)) {
- // do some stuff...
+ if (x.info(info_flags::integer) && !x.info(info_flags::positive))
+ throw (std::domain_error("psi_eval(): simple pole"));
+ if (x.info(info_flags::positive)) {
+ // psi(n) -> 1 + 1/2 +...+ 1/(n-1) - EulerGamma
+ if (x.info(info_flags::integer)) {
+ numeric rat(0);
+ for (numeric i(ex_to_numeric(x)-numONE()); i.is_positive(); --i)
+ rat += i.inverse();
+ return rat-EulerGamma;
+ }
+ // psi((2m+1)/2) -> 2/(2m+1) + 2/2m +...+ 2/1 - EulerGamma - 2log(2)
+ if ((exTWO()*x).info(info_flags::integer)) {
+ numeric rat(0);
+ for (numeric i((ex_to_numeric(x)-numONE())*numTWO()); i.is_positive(); i-=numTWO())
+ rat += numTWO()*i.inverse();
+ return rat-EulerGamma-exTWO()*log(exTWO());
+ }
+ if (x.compare(exONE())==1) {
+ // should call numeric, since >1
+ }
+ }
}
return psi(x).hold();
-}
+}
static ex psi1_evalf(ex const & x)
{
return psi(exONE(), x);
}
-static ex psi1_series(ex const & x, symbol const & s, ex const & point, int order)
-{
- throw(std::logic_error("Nobody told me how to series expand the psi function. :-("));
-}
-
-unsigned function_index_psi1 = function::register_new("psi", psi1_eval, psi1_evalf, psi1_diff, psi1_series);
+const unsigned function_index_psi1 = function::register_new("psi", psi1_eval, psi1_evalf, psi1_diff, NULL);
//////////
// Psi-functions (aka polygamma-functions) psi(0,x)==psi(x)
{
// psi(0,x) -> psi(x)
if (n.is_zero())
- return psi(x).hold();
+ return psi(x);
if (n.info(info_flags::numeric) && x.info(info_flags::numeric)) {
// do some stuff...
}
return psi(n+1, x);
}
-static ex psi2_series(ex const & n, ex const & x, symbol const & s, ex const & point, int order)
-{
- throw(std::logic_error("Nobody told me how to series expand the psi functions. :-("));
-}
-
-unsigned function_index_psi2 = function::register_new("psi", psi2_eval, psi2_evalf, psi2_diff, psi2_series);
+const unsigned function_index_psi2 = function::register_new("psi", psi2_eval, psi2_evalf, psi2_diff, NULL);
} // namespace GiNaC
if (y.is_integer()) {
if (y.is_zero())
return -exHALF();
- if (!x.compare(exONE()))
+ if (x.is_equal(exONE()))
throw(std::domain_error("zeta(1): infinity"));
if (x.info(info_flags::posint)) {
if (x.info(info_flags::odd))
inline numeric denom(numeric const & x)
{ return x.denom(); }
-ex IEvalf(void);
+// numeric evaluation functions for class constant objects:
+
ex PiEvalf(void);
ex EulerGammaEvalf(void);
ex CatalanEvalf(void);
numeric operator-(numeric const & lh)
{
- return (numeric(-1)*lh);
+ return numMINUSONE()*lh;
+}
+
+/** Numeric prefix increment. Adds 1 and returns incremented number. */
+numeric& operator++(numeric & rh)
+{
+ rh = rh+numONE();
+ return rh;
+}
+
+/** Numeric prefix decrement. Subtracts 1 and returns decremented number. */
+numeric& operator--(numeric & rh)
+{
+ rh = rh-numONE();
+ return rh;
+}
+
+/** Numeric postfix increment. Returns the number and leaves the original
+ * incremented by 1. */
+numeric operator++(numeric & lh, int)
+{
+ numeric tmp = lh;
+ lh = lh+numONE();
+ return tmp;
+}
+
+/** Numeric Postfix decrement. Returns the number and leaves the original
+ * decremented by 1. */
+numeric operator--(numeric & lh, int)
+{
+ numeric tmp = lh;
+ lh = lh-numONE();
+ return tmp;
}
// binary relational operators ex with ex
numeric operator+(numeric const & lh);
numeric operator-(numeric const & lh);
+numeric& operator++(numeric & rh);
+numeric& operator--(numeric & rh);
+numeric operator++(numeric & lh, int);
+numeric operator--(numeric & lh, int);
// binary relational operators ex with ex
relational operator==(ex const & lh, ex const & rh);