#include <math.h>

typedef unsigned char byte;
typedef unsigned short word;

typedef byte byte4 __attribute__ ((ext_vector_type(4)));
typedef byte byte8 __attribute__ ((ext_vector_type(8)));
typedef word word4 __attribute__ ((ext_vector_type(4)));
typedef word word8 __attribute__ ((ext_vector_type(8)));
typedef int int4 __attribute__ ((ext_vector_type(4)));
typedef float float4 __attribute__ ((ext_vector_type(4)));

#include  <immintrin.h>

static inline float4 rsqrt4(const float4 a) {
    return _mm_rsqrt_ps(a);
}

// целочисленное 16-битное умножение 4 пикселей на 4 коэфф-та
static inline void mulbytes(byte8 * Src, byte8 * Dst, float4 k) {
   float4 MinV = {0.0, 0.0, 0.0, 0.0};
   float4 MaxV = {255.0, 255.0, 255.0, 255.0};

   float4 kLim = _mm_max_ps(_mm_min_ps(k, MaxV), MinV);
   word4 kw =  __builtin_convertvector(kLim, word4);
   // для r,g,b,a коэфф. одинаковый, раскладываем word4 на 2*word8
   word8 k1 = (word8){kw[0], kw[0], kw[0], kw[0], kw[1], kw[1], kw[1], kw[1]};
   word8 k2 = (word8){kw[2], kw[2], kw[2], kw[2], kw[3], kw[3], kw[3], kw[3]};
   word8 w1 =  __builtin_convertvector(Src[0], word8); // byte8 -> word8
   word8 w2 =  __builtin_convertvector(Src[1], word8);
   Dst[0] = __builtin_convertvector(w1 * k1 >> 8, byte8); // умножение
   Dst[1] = __builtin_convertvector(w2 * k2 >> 8, byte8); // и word8 -> byte8
}

// функция считает 1 строку
void DotBump4(byte4 * Src, byte * NMap, byte4 * Dst, float * LVec, int NMapStride, int Cnt)
{
    // x/y/z в карте нормалей заранее разложены в разные массивы 
    // и преобразованы во float, можно оставить в байтах, будет -10..20%
    float4 * nx = (float4*) NMap;
    float4 * ny = (float4*) (NMap + NMapStride);
    float4 * nz = (float4*) (NMap + NMapStride*2);

    float lxb = LVec[0], ly = LVec[1], lz = LVec[2];
    float4 inc = (float4) {0, 1, 2, 3};

    // считаем по 4 пикселя
    for (int i = 0; i < (Cnt / 4); i++)
    {
        float4 fnx = nx[i & 63]; // & зацикливает тестуру по размеру 256
        float4 fny = ny[i & 63];
        float4 fnz = nz[i & 63];

        // X источника света для 4-х элементов, Y,Z постоянны
        float4 lx = lxb - ((float)(i * 4) + inc);
        
        // dot product (косинус угла)
        float4 d = (fnx * lx) + (fny * ly) + (fnz * lz);
        // квадрат расстояния
        float4 dx = fnx - lx, dy = fny - ly, dz = fnz - lz;
        float4 d2 = (dx*dx + dy*dy + dz*dz);

        // свет = dot / дистанция
        d = d * rsqrt4(d2);
        // можно отрезать элементы, в которых 0 - но вроде работает и так
        /*
        float4 d1 = d * rsqrt4(d2);
        int4 c = ((d2 > 0) && (d > 0));
        d = (float4) ( ((int4)d1 & c) | ((int4)d & !c) );
        */
        // в d получили 4 коэфф-та света
        // умножаем на них 4 rgba пикселя 
        mulbytes((byte8*)&Src[(i*4) & 255], (byte8*)Dst, d);
        Dst+=4;
    }
}