From: Chris Dams <Chris.Dams@mi.infn.it>
Date: Thu, 9 Mar 2006 19:13:11 +0000 (+0000)
Subject: Add step function to GiNaCs built-in functions.
X-Git-Tag: release_1-4-0~107
X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=commitdiff_plain;h=5ea3a700d1346284a985872d8efc5695130ba487;ds=sidebyside

Add step function to GiNaCs built-in functions.
---

diff --git a/doc/tutorial/ginac.texi b/doc/tutorial/ginac.texi
index 92d8fca3..d907256d 100644
--- a/doc/tutorial/ginac.texi
+++ b/doc/tutorial/ginac.texi
@@ -1393,6 +1393,8 @@ evaluated immediately:
 @cindex @code{imag()}
 @item @code{csgn(z)}
 @tab complex sign (returns an @code{int})
+@item @code{step(x)}
+@tab step function (returns a @code{numeric})
 @item @code{numer(z)}
 @tab numerator of rational or complex rational number
 @item @code{denom(z)}
@@ -5698,12 +5700,15 @@ GiNaC contains the following predefined mathematical functions:
 @item @code{abs(x)}
 @tab absolute value
 @cindex @code{abs()}
+@item @code{step(x)}
+@tab step function
+@cindex @code{step()}
 @item @code{csgn(x)}
 @tab complex sign
 @cindex @code{conjugate()}
 @item @code{conjugate(x)}
 @tab complex conjugation
-@cindex @code{csgn()}
+@cindex @code{conjugate()}
 @item @code{sqrt(x)}
 @tab square root (not a GiNaC function, rather an alias for @code{pow(x, numeric(1, 2))})
 @cindex @code{sqrt()}
diff --git a/ginac/inifcns.cpp b/ginac/inifcns.cpp
index 9a488fe9..830b1bff 100644
--- a/ginac/inifcns.cpp
+++ b/ginac/inifcns.cpp
@@ -123,6 +123,82 @@ REGISTER_FUNCTION(abs, eval_func(abs_eval).
                        conjugate_func(abs_conjugate).
                        power_func(abs_power));
 
+//////////
+// Step function
+//////////
+
+static ex step_evalf(const ex & arg)
+{
+	if (is_exactly_a<numeric>(arg))
+		return step(ex_to<numeric>(arg));
+	
+	return step(arg).hold();
+}
+
+static ex step_eval(const ex & arg)
+{
+	if (is_exactly_a<numeric>(arg))
+		return step(ex_to<numeric>(arg));
+	
+	else if (is_exactly_a<mul>(arg) &&
+	         is_exactly_a<numeric>(arg.op(arg.nops()-1))) {
+		numeric oc = ex_to<numeric>(arg.op(arg.nops()-1));
+		if (oc.is_real()) {
+			if (oc > 0)
+				// step(42*x) -> step(x)
+				return step(arg/oc).hold();
+			else
+				// step(-42*x) -> step(-x)
+				return step(-arg/oc).hold();
+		}
+		if (oc.real().is_zero()) {
+			if (oc.imag() > 0)
+				// step(42*I*x) -> step(I*x)
+				return step(I*arg/oc).hold();
+			else
+				// step(-42*I*x) -> step(-I*x)
+				return step(-I*arg/oc).hold();
+		}
+	}
+	
+	return step(arg).hold();
+}
+
+static ex step_series(const ex & arg,
+                      const relational & rel,
+                      int order,
+                      unsigned options)
+{
+	const ex arg_pt = arg.subs(rel, subs_options::no_pattern);
+	if (arg_pt.info(info_flags::numeric)
+	    && ex_to<numeric>(arg_pt).real().is_zero()
+	    && !(options & series_options::suppress_branchcut))
+		throw (std::domain_error("step_series(): on imaginary axis"));
+	
+	epvector seq;
+	seq.push_back(expair(step(arg_pt), _ex0));
+	return pseries(rel,seq);
+}
+
+static ex step_power(const ex & arg, const ex & exp)
+{
+	if (exp.info(info_flags::positive))
+		return step(arg);
+	
+	return power(step(arg), exp).hold();
+}
+
+static ex step_conjugate(const ex& arg)
+{
+	return step(arg);
+}
+
+REGISTER_FUNCTION(step, eval_func(step_eval).
+                        evalf_func(step_evalf).
+                        series_func(step_series).
+                        conjugate_func(step_conjugate).
+                        power_func(step_power));
+
 //////////
 // Complex sign
 //////////
diff --git a/ginac/inifcns.h b/ginac/inifcns.h
index 05043819..0e13a040 100644
--- a/ginac/inifcns.h
+++ b/ginac/inifcns.h
@@ -35,6 +35,9 @@ DECLARE_FUNCTION_1P(conjugate_function)
 /** Absolute value. */
 DECLARE_FUNCTION_1P(abs)
 	
+/** Step function. */
+DECLARE_FUNCTION_1P(step)
+
 /** Complex sign. */
 DECLARE_FUNCTION_1P(csgn)
 
diff --git a/ginac/numeric.cpp b/ginac/numeric.cpp
index df525d60..b8170e71 100644
--- a/ginac/numeric.cpp
+++ b/ginac/numeric.cpp
@@ -927,6 +927,18 @@ const numeric numeric::inverse() const
 	return numeric(cln::recip(value));
 }
 
+/** Return the step function of a numeric. The imaginary part of it is
+ *  ignored because the step function is generally considered real but
+ *  a numeric may develop a small imaginary part due to rounding errors.
+ */
+numeric numeric::step() const
+{	cln::cl_R r = cln::realpart(value);
+	if(cln::zerop(r))
+		return numeric(1,2);
+	if(cln::plusp(r))
+		return 1;
+	return 0;
+}
 
 /** Return the complex half-plane (left or right) in which the number lies.
  *  csgn(x)==0 for x==0, csgn(x)==1 for Re(x)>0 or Re(x)=0 and Im(x)>0,
diff --git a/ginac/numeric.h b/ginac/numeric.h
index 0456b81d..f0e3362f 100644
--- a/ginac/numeric.h
+++ b/ginac/numeric.h
@@ -148,6 +148,7 @@ public:
 	const numeric & operator=(double d);
 	const numeric & operator=(const char *s);
 	const numeric inverse() const;
+	numeric step() const;
 	int csgn() const;
 	int compare(const numeric &other) const;
 	bool is_equal(const numeric &other) const;
@@ -250,6 +251,9 @@ inline const numeric pow(const numeric &x, const numeric &y)
 inline const numeric inverse(const numeric &x)
 { return x.inverse(); }
 
+inline numeric step(const numeric &x)
+{ return x.step(); }
+
 inline int csgn(const numeric &x)
 { return x.csgn(); }