how do you fill the image? maybe you index it by x+y*height instead of x+y*width or something.

Yea, you are right. I switched width and height in my for-loops. That's kinda embarrassing.

But I think there is still something wrong with the rendering because this looks just wrong: https://picload.org/view/dapagoil/test.png.html

My intersection code:

double intersectSphere(const struct Sphere s, const struct Ray r) {
    double a = dot(r.dir, sub(r.origin, s.center));
    double delta1 =  a * a;
    double delta2 = length(sub(r.origin,s.center)) * length(sub(r.origin,s.center));
    double rsqr = s.radius * s.radius;

    double delta = delta1 - delta2 + rsqr;

    if(delta < 0)
        return -1.0;
    else {
        double left = -a;
        double d1 = left + sqrt(delta);
        double d2 = left - sqrt(delta);
        return fmin(d1,d2);
    }
};

(from here: https://en.wikipedia.org/wiki/Line–sphere_intersection)

I guess

fmin(d1,d2);

could be wrong?

Hard-code one specific ray, that hits something, and step in your debugger through the code. On CPU it's quite easy to figure out where the expected value diverges from the one generated by your source

fmin(d1,d2) sounds ok to me.

I don't understand how this could be helpfull? I would apply the same algorithms I implemented in my code, so there would be no difference in my debugged ray and the one the program calculates, right?

Maybe there is something wrong with the diffuse reflection, there is no way the lower sphere could be so bright on the bottom with the light being above it.

https://picload.org/view/dapopdwr/test.png.html

Edit: I don't know why but somehow removing this two lines of code:

        if(dot(rnd,nrm) < 0.0): #wrong hemisphere
            rnd = mult(rnd,-1.0);

fixed the lighting problem:

https://picload.org/view/dapopwaw/test.png.html

But now I am wondering why this results lacks anti-aliasing even though I use the halton sequence:

    const double fov = 160.0 * M_PI / 180.0;
    const double zdir = 1.0 / tan(fov);

    const double aspect = (double)h / (double)w;

    const double jitterX = (halton(2, samples)) - 0.5;
    const double jitterY = (halton(3, samples)) - 0.5;

    const double xH = x + jitterX;
    const double yH = y + jitterY;

    const double xdir = (xH / (double)w) * 2.0 - 1.0;
    const double ydir = ((yH/ (double)h) * 2.0 - 1.0) * aspect;

    const struct Point dir = norm((struct Point){.x = xdir, .y = ydir, .z = zdir});

    return (struct Ray){.origin = c.pos, .dir = dir};

double halton(int base, int ix) {
    double r = 0;
    double f = 1.0 / base;
    int i = ix;
    while(i > 0) {
        r += f * (i % base);
        f /= base;
        i = (int)floor((double)i / base);
    }

    return r;
}