Asm.js практика

Этим прохладным днём я искал алгоритмы и реализации вычисления числа пи. Алгоритмов нашлось какое-то несметное множество, но тут нашёлся пост с описанием алгоритма и его реализацией на си.
Алгоритм подкупает своей скоростью, хоть и выдаёт hex представление, но так уж вышло, что мне нужен был вариант на js. Моментальная, практически, переработка на обычный js показала очень плохую статистику, работа при подсчёте 1000000-ого знака заняла… 48 секунд (4ГГц FF).
О том, как возился с asmjs и каких камней повстречал можно узнать под катом.

Для нетерпеливых, результат на гитхабе.

После беглого просмотра стало понятно, что нам не нужно выносить работу со всей генерацией в модуль asm, а используются только две функции: expm и series. Т.к. expm вызывается внутри series, то нам из модуля следует экспортировать только series.

double series (int m, int id)

/*  This routine evaluates the series  sum_k 16^(id-k)/(8*k+m) 
    using the modular exponentiation technique. */

{
  int k;
  double ak, eps, p, s, t;
  double expm (double x, double y);
#define eps 1e-17

  s = 0.;

/*  Sum the series up to id. */

  for (k = 0; k < id; k++){
    ak = 8 * k + m;
    p = id - k;
    t = expm (p, ak);
    s = s + t / ak;
    s = s - (int) s;
  }

/*  Compute a few terms where k >= id. */

  for (k = id; k <= id + 100; k++){
    ak = 8 * k + m;
    t = pow (16., (double) (id - k)) / ak;
    if (t < eps) break;
    s = s + t;
    s = s - (int) s;
  }
  return s;
}

double expm (double p, double ak)

/*  expm = 16^p mod ak.  This routine uses the left-to-right binary 
    exponentiation scheme. */

{
  int i, j;
  double p1, pt, r;
#define ntp 25
  static double tp[ntp];
  static int tp1 = 0;

/*  If this is the first call to expm, fill the power of two table tp. */

  if (tp1 == 0) {
    tp1 = 1;
    tp[0] = 1.;

    for (i = 1; i < ntp; i++) tp[i] = 2. * tp[i-1];
  }

  if (ak == 1.) return 0.;

/*  Find the greatest power of two less than or equal to p. */

  for (i = 0; i < ntp; i++) if (tp[i] > p) break;

  pt = tp[i-1];
  p1 = p;
  r = 1.;

/*  Perform binary exponentiation algorithm modulo ak. */

  for (j = 1; j <= i; j++){
    if (p1 >= pt){
      r = 16. * r;
      r = r - (int) (r / ak) * ak;
      p1 = p1 - pt;
    }
    pt = 0.5 * pt;
    if (pt >= 1.){
      r = r * r;
      r = r - (int) (r / ak) * ak;
    }
  }

  return r;
}

Базовый шаблон.

По сути мы создаём чёрный ящик с некоторым интерфейсом, поэтому всё, что мы можем, так передать что-либо в модуль и получить из него на выходе набор методом и/или значений.
В просмотренных мной кодах устоялась конструкция вида:

(function (window) {
  "use strict";
  // переменные
  var buffer = new ArrayBuffer(BUFFER_SIZE); // буфер для работы представлений типизированного массива

  var functionNameOrStuctureName = (function (stdlib, env, heap) {
    "use asm";
    // переменные
    // тело модуля

    return {
      methodNameExport: methodNameInModule,
      methodName2Export: methodName2InModule,
    }; // или просто return methodNameInModule
  })(
    { // классы и объекты (stdlib)
      Uint8Array:  Uint8Array,
      Int8Array:   Int8Array,
      Uint16Array: Uint16Array,
      Int16Array:  Int16Array,
      Uint32Array: Uint32Array,
      Int32Array:  Int32Array,
      Float32Array:Float32Array,
      Float64Array:Float64Array,
      Math:        Math
    },
    { // переменные (env)
      NTP:NTP
    },
    buffer // и буфер, крайне немаловажен, при чём размер > 4096 и равен степени дв0йки
  );
})(window); // wrapper end

