Yesterday I ran into smallpt, a tiny raytracer in 100LOC in C++.
It turns out it works as a nice cross-lang CPU benchmark! So I ported it to crystal and ran some benchmarks comparing with other implementations and … well, the result is README.md
Tl;DR: crystal gives you readable code with practically the same performance as C++ (as expected) and multithreaded performance scales linearly (this is an ideal workload TBH)
I ran your code on my laptop, here are the times (secs).
System:
Linux 7.0.0, Tuxedo OS
Tuxedo Gemini Gen3, AMD Ryzen 9, 7945HX, 5.4GHz, 16C|32T
$ ./smallpt-raytracer xyz
$ SMALLPT_WORKERS=1 ./smallpt-raytracer xyz
pixels | 32 threads | 1 thread |
--------------|-----------------|-------------------|
500 / 512 | 19.66 / 19.69 | 312/65 / 319.66 |
--------------|-----------------|-------------------|
1000 / 1024 | 39.08 / 40.04 | 623.17 / 639.09 |
--------------|-----------------|-------------------|
2000 / 2048 | 78.41 / 80.22 | 1245.02 / 1269.47 |
--------------|-----------------|-------------------|
I got|ran the C++ code (2008) to compare with Crystal on my system.
https://www.kevinbeason.com/smallpt/smallpt.tar.gz
I have gcc 13.3.0 on my system.
smallpt.ccp compiled with one warning, but ran.
$ g++ -O3 -march=native -fopenmp smallpt.cpp -o smallpt-cpp
smallpt.cpp: In function ‘int main(int, char**)’:
smallpt.cpp:82:44: warning: narrowing conversion of ‘((y * y) * y)’ from ‘int’ to ‘short unsigned int’ [-Wnarrowing]
82 | for (unsigned short x=0, Xi[3]={0,0,y*y*y}; x<w; x++) // Loop cols
The code for smallpt4k.cpp had the same warning, but also an error and wouldn’t compile.
I didn’t try to fix the error in the code so it would compile, but maybe someone else can.
$ g++ -O3 -march=native -fopenmp smallpt4k.cpp -o smallpt4k-cpp
smallpt4k.cpp: In function ‘void _start()’:
smallpt4k.cpp:83:44: warning: narrowing conversion of ‘((y * y) * y)’ from ‘int’ to ‘short unsigned int’ [-Wnarrowing]
83 | for (unsigned short x=0, Xi[3]={0,0,y*y*y}; x<w; x++) // Loop cols
| ~~~^~
/usr/bin/ld: /tmp/cc96YXat.o: in function `_start':
smallpt4k.cpp:(.text+0x2870): multiple definition of `_start'; /usr/lib/gcc/x86_64-linux-gnu/13/../../../x86_64-linux-gnu/Scrt1.o:(.text+0x0): first defined here
/usr/bin/ld: /usr/lib/gcc/x86_64-linux-gnu/13/../../../x86_64-linux-gnu/Scrt1.o: in function `_start':
(.text+0x1b): undefined reference to `main'
collect2: error: ld returned 1 exit status
The C++ code is a little faster than Crystal, here 33s faster for 5000 spp.
$ time ./smallpt-cr 5000
Rendering took 203.051049641s
./smallpt-cr 5000 6357.76s user 6.08s system 3132% cpu 3:23.14 total
$ time ./smallpt-cpp 5000
Rendering (5000 spp) 100.00%./smallpt-cpp 5000 5320.13s user 5.11s system 3129% cpu 2:50.18 total
Here are the C++ source files.
smallpt.cpp
#include <math.h> // smallpt, a Path Tracer by Kevin Beason, 2008
#include <stdlib.h> // Make : g++ -O3 -fopenmp smallpt.cpp -o smallpt
#include <stdio.h> // Remove "-fopenmp" for g++ version < 4.2
struct Vec { // Usage: time ./smallpt 5000 && xv image.ppm
double x, y, z; // position, also color (r,g,b)
Vec(double x_=0, double y_=0, double z_=0){ x=x_; y=y_; z=z_; }
Vec operator+(const Vec &b) const { return Vec(x+b.x,y+b.y,z+b.z); }
Vec operator-(const Vec &b) const { return Vec(x-b.x,y-b.y,z-b.z); }
Vec operator*(double b) const { return Vec(x*b,y*b,z*b); }
Vec mult(const Vec &b) const { return Vec(x*b.x,y*b.y,z*b.z); }
Vec& norm(){ return *this = *this * (1/sqrt(x*x+y*y+z*z)); }
double dot(const Vec &b) const { return x*b.x+y*b.y+z*b.z; } // cross:
Vec operator%(Vec&b){return Vec(y*b.z-z*b.y,z*b.x-x*b.z,x*b.y-y*b.x);}
};
struct Ray { Vec o, d; Ray(Vec o_, Vec d_) : o(o_), d(d_) {} };
enum Refl_t { DIFF, SPEC, REFR }; // material types, used in radiance()
struct Sphere {
double rad; // radius
Vec p, e, c; // position, emission, color
Refl_t refl; // reflection type (DIFFuse, SPECular, REFRactive)
Sphere(double rad_, Vec p_, Vec e_, Vec c_, Refl_t refl_):
rad(rad_), p(p_), e(e_), c(c_), refl(refl_) {}
double intersect(const Ray &r) const { // returns distance, 0 if nohit
Vec op = p-r.o; // Solve t^2*d.d + 2*t*(o-p).d + (o-p).(o-p)-R^2 = 0
double t, eps=1e-4, b=op.dot(r.d), det=b*b-op.dot(op)+rad*rad;
if (det<0) return 0; else det=sqrt(det);
return (t=b-det)>eps ? t : ((t=b+det)>eps ? t : 0);
}
};
Sphere spheres[] = {//Scene: radius, position, emission, color, material
Sphere(1e5, Vec( 1e5+1,40.8,81.6), Vec(),Vec(.75,.25,.25),DIFF),//Left
Sphere(1e5, Vec(-1e5+99,40.8,81.6),Vec(),Vec(.25,.25,.75),DIFF),//Rght
Sphere(1e5, Vec(50,40.8, 1e5), Vec(),Vec(.75,.75,.75),DIFF),//Back
Sphere(1e5, Vec(50,40.8,-1e5+170), Vec(),Vec(), DIFF),//Frnt
Sphere(1e5, Vec(50, 1e5, 81.6), Vec(),Vec(.75,.75,.75),DIFF),//Botm
Sphere(1e5, Vec(50,-1e5+81.6,81.6),Vec(),Vec(.75,.75,.75),DIFF),//Top
Sphere(16.5,Vec(27,16.5,47), Vec(),Vec(1,1,1)*.999, SPEC),//Mirr
Sphere(16.5,Vec(73,16.5,78), Vec(),Vec(1,1,1)*.999, REFR),//Glas
Sphere(600, Vec(50,681.6-.27,81.6),Vec(12,12,12), Vec(), DIFF) //Lite
};
inline double clamp(double x){ return x<0 ? 0 : x>1 ? 1 : x; }
inline int toInt(double x){ return int(pow(clamp(x),1/2.2)*255+.5); }
inline bool intersect(const Ray &r, double &t, int &id){
double n=sizeof(spheres)/sizeof(Sphere), d, inf=t=1e20;
for(int i=int(n);i--;) if((d=spheres[i].intersect(r))&&d<t){t=d;id=i;}
return t<inf;
}
Vec radiance(const Ray &r, int depth, unsigned short *Xi){
double t; // distance to intersection
int id=0; // id of intersected object
if (!intersect(r, t, id)) return Vec(); // if miss, return black
const Sphere &obj = spheres[id]; // the hit object
Vec x=r.o+r.d*t, n=(x-obj.p).norm(), nl=n.dot(r.d)<0?n:n*-1, f=obj.c;
double p = f.x>f.y && f.x>f.z ? f.x : f.y>f.z ? f.y : f.z; // max refl
if (++depth>5) if (erand48(Xi)<p) f=f*(1/p); else return obj.e; //R.R.
if (obj.refl == DIFF){ // Ideal DIFFUSE reflection
double r1=2*M_PI*erand48(Xi), r2=erand48(Xi), r2s=sqrt(r2);
Vec w=nl, u=((fabs(w.x)>.1?Vec(0,1):Vec(1))%w).norm(), v=w%u;
Vec d = (u*cos(r1)*r2s + v*sin(r1)*r2s + w*sqrt(1-r2)).norm();
return obj.e + f.mult(radiance(Ray(x,d),depth,Xi));
} else if (obj.refl == SPEC) // Ideal SPECULAR reflection
return obj.e + f.mult(radiance(Ray(x,r.d-n*2*n.dot(r.d)),depth,Xi));
Ray reflRay(x, r.d-n*2*n.dot(r.d)); // Ideal dielectric REFRACTION
bool into = n.dot(nl)>0; // Ray from outside going in?
double nc=1, nt=1.5, nnt=into?nc/nt:nt/nc, ddn=r.d.dot(nl), cos2t;
if ((cos2t=1-nnt*nnt*(1-ddn*ddn))<0) // Total internal reflection
return obj.e + f.mult(radiance(reflRay,depth,Xi));
Vec tdir = (r.d*nnt - n*((into?1:-1)*(ddn*nnt+sqrt(cos2t)))).norm();
double a=nt-nc, b=nt+nc, R0=a*a/(b*b), c = 1-(into?-ddn:tdir.dot(n));
double Re=R0+(1-R0)*c*c*c*c*c,Tr=1-Re,P=.25+.5*Re,RP=Re/P,TP=Tr/(1-P);
return obj.e + f.mult(depth>2 ? (erand48(Xi)<P ? // Russian roulette
radiance(reflRay,depth,Xi)*RP:radiance(Ray(x,tdir),depth,Xi)*TP) :
radiance(reflRay,depth,Xi)*Re+radiance(Ray(x,tdir),depth,Xi)*Tr);
}
int main(int argc, char *argv[]){
int w=1024, h=768, samps = argc==2 ? atoi(argv[1])/4 : 1; // # samples
Ray cam(Vec(50,52,295.6), Vec(0,-0.042612,-1).norm()); // cam pos, dir
Vec cx=Vec(w*.5135/h), cy=(cx%cam.d).norm()*.5135, r, *c=new Vec[w*h];
#pragma omp parallel for schedule(dynamic, 1) private(r) // OpenMP
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;
}
}
FILE *f = fopen("image.ppm", "w"); // Write image to PPM file.
fprintf(f, "P3\n%d %d\n%d\n", w, h, 255);
for (int i=0; i<w*h; i++)
fprintf(f,"%d %d %d ", toInt(c[i].x), toInt(c[i].y), toInt(c[i].z));
}
smallpt4k.cpp
#include <math.h> // smallpt, a Path Tracer by Kevin Beason, 2008
#include <stdlib.h> // Make : g++ -O3 -fopenmp smallpt4k.cpp -o smallpt4k
#include <stdio.h> // Remove "-fopenmp" for g++ version < 4.2
struct Vec { // Usage: time ./smallpt4k && xv image.ppm
double x, y, z; // position, also color (r,g,b)
Vec(double x_=0, double y_=0, double z_=0){ x=x_; y=y_; z=z_; }
Vec operator+(const Vec &b) const { return Vec(x+b.x,y+b.y,z+b.z); }
Vec operator-(const Vec &b) const { return Vec(x-b.x,y-b.y,z-b.z); }
Vec operator*(double b) const { return Vec(x*b,y*b,z*b); }
Vec mult(const Vec &b) const { return Vec(x*b.x,y*b.y,z*b.z); }
Vec& norm(){ return *this = *this * (1/sqrt(x*x+y*y+z*z)); }
double dot(const Vec &b) const { return x*b.x+y*b.y+z*b.z; } // cross:
Vec operator%(Vec&b){return Vec(y*b.z-z*b.y,z*b.x-x*b.z,x*b.y-y*b.x);}
};
struct Ray { Vec o, d; Ray(Vec o_, Vec d_) : o(o_), d(d_) {} };
enum Refl_t { DIFF, SPEC, REFR }; // material types, used in radiance()
struct Sphere {
double rad; // radius
Vec p, e, c; // position, emission, color
Refl_t refl; // reflection type (DIFFuse, SPECular, REFRactive)
Sphere(){}
Sphere(double rad_, Vec p_, Vec e_, Vec c_, Refl_t refl_):
rad(rad_), p(p_), e(e_), c(c_), refl(refl_) {}
double intersect(const Ray &r) const { // returns distance, 0 if nohit
Vec op = p-r.o; // Solve t^2*d.d + 2*t*(o-p).d + (o-p).(o-p)-R^2 = 0
double t, eps=1e-4, b=op.dot(r.d), det=b*b-op.dot(op)+rad*rad;
if (det<0) return 0; else det=sqrt(det);
return (t=b-det)>eps ? t : ((t=b+det)>eps ? t : 0);
}
};
Sphere spheres[9]; //Scene: radius, position, emission, color, material
inline double clamp(double x){ return x<0 ? 0 : x>1 ? 1 : x; }
inline int toInt(double x){ return int(pow(clamp(x),1/2.2)*255+.5); }
inline bool intersect(const Ray &r, double &t, int &id){
double n=sizeof(spheres)/sizeof(Sphere), d, inf=t=1e20;
for(int i=int(n);i--;) if((d=spheres[i].intersect(r))&&d<t){t=d;id=i;}
return t<inf;
}
Vec radiance(const Ray &r, int depth, unsigned short *Xi){
double t; // distance to intersection
int id=0; // id of intersected object
if (!intersect(r, t, id)) return Vec(); // if miss, return black
const Sphere &obj = spheres[id]; // the hit object
Vec x=r.o+r.d*t, n=(x-obj.p).norm(), nl=n.dot(r.d)<0?n:n*-1, f=obj.c;
double p = f.x>f.y && f.x>f.z ? f.x : f.y>f.z ? f.y : f.z; // max refl
if (++depth>5) if (erand48(Xi)<p) f=f*(1/p); else return obj.e; //R.R.
if (obj.refl == DIFF){ // Ideal DIFFUSE reflection
double r1=2*M_PI*erand48(Xi), r2=erand48(Xi), r2s=sqrt(r2);
Vec w=nl, u=((fabs(w.x)>.1?Vec(0,1):Vec(1))%w).norm(), v=w%u;
Vec d = (u*cos(r1)*r2s + v*sin(r1)*r2s + w*sqrt(1-r2)).norm();
return obj.e + f.mult(radiance(Ray(x,d),depth,Xi));
} else if (obj.refl == SPEC) // Ideal SPECULAR reflection
return obj.e + f.mult(radiance(Ray(x,r.d-n*2*n.dot(r.d)),depth,Xi));
Ray reflRay(x, r.d-n*2*n.dot(r.d)); // Ideal dielectric REFRACTION
bool into = n.dot(nl)>0; // Ray from outside going in?
double nc=1, nt=1.5, nnt=into?nc/nt:nt/nc, ddn=r.d.dot(nl), cos2t;
if ((cos2t=1-nnt*nnt*(1-ddn*ddn))<0) // Total internal reflection
return obj.e + f.mult(radiance(reflRay,depth,Xi));
Vec tdir = (r.d*nnt - n*((into?1:-1)*(ddn*nnt+sqrt(cos2t)))).norm();
double a=nt-nc, b=nt+nc, R0=a*a/(b*b), c = 1-(into?-ddn:tdir.dot(n));
double Re=R0+(1-R0)*c*c*c*c*c,Tr=1-Re,P=.25+.5*Re,RP=Re/P,TP=Tr/(1-P);
return obj.e + f.mult(depth>2 ? (erand48(Xi)<P ? // Russian roulette
radiance(reflRay,depth,Xi)*RP:radiance(Ray(x,tdir),depth,Xi)*TP) :
radiance(reflRay,depth,Xi)*Re+radiance(Ray(x,tdir),depth,Xi)*Tr);
}
extern "C" void _start() {
spheres[0]=Sphere(1e5, Vec( 1e5+1,40.8,81.6), Vec(),Vec(.75,.25,.25),DIFF);//Left
spheres[1]=Sphere(1e5, Vec(-1e5+99,40.8,81.6),Vec(),Vec(.25,.25,.75),DIFF);//Rght
spheres[2]=Sphere(1e5, Vec(50,40.8, 1e5), Vec(),Vec(.75,.75,.75),DIFF);//Back
spheres[3]=Sphere(1e5, Vec(50,40.8,-1e5+170), Vec(),Vec(), DIFF);//Frnt
spheres[4]=Sphere(1e5, Vec(50, 1e5, 81.6), Vec(),Vec(.75,.75,.75),DIFF);//Botm
spheres[5]=Sphere(1e5, Vec(50,-1e5+81.6,81.6),Vec(),Vec(.75,.75,.75),DIFF);//Top
spheres[6]=Sphere(16.5,Vec(27,16.5,47), Vec(),Vec(1,1,1)*.999, SPEC);//Mirr
spheres[7]=Sphere(16.5,Vec(73,16.5,78), Vec(),Vec(1,1,1)*.999, REFR);//Glas
spheres[8]=Sphere(600, Vec(50,681.6-.27,81.6),Vec(12,12,12), Vec(), DIFF);//Lite
int w=1024, h=768, samps = 5000/4; // # samples
Ray cam(Vec(50,52,295.6), Vec(0,-0.042612,-1).norm()); // cam pos, dir
Vec cx=Vec(w*.5135/h), cy=(cx%cam.d).norm()*.5135, r,
*c=(Vec*)malloc(sizeof(Vec)*w*h);
#pragma omp parallel for schedule(dynamic, 1) private(r) // OpenMP
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;
}
}
FILE *f = fopen("image.ppm", "w"); // Write image to PPM file.
fprintf(f, "P3\n%d %d\n%d\n", w, h, 255);
for (int i=0; i<w*h; i++)
fprintf(f,"%d %d %d ", toInt(c[i].x), toInt(c[i].y), toInt(c[i].z));
fclose(f);
exit(0);
}
Very nice!
We can optimize MT further: the fibers are fully CPU-bound, and don’t need to communicate aside from the eventual “I’m done” (and a few event-loop writes). Instead of relying on auto-scaling to eventually scale the fibers across to N threads (aka slow-parallelism) we can start N isolated fibers that immediately start rendering.
It ain’t more complex:
spawn(name, &) with Fiber::ExecutionContext::Isolated.new(name, &).On my Intel 14700k, Ubuntu 24.04 and rusage:
Parallel (1.26s WALL, 1733% CPU):
$ crystal build src/smallpt.cr -o smallpt-parallel --mcpu=native --release -Dgc_none
$ rusage time smallpt-parallel 32
Rendering took 1.26085764s
23.16user 0.03system 0:01.33elapsed 1733%CPU (0avgtext+0avgdata 33236maxresident)k
0inputs+15704outputs (0major+8435minor)pagefaults 0swaps
took 1,340,655µs wall time
ballooned to 33,236kb in size
needed 23,199,017us cpu (0% kernel)
caused 8,521 page faults (100% memcpy)
1,632 context switches (11% consensual)
performed 0 read and 15,704 write i/o operations
Isolated (0.91s, 2510% CPU):
$ crystal build src/smallpt.cr -o smallpt-parallel --mcpu=native --release -Dgc_none
$ rusage time smallpt-isolated 32
Rendering took 0.914614525s
24.76user 0.03system 0:00.98elapsed 2510%CPU (0avgtext+0avgdata 30072maxresident)k
0inputs+15704outputs (0major+8332minor)pagefaults 0swaps
took 989,117µs wall time
ballooned to 30,072kb in size
needed 24,801,556us cpu (0% kernel)
caused 8,420 page faults (100% memcpy)
2,277 context switches (8% consensual)
performed 0 read and 15,704 write i/o operations
GCC 13, OpenMP 0.7.3 (1.10s WALL; 2522% CPU):
$ c++ -O3 -march=native -fopenmp smallpt.cpp -o smallpt-cpp
$ rusage time smallpt-isolated 32
27.76user 0.03system 0:01.10elapsed 2522%CPU (0avgtext+0avgdata 35476maxresident)k
0inputs+15680outputs (0major+8809minor)pagefaults 0swaps
took 1,104,450µs wall time
ballooned to 35,476kb in size
needed 27,803,973us cpu (0% kernel)
caused 8,896 page faults (100% memcpy)
2,768 context switches (4% consensual)
performed 0 read and 15,680 write i/o operations
NOTE: I’m using a low sample (32) on purpose. With larger samples, for example 512, there is no noticeable difference between Isolated and Parallel anymore (maybe 1% at most, might be a fluke). Though GCC and OpenMP stay significantly slower (~10%).
That sounds like there should perhaps be an option to start fully scaled out on the things..
Can you show the different code.
Hm, parallel doesn’t start all the threads because they’d spin and stop until fibers are spawned, but we should be able to manually scale up, yes 
Scaling down is probably nonsensical: it already happens when there’s nothing to do, and otherwise the monitor thread would re-scale up anyway.
@jzakiya I already gave the code. Read my previous post: it’s one line to remove and one line to edit.
I simplified|shortened the code, and used both Parallel and Isolated contexts, but I
found no appreciable performance difference beteen them.
workers = ENV["SMALLPT_WORKERS"]?.try(&.to_i) || System.cpu_count
Fiber::ExecutionContext.default.resize(workers)
next_row = Atomic(Int32).new(0)
done_rows = Atomic(Int32).new(0)
wg = WaitGroup.new(workers)
workers.times do
spawn do
while (y = next_row.add(1)) < height
....
end
wg.done
end end
wg.wait
---------------------------------------------------------
workers = ENV["SMALLPT_WORKERS"]?.try(&.to_i) || System.cpu_count
next_row = Atomic(Int32).new(0)
done_rows = Atomic(Int32).new(0)
wg = WaitGroup.new(workers)
workers.times do
Fiber::ExecutionContext::Isolated.new("workers") do
while (y = next_row.add(1)) < height
....
end
wg.done
end end
wg.wait
What did make a HUGE difference (~20% faster) was turning off gc with -Dgc_none.
crystal build --release --mcpu=native -Dgc_none smallpt-cr.cr (-o xxxxx)
For input 5000 spp, with gc is ~200s, with Dgc_none it went down to ~162-5.
This means Crystal is now faster than C++, ~2m40s vs ~2m50s, for 5000 spp.
I didn’t realize you could turn off garbage collection like this, so this is new knowledge.
As Crystal is evolving rapidly in this area, it would be very useful to update the
performance tutorials on what|how users can use compiler directives, etc, to optimize for speed.
In fact, this code provides a nice non-trivial example of various ways to increase speed.
Here’s my simplified|shortened code, without the line comments. I did $ diff image1.ppm image2.ppm to verify the outputs.
# Optimal compile as:
# crystal build --release --mcpu=native -Dgc_none smallpt-cr.cr (-o xxxxx)
require "wait_group"
record(Vec, x : Float64 = 0.0, y : Float64 = 0.0, z : Float64 = 0.0) do
def +(other : Vec) Vec.new(x + other.x, y + other.y, z + other.z) end
def -(other : Vec) Vec.new(x - other.x, y - other.y, z - other.z) end
def *(other : Float64) Vec.new(x * other, y * other, z * other) end
def /(other : Float64) a = 1.0/other; Vec.new(x * a, y * a , z * a) end
def mult(other : Vec) Vec.new(x * other.x, y * other.y, z * other.z) end
def %(other : Vec) Vec.new(y * other.z - z * other.y, z * other.x - x * other.z, x * other.y - y * other.x) end
def dot(other : Vec) x * other.x + y * other.y + z * other.z end
def norm; self / Math.sqrt(x * x + y * y + z * z) end
end
record(Ray, o : Vec, d : Vec)
enum ReflT
Diffuse
Specular
Refractive
end
struct Sphere
getter radius : Float64, position : Vec, emission : Vec, color : Vec, reflection : ReflT
def initialize(@radius, @position, @emission, @color, @reflection) end
def intersect(ray : Ray)
op = position - ray.o
epsilon, b = 1e-4, op.dot(ray.d)
det = b * b - op.dot(op) + radius * radius
return 0.0 if det < 0
det = Math.sqrt(det)
t = b - det
return t if t > epsilon
t = b + det
t > epsilon ? t : 0.0
end end
SPHERES = [
Sphere.new(1e5, Vec.new(1e5 + 1, 40.8, 81.6), Vec.new, Vec.new(0.75, 0.25, 0.25), ReflT::Diffuse), # Left
Sphere.new(1e5, Vec.new(-1e5 + 99, 40.8, 81.6), Vec.new, Vec.new(0.25, 0.25, 0.75), ReflT::Diffuse), # Right
Sphere.new(1e5, Vec.new(50, 40.8, 1e5), Vec.new, Vec.new(0.75, 0.75, 0.75), ReflT::Diffuse), # Back
Sphere.new(1e5, Vec.new(50, 40.8, -1e5 + 170), Vec.new, Vec.new, ReflT::Diffuse), # Front
Sphere.new(1e5, Vec.new(50, 1e5, 81.6), Vec.new, Vec.new(0.75, 0.75, 0.75), ReflT::Diffuse), # Bottom
Sphere.new(1e5, Vec.new(50, -1e5 + 81.6, 81.6), Vec.new, Vec.new(0.75, 0.75, 0.75), ReflT::Diffuse), # Top
Sphere.new(16.5, Vec.new(27, 16.5, 47), Vec.new, Vec.new(1, 1, 1) * 0.999, ReflT::Specular), # Mirror
Sphere.new(16.5, Vec.new(73, 16.5, 78), Vec.new, Vec.new(1, 1, 1) * 0.999, ReflT::Refractive), # Glass
Sphere.new(600, Vec.new(50, 681.6 - 0.27, 81.6), Vec.new(12, 12, 12), Vec.new, ReflT::Diffuse), # Light
]
def clamp(value : Float64) value < 0 ? 0.0 : value > 1 ? 1.0 : value end
def to_int(value : Float64) (clamp(value) ** (1 / 2.2) * 255 + 0.5).to_i end
def intersect(ray : Ray)
t, id = 1e20, 0
SPHERES.size.downto(1) do |index|
distance = SPHERES[index - 1].intersect(ray)
next unless distance != 0.0 && distance < t
t, id = distance, index - 1
end
{t < 1e20, t, id}
end
TWO_POW_32_INV = 1.0 / 4294967296.0
def next_f(rng : Random::PCG32) : Float64; rng.next_u * TWO_POW_32_INV end
def radiance(ray : Ray, depth : Int32, rng : Random::PCG32) : Vec
hit, t, sphere_index = intersect(ray)
return Vec.new unless hit
object = SPHERES[sphere_index]
x = ray.o + ray.d * t # hit point
n = (x - object.position).norm # geometric normal
nl = n.dot(ray.d) < 0 ? n : n * -1.0 # normal facing the incoming ray
f = object.color
p = f.x > f.y && f.x > f.z ? f.x : f.y > f.z ? f.y : f.z
depth &+= 1
next_f(rng) < p ? (f = f / p) : return object.emission if depth > 5
case object.reflection
in .diffuse?
r1, r2 = 2 * Math::PI * next_f(rng), next_f(rng)
r2s = Math.sqrt(r2)
w = nl
u = ((w.x.abs > 0.1 ? Vec.new(0, 1, 0) : Vec.new(1, 0, 0)) % w).norm
v = w % u
d = (u * (Math.cos(r1) * r2s) + v * (Math.sin(r1) * r2s) + w * Math.sqrt(1 - r2)).norm
object.emission + f.mult(radiance(Ray.new(x, d), depth, rng))
in .specular?
object.emission + f.mult(radiance(Ray.new(x, ray.d - n * (2 * n.dot(ray.d))), depth, rng))
in .refractive?
radiance_refractive(object, ray, x, n, nl, f, depth, rng)
end end
def radiance_refractive(object : Sphere, ray : Ray, x : Vec, n : Vec, nl : Vec,
f : Vec, depth : Int32, rng : Random::PCG32) : Vec
refl_ray = Ray.new(x, ray.d - n * (2 * n.dot(ray.d)))
into = n.dot(nl) > 0 # entering or exiting the glass?
nc, nt = 1.0, 1.5 # index of refraction of air|glass
nnt = into ? nc / nt : nt / nc
ddn = ray.d.dot(nl)
cos2t = 1 - nnt * nnt * (1 - ddn * ddn)
return object.emission + f.mult(radiance(refl_ray, depth, rng)) if cos2t < 0
tdir = (ray.d * nnt - n * ((into ? 1.0 : -1.0) * (ddn * nnt + Math.sqrt(cos2t)))).norm
a, b = nt - nc, nt + nc
r0 = a * a / (b * b) # Fresnel reflectance at normal incidence
c = 1 - (into ? -ddn : tdir.dot(n))
re = r0 + (1 - r0) * c * c * c * c * c # Schlick's approximation
tr, prob = 1 - re, 0.25 + 0.5 * re
rp, tp = re / prob, tr / (1 - prob)
object.emission +
f.mult(
if depth > 2
next_f(rng) < prob ? radiance(refl_ray, depth, rng) * rp : radiance(Ray.new(x, tdir), depth, rng) * tp
else
radiance(refl_ray, depth, rng) * re + radiance(Ray.new(x, tdir), depth, rng) * tr
end
)
end
width, height = 1024, 768
samples : Int32 = ARGV.size == 1 ? (ARGV[0].to_i // 4) : 1
camera = Ray.new(Vec.new(50, 52, 295.6), Vec.new(0, -0.042612, -1).norm)
cx = Vec.new(width * 0.5135 / height)
cy = (cx % camera.d).norm * 0.5135
canvas = Slice(Vec).new(width * height, Vec.new)
def render_row(y : Int32, width : Int32, height : Int32, samples : Int32,
camera : Ray, cx : Vec, cy : Vec, canvas : Slice(Vec))
rng = Random::PCG32.new(UInt64.new(y * y * y))
width.times do |x|
2.times do |subpixel_y|
index = (height - y - 1) * width + x
2.times do |subpixel_x|
accumulated = Vec.new
samples.times do
r1, r2 = 2 * next_f(rng), 2 * next_f(rng)
dx = r1 < 1 ? Math.sqrt(r1) - 1 : 1 - Math.sqrt(2 - r1)
dy = r2 < 1 ? Math.sqrt(r2) - 1 : 1 - Math.sqrt(2 - r2)
direction = cx * (((subpixel_x + 0.5 + dx) / 2 + x) / width - 0.5) +
cy * (((subpixel_y + 0.5 + dy) / 2 + y) / height - 0.5) + camera.d
accumulated += radiance(Ray.new(camera.o + direction * 140, direction.norm), 0, rng) / samples.to_f
end
canvas[index] += Vec.new(clamp(accumulated.x), clamp(accumulated.y), clamp(accumulated.z)) / 4.0
end end end end
start_time = Time.instant
workers = ENV["SMALLPT_WORKERS"]?.try(&.to_i) || System.cpu_count
next_row = Atomic(Int32).new(0)
done_rows = Atomic(Int32).new(0)
wg = WaitGroup.new(workers)
workers.times do
Fiber::ExecutionContext::Isolated.new("workers") do
while (y = next_row.add(1)) < height
render_row(y, width, height, samples, camera, cx, cy, canvas)
done = done_rows.add(1) &+ 1
STDERR.printf("\rRendering (%d spp) %5.2f%%", samples * 4, 100.0 * done / height) if done.divisible_by?(16)
end
wg.done
end end
wg.wait
elapsed = Time.instant - start_time
STDERR.puts "\rRendering took #{elapsed.total_milliseconds / 1000.0}s "
def build_string(canvas : Slice(Vec), width : Int32, height : Int32)
String.build do |io|
io << "P3\n#{width} #{height}\n255\n"
canvas.each { |pixel| io << to_int(pixel.x) << ' ' << to_int(pixel.y) << ' ' << to_int(pixel.z) << ' ' }
end end
File.write("image.ppm", build_string(canvas, width, height))
I said “impact is only for small samples like 32”. At 512 it’s already invisible, so 5000 is invisible.
Here, parallel processing can also be done using thread pools.
Below are 3 examples for doing it, listed from most to least amount of code.
These implementations are also more similar to the C++ implementation.
Compiled as: crystal build --release --mcpu=native -Dgc_none smallpt-xyz.cr
Also changed: done = done_rows.add(1) + 1 to just done = done_rows.add(1).
The last version (least code) was definitively fastest for all inputs (tested up to 10000 spp).
So for this task, I like the last implementation best, as it’s shortest code is intuitively easiest
to understand, it’s the fastest, and only uses a smidgen more (~200MB) memory than the others.
Used $ diff image-xyw.ppm image-abc.ppm to verify all images were the same for given inputs.
next_row = Atomic(Int32).new(0)
done_rows = Atomic(Int32).new(0)
wg = WaitGroup.new(height)
Fiber::ExecutionContext.default.resize(workers)
height.times do
spawn do
y = next_row.add(1)
render_row(y, width, height, samples, camera, cx, cy, canvas)
done = done_rows.add(1)
STDERR.printf("\rRendering (%d spp) %5.2f%%", samples * 4, 100.0 * done / height) if done.divisible_by?(16)
wg.done
end end
wg.wait
------------------------------------------------------------------------------
done_rows = Atomic(Int32).new(0)
wg = WaitGroup.new(height)
Fiber::ExecutionContext.default.resize(workers)
height.times do |y|
spawn do
render_row(y, width, height, samples, camera, cx, cy, canvas)
done = done_rows.add(1)
STDERR.printf("\rRendering (%d spp) %5.2f%%", samples * 4, 100.0 * done / height) if done.divisible_by?(16)
wg.done
end end
wg.wait
-------------------------------------------------------------------------------
done_rows = Atomic(Int32).new(0)
wg = WaitGroup.new(height)
height.times do |y|
Fiber::ExecutionContext::Isolated.new("workers") do
render_row(y, width, height, samples, camera, cx, cy, canvas)
done = done_rows.add(1)
STDERR.printf("\rRendering (%d spp) %5.2f%%", samples * 4, 100.0 * done / height) if done.divisible_by?(16)
wg.done
end end
wg.wait
Hi, You are not added the github link into the blog, right?
Are you asking if my version is on the github link? (No it isn’t)
Hey @ralsina, I |ran your go version on my system and their diffd images give this error:
\ No newline at end of file
All the C++|Crystal images for a given input are the same between them.
Also the compiler flags docs for Crystal didn’t list -Dgc_none. This not only significantly increases its speed, but on my system the binary went from 8685302 to 678192 bytes (no gc code).
The easiest way do a parallel Rust version is to use the rayon crate, which has become standard practice for Rust, like using OpenMP is for C++. Doing something like height.par_iter() will make its code look|perform similar as C++|Crystal.
I think -Dgc_none is kinda cheating?
Why? Here’s what you say on your github page.
A Crystal port of smallpt, Kevin Beason's famous ~100-line unbiased path tracer, written to explore the performance gap between optimized C++ and optimized Crystal.
Until @ysbaddaden showed this I didn’t know anything about it. So this is obviously a valid option, for just this purpose. Exposing all the possibilities to optimize Crystal code can only be in the best interest of its users, the community, and the future of the language. Why use gc if your code doesn’t need it? Some might say using OpenMP is cheating for C++ because is a total 3rd party application to parallelize it, while EC is part of the Crystal standard environment.
I say, if your intent is to showcase how to optimize the use of Crystal then do it all the way. 
Compiling without garbage collection puts a hard upper bound on how much work the process can do, making it unfit for the use cases ray tracers are typically used in. In order for a microbenchmark to be useful, it must represent how it’s used in a real-world program. Cutting out GC here optimizes for the benchmark itself.
Compiling without garbage collection puts a hard upper bound on how much work the process can do, making it unfit for the use cases ray tracers are typically used in. In order for a microbenchmark to be useful, it must represent how it’s used in a real-world program.
Huh???
Cutting out GC here optimizes for the benchmark itself.
Exactly!!. Which is the stated purpose of this exercise.
Also, the record shows the gitbub page also states results using different compiler flags.
Native-CPU code generation was attempted for all four: it only helped Crystal (--mcpu=native, ~18% faster) and C++ (-march=native). For Rust, RUSTFLAGS="-C target-cpu=native" changed nothing measurable; for Go, neither GOAMD64=v3 nor GOGC=off helped (GOAMD64=v4 requires AVX-512 and won't run on the test machine at all).
This shows @ralsina tried using GOGC=false, but he said it made no difference for go.
It’s logically inconsistent to try turning off garbage collection for go but not for Crystal.
Because Rust|C++ don’t use gc, Crystal is at a disadvantage using it. So using -Dgc_none here creates a fairer apples-to-apples comparison with them.
So @ralsina states he used various compiler flags for all the languages to optimize speed.
And @ysbaddaden showing compiling with -Dgc_none is consistent with the stated goal to create a (most) optimized Crystal version. Thus, you should be directing your ire at him, not me. I’m merely advocating its use here, which I thank him for revealing. Of course, people are free to use|not use it as they feel. But you can’t consciously not use something you didn’t know you could.
Thus, if turning off gc for go was acceptable, then so too for Crystal; and its benefits (speed|binary reduction) should be listed as well.
Has another drama started? I do not see any ire at all in jgaskins’s comment.
I think this comes from a difference in what we want from programming. From his usual comments, jgaskins seems interested in software that runs on servers for a long time and handles a large throughput efficiently. A benchmark is only a tool for that purpose.
People like jzakiya and kojix2, on the other hand, enjoy programming more as a hobby. jzakiya values speed. I prefer programs that process large amounts of data with little memory. The former world is much larger. The latter has little economic sustainability, so we do not see it very often.
Another approach is to put 256 GB of memory in a desktop PC and treat the computation as a short race: it is enough if the program finishes before the memory fills up. This approach is not without value. But we should understand why it does not receive much support from industry.
(Translated by ChatGPT, as usual.)
That is indeed my usual perspective. And in this case, ray tracers typically execute on every render loop inside long-running processes, in video games. GC would be needed for this.
FWIW the raytracer itself should be essentially allocation-free. Which should of course be expected from a performance-concious algorithm that’s intended to run in a tight loop.
The only allocations controled by the algorithm are for the scene, the canvas and the PPM format.
All these can be explicitly managed and don’t need garbage collection.
The primary work for the GC is one-off runtime setup.
Thats external overhead for the raytracing algorithm.
And nothing there gets ever collected, yet the GC doesn’t know that.
Is it possible then to use GC.disable|GC.enable in the code to fine-tune this behavior?