Initial SIMD demos life and earth for PNaCl.

native_client_sdk/src/examples/demo/earth_simd/earth.cc

 // Copyright (c) 2013 The Chromium Authors. All rights reserved.
 // Use of this source code is governed by a BSD-style license that can be
 // found in the LICENSE file.
 
 #include <assert.h>
 #include <math.h>
 #include <ppapi/c/ppb_input_event.h>
 #include <ppapi/cpp/input_event.h>
 #include <ppapi/cpp/var.h>
 #include <ppapi/cpp/var_array.h>
@@ skipping 12 matching lines @@
 #include "ppapi_simple/ps.h"
 #include "ppapi_simple/ps_context_2d.h"
 #include "ppapi_simple/ps_event.h"
 #include "ppapi_simple/ps_interface.h"
 #include "ppapi_simple/ps_main.h"
 #include "sdk_util/macros.h"
 #include "sdk_util/thread_pool.h"
 
 using namespace sdk_util;  // For sdk_util::ThreadPool
 
+#define INLINE inline __attribute__((always_inline))
+
+// 128 bit SIMD vector types
+typedef uint8_t u8x16_t __attribute__ ((vector_size (16)));
+typedef int32_t i32x4_t __attribute__ ((vector_size (16)));
+typedef uint32_t u32x4_t __attribute__ ((vector_size (16)));
+typedef float f32x4_t __attribute__ ((vector_size (16)));
+
 // Global properties used to setup Earth demo.
 namespace {
 const float kPI = M_PI;
 const float kTwoPI = kPI * 2.0f;
 const float kOneOverPI = 1.0f / kPI;
 const float kOneOver2PI = 1.0f / kTwoPI;
-const float kOneOver255 = 1.0f / 255.0f;
 const int kArcCosineTableSize = 4096;
 const int kFramesToBenchmark = 100;
 const float kZoomMin = 1.0f;
 const float kZoomMax = 50.0f;
 const float kWheelSpeed = 2.0f;
 const float kLightMin = 0.0f;
 const float kLightMax = 2.0f;
 
 // Timer helper for benchmarking.  Returns seconds elapsed since program start,
 // as a double.
 timeval start_tv;
 int start_tv_retv = gettimeofday(&start_tv, NULL);
 
 inline double getseconds() {
   const double usec_to_sec = 0.000001;
   timeval tv;
   if ((0 == start_tv_retv) && (0 == gettimeofday(&tv, NULL)))
     return (tv.tv_sec - start_tv.tv_sec) + tv.tv_usec * usec_to_sec;
   return 0.0;
 }
 
-// RGBA helper functions, used for extracting color from RGBA source image.
-inline float ExtractR(uint32_t c) {
-  return static_cast<float>(c & 0xFF) * kOneOver255;
+// SIMD Vector helper functions.
+INLINE f32x4_t min(f32x4_t a, f32x4_t b) {
+  i32x4_t m = a < b;
+  return (f32x4_t)(((i32x4_t)a & m) | ((i32x4_t)b & ~m));
 }
 
-inline float ExtractG(uint32_t c) {
-  return static_cast<float>((c & 0xFF00) >> 8) * kOneOver255;
+INLINE f32x4_t max(f32x4_t a, f32x4_t b) {
+  i32x4_t m = a > b;
+  return (f32x4_t)(((i32x4_t)a & m) | ((i32x4_t)b & ~m));
 }
 
-inline float ExtractB(uint32_t c) {
-  return static_cast<float>((c & 0xFF0000) >> 16) * kOneOver255;
+INLINE float dot3(f32x4_t a, f32x4_t b) {
+  f32x4_t c = a * b;
+  return c[0] + c[1] + c[2];
+}
+
+INLINE f32x4_t broadcast(float x) {
+  f32x4_t r = {x, x, x, x};
+  return r;
+}
+
+// SIMD RGBA helper functions, used for extracting color from RGBA source image.
+INLINE f32x4_t ExtractRGBA(uint32_t c) {
+  const f32x4_t kOneOver255 = broadcast(1.0f / 255.0f);
+  const i32x4_t kZero = {0, 0, 0, 0};
+  i32x4_t v = {c, c, c, c};
+  // zero extend packed color into 32x4 integer vector
+  v = (i32x4_t)__builtin_shufflevector((u8x16_t)v, (u8x16_t)kZero,
+      0, 16, 16, 16, 1, 16, 16, 16, 2, 16, 16, 16, 3, 16, 16, 16);
+  // convert color values to float, range 0..1
+  f32x4_t f = __builtin_convertvector(v, f32x4_t) * kOneOver255;
+  return f;
+}
+
+// SIMD BGRA helper function, for constructing a pixel for a BGRA buffer.
+INLINE uint32_t PackBGRA(f32x4_t f) {
+  const f32x4_t kZero = broadcast(0.0f);
+  const f32x4_t kHalf = broadcast(0.5f);
+  const f32x4_t k255 = broadcast(255.0f);
+  f = max(f, kZero);
+  f = f * k255 + kHalf;
+  f = min(f, k255);
+  i32x4_t i = __builtin_convertvector(f, i32x4_t);
+  u32x4_t p = (u32x4_t)__builtin_shufflevector((u8x16_t)i, (u8x16_t)i,
+      8, 4, 0, 12, 8, 4, 0, 12, 8, 4, 0, 12, 8, 4, 0, 12);
+  return p[0];
 }
 
 // BGRA helper function, for constructing a pixel for a BGRA buffer.
-inline uint32_t MakeBGRA(uint32_t b, uint32_t g, uint32_t r, uint32_t a) {
+INLINE uint32_t MakeBGRA(uint32_t b, uint32_t g, uint32_t r, uint32_t a) {
   return (((a) << 24) | ((r) << 16) | ((g) << 8) | (b));
 }
 
 // simple container for earth texture
 struct Texture {
   int width, height;
   uint32_t* pixels;
   Texture(int w, int h) : width(w), height(h) {
     pixels = new uint32_t[w * h];
     memset(pixels, 0, sizeof(uint32_t) * w * h);
@@ skipping 20 matching lines @@
   // build a slightly larger table to allow for numeric imprecision
   for (int i = 0; i < (kArcCosineTableSize + 2); ++i) {
     float f = static_cast<float>(i) / kArcCosineTableSize;
     f = f * 2.0f - 1.0f;
     table[i] = acos(f);
   }
 }
 
 // looks up acos(f) using a table and lerping between entries
 // (it is expected that input f is between -1 and 1)
-float ArcCosine::TableLerp(float f) {
+INLINE float ArcCosine::TableLerp(float f) {
   float x = (f + 1.0f) * 0.5f;
   x = x * kArcCosineTableSize;
   int ix = static_cast<int>(x);
   float fx = static_cast<float>(ix);
   float dx = x - fx;
   float af = table[ix];
   float af2 = table[ix + 1];
   return af + (af2 - af) * dx;
 }
 
 // Helper functions for quick but approximate sqrt.
 union Convert {
   float f;
   int i;
   Convert(int x) { i = x; }
   Convert(float x) { f = x; }
   int AsInt() { return i; }
   float AsFloat() { return f; }
 };
 
-inline const int AsInteger(const float f) {
+INLINE const int AsInteger(const float f) {
   Convert u(f);
   return u.AsInt();
 }
 
-inline const float AsFloat(const int i) {
+INLINE const float AsFloat(const int i) {
   Convert u(i);
   return u.AsFloat();
 }
 
 const long int kOneAsInteger = AsInteger(1.0f);
 
-inline float inline_quick_sqrt(float x) {
+INLINE float inline_quick_sqrt(float x) {
   int i;
   i = (AsInteger(x) >> 1) + (kOneAsInteger >> 1);
   return AsFloat(i);
 }
 
-inline float inline_sqrt(float x) {
+INLINE float inline_sqrt(float x) {
   float y;
   y = inline_quick_sqrt(x);
   y = (y * y + x) / (2.0f * y);
   y = (y * y + x) / (2.0f * y);
   return y;
 }
 
-// takes a -0..1+ color, clamps it to 0..1 and maps it to 0..255 integer
-inline uint32_t Clamp255(float x) {
-  if (x < 0.0f) {
-    x = 0.0f;
-  } else if (x > 1.0f) {
-    x = 1.0f;
-  }
-  return static_cast<uint32_t>(x * 255.0f);
-}
 }  // namespace
 
 
 // The main object that runs the Earth demo.
 class Planet {
  public:
   Planet();
   virtual ~Planet();
   // Runs a tick of the simulations, update 2D output.
   void Update();
@@ skipping 187 matching lines @@
   *w = ps_context_->width;
   *y = r;
   *h = 1;
 }
 
 
 inline uint32_t* Planet::wGetAddr(int x, int y) {
   return ps_context_->data + x + y * ps_context_->stride / sizeof(uint32_t);
 }
 
-// This is the meat of the ray tracer.  Given a pixel span (x0, x1) on
+// This is the inner loop of the ray tracer.  Given a pixel span (x0, x1) on
 // scanline y, shoot rays into the scene and render what they hit.  Use
-// scanline coherence to do a few optimizations
+// scanline coherence to do a few optimizations.
+// This version uses portable SIMD 4 element single precision floating point
+// vectors to perform many of the calculations, and builds only on PNaCl.
 void Planet::wRenderPixelSpan(int x0, int x1, int y) {
   if (!base_tex_ || !night_tex_)
     return;
-  const int kColorBlack = MakeBGRA(0, 0, 0, 0xFF);
+  const uint32_t kColorBlack = MakeBGRA(0, 0, 0, 0xFF);
+  const uint32_t kSolidAlpha = MakeBGRA(0, 0, 0, 0xFF);
+  const f32x4_t kOne = {1.0f, 1.0f, 1.0f, 1.0f};
+  const f32x4_t diffuse = {diffuse_r_, diffuse_g_, diffuse_b_, 0.0f};
+  const f32x4_t ambient = {ambient_r_, ambient_g_, ambient_b_, 0.0f};
+  const f32x4_t light_pos = {light_x_, light_y_, light_z_, 1.0f};
+  const f32x4_t planet_pos = {planet_x_, planet_y_, planet_z_, 1.0f};
+  const f32x4_t planet_one_over_radius = broadcast(planet_one_over_radius_);
+  const f32x4_t planet_equator = {
+      planet_equator_x_, planet_equator_y_, planet_equator_z_, 0.0f};
+  const f32x4_t planet_pole = {
+      planet_pole_x_, planet_pole_y_, planet_pole_z_, 1.0f};
+  const f32x4_t planet_pole_x_equator = {
+      planet_pole_x_equator_x_, planet_pole_x_equator_y_,
+      planet_pole_x_equator_z_, 0.0f};
+
   float width = ps_context_->width;
   float height = ps_context_->height;
   float min_dim = width < height ? width : height;
   float offset_x = width < height ? 0 : (width - min_dim) * 0.5f;
   float offset_y = width < height ? (height - min_dim) * 0.5f : 0;
   float y0 = eye_y_;
   float z0 = eye_z_;
   float y1 = (static_cast<float>(y - offset_y) / min_dim) * 2.0f - 1.0f;
   float z1 = 0.0f;
   float dy = (y1 - y0);
@@ skipping 20 matching lines @@
     // calculate discriminant
     float disc = b * b - 4.0f * a * c;
 
     // Did ray hit the sphere?
     if (disc < 0.0f) {
       *pixels = kColorBlack;
       ++pixels;
       continue;
     }
 
-    // calc parametric t value
+    f32x4_t delta = {dx, dy, dz, 1.0f};
+    f32x4_t base = {x0, y0, z0, 1.0f};
+
+    // Calc parametric t value.
     float t = (-b - inline_sqrt(disc)) / (2.0f * a);
-    float px = x0 + t * dx;
-    float py = y0 + t * dy;
-    float pz = z0 + t * dz;
-    float nx = (px - planet_x_) * planet_one_over_radius_;
-    float ny = (py - planet_y_) * planet_one_over_radius_;
-    float nz = (pz - planet_z_) * planet_one_over_radius_;
+
+    f32x4_t pos = base + broadcast(t) * delta;
+    f32x4_t normal = (pos - planet_pos) * planet_one_over_radius;
 
     // Misc raytrace calculations.
-    float Lx = (light_x_ - px);
-    float Ly = (light_y_ - py);
-    float Lz = (light_z_ - pz);
-    float Lq = 1.0f / inline_quick_sqrt(Lx * Lx + Ly * Ly + Lz * Lz);
-    Lx *= Lq;
-    Ly *= Lq;
-    Lz *= Lq;
-    float d = (Lx * nx + Ly * ny + Lz * nz);
-    float pr = (diffuse_r_ * d) + ambient_r_;
-    float pg = (diffuse_g_ * d) + ambient_g_;
-    float pb = (diffuse_b_ * d) + ambient_b_;
-    float ds = -(nx * planet_pole_x_ +
-                 ny * planet_pole_y_ +
-                 nz * planet_pole_z_);
+    f32x4_t L = light_pos - pos;
+    float Lq = 1.0f / inline_quick_sqrt(dot3(L, L));
+    L = L * broadcast(Lq);
+    float d = dot3(L, normal);
+    f32x4_t p = diffuse * broadcast(d) + ambient;
+    float ds = -dot3(normal, planet_pole);
     float ang = acos_.TableLerp(ds);
     float v = ang * kOneOverPI;
-    float dp = planet_equator_x_ * nx +
-               planet_equator_y_ * ny +
-               planet_equator_z_ * nz;
-    float w = dp / sin(ang);
+    float dp = dot3(planet_equator, normal);
+    float w = dp / sinf(ang);
     if (w > 1.0f) w = 1.0f;
     if (w < -1.0f) w = -1.0f;
     float th = acos_.TableLerp(w) * kOneOver2PI;
-    float dps = planet_pole_x_equator_x_ * nx +
-                planet_pole_x_equator_y_ * ny +
-                planet_pole_x_equator_z_ * nz;
+    float dps = dot3(planet_pole_x_equator, normal);
     float u;
     if (dps < 0.0f)
       u = th;
     else
       u = 1.0f - th;
 
     // Look up daylight texel.
     int tx = static_cast<int>(u * base_tex_->width);
     int ty = static_cast<int>(v * base_tex_->height);
     int offset = tx + ty * base_tex_->width;
     uint32_t base_texel = base_tex_->pixels[offset];
-    float tr = ExtractR(base_texel);
-    float tg = ExtractG(base_texel);
-    float tb = ExtractB(base_texel);
-
-    float ipr = 1.0f - pr;

[binji, 2014/05/15 18:58:37]

You don't clamp this anymore?

[nfullagar, 2014/05/15 21:25:13]

good catch! (Its only slightly visual when adjusting the light way up.)
On 2014/05/15 18:58:37, binji wrote:

-    if (ipr < 0.0f) ipr = 0.0f;
-    float ipg = 1.0f - pg;
-    if (ipg < 0.0f) ipg = 0.0f;
-    float ipb = 1.0f - pb;
-    if (ipb < 0.0f) ipb = 0.0f;
+    f32x4_t dc = ExtractRGBA(base_texel);
 
     // Look up night texel.
     int nix = static_cast<int>(u * night_tex_->width);
     int niy = static_cast<int>(v * night_tex_->height);
     int noffset = nix + niy * night_tex_->width;
     uint32_t night_texel = night_tex_->pixels[noffset];
-    float nr = ExtractR(night_texel);
-    float ng = ExtractG(night_texel);
-    float nb = ExtractB(night_texel);
+    f32x4_t nc = ExtractRGBA(night_texel);
 
-    // Final color value is lerp between day and night texels.
-    unsigned int ir = Clamp255(pr * tr + nr * ipr);
-    unsigned int ig = Clamp255(pg * tg + ng * ipg);
-    unsigned int ib = Clamp255(pb * tb + nb * ipb);
+    // Blend between daylight (dc) and nighttime (nc) color.
+    f32x4_t fc = dc * p + nc * (kOne - p);
+    uint32_t color = PackBGRA(fc);
 
-    unsigned int color = MakeBGRA(ib, ig, ir, 0xFF);
-
-    *pixels = color;
+    *pixels = color | kSolidAlpha;
     ++pixels;
   }
 }
 
 // Renders a rectangular area of the screen, scan line at a time
 void Planet::wRenderRect(int x, int y, int w, int h) {
   for (int j = y; j < (y + h); ++j) {
     this->wRenderPixelSpan(x, x + w - 1, j);
   }
 }
@@ skipping 302 matching lines @@
     // Do simulation, render and present.
     earth.Update();
   }
 
   return 0;
 }
 
 // Register the function to call once the Instance Object is initialized.
 // see: pappi_simple/ps_main.h
 PPAPI_SIMPLE_REGISTER_MAIN(example_main);