Add step function to GiNaCs built-in functions.

diff --git a/doc/tutorial/ginac.texi b/doc/tutorial/ginac.texi
index 92d8fca38933bdfd9ff8d307849ead1139a071bc..d907256dca46dc8e597c2b53e649b8cee3d12a90 100644 (file)
--- 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})
 @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)}
 @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{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
 @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()}
 @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 9a488fe9219b383ff2111a6bc0a808bf14d67a6e..830b1bff41ff556816e6c398add542179ad6f7e2 100644 (file)
--- 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));
 
                        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
 //////////
 //////////
 // Complex sign
 //////////

diff --git a/ginac/inifcns.h b/ginac/inifcns.h
index 05043819e9deeceb6a42d54cc03c5140238399e9..0e13a040f1b42d85644e4525fd227ca11eda9c4c 100644 (file)
--- a/ginac/inifcns.h
+++ b/ginac/inifcns.h
@@ -35,6 +35,9 @@ DECLARE_FUNCTION_1P(conjugate_function)
 /** Absolute value. */
 DECLARE_FUNCTION_1P(abs)
        
 /** Absolute value. */
 DECLARE_FUNCTION_1P(abs)
        

+/** Step function. */
+DECLARE_FUNCTION_1P(step)
+

 /** Complex sign. */
 DECLARE_FUNCTION_1P(csgn)
 
 /** Complex sign. */
 DECLARE_FUNCTION_1P(csgn)
 

diff --git a/ginac/numeric.cpp b/ginac/numeric.cpp
index df525d6025eaf4221cbcb3c2fda145068eede1f6..b8170e716ff3220685cb490d1a9ccf6473c23947 100644 (file)
--- a/ginac/numeric.cpp
+++ b/ginac/numeric.cpp
@@ -927,6 +927,18 @@ const numeric numeric::inverse() const
        return numeric(cln::recip(value));
 }
 
        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,
 
 /** 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 0456b81de8af5d4877b793096d6acadae66c5d69..f0e3362f8965aec7bcd5aaf199e6d65f1f3362d7 100644 (file)
--- 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;
        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;
        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 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(); }
 
 inline int csgn(const numeric &x)
 { return x.csgn(); }