#include <cmath>
#include <vector>
#include <string>
#include <limits>
#include <iomanip>
#include <iostream>
// glad.h is generated by the following command:
// glad --out-path=. --generator=c --omit-khrplatform --api="gl=3.3" --profile=core --extensions=
#include "glad/glad.h"
#include <GL/freeglut.h>
#include <glm/glm.hpp>
using glm::vec4;

GLuint vao, vbo;
GLuint texFBO;
GLuint inputTextures[2];
GLuint program;
GLuint fbo;
int width=1, height=2;

void printShaderOutput(int texW, int texH)
{
    glActiveTexture(GL_TEXTURE0);
    glBindTexture(GL_TEXTURE_2D, texFBO);

    std::vector<vec4> data(texW*texH);
    glGetTexImage(GL_TEXTURE_2D, 0, GL_RGBA, GL_FLOAT, data.data());
    std::cout << "a,b,sum,product,relError(sum),relError(product),note\n";
    for(int i=0;i<width;++i)
    {
        const auto a=double(data[i+width*0].x)+double(data[i+width*0].y);
        const auto b=double(data[i+width*0].z)+double(data[i+width*0].w);
        const auto sum=double(data[i+width*1].x)+double(data[i+width*1].y);
        const auto product=double(data[i+width*1].z)+double(data[i+width*1].w);
        const auto trueSum=a+b, trueProduct=a*b;
        const auto sumErr=(sum-trueSum)/trueSum;
        const auto prodErr=(product-trueProduct)/trueProduct;
        std::cout << std::setprecision(std::numeric_limits<double>::max_digits10)
                  << a << ',' << b << ','
                  << sum << ',' << product << ','
                  << std::setprecision(3)
                  << sumErr << ',' << prodErr << ','
                  << (std::abs(sumErr)>1e-14 || std::abs(prodErr)>1e-14 ? "WARN" : "OK")
                  << '\n';
    }
    std::cout.flush();
}

GLuint makeShader(GLenum type, std::string const& srcStr)
{
    const auto shader=glCreateShader(type);
    const GLint srcLen=srcStr.size();
    const GLchar*const src=srcStr.c_str();
    glShaderSource(shader, 1, &src, &srcLen);
    glCompileShader(shader);
    GLint status=-1;
    glGetShaderiv(shader, GL_COMPILE_STATUS, &status);
    assert(glGetError()==GL_NO_ERROR);
    assert(status);
    return shader;
}

void loadShaders()
{
    program=glCreateProgram();

    const auto vertexShader=makeShader(GL_VERTEX_SHADER, 1+R"(
#version 330
in vec4 vertex;
void main() { gl_Position=vertex; }
)");
    glAttachShader(program, vertexShader);

    const auto fragmentShader=makeShader(GL_FRAGMENT_SHADER, 1+R"(
#version 330
#extension GL_ARB_gpu_shader5 : require

// Emulation based on Fortran-90 double-single package. See http://crd-legacy.lbl.gov/~dhbailey/mpdist/
// Add: res = ds_add(a, b) => res = a + b
vec2 ds_add(vec2 dsa, vec2 dsb)
{
    precise float t1 = dsa.x + dsb.x;
    precise float e = t1 - dsa.x;
    precise float t2 = ((dsb.x - e) + (dsa.x - (t1 - e))) + dsa.y + dsb.y;

    precise vec2 dsc;
    dsc.x = t1 + t2;
    dsc.y = t2 - (dsc.x - t1);
    return dsc;
}

// Multiply: res = ds_mul(a, b) => res = a * b
vec2 ds_mul(vec2 dsa, vec2 dsb)
{
    const float split = 8193.;

    precise float cona = dsa.x * split;
    precise float conb = dsb.x * split;
    precise float a1 = cona - (cona - dsa.x);
    precise float b1 = conb - (conb - dsb.x);
    precise float a2 = dsa.x - a1;
    precise float b2 = dsb.x - b1;

    precise float c11 = dsa.x * dsb.x;
    precise float c21 = a2 * b2 + (a2 * b1 + (a1 * b2 + (a1 * b1 - c11)));

    precise float c2 = dsa.x * dsb.y + dsa.y * dsb.x;

    precise float t1 = c11 + c2;
    precise float e = t1 - c11;
    precise float t2 = dsa.y * dsb.y + ((c2 - e) + (c11 - (t1 - e))) + c21;

    precise vec2 dsc;
    dsc.x = t1 + t2;
    dsc.y = t2 - (dsc.x - t1);

    return dsc;
}

uniform sampler1D opA, opB;
uniform float rtWidth;
out vec4 color;

void main()
{
    float texCoord=gl_FragCoord.x/rtWidth;
    vec2 a=texture(opA,texCoord).xy;
    vec2 b=texture(opB,texCoord).xy;
    if(gl_FragCoord.y<1)   // first row
        color=vec4(a,b);
    else if(gl_FragCoord.y<2)   // second row
        color=vec4(ds_add(a,b),ds_mul(a,b));
}

)");
    glAttachShader(program, fragmentShader);

    glLinkProgram(program);
    GLint status=0;
    glGetProgramiv(program, GL_LINK_STATUS, &status);
    assert(glGetError()==GL_NO_ERROR);
    assert(status);

    glDetachShader(program, fragmentShader);
    glDeleteShader(fragmentShader);

    glDetachShader(program, vertexShader);
    glDeleteShader(vertexShader);
}