Бывает вылетает ошибка вида «TypeError: asm.js type error: asm.js module must end with a return export statement», удостоверьтесь, что модуль возвращает что-либо. Если же возвращает как положено, то следует убедиться в правильности деклараций переменных. У меня была ошибка, когда после декларации я пытался что-то ещё делать с одной интовой переменной.
Надеюсь базовые вещи уже успели узнать, но всё же:

  function f1(i, d) {
    i = i | 0; // integer заявляем, что i  целочисленная переменная
    d = +d; // double заявляем, что d  переменная с плавающей точкой
    var char = new Uint8Array(heap); // строки, в данной статье рассмотрены не будут
    var i2 = 0; // объявляем целочисленную переменную (на самом деле она сейчас fixnum)
    var d2 = 0.0; // объявляем переменную с плавающей точкой
    i2 = i2 | 0; // конвертируем fixnum в integer
    return +d2; // функция имеет тип double  и может возвращать только double
  }

Подводный камень №1: переменные и функции модуля

Я привык сначала декларировать все переменные, а уже потом присваивать им значения. В asm.js всё несколько хитрее: сначала мы декларируем все переменные, которые используем через замыкания, при чём функции математической библиотеки тоже надо описать здесь, а не вызывать ниже.
Вот что вышло:

   "use asm"; 
    // об stdlib выше
    var floor = stdlib.Math.floor; // некоторые инициализируют тут все ссылки на функции, но я ленив
    var pow = stdlib.Math.pow; // поэтому тут только те, что использовал ниже
    var tp1 = 0; // этот флаг используется ниже
    var tp = new stdlib.Float64Array(heap); // представление типизированного массива
    var ntp = env.NTP | 0; // для работы с глобальной переменной используем env, о нём выше

Функции же нельзя создавать предварительно инициализируя переменную и присваивая ей функцию. Поэтому следом идут функции.

    function expm(p, ak) {
        // тело
    }

Подводный камень №2: переменные в функциях модуля.

А вот в теле функций модуля сначала следует указать тип переменных, принимаемых на входе, потом инициализировать внутренние переменные, при чём если это int, то var i = 0, если double, то var d = 0.0, если массив, то через new и тип мссива с передачей в него heap, а уже после инициализации интов советую их «доинициализировать» путём присвоения вида i = i|0. К слову: инициализация переменных заранее в стиле си не обязательна. Числа вида 1e-17 выдают ошибку выхода за границы, используйте 0.00000000000000001
На выходе:

    function expm(p, ak) {
      p = +p;
      ak = +ak;
      var i = 0;
      var pt = 0.0;
      // ....
      i = i | 0;
      ntp = ntp | 0;
      // тело
    }

Подводный камень №3: сравнения и итераторы

Скажу честно, тут я смеялся долго, но int и int сравнению не подлежат. Если вы сделаете что-то вроде i == k или i < l, то вывалится компиляция с ошибкой вида
«TypeError: asm.js type error: arguments to a comparison must both be signed, unsigned or doubles, int and int are given».
Ещё немного смеху добавило сравнение int и целого числа (i ==0)
«TypeError: asm.js type error: arguments to a comparison must both be signed, unsigned or doubles, fixnum and int are given».
Только у чисел с плавающей точкой всё хорошо (например pt == 1.0).
В итоге, если вы хотите всё-таки сравнить int с другим int. надо продекларировать, конструкцию вида (i | 0) < (ntp | 0).
Что касается итераторов, то тут просто «прелесть»: вместо всем нам привычного i++ мы имеем i = (i | 0) + 1 | 0.
Результат:

      if ((tp1 | 0) == 0) {
        // vars
        for (i = 1; (i | 0) < (ntp | 0); i = (i | 0) + 1 | 0) {
          // surprise, read more)
        }
      }

Подводный камень №4: Математические и прочие внешние функции.

Тут всё просто, если вы хотите использовать floor, sin и т.п., то вам нужно декларировать их в начале модуля (сразу после «use asm»). Если в функции написать stdlib.Math.floor, у меня выкидывал ошибку возвращаемого типа. Видимо из-за обращения к свойствам объекта.

Подводный камень №5: Буферы.

А вот тут всё очень и очень интересно. Что мы делаем, когда хотим получить/установать значение из/в массива arr с индексом i? arr[i]. Допустим мы так и поступим, но тогда мы получим ошибку вида.

          arr[i] = +1;

«TypeError: asm.js type error: index expression isn't shifted; must be an Int8/Uint8 access»
Нам тонко намеают, что должен быть сдвиг. У одного гуру я нашёл сдвиг на 2 вправо.

          arr[i >> 2] = +1;

«TypeError: asm.js type error: shift amount must be 3»
Нам как бы тонко намекают, что он должно быть трём.

          arr[i << 3 >> 3] = +1;

