Should I convert FORTRAN code to C?
kenny at m.cs.uiuc.edu
kenny at m.cs.uiuc.edu
Wed Jul 6 13:39:00 AEST 1988
On the discussion of whether or not `recursion is necessary,' I have
at least a nickel's worth of comment. The example shown is frightful!
In article <817 at garth.UUCP> someone writes (correctly):
>However convenient it might be, recursion is not strictly necessary in
>practical programs....
/* Written 9:15 am Jun 29, 1988 by shenkin at cubsun.BIO.COLUMBIA.EDU in m.cs.uiuc.edu:comp.lang.c */
Peter Shenkin <shenkin at cubsun.bio.columbia.edu> replies:
>I used to agree, and most people agree. Most textbook examples of recursion
>-- eg, calculation of factorial(n) -- would be much clearer and faster, in
>practice if coded using simple iteration. But I've come across a class of
>problem which can only be generally handled through recursion, I think. I'd
>be intrigued if someone can come up with a Fortran solution. (Mail it to me,
>and I'll post a summary to comp.lang.fortran.)
>This class involves, conceptually, nested DO-loops, when the depth of the nest
>is unknown at compile time. For example, the following C program generates
>grid points within a hypercube, given the number of degrees of freedom
>(dimensionality) of the space, and the number of gridpoints (number of levels)
>along each coordinate. Turns out I need to do more complicated versions of
>this sort of thing all the time in numerical work.
>NB: This *can* be coded in Fortran by setting up the maximum number
>of needed DO-loops and escaping early, if need be. But this is
>unwieldy, at best , unless this maximum is small. And here (if one
>malloc()'s -- see my note on MAXDF) a maximum need not be specified.
That comment is utter poppycock. Recursion can *always* be
eliminated, if necessary by setting up an auxiliary stack. I shall
demonstrate the techniques with your sample program.
I'm posting this, rather than emailing it, because I think it also
gives a good example of stepwise refinement of a program. I give ten
versions of the same program, using the same algorithm, each of which
attempts to be an improvement on its predecessors. These discussions
should give some idea of how I, at least, attack this type of problem.
(The `>'s are omitted in the example program)
PROGRAM 0. The example program, as originally posted.
/*****************************************************************************/
/* grid.c: generate grid points in a hypercube */
/* compile by saying, "cc -o grid grid.c"
/* invoke by saying, eg, "grid 4 2", to generate the vertices, only, of
* a 4-D hypercube.
*/
#define MAXDF 100
/* max number of degrees of freedom -- not really needed,
* since we could malloc() value[] in grdgen() */
main(argc,argv)
int argc;
char *argv[];
{
int df,nlev; /* number of deg. of freedom, number levels per df */
df= atoi(argv[1]);
nlev= atoi(argv[2]);
printf("df=%d; nlev=%d\n",df,nlev);
grdgen(df,nlev);
}
/*****************************************************************************/
grdgen(df,nlev)
int df; /* no. of degrees of freedom */
int nlev; /* no. of levels each df */
{
static int begin=1; /* ==0 if func was called by itself */
static int totaldf; /* total no. of degrees of freedom (1st call)*/
static int value[MAXDF]; /* to hold current values of all df's */
int i,j;
if(begin)
totaldf=df;
if(--df ==0) { /* iterate through last df and print values */
for(i=0; i<nlev; ++i) {
value[df] = i;
for(j=0; j<totaldf; ++j) {
printf("%d ",value[j]); }
printf("\n");
}
begin=1;
}
else
for(i=0; i<nlev; ++i) {
value[df] = i;
begin = 0;
grdgen(df,nlev);
}
}
Not only can recursion be eliminated in this type of program; the
recursion can be eliminated mechanically; a *really* good optimizing
compiler could conceivably manage it.
Let's first recode your example to a `purer' recursive form,
eliminating the innermost _for_ loop. That change gives us simpler
code to work with:
PROGRAM 1
[main program is the same as before]
static int nlev;
static int totaldf;
static int value [MAXDF]; /* to hold current values of all df's */
grdgen_aux (df) {
int i;
if (df < 0) { /* Print a grid point */
register int j;
for (j = 0; j < totaldf - 1; ++j)
printf ("%d ", value [j]);
printf ("%d\n", value [j]);
}
else { /* Generate grid points with df degrees */
/* of freedom */
for (i = 0; i < nlev; ++i) {
value [df] = i;
grdgen_aux (df - 1);
}
}
}
grdgen(df,nl)
int df; /* no. of degrees of freedom */
int nl; /* no. of levels each df */
{
totaldf=df;
nlev = nl;
grdgen_aux (df - 1);
}
Here, all the recursion is handled in the `grdgen_aux' procedure. The
useless variable `begin' is removed, and the control structure has
been simplified significantly.
Now let's work on eliminating the recursion in `grdgen_aux.' We note
that the only two automatic variables that are destroyed on recursion
are _df_ and _i_. We also note that _df_ is always decremented by 1,
and that _i_ can be recovered from the _value_ array. (Were this not
the case, we could maintain an auxiliary stack of values.) We can
eliminate the recursion by replacing the recursive call by a _goto_ to
the function entry (If you think that _goto_ is a dirty word, take
heart; I'm going to take the _gotos_ out again!), and the function
return with a test for the top-level invocation and a _goto_ to the
return point, giving the following version of the program.
PROGRAM 2.
[main program is still the same]
static int nlev;
static int totaldf;
static int value [MAXDF]; /* to hold current values of all df's */
grdgen_aux (df) {
int i;
start:
if (df < 0) { /* Print a grid point */
register int j;
for (j = 0; j < totaldf - 1; ++j)
printf ("%d ", value [j]);
printf ("%d\n", value [j]);
}
else { /* Generate grid points with df degrees */
/* of freedom */
for (i = 0; i < nlev; ++i) {
value [df--] = i;
goto start;
finish:
i = value [++df];
}
}
if (df < totaldf - 1) goto finish;
}
grdgen(df,nl)
int df; /* no. of degrees of freedom */
int nl; /* no. of levels each df */
{
totaldf=df;
nlev = nl;
grdgen_aux (df - 1);
}
Now our problem is getting rid of the _goto_ statements, and
untangling this bowl of spaghetti back into something readable.
There are a couple of ways to attack this.
METHOD I. Use a finite-state machine.
Here, we take some of the state information out of the program
counter, and put it in an auxiliary variable. We use a _switch_
statement on this auxiliary variable to control program flow. Fortran
hackers will recognize that a computed or assigned GO TO will do the
switch quite nicely.
PROGRAM 3.
[main program is still the same]
static int nlev;
static int totaldf;
static int value [MAXDF]; /* to hold current values of all df's */
grdgen_aux (df) {
int i;
enum { EXIT, WORK, PUSH, POP } state;
state = WORK;
for (;;) {
switch (state) {
default:
return;
case WORK:
if (df >= 0) {
i = 0;
state = PUSH;
}
else {
register int j;
for (j = 0; j < (totaldf - 1); ++j)
printf ("%d ", value [j]);
printf ("%d\n", value [j]);
state = POP;
}
break;
case PUSH:
if (i >= nlev) state = POP;
else {
value [df--] = i;
state = WORK;
}
break;
case POP:
if (df >= (totaldf - 1))
state = EXIT;
else {
i = value [++df] + 1;
state = PUSH;
}
break;
}
}
}
grdgen(df,nl)
int df; /* no. of degrees of freedom */
int nl; /* no. of levels each df */
{
totaldf=df;
nlev = nl;
grdgen_aux (df - 1);
}
Folding the `grdgen_aux' procedure back into `grdgen', now that the
former is no longer recursive, and performing a few trivial
optimizations, gives:
PROGRAM 4.
[main program is still the same]
grdgen (df, nlev)
register int df; /* no. of degrees of freedom */
register int nlev; /* no. of levels each df */
{
register int i, j;
register int totaldfm1;
static int value [MAXDF]; /* to hold current values of all df's */
register enum { EXIT = 0, WORK, PUSH, POP } state;
totaldfm1 = df - 1;
state = WORK;
while (state) {
switch (state) {
case WORK:
if (df >= 0) {
i = 0;
state = PUSH;
}
else {
register int j;
for (j = 0; j < totaldfm1; ++j)
printf ("%d ", value [j]);
printf ("%d\n", value [j]);
state = POP;
}
break;
case PUSH:
if (i >= nlev) state = POP;
else {
value [df--] = i;
state = WORK;
}
break;
case POP:
if (df >= totaldfm1)
state = EXIT;
else {
i = value [++df] + 1;
state = PUSH;
}
break;
}
}
}
Pascal devotees will note that Program 4 doesn't require any
non-Pascal control structures (although the _switch_ has to be
replaced with a chain of if ... else if ... else ... constructs.
That argument is probably the only really good argument in favor of
the state-variable technique.
METHOD II. Code copying.
If code fragments in the `grdgen_aux' procedure are duplicated in a
disciplined fashion, the code can indeed be unscrambled to a
_goto_-less version. Note that any optimizer that performs interval
analysis, node listings, or T1-T2 path compression will need to do the
same sort of code copying as part of the optimization phase; I wouldn't
be surprised if the object code generated by a good optimizing
compiler for Program 2 resembles that generated for one of the
following.
PROGRAM 5.
[main program is still the same]
static int nlev;
static int totaldf;
static int value [MAXDF]; /* to hold current values of all df's */
grdgen_aux (df) {
int i;
for (;;) {
for (;;) {
if (df < 0) { /* Print a node */
register int j;
for (j = 0; j < (totaldf - 1); ++j)
printf ("%d ", value [j]);
printf ("%d\n", value [j]);
break;
}
i = 0;
if (i >= nlev) break; /* This probably can't happen */
value [df--] = i;
}
for (;;) { /* Develop another node */
if (df >= (totaldf - 1)) return;
i = value [++df] + 1;
if (i < nlev) {
value [df--] = i;
break;
}
}
}
}
grdgen(df,nl)
int df; /* no. of degrees of freedom */
int nl; /* no. of levels each df */
{
totaldf=df;
nlev = nl;
grdgen_aux (df - 1);
}
Program 5 was generated by a reasonable mechanical technique, given
Program 2. We see that the control structures are needlessly complex;
in particular, the statement marked, /*this probably can't happen*/,
is executed only if the number of levels is zero or fewer. Since the
program generates null output in that case, it probably pays to check
nlev in advance and simplify the control structures, yielding the
following.
PROGRAM 6.
[guess what? I haven't changed the main program.]
static int nlev;
static int totaldfm1;
static int value [MAXDF]; /* to hold current values of all df's */
grdgen_aux (df) {
int i;
for (;;) {
register int j;
while (df >= 0) {
i = 0;
value [df--] = i;
}
for (j = 0; j < totaldfm1; ++j) /* Print a node */
printf ("%d ", value [j]);
printf ("%d\n", value [j]);
for (;;) { /* Develop another node */
if (df >= totaldfm1) return;
i = value [++df] + 1;
if (i < nlev) {
value [df--] = i;
break;
}
}
}
}
grdgen(df,nl)
int df; /* no. of degrees of freedom */
int nl; /* no. of levels each df */
{
if (nl <= 0) return; /* Sanity check */
totaldfm1 = df - 1;
nlev = nl;
grdgen_aux (totaldfm1);
}
Once again, we can move `grdgen_aux' in line, and provide trivial
optimizations. The final version then looks like:
PROGRAM 7.
[main program still hasn't changed]
grdgen (df, nlev)
register int df; /* no. of degrees of freedom */
register int nlev; /* no. of levels each df */
{
register int i;
register int j;
register int totaldfm1 = df - 1;
static int value [MAXDF]; /* to hold current values of all df's */
for (;;) {
while (df >= 0) { /* Set up initial values */
i = 0;
value [df--] = i;
}
for (j = 0; j < totaldfm1; ++j) /* Print a node */
printf ("%d ", value [j]);
printf ("%d\n", value [j]);
for (;;) { /* Develop another node */
if (df >= totaldfm1) return;
i = value [++df] + 1;
if (i < nlev) {
value [df--] = i;
break;
}
}
}
}
There is another choice for code to be copied: we can, instead of
Program 5, come up with:
PROGRAM 8.
[same main program]
static int nlev;
static int totaldf;
static int value [MAXDF]; /* to hold current values of all df's */
grdgen_aux (df) {
int i;
for (;;) {
if (df >= 0) i = 0; /* Initialize a level */
else { /* Print a node */
register int j;
for (j = 0; j < (totaldf - 1); ++j)
printf ("%d ", value [j]);
printf ("%d\n", value [j]);
if (df >= totaldf - 1) return; /* Probably can't happen */
i = value [++df] + 1; /* Advance to next node */
}
for (;;) { /* Finish a level */
if (i < nlev) {
value [df--] = i;
break;
}
if (df >= (totaldf - 1)) return;
i = value [++df] + 1; /* Advance to next node */
}
}
}
grdgen(df,nl)
int df; /* no. of degrees of freedom */
int nl; /* no. of levels each df */
{
totaldf=df;
nlev = nl;
grdgen_aux (df - 1);
}
Once again, there are branches that can't be taken, this time unless
the number of degrees of freedom is <=0, and there are needlessly
complex control structures. We can clean those up as before:
PROGRAM 9.
[same main program]
static int nlev;
static int totaldfm1;
static int value [MAXDF]; /* to hold current values of all df's */
grdgen_aux (df) {
int i;
for (;;) {
if (df >= 0) i = 0; /* Initialize a level */
else { /* Print a node */
register int j;
for (j = 0; j < totaldfm1; ++j)
printf ("%d ", value [j]);
printf ("%d\n", value [j]);
i = value [++df] + 1; /* Advance to next node */
}
while (i >= nlev) {
if (df >= totaldfm1) return;
i = value [++df] + 1; /* Advance to next node */
}
value [df--] = i;
}
}
grdgen(df,nl)
int df; /* no. of degrees of freedom */
int nl; /* no. of levels each df */
{
if (df <= 0) return; /* Sanity check */
totaldfm1 = df - 1;
nlev = nl;
grdgen_aux (df - 1);
}
And, just as with Program 7, we bring _grdgen_aux_ inline, and do some
more trivial optimization:
PROGRAM 10.
[still the same main program]
grdgen (df, nlev)
int df; /* no. of degrees of freedom */
register int nlev; /* no. of levels each df */
{
register int totaldfm1;
register int i;
register int j;
static int value [MAXDF]; /* to hold current values of all df's */
if (df <= 0) return; /* Sanity check */
totaldfm1 = df - 1;
for (;;) {
if (df >= 0) i = 0; /* Initialize a level */
else { /* Print a node */
for (j = 0; j < totaldfm1; ++j)
printf ("%d ", value [j]);
printf ("%d\n", value [j]);
++i; /* Advance to next node */
df = 0;
}
while (i >= nlev) {
if (df >= totaldfm1) return;
i = value [++df] + 1; /* Advance to next node */
}
value [df--] = i;
}
}
CONCLUSIONS.
We have compe up with three `final versions:' Programs 4, 7,
and 10. How do we choose among them.
Program 4 is likely to compile to the smallest amount of code,
as none of the program text has been duplicated. It also does not use
any control structures that can't be coded in Pascal without a GO TO;
this is useful in places where Pascal is the language of choice and GO
TOs are verboten.
The disadvantage of Program 4 is that the control structure
remains somewhat unclear; the `spaghetti' _gotos_ have simply been
replaced with assignments to a state variable, and the program state
still shifts in a way that's difficult to describe.
Programs 7 and 10 don't have that problem quite as much, since
the control flows from top to bottom. In my opinion, Program 7 wins
somewhat on readability, since its action is divided into several
distinct phases, that correspond to intuitive actions like `begin a
new level.' Program 10 has a somewhat simpler control structure,
though, and some may prefer it.
If performance is the object (and performance may often be
improved substantially by eliminating recursion), then I'd code all
three, instrument them, and use the fastest. There is no substitute
for the actual data from an instrumented program if you're after
speed.
Kevin Kenny UUCP: {ihnp4,pur-ee,convex}!uiucdcs!kenny
Department of Computer Science ARPA Internet or CSNet: kenny at CS.UIUC.EDU
University of Illinois BITNET: kenny at UIUCDCS.BITNET
1304 W. Springfield Ave.
Urbana, Illinois, 61801 Voice: (217) 333-6680
More information about the Comp.lang.c
mailing list