Kevin Beason's smallpt estimates the pixel radiance by accumulating $2\times2$ subpixel radiance estimates using the following expression:

c[i] = c[i] + Vec(clamp(r.x),clamp(r.y),clamp(r.z))*.25;

Here, c[i] represents the radiance estimate of the pixel with flattened index i, r represents the radiance estimate of one of its four subpixels. As one notices a clamp operation is used to clamp the radiance estimate of each subpixel separately to the $[0,1]^3$ range before accumulation. But this seems wrong from an unbiased path tracing perspective, since this effectively introduces bias to the radiance estimate of each pixel. So I wonder if I am missing something (e.g., clever trick at the cost of bias), given that Kevin Beason's smallpt is already around and known for nearly ten years?

Apparently, Kevin Beason also added a presentation by David Cline explaining every line. The aforementioned expression is explained as follows:

Add the gamma-corrected subpixel color estimate to the Pixel color c[i].

This seems, however, completely unrelated since no gamma correction is involved yet at this stage in the code. And since exponential functions are non-linear transformations.

Reference code:

for (int y=0; y<h; y++){                       // Loop over image rows 
    fprintf(stderr,"\rRendering (%d spp) %5.2f%%",samps*4,100.*y/(h-1)); 
    for (unsigned short x=0, Xi[3]={0,0,y*y*y}; x<w; x++)   // Loop cols 
        for (int sy=0, i=(h-y-1)*w+x; sy<2; sy++)     // 2x2 subpixel rows 
            for (int sx=0; sx<2; sx++, r=Vec()){        // 2x2 subpixel cols 
                for (int s=0; s<samps; s++){ 
                    double r1=2*erand48(Xi), dx=r1<1 ? sqrt(r1)-1: 1-sqrt(2-r1); 
                    double r2=2*erand48(Xi), dy=r2<1 ? sqrt(r2)-1: 1-sqrt(2-r2); 
                    Vec d = cx*( ( (sx+.5 + dx)/2 + x)/w - .5) + 
                            cy*( ( (sy+.5 + dy)/2 + y)/h - .5) + cam.d; 
                    r = r + radiance(Ray(cam.o+d*140,d.norm()),0,Xi)*(1./samps); 
                } // Camera rays are pushed ^^^^^ forward to start in interior 
                c[i] = c[i] + Vec(clamp(r.x),clamp(r.y),clamp(r.z))*.25; 
            } 
}