/*
 * Copyright 1993-2015 NVIDIA Corporation.  All rights reserved.
 *
 * Please refer to the NVIDIA end user license agreement (EULA) associated
 * with this source code for terms and conditions that govern your use of
 * this software. Any use, reproduction, disclosure, or distribution of
 * this software and related documentation outside the terms of the EULA
 * is strictly prohibited.
 *
 */

#include "Mandelbrot_gold.h"

#define ABS(n) ((n) < 0 ? -(n) : (n))

/* dfloat class declaration */
class dfloat {
private:
    float val[2];
public:
    dfloat() { val[0] = val[1] = 0; }
    dfloat(float a,float b) { val[0] = a; val[1] = b; }
    dfloat(double b);
    inline float operator[](unsigned idx) const { return val[idx]; }
};
inline dfloat operator+(const dfloat &dsa, const dfloat &dsb);
inline dfloat operator-(const dfloat &dsa, const dfloat &dsb);
inline dfloat operator*(const dfloat &dsa, const dfloat &dsb);
inline int operator<(const dfloat &a, float b) {return a[0]<b;}



// The core Mandelbrot calculation function template
template<class T>
inline int CalcMandelbrot(const T xPos, const T yPos, const T xJParam, const T yJParam, const int crunch,
                          const bool isJulia)
{
    T x, y, xx, yy, xC, yC;
    int i = crunch;

    if (isJulia)
    {
        xC = xJParam;
        yC = yJParam;
        y = yPos;
        x = xPos;
        yy = y * y;
        xx = x * x;
    }
    else
    {
        xC = xPos;
        yC = yPos;
        x = y = 0;
        xx = yy = 0;
    }

    while (--i && (xx + yy < 4.0f))
    {
    	y = x * y +  x * y + yC ;
        x = xx - yy + xC ;
        yy = y * y;
        xx = x * x;
    }

    return i;
} // CalcMandelbrot


inline void updatePixel(uchar4 &dst, const uchar4 &color, int frame)
{
    int frame1 = frame + 1;
    int frame2 = frame1 / 2;
    dst.x = (dst.x * frame + color.x + frame2) / frame1;
    dst.y = (dst.y * frame + color.y + frame2) / frame1;
    dst.z = (dst.z * frame + color.z + frame2) / frame1;

}

inline void setColor(uchar4 &dst, const uchar4 &colors, int &m, const int animationFrame)
{
    if (m == 0)
    {
        dst.x = 0;
        dst.y = 0;
        dst.z = 0;
        return;
    }

    m += animationFrame;
    dst.x = m * colors.x;
    dst.y = m * colors.y;
    dst.z = m * colors.z;
}

template<class T,class T_>
void runMandelbrotGold0(uchar4 *dst, const int imageW, const int imageH, const int crunch, const  T xOff, const T yOff,
        const T xJParam, const T yJParam, const T scale, const uchar4 colors, const int frame,
        const int animationFrame, const bool isJulia)
{
    for (int iy = 0; iy < imageH; iy++)
        for (int ix = 0; ix < imageW; ix++)
        {
            // Calculate the location
            const T_ xPos = (T)ix * scale + xOff;
            const T_ yPos = (T)iy * scale + yOff;

            // Calculate the Mandelbrot index for the current location
            int m = CalcMandelbrot<T_>(xPos, yPos, xJParam, yJParam, crunch, isJulia);
            m = m > 0 ? crunch - m : 0;

            // Convert the Mandelbrot index into a color
            uchar4 color;

            setColor (color, colors, m, animationFrame);

            // Output the pixel
            int pixel = imageW * iy + ix;

            if (frame == 0)
            {
                color.w = 0;
                dst[pixel] = color;
            }
            else
                updatePixel(dst[pixel], color, frame);
        }

} //runMandelbrotGold0_


// Determine if two pixel colors are within tolerance
inline int CheckColors(const uchar4 &color0, const uchar4 &color1)
{
    int x = color1.x - color0.x;
    int y = color1.y - color0.y;
    int z = color1.z - color0.z;
    return (ABS(x) > 10) || (ABS(y) > 10) || (ABS(z) > 10);
} // CheckColors


