author Chris Dams Thu, 9 Mar 2006 19:13:11 +0000 (19:13 +0000) committer Chris Dams Thu, 9 Mar 2006 19:13:11 +0000 (19:13 +0000)
 doc/tutorial/ginac.texi patch | blob | history ginac/inifcns.cpp patch | blob | history ginac/inifcns.h patch | blob | history ginac/numeric.cpp patch | blob | history ginac/numeric.h patch | blob | history

@@ -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()}
@@ -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
//////////
index 05043819e9deeceb6a42d54cc03c5140238399e9..0e13a040f1b42d85644e4525fd227ca11eda9c4c 100644 (file)
@@ -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)

index df525d6025eaf4221cbcb3c2fda145068eede1f6..b8170e716ff3220685cb490d1a9ccf6473c23947 100644 (file)
@@ -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,
@@ -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(); }