The easiest and most instructive way to start extending GiNaC is probably to create your own symbolic functions. These are implemented with the help of two preprocessor macros:
DECLARE_FUNCTION_<n>P(<name>)
REGISTER_FUNCTION(<name>, <options>)
The DECLARE_FUNCTION macro will usually appear in a header file. It
declares a C++ function with the given ‘name’ that takes exactly ‘n’
parameters of type ex and returns a newly constructed GiNaC
function object that represents your function.
The REGISTER_FUNCTION macro implements the function. It must be passed
the same ‘name’ as the respective DECLARE_FUNCTION macro, and a
set of options that associate the symbolic function with C++ functions you
provide to implement the various methods such as evaluation, derivative,
series expansion etc. They also describe additional attributes the function
might have, such as symmetry and commutation properties, and a name for
LaTeX output. Multiple options are separated by the member access operator
‘.’ and can be given in an arbitrary order.
(By the way: in case you are worrying about all the macros above we can assure you that functions are GiNaC's most macro-intense classes. We have done our best to avoid macros where we can.)
Here is an example for the implementation of a function with two arguments that is not further evaluated:
DECLARE_FUNCTION_2P(myfcn)
REGISTER_FUNCTION(myfcn, dummy())
Any code that has seen the DECLARE_FUNCTION line can use myfcn()
in algebraic expressions:
{
...
symbol x("x");
ex e = 2*myfcn(42, 1+3*x) - x;
cout << e << endl;
// prints '2*myfcn(42,1+3*x)-x'
...
}
The dummy() option in the REGISTER_FUNCTION line signifies
"no options". A function with no options specified merely acts as a kind of
container for its arguments. It is a pure "dummy" function with no associated
logic (which is, however, sometimes perfectly sufficient).
Let's now have a look at the implementation of GiNaC's cosine function for an example of how to make an "intelligent" function.
The GiNaC header file inifcns.h contains the line
DECLARE_FUNCTION_1P(cos)
which declares to all programs using GiNaC that there is a function ‘cos’
that takes one ex as an argument. This is all they need to know to use
this function in expressions.
The implementation of the cosine function is in inifcns_trans.cpp. Here
is its REGISTER_FUNCTION line:
REGISTER_FUNCTION(cos, eval_func(cos_eval).
evalf_func(cos_evalf).
derivative_func(cos_deriv).
latex_name("\\cos"));
There are four options defined for the cosine function. One of them
(latex_name) gives the function a proper name for LaTeX output; the
other three indicate the C++ functions in which the "brains" of the cosine
function are defined.
The eval_func() option specifies the C++ function that implements
the eval() method, GiNaC's anonymous evaluator. This function takes
the same number of arguments as the associated symbolic function (one in this
case) and returns the (possibly transformed or in some way simplified)
symbolically evaluated function (See Automatic evaluation, for a description
of the automatic evaluation process). If no (further) evaluation is to take
place, the eval_func() function must return the original function
with .hold(), to avoid a potential infinite recursion. If your
symbolic functions produce a segmentation fault or stack overflow when
using them in expressions, you are probably missing a .hold()
somewhere.
The eval_func() function for the cosine looks something like this
(actually, it doesn't look like this at all, but it should give you an idea
what is going on):
static ex cos_eval(const ex & x)
{
if ("x is a multiple of 2*Pi")
return 1;
else if ("x is a multiple of Pi")
return -1;
else if ("x is a multiple of Pi/2")
return 0;
// more rules...
else if ("x has the form 'acos(y)'")
return y;
else if ("x has the form 'asin(y)'")
return sqrt(1-y^2);
// more rules...
else
return cos(x).hold();
}
This function is called every time the cosine is used in a symbolic expression:
{
...
e = cos(Pi);
// this calls cos_eval(Pi), and inserts its return value into
// the actual expression
cout << e << endl;
// prints '-1'
...
}
In this way, cos(4*Pi) automatically becomes 1,
cos(asin(a+b)) becomes sqrt(1-(a+b)^2), etc. If no reasonable
symbolic transformation can be done, the unmodified function is returned
with .hold().
GiNaC doesn't automatically transform cos(2) to ‘-0.416146...’.
The user has to call evalf() for that. This is implemented in a
different function:
static ex cos_evalf(const ex & x)
{
if (is_a<numeric>(x))
return cos(ex_to<numeric>(x));
else
return cos(x).hold();
}
Since we are lazy we defer the problem of numeric evaluation to somebody else,
in this case the cos() function for numeric objects, which in
turn hands it over to the cos() function in CLN. The .hold()
isn't really needed here, but reminds us that the corresponding eval()
function would require it in this place.
Differentiation will surely turn up and so we need to tell cos
what its first derivative is (higher derivatives, .diff(x,3) for
instance, are then handled automatically by basic::diff and
ex::diff):
static ex cos_deriv(const ex & x, unsigned diff_param)
{
return -sin(x);
}
The second parameter is obligatory but uninteresting at this point. It specifies which parameter to differentiate in a partial derivative in case the function has more than one parameter, and its main application is for correct handling of the chain rule.
An implementation of the series expansion is not needed for cos() as
it doesn't have any poles and GiNaC can do Taylor expansion by itself (as
long as it knows what the derivative of cos() is). tan(), on
the other hand, does have poles and may need to do Laurent expansion:
static ex tan_series(const ex & x, const relational & rel,
int order, unsigned options)
{
// Find the actual expansion point
const ex x_pt = x.subs(rel);
if ("x_pt is not an odd multiple of Pi/2")
throw do_taylor(); // tell function::series() to do Taylor expansion
// On a pole, expand sin()/cos()
return (sin(x)/cos(x)).series(rel, order+2, options);
}
The series() implementation of a function must return a
pseries object, otherwise your code will crash.
GiNaC functions understand several more options which are always
specified as .option(params). None of them are required, but you
need to specify at least one option to REGISTER_FUNCTION(). There
is a do-nothing option called dummy() which you can use to define
functions without any special options.
eval_func(<C++ function>)
evalf_func(<C++ function>)
derivative_func(<C++ function>)
series_func(<C++ function>)
conjugate_func(<C++ function>)
These specify the C++ functions that implement symbolic evaluation,
numeric evaluation, partial derivatives, and series expansion, respectively.
They correspond to the GiNaC methods eval(), evalf(),
diff() and series().
The eval_func() function needs to use .hold() if no further
automatic evaluation is desired or possible.
If no series_func() is given, GiNaC defaults to simple Taylor
expansion, which is correct if there are no poles involved. If the function
has poles in the complex plane, the series_func() needs to check
whether the expansion point is on a pole and fall back to Taylor expansion
if it isn't. Otherwise, the pole usually needs to be regularized by some
suitable transformation.
latex_name(const string & n)
specifies the LaTeX code that represents the name of the function in LaTeX
output. The default is to put the function name in an \mbox{}.
do_not_evalf_params()
This tells evalf() to not recursively evaluate the parameters of the
function before calling the evalf_func().
set_return_type(unsigned return_type, const return_type_t * return_type_tinfo)
This allows you to explicitly specify the commutation properties of the
function (See Non-commutative objects, for an explanation of
(non)commutativity in GiNaC). For example, with an object of type
return_type_t created like
return_type_t my_type = make_return_type_t<matrix>();
you can use set_return_type(return_types::noncommutative, &my_type) to
make GiNaC treat your function like a matrix. By default, functions inherit the
commutation properties of their first argument. The utilized template function
make_return_type_t<>()
template<typename T> inline return_type_t make_return_type_t(const unsigned rl = 0)
can also be called with an argument specifying the representation label of the non-commutative function (see section on dirac gamma matrices for more details).
set_symmetry(const symmetry & s)
specifies the symmetry properties of the function with respect to its arguments. See Indexed objects, for an explanation of symmetry specifications. GiNaC will automatically rearrange the arguments of symmetric functions into a canonical order.
Sometimes you may want to have finer control over how functions are
displayed in the output. For example, the abs() function prints
itself as ‘abs(x)’ in the default output format, but as ‘|x|’
in LaTeX mode, and fabs(x) in C source output. This is achieved
with the
print_func<C>(<C++ function>)
option which is explained in the next section.
The DECLARE_FUNCTION and REGISTER_FUNCTION macros define
functions with a fixed number of arguments. Sometimes, though, you may need
to have a function that accepts a variable number of expressions. One way to
accomplish this is to pass variable-length lists as arguments. The
Li() function uses this method for multiple polylogarithms.
It is also possible to define functions that accept a different number of
parameters under the same function name, such as the psi() function
which can be called either as psi(z) (the digamma function) or as
psi(n, z) (polygamma functions). These are actually two different
functions in GiNaC that, however, have the same name. Defining such
functions is not possible with the macros but requires manually fiddling
with GiNaC internals. If you are interested, please consult the GiNaC source
code for the psi() function (inifcns.h and
inifcns_gamma.cpp).