How to make my function work for arbitrary dimensionality in C

Question

I have a math function that I would like to use in 1d, 2d, and 3d spaces. I wrote the code to work for the 3-d case and I can easily modify it for the other cases, however I was wondering if there is a way I can write this so that it works with any arbitrary dimensionality. Basically for every dimension there is another nested loop and I don't know how to make multiple nested loops without actually writing them out in the code.

void gen_kernel_Gk(double *Gk, double L, double b, int np, int Z, int m)
{
    // Generates Kernel Function G in k-space
    // Gk is the array of the function values, it has Z * m number of values
    // L is cuttoff distance
    // b is monomer size
    // np is the degree of polymerization
    // Z is the number of modes
    // m is the dimensionality of the real space

    int i, j, k, l, index, Zhalf;
    double pref, kvec[m], k2;

    pref = pow( 2 * PI / (double)L, 2) * pow( b, 2) * (double)np / 6.0;

    Zhalf = Z / 2;

    for( i = 0; i < Z; ++i)
    {
        if( i < Zhalf)
        {
            kvec[0] = (double) pow( i, 2);
        }
        else
        {
            kvec[0] = (double) pow ( Z - i ,2);
        }

        for( j = 0; j < Z; ++j)
        {
            if( i < Zhalf)
            {
                kvec[1] = (double) pow( j, 2);
            }
            else
            {
                kvec[1] = (double) pow ( Z - j ,2);
            }

            for( k = 0 ; k < Z ; ++k)
            {
                if( k < Zhalf)
                {
                    kvec[2] = (double) pow( k, 2);
                }
                else
                {
                    kvec[2] = (double) pow ( Z - k ,2);
                }

                for( k2 = 0, l = 0; l < m; ++l)
                {
                    k2 += kvec[l];
                }

                index = i * pow( Z, 2) + j * Z + k;
                Gk[index] = 2.0 * ( exp( -k2) + k2 - 1.0) / pow( k2, 2);
            } 
        }
    }
    return;
}

Answer 1

If you want an arbitrary number of nested for loops, you can go for a solution that uses recursion.

Exactly how to do it depends on your specific task. The example below isn't based exactly on your code but it's a more general code example to show the principle. Hopefully you can see the idea and apply it to your specific task.

The basic idea is that the recursive function needs arguments that describe stuff like the number of dimensions and the size of each dimension (in an array) and the current level/index. For solving a specific task, you probably need to add additional arguments.

Like:

#include <stdio.h>
#include <assert.h>

void opr(const unsigned dimensions, const unsigned * dimsize, 
         unsigned curdim, unsigned * curindex)
{
  assert(curdim < dimensions);
  for (unsigned u = 0; u < dimsize[curdim]; ++u)
  {
    // Save the current index
    curindex[curdim] = u;

    // Check if the last dimension has been reached
    if (curdim == (dimensions - 1))
    {
      // yes, last dimension has been reached - now do the operation you want...
      // Here the code just do a simple print
      printf("Handle array element at ");
      for (unsigned t = 0; t < dimensions; ++t)
      {
        printf("[%u]", curindex[t]);
      }
      printf("\n");
    }
    else
    {
      // no, last dimension has NOT been reached
      // so do a recursive call for the next dimension
      opr(dimensions, dimsize, curdim + 1, curindex);
                    // notice: ^^^^^^^^^^ increase current dimension 
    }
  }
}

int main(void)
{
  unsigned dimsizeA[] = {2, 3};  // [2][3]
  unsigned dimensionsA = sizeof dimsizeA / sizeof dimsizeA[0];
  unsigned curindexA[] = {0, 0};
  opr(dimensionsA, dimsizeA, 0, curindexA);

  printf("-----------------------------------\n");

  unsigned dimsizeB[] = {2, 3, 3, 2};  // [2][3][3][2]
  unsigned dimensionsB = sizeof dimsizeB / sizeof dimsizeB[0];
  unsigned curindexB[] = {0, 0, 0, 0};
  opr(dimensionsB, dimsizeB, 0, curindexB);

  return 0;
}

Output:

Handle array element at [0][0]
Handle array element at [0][1]
Handle array element at [0][2]
Handle array element at [1][0]
Handle array element at [1][1]
Handle array element at [1][2]
-----------------------------------
Handle array element at [0][0][0][0]
Handle array element at [0][0][0][1]
Handle array element at [0][0][1][0]
Handle array element at [0][0][1][1]
Handle array element at [0][0][2][0]
Handle array element at [0][0][2][1]
Handle array element at [0][1][0][0]
Handle array element at [0][1][0][1]
Handle array element at [0][1][1][0]
Handle array element at [0][1][1][1]
Handle array element at [0][1][2][0]
Handle array element at [0][1][2][1]
Handle array element at [0][2][0][0]
Handle array element at [0][2][0][1]
Handle array element at [0][2][1][0]
Handle array element at [0][2][1][1]
Handle array element at [0][2][2][0]
Handle array element at [0][2][2][1]
Handle array element at [1][0][0][0]
Handle array element at [1][0][0][1]
Handle array element at [1][0][1][0]
Handle array element at [1][0][1][1]
Handle array element at [1][0][2][0]
Handle array element at [1][0][2][1]
Handle array element at [1][1][0][0]
Handle array element at [1][1][0][1]
Handle array element at [1][1][1][0]
Handle array element at [1][1][1][1]
Handle array element at [1][1][2][0]
Handle array element at [1][1][2][1]
Handle array element at [1][2][0][0]
Handle array element at [1][2][0][1]
Handle array element at [1][2][1][0]
Handle array element at [1][2][1][1]
Handle array element at [1][2][2][0]
Handle array element at [1][2][2][1]