NRook.h

Go to the documentation of this file.

#ifndef _ZFXMATH_INCLUDE_NROOK_H_
#define _ZFXMATH_INCLUDE_NROOK_H_

namespace ZFXMath
{

namespace SphericalHarmonics
{

    template<class PrecisionType>
    SphericalHarmonics::TSample<PrecisionType>* GenerateSamples( const int nb_samples )
    {
        // Generate nrook cells
        int* cells = new int[nb_samples];
        for (int i = 0; i < nb_samples; i++)
            cells[i] = i;
        MultiStageNRooks<PrecisionType>(nb_samples, cells);

        // Allocate the samples list
        SphericalHarmonics::TSample<PrecisionType>* samples = new SphericalHarmonics::TSample<PrecisionType>[nb_samples];

        for (int i = 0; i < nb_samples; i++)
        {
            // Generate uniform random sample
            PrecisionType x = (i + drand48<PrecisionType>()) / nb_samples;
            PrecisionType y = (cells[i] + drand48<PrecisionType>()) / nb_samples;
            
            // Generate spherical/cartesian co-ordinates
            samples[i] = SphericalHarmonics::TSample<PrecisionType>( PrecisionType( 2 * acos(sqrt(1 - x)) ), PrecisionType( 2 * PI * y ) );
        }

        delete [] cells;

        return samples;
    }

    template<class PrecisionType>
    PrecisionType drand48()
    {
        return ((PrecisionType)rand() / RAND_MAX);
    }

    int lrand48()
    {
        return (rand());
    }

    template<class PrecisionType>
    void MultiStageNRooks(const int size, int* cells)
    {
        if (size == 1)
            return;

        int size1 = size >> 1;
        int size2 = size >> 1;

        if (size & 1)
        {
            if (drand48<PrecisionType>() > 0.5)
                size1++;
            else
                size2++;
        }

        int* upper_cells = new int[size1];
        int* lower_cells = new int[size2];

        for (int i = 0, j = 0; i < size - 1; i += 2, j++)
        {
            if (lrand48() & 1)
            {
                upper_cells[j] = cells[i];
                lower_cells[j] = cells[i + 1];
            }
            else
            {
                upper_cells[j] = cells[i + 1];
                lower_cells[j] = cells[i];
            }
        }

        if (size1 != size2)
        {
            if (size1 > size2)
                upper_cells[j] = cells[i];
            else
                lower_cells[j] = cells[i];
        }

        MultiStageNRooks<PrecisionType>(size1, upper_cells);
        memcpy(cells, upper_cells, size1 * sizeof(int));
        delete [] upper_cells;

        MultiStageNRooks<PrecisionType>(size2, lower_cells);
        memcpy(cells + size1, lower_cells, size2 * sizeof(int));
        delete [] lower_cells;
    }

}; // namespace SphericalHarmonics

}; // namespace ZFXMath

#endif //_ZFXMATH_INCLUDE_NROOK_H__

---