Asm.js практика

от автора

image
Этим прохладным днём я искал алгоритмы и реализации вычисления числа пи. Алгоритмов нашлось какое-то несметное множество, но тут нашёлся пост с описанием алгоритма и его реализацией на си.
Алгоритм подкупает своей скоростью, хоть и выдаёт 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.

А именно

time ./pi 10000000
>> real 0m25.345s
>> user 0m25.176s
>> sys 0m0.019s

console.time(‘s’);
pi(1000000);
console.timeEnd(‘s’);
>> s: 1868.5мс

И ещё:
time ./pi 100000000
>> real 0m25.345s
>> user 0m25.176s
>> sys 0m0.019s

console.time(‘s’);
pi(10000000);
console.timeEnd(‘s’);
>> s: 22152.83мс

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

ссылка на оригинал статьи http://habrahabr.ru/post/193642/


Комментарии

Добавить комментарий

Ваш адрес email не будет опубликован. Обязательные поля помечены *