template<class T,class T_>
void runMandelbrotGold1(uchar4 *dst, const int imageW, const int imageH, const int crunch, const T xOff, const T yOff,
                          const T xJParam, const T yJParam, const T scale, const uchar4 colors, const int frame,
                          const int animationFrame, const bool isJulia)
{
    for (int iy = 0; iy < imageH; iy++)
        for (int ix = 0; ix < imageW; ix++)
        {
            // Get the current pixel color
            int pixel = imageW * iy + ix;
            uchar4 pixelColor = dst[pixel];
            int count = 0;

            // Search for pixels out of tolerance surrounding the current pixel
            if (ix > 0)
                count += CheckColors(pixelColor, dst[pixel - 1]);

            if (ix + 1 < imageW)
                count += CheckColors(pixelColor, dst[pixel + 1]);

            if (iy > 0)
                count += CheckColors(pixelColor, dst[pixel - imageW]);

            if (iy + 1 < imageH)
                count += CheckColors(pixelColor, dst[pixel + imageW]);

            if (count)
            {
                // Calculate the location
                const T_ xPos = (T)ix * scale + xOff;
                const T_ yPos = (T)iy * scale + yOff;

                // Calculate the Mandelbrot index for the current location
                int m = CalcMandelbrot<T_>(xPos, yPos, xJParam, yJParam, crunch, isJulia);
                m = m > 0 ? crunch - m : 0;

                // Convert the Mandelbrot index into a color
                uchar4 color;

                setColor (color, colors, m, animationFrame);

                // Output the pixel
                updatePixel(dst[pixel], color, frame);
            }
        }
} // RunMandelbrotGold1_



/* Implementation of exported functions */
void RunMandelbrotGold1(uchar4 *dst, const int imageW, const int imageH, const int crunch, const float xOff, const float yOff,
                        const float xJParam, const float yJParam, const float scale, const uchar4 colors, const int frame,
                        const int animationFrame, const bool isJulia)
{
    runMandelbrotGold1<float,float> (dst, imageW, imageH, crunch,
                                     xOff, yOff, xJParam, yJParam, scale, colors, frame,
                                     animationFrame, isJulia);

} // RunMandelbrotGold1

void RunMandelbrotDSGold1(uchar4 *dst, const int imageW, const int imageH, const int crunch, const double xOff, const double yOff,
                          const double xJParam, const double yJParam, const double scale, const uchar4 colors, const int frame,
                          const int animationFrame, const bool isJulia)
{
    runMandelbrotGold1<double,dfloat> (dst, imageW, imageH, crunch,
                                       xOff, yOff, xJParam, yJParam, scale, colors, frame,
                                       animationFrame, isJulia);

} // RunMandelbrotDSGold1

void RunMandelbrotGold0(uchar4 *dst, const int imageW, const int imageH, const int crunch, const float xOff, const float yOff,
                        const float xJParam, const float yJParam, const float scale, const uchar4 colors, const int frame,
                        const int animationFrame, const bool isJulia)
{
    runMandelbrotGold0<float,float> (dst, imageW, imageH, crunch,
                                     xOff, yOff, xJParam, yJParam, scale, colors, frame,
                                     animationFrame, isJulia);
} // RunMandelbrotGold0

void RunMandelbrotDSGold0(uchar4 *dst, const int imageW, const int imageH, const int crunch, const double xOff, const double yOff,
                          const double xJParam, const double yJParam, const double scale, const uchar4 colors, const int frame,
                          const int animationFrame, const bool isJulia)
{
    runMandelbrotGold0<double,dfloat> (dst, imageW, imageH, crunch,
                                       xOff, yOff, xJParam, yJParam, scale, colors, frame,
                                       animationFrame, isJulia);
} // RunMandelbrotDSGold0


/*dfloat operations implementation */

/* Construct a DS number equal to a double precision floating point number b*/
dfloat::dfloat(double b)
{
    val[0] = (float)b;
    val[1] = (float)(b-val[0]);
}

inline dfloat operator+(const dfloat &dsa, const dfloat &dsb)
{
    // Compute dsa + dsb using Knuth's trick.
    float t1 = dsa[0] + dsb[0];
    float e = t1 - dsa[0];
    float t2 = ((dsb[0] - e) + (dsa[0] - (t1 - e))) + dsa[1] + dsb[1];

    // The result is t1 + t2, after normalization.
    e = t1 + t2;
    return dfloat(e, t2- ( e - t1));
}

inline dfloat operator-(const dfloat &dsa, const dfloat &dsb)
{
    // Compute dsa - dsb using Knuth's trick.
    float t1 = dsa[0] - dsb[0];
    float e = t1 - dsa[0];
    float t2 = ((-dsb[0] - e) + (dsa[0] - (t1 - e))) + dsa[1] - dsb[1];

    // The result is t1 + t2, after normalization.
    e = t1 + t2;
    return dfloat(e, t2 - (e - t1));
}

inline dfloat operator*(const dfloat &dsa, const dfloat &dsb)
{
    // This splits dsa(1) and dsb(1) into high-order and low-order words.
    float c11 = dsa[0] * dsb[0];
    float c21 = dsa[0] * dsb[0] - c11;

    // Compute dsa[0] * dsb[1] + dsa[1] * dsb[0] (only high-order word is needed).
    float c2 = dsa[0] * dsb[1] + dsa[1] * dsb[0];

    // Compute (c11, c21) + c2 using Knuth's trick, also adding low-order product.
    float t1 = c11 + c2;
    float e = t1 - c11;
    float t2 = ((c2 - e) + (c11 - (t1 - e))) + c21 + dsa[1] * dsb[1];

    // The result is t1 + t2, after normalization.
    e = t1 + t2;
    return dfloat (e, t2 - (e - t1));
}