void setupBuffers()
{
    glGenVertexArrays(1, &vao);
    glBindVertexArray(vao);
    glGenBuffers(1, &vbo);
    glBindBuffer(GL_ARRAY_BUFFER, vbo);
    const GLfloat vertices[]=
    {
        -1, -1,
         1, -1,
        -1,  1,
         1,  1,
    };
    glBufferData(GL_ARRAY_BUFFER, sizeof vertices, vertices, GL_STATIC_DRAW);
    constexpr GLuint attribIndex=0;
    constexpr int coordsPerVertex=2;
    glVertexAttribPointer(attribIndex, coordsPerVertex, GL_FLOAT, false, 0, 0);
    glEnableVertexAttribArray(attribIndex);
    glBindVertexArray(0);
}

void setupRenderTarget()
{
    glGenTextures(1, &texFBO);
    glGenFramebuffers(1,&fbo);
}

void setupInputTextures()
{
    glGenTextures(2, inputTextures);
    {
        glBindTexture(GL_TEXTURE_1D, inputTextures[0]);
        glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
        glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
#define SPLIT_DOUBLE_TO_FLOATS(x) GLfloat(x),GLfloat(x-GLfloat(x))
        static constexpr GLfloat data[]=
        {
            SPLIT_DOUBLE_TO_FLOATS(3.1415926535897932),
        };
        glTexImage1D(GL_TEXTURE_1D, 0, GL_RGBA32F, sizeof data/sizeof*data/2, 0, GL_RG, GL_FLOAT, data);
    }
    {
        glBindTexture(GL_TEXTURE_1D, inputTextures[1]);
        glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
        glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
        static constexpr GLfloat data[]=
        {
            SPLIT_DOUBLE_TO_FLOATS(2.7182818284590452),
        };
        glTexImage1D(GL_TEXTURE_1D, 0, GL_RGBA32F, sizeof data/sizeof*data/2, 0, GL_RG, GL_FLOAT, data);
    }
}

bool init()
{
    if(!gladLoadGL())
    {
        std::cerr << "Failed to initialize GLAD\n";
        return false;
    }
    if(!GLAD_GL_VERSION_3_3)
    {
        std::cerr << "OpenGL 3.3 not supported\n";
        return false;
    }

    setupRenderTarget();
    setupInputTextures();
    loadShaders();
    setupBuffers();

    glViewport(0,0,width,height);

    glBindTexture(GL_TEXTURE_2D,texFBO);
    glTexImage2D(GL_TEXTURE_2D,0,GL_RGBA32F,width,height,0,GL_RGBA,GL_UNSIGNED_BYTE,nullptr);
    glBindTexture(GL_TEXTURE_2D,0);
    glBindFramebuffer(GL_FRAMEBUFFER,fbo);
    glFramebufferTexture2D(GL_FRAMEBUFFER,GL_COLOR_ATTACHMENT0,GL_TEXTURE_2D,texFBO,0);
    const auto status=glCheckFramebufferStatus(GL_FRAMEBUFFER);
    assert(status==GL_FRAMEBUFFER_COMPLETE);
    glBindFramebuffer(GL_FRAMEBUFFER,0);

    return true;
}

void display()
{
    const static bool inited=init();
    if(!inited) std::exit(1);

    glBindFramebuffer(GL_FRAMEBUFFER,fbo);

    glUseProgram(program);
    glActiveTexture(GL_TEXTURE0);
    glBindTexture(GL_TEXTURE_1D, inputTextures[0]);
    glActiveTexture(GL_TEXTURE1);
    glBindTexture(GL_TEXTURE_1D, inputTextures[1]);
    glUniform1i(glGetUniformLocation(program,"opA"),0);
    glUniform1i(glGetUniformLocation(program,"opB"),1);
    glUniform1f(glGetUniformLocation(program,"rtWidth"),width);

    glBindVertexArray(vao);
    glDrawArrays(GL_TRIANGLE_STRIP, 0, 4);
    glBindVertexArray(0);

    printShaderOutput(width, height);
    std::exit(0);

    glFinish();
}

int main(int argc, char** argv)
{
    glutInitContextVersion(3,3);
    glutInitContextProfile(GLUT_CORE_PROFILE);
    glutInit(&argc, argv);

    glutInitDisplayMode(GLUT_RGB);

    glutInitWindowSize(width, height);
    glutCreateWindow("Test");
    glutDisplayFunc(display);

    glutMainLoop();
}