Выходит в итоге. Сдвигом в лево мы скомпенсировали сдвиг в право. Вроде всё тоже самое, а работает.

/**
 * Created with JetBrains WebStorm.
 * User: louter
 * Date: 12.09.13
 * Time: 17:49
 */
(function (window) {
  "use strict";
  var ihex;
  var NTP = 25;
  var buffer = new ArrayBuffer(1024 * 1024 * 8);

  var series = (function (stdlib, env, heap) {
    "use asm";
    var floor = stdlib.Math.floor;
    var pow = stdlib.Math.pow;
    var tp1 = 0;
    var tp = new stdlib.Float64Array(heap);
    var ntp = env.NTP | 0;

    function expm(p, ak) {
      p = +p;
      ak = +ak;
      var i = 0;
      var j = 0;
      var p1 = 0.0;
      var pt = 0.0;
      var r = 0.0; // float as double
      i = i | 0;
      j = j | 0;
      ntp = ntp | 0;

      if ((tp1 | 0) == 0) {
        tp1 = 1 | 0;
        tp[0] = +1;
        for (i = 1; (i | 0) < (ntp | 0); i = (i | 0) + 1 | 0) {
          tp[(i << 3) >> 2] = +(+2 * tp[((i - 1) << 3) >> 3]);
        }
      }
      if (ak == 1.0) {
        return +0;
      }
      for (i = 0; (i | 0) < (ntp | 0); i = (i | 0) + 1 | 0) {
        if (+tp[i << 3 >> 3] > p) {
          break;
        }
      }
      pt = +tp[(i - 1) << 3 >> 3];
      p1 = +p;
      r = +1;
      for (j = 0; (j | 0) <= (i | 0); j = (j | 0) + 1 | 0) {
        if (p1 >= pt) {
          r = +16 * r;
          r = r - (+(floor(r / ak))) * ak;
          p1 = p1 - pt;
        }
        pt = 0.5 * pt;
        if (pt >= 1.) {
          r = r * r;
          r = r - (+(floor(r / ak))) * ak;
        }
      }

      return +r;
    }

    function series(m, id) {
      m = m | 0;
      id = id | 0;
      var k = 0;
      var ak = 0.0;
      var eps = 0.0;
      var p = 0.0;
      var s = 0.0;
      var t = 0.0;
      eps = 0.00000000000000001;
      k = 0 | 0;
      for (k; (k | 0) < (id | 0); k = (k | 0) + 1 | 0) {
        ak = +8 * (+(k | 0)) + (+(m | 0));
        p = (+(id | 0) - +(k | 0));
        t = +expm(p, ak);
        s = s + t / ak;
        s = s - floor(s);
      }
      for (k = (id | 0); (k | 0) <= ((id + 100) | 0); k = (k | 0) + 1 | 0) {
        ak = +8 * (+(k | 0)) + (+(m | 0));
        t = pow(+16, +(id | 0) - (+(k | 0))) / +ak;
        if (t < eps) break;
        s = s + t;
        s = s - floor(s);
      }
      return +s;
    }

    return series;
  })(
    {
      Uint8Array:  Uint8Array,
      Int8Array:   Int8Array,
      Uint16Array: Uint16Array,
      Int16Array:  Int16Array,
      Uint32Array: Uint32Array,
      Int32Array:  Int32Array,
      Float32Array:Float32Array,
      Float64Array:Float64Array,
      Math:        Math
    },
    {
      NTP:NTP
    },
    buffer
  );
  ihex = function (x, nhx, chx) {
    var i, y, hx = "0123456789ABCDEF";

    y = Math.abs(x);
    for (i = 0; i < nhx; i++) {
      y = 16. * (y - (y | 0));
      chx[i] = hx[y | 0];
    }
  };
  window.pi = function (id) {
    var pid, s1, s2, s3, s4
      , hex = [];
    s1 = series(1, id);
    s2 = series(4, id);
    s3 = series(5, id);
    s4 = series(6, id);
    pid = 4 * s1 - 2 * s2 - s3 - s4;
    pid = pid - (pid | 0) + 1;
    ihex(pid, 16, hex);

    return {
      hex:     hex.join('').substr(0, 10),
      fraction:pid
    };
  };
})(window);

P.S. Я уж не знаю, кто или что не правы, но(!) почему-то скомпилированная программа выдала результаты хуже, чем asm.js.

Не верите? Проверьте. Исходник программы в вышеописанном посте. Исходники мои так же выше указаны.

Автор: Louter

Источник