C++ Implementation (by Peter Sturdza)

The implementation of this method in C/C++ algorithms with a C++ compiler is facilitated by the fact that the intrinsic functions and operators can be redefined and overloaded. It is easier to implement than in Fortran, and the procedure is as follows:

• Make sure you have a reasonably up to date C++ compiler (1998 ANSI/ISO).
• Download the complexify.h C++ program file and #include it into your source code. It isn’t really a header file, so if you have many source files, then you may want to split complexify.h into a real header file to be included everywhere, and a complexify.cc program file.
• Change all variable declarations from double to cplx in your program.
• Fix the input/output methods. For ascii input/output, it is reasonably straightforward. You have several choices:
  • Manually control the real and imaginary parts, as in:

    #include <stdio.h>
    #include "complexify.h"
    int main() {
    cplx age = cplx(45.3,0.01);
    printf("I am %g years old,\n",real(age));
    printf("and am loosing %g strands per day.\n",imag(age));
    }
    

  • Use the C++ iostream methods:

    #include <fstream.h>
    #include "complexify.h"
    int main() {
    char inputfile[10]="in.dat";
    cplx age;
    ifstream input(inputfile,ios::in);
    input >> age;
    input.close();
    cout << age << endl;
    }
    

    This way, the input file can contain a complex valued number, (45.3,0.01), or a real number, 45.3. The » operator sets the imaginary part to zero in the latter case.

  • A combination of the above: It is convenient to use the C++ stream library for ascii input because it automatically deals with real and complex numbers. In fact, input files for real-valued codes can be used unmodified after the code is complexified. Only the variables that need imaginary parts have to then be changed. For output, however, it is often nice to separate the real and imaginary parts (and only print imaginary parts for certain variables of interest – or print the derivative itself).
• Compile the new code and run it.
• Alternatively, a simple form of algorithmic differentiation can be used instead of the complex-step approach. The only difference is that the complexify.h file is replaced with derivify.h, and occurrences of double must now be replaced with surreal. With this method the “step size” no longer has to be small.

C++ Template Implementation

Using the template facility in C++, one program can contain either the real-valued, complex-step or algorithmic differentiation codes (and single and double precision versions of each, if you really felt like it). This is not simple to implement, but is very elegant, and once done, the user can decide which version to execute at run time – they are all compiled into one big executable. The resulting source code, however, remains compact and easy to maintain. Memory usage and speed of execution are not noticeably affected when templates are used in this way.

Below is an example that shows how templates can be used to complexify or derivify a code (or both). Only the main() routine declares the type of the variable age. The rest of the program is independent of that choice by using the placeholder type D.

The primary advantage is version control: when the program is modified, the real, complex-step and algorithmic differentiation versions remain identical. This also means that debugging the real-valued version automatically debugs the other versions of the program as well.

The disadvantage is that the C++ template syntax is complicated. Also, if the original code is written in C, then it must be made “object oriented.” But that is not too bad. You can cheat, as was done here, by creating a template class that just contains all of the original subroutines and the original global variables. The constructor of this new class is the main() in the original C program, and the new main() routine decides whether to run the real-valued, complex-step or algorithmic differentiation version. There may be a better way, but this general idea should work fine.

Starting from the following non-object-oriented and real-valued program:

#include <stdio.h>
#include <iostream.h>

void printage(double age)
{
printf("I am %g years old,\n",age);
}

double getage()
{
double age;
cout << "Type in age: " << endl;
cin >> age;
return age;
}

int main(void)
{
printage(getage());
}

The application of templates allows the program to be complexified and derivified like so: 
#include <stdio.h>
#include <iostream.h>
#include <string.h>
#include "complexify.h"
#include "derivify.h"

template<class D>
class Age {
public:
Age() {printage(getage());}
void printage(D);
D getage();
};

template<class D>
void Age<D>::printage(D age)
{
printf("I am %g years old,\n",real(age));
printf("and am losing %g strands per day.\n",imag(age));
}

template<class D>
D Age<D>::getage()
{
D age;
cout << "Type in age, 'r', for real," << endl;
cout << " and '(r,i)', for complex or surreal: ";
cin >> age;
return age;
}

int main(int argc, char **argv)
{
if(argc!=2) {
cerr << "usage: " << argv[0] << " real" << endl;
cerr << " or: " << argv[0] << " complex" << endl;
cerr << " or: " << argv[0] << " surreal" << endl;
}
else {
if(!strcmp(argv[1],"real")) Age<double> a;
else if(!strcmp(argv[1],"complex")) Age<cplx> a;
else if(!strcmp(argv[1],"surreal")) Age<surreal> a;
}
return 0;
}

C Implementation

The 1999 ANSI version of C does have a built-in complex type, but there is no clean way to overload operators and functions. It is, however, still possible to implement the method by writing the new function definitions under different names and changing the function calls to those names. For example, you would write a function called complex_abs() and substitute all the calls to abs() with this one. You also need to change the comparison operations to compare only the real parts of the arguments.