[WIP] Solution to http://stackoverflow.com/questions/7304511/partial-derivatives

main.cpp

#include <cstring>
#include <numeric>
#include <functional>

/**
 * Hackers Delight: http://www.hackersdelight.org/HDcode/nlz.c.txt
 */
int pop(unsigned x) {
   x = x - ((x >> 1) & 0x55555555);
   x = (x & 0x33333333) + ((x >> 2) & 0x33333333);
   x = (x + (x >> 4)) & 0x0F0F0F0F;
   x = x + (x << 8);
   x = x + (x << 16);
   return x >> 24;
}

/**
 * Counts the number of leading zeros of an unsigned int
 * Hackers Delight: http://www.hackersdelight.org/HDcode/nlz.c.txt
 */
int nlz(unsigned x) {
   int pop(unsigned x);

   x = x | (x >> 1);
   x = x | (x >> 2);
   x = x | (x >> 4);
   x = x | (x >> 8);
   x = x | (x >>16);
   return pop(~x);
}

// This is set globally for purposes of demonstration only
// the actual code has most of these methods implemented
// in classes and such
size_t _sizeOfElements = 9;
size_t _sizeOfEachDimension[] = {3, 3};
size_t _numberOfDimensions = 2;

size_t indexFromPoint(const size_t* pnt)
{
  size_t index = pnt[0];
  size_t lastDimMemSize = 1;
  size_t lastDimLength = _sizeOfEachDimension[0];
  //if(pnt[0] >= _sizeOfEachDimension[0])
  //  throw std::out_of_range("Calculated index will be out of range.");

  // [1, n-1], case 0 was done above.
  for (size_t dim = 1; dim < _numberOfDimensions; dim++)
  {
    //if(pnt[dim] >= _sizeOfEachDimension[dim])
    //  throw std::out_of_range("Calculated index will be out of range.");
    size_t thisDimMemSize = lastDimMemSize * lastDimLength;
    index += thisDimMemSize * pnt[dim];
    lastDimMemSize = thisDimMemSize;
    lastDimLength = _sizeOfEachDimension[dim];
  }

  return index;
}

size_t* pointFromIndex(size_t index)
{
  size_t* point = new size_t[_numberOfDimensions];
  memset(point, 0, sizeof(size_t) * _numberOfDimensions);

  //if(index >= _indexUpperBound)
  //  throw std::out_of_range("Index is too large to be within bounds.");

  for(size_t dim = _numberOfDimensions; dim-- > 0;)
  {
    size_t memSizeCurrentDim = std::accumulate(_sizeOfEachDimension,
                                               _sizeOfEachDimension + dim, 1,
                                               std::multiplies<size_t>());
    if(dim == 0)
      memSizeCurrentDim = 1;

    point[dim] = (index / memSizeCurrentDim);
    index = index % memSizeCurrentDim;
  }

  return point;
}

size_t* pointBefore(const size_t* pnt, size_t dimensionFlags)
{
  size_t dimension = (sizeof(size_t) * 8) - nlz(dimensionFlags) - 1;
  size_t* newPoint = new size_t[_numberOfDimensions];
  memcpy(newPoint, pnt, sizeof(size_t) * _numberOfDimensions);
  if(newPoint[dimension] != 0)
    newPoint[dimension]--;
  return newPoint;
}

size_t* pointAfter(const size_t* pnt, size_t dimensionFlags)
{
  size_t dimension = (sizeof(size_t) * 8) - nlz(dimensionFlags) - 1;
  size_t* newPoint = new size_t[_numberOfDimensions];
  memcpy(newPoint, pnt, sizeof(size_t) * _numberOfDimensions);
  if(newPoint[dimension] != _sizeOfEachDimension[dimension] - 1)
    newPoint[dimension]++;
  return newPoint;
}

double partialDerivative(double f1, double f2)
{
  double h = 1.0; // The one is only valid for this limited case
  return (f1 - f2) / (h * 2);
}

double partialDerivative(const double* points, const size_t* pnt, size_t dimension)
{
  size_t* beforePnt = pointBefore(pnt, dimension);
  size_t beforeIndex = indexFromPoint(beforePnt);
  double beforeValue = points[beforeIndex];
  delete beforePnt;

  size_t* afterPnt = pointAfter(pnt, dimension);
  size_t afterIndex = indexFromPoint(afterPnt);
  double afterValue = points[afterIndex];
  delete afterPnt;

  return partialDerivative(afterValue, beforeValue);
}

size_t topBit(size_t x)
{
  size_t numberOfLeadingZeros = nlz(x);
  // Moves the 1 bit to the same position it was in x,
  // while keeping all other bits 0
  size_t result = 0x80000000 >> numberOfLeadingZeros;
  return result;
}

/**
 * Calculates a partial derivative for a given point by
 * using the following bit masking table.
 *
 * w z y x |  Der. | Derivative Selector (derivSelector)
 * --------|-------|------------------------------------
 * 0 0 0 1 | dx    | 1
 * 0 0 1 0 | dy    | 2
 * 0 0 1 1 | dyx   | 3
 * 0 1 0 0 | dz    | 4
 * 0 1 0 1 | dzx   | 5
 * 0 1 1 0 | dzy   | 6
 * 0 1 1 1 | dzyx  | 7
 * 1 0 0 0 | dw    | 8
 * 1 0 0 1 | dwx   | 9
 * 1 0 1 0 | dwy   | 10
 * 1 0 1 1 | dwyx  | 11
 * 1 1 0 0 | dwz   | 12
 * 1 1 0 1 | dwzx  | 13
 * 1 1 1 0 | dwzy  | 14
 * 1 1 1 1 | dwzyx | 15
 * --------|-------|------------------------------------
 */
double partial_der(size_t* point, size_t derivSelector)
{
  // 32 - number of leading zeros
  size_t neighborSelector = topBit(derivSelector);
  size_t valueSelector = derivSelector ^ neighborSelector;

  size_t* prevPoint = pointBefore(point, neighborSelector);
  size_t* nextPoint = pointAfter(point, neighborSelector);

  if(valueSelector == 0)
  {
    size_t nextIndex = indexFromPoint(nextPoint);
    size_t prevIndex = indexFromPoint(prevPoint);
    delete prevPoint;
    delete nextPoint;
    return partialDerivative(depVariables[nextIndex], depVariables[prevIndex]);
  }
  else
  {
    double prevPartial = partial_der(prevPoint, valueSelector);
    double nextPartial = partial_der(nextPoint, valueSelector);
    delete prevPoint;
    delete nextPoint;
    return partialDerivative(nextPartial, prevPartial);
  }
}

void main()
{
  // This is the same as the column f(x, y)
  double depVariables[] = {4, 3, 0, 3, 2, -1, 0, -1, -4};

  // This is where the resulting values should be stored
  double dx[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
  double dy[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
  double dyx[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};

  // Create an array of them so its easy to address them all
  double* derivatives[] = {depVariables, dx, dy, dyx};

  size_t* point = new size_t[_numberOfDimensions];
  memset(point, 0, sizeof(size_t) * _numberOfDimensions);

  // 2^2 = 4
  for(size_t der = 0; der < 4 - 1; der++)
  {
    double* data = derivatives[der];
    for(size_t dim = 0; dim < _numberOfDimensions; dim++)
    {
      // loop to 9 since that is the size of depVariables
      for(size_t index = 0; index < 9; index++)
      {
        size_t* point = pointFromIndex(index);
        derivatives[dim + 1][index] = partialDerivative(data, point, dim);
        delete point;
      }
    }
  }

  delete [] point;
}