#pragma once

#include <algorithm>
#include <cmath>
#include <cstdint>

#include "core/image.hpp"

namespace effector {

// Reads the colour at a fractional position, mixing the four neighbouring
// pixels; positions off the image take the nearest edge pixel. The weights are
// cut to 1/256 steps so the mix itself is integer arithmetic.
inline void sampleBilinear(const Image& image, double x, double y, std::uint8_t out[3]) {
    const int maxX = image.width() - 1;
    const int maxY = image.height() - 1;
    x = std::clamp(x, 0.0, static_cast<double>(maxX));
    y = std::clamp(y, 0.0, static_cast<double>(maxY));

    const int x0 = static_cast<int>(x);
    const int y0 = static_cast<int>(y);
    const int x1 = std::min(x0 + 1, maxX);
    const int y1 = std::min(y0 + 1, maxY);
    const int fx = static_cast<int>((x - x0) * 256.0);
    const int fy = static_cast<int>((y - y0) * 256.0);

    const std::uint8_t* a = image.row(y0) + static_cast<std::size_t>(x0) * 3;
    const std::uint8_t* b = image.row(y0) + static_cast<std::size_t>(x1) * 3;
    const std::uint8_t* c = image.row(y1) + static_cast<std::size_t>(x0) * 3;
    const std::uint8_t* d = image.row(y1) + static_cast<std::size_t>(x1) * 3;
    for (int ch = 0; ch < 3; ++ch) {
        const int top = a[ch] * (256 - fx) + b[ch] * fx;
        const int bottom = c[ch] * (256 - fx) + d[ch] * fx;
        out[ch] = static_cast<std::uint8_t>((top * (256 - fy) + bottom * fy + 32768) >> 16);
    }
}

// Sine and cosine of an angle in degrees from nothing but additions and
// multiplications, so a rotation lands on the same pixels whatever the
// machine's libm makes of std::sin. The angle is folded into one octant and a
// short Taylor series does the rest, far below a pixel's worth of error.
inline void sinCosDegrees(double degrees, double& sine, double& cosine) {
    const double turns = degrees / 90.0;
    const double quadrant = std::floor(turns + 0.5);
    const double x = (turns - quadrant) * 1.5707963267948966;  // [-pi/4, pi/4]
    const double x2 = x * x;
    const double s =
        x * (1.0 - x2 / 6.0 * (1.0 - x2 / 20.0 * (1.0 - x2 / 42.0 * (1.0 - x2 / 72.0))));
    const double c =
        1.0 - x2 / 2.0 * (1.0 - x2 / 12.0 * (1.0 - x2 / 30.0 * (1.0 - x2 / 56.0 * (1.0 - x2 / 90.0))));
    switch (static_cast<int>(quadrant - 4.0 * std::floor(quadrant / 4.0))) {
        case 0: sine = s; cosine = c; break;
        case 1: sine = c; cosine = -s; break;
        case 2: sine = -s; cosine = -c; break;
        default: sine = -c; cosine = s; break;
    }
}

// A uniform double in [0, 1) from the top 53 bits of one rng draw.
inline double unitDouble(std::uint64_t bits) {
    return static_cast<double>(bits >> 11) * (1.0 / 9007199254740992.0);
}

// Whether an event with the given probability happens, from one rng draw.
inline bool chance(std::uint64_t bits, double probability) {
    return unitDouble(bits) < probability;
}

}  // namespace effector