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