UPD 26.10.2017
Вывод алгоритма
Известно, что:
Квадрат суммы нескольких чисел - это сумма их квадратов, плюс удвоенная сумма всех их попарных произведений
Квадрат бита (0 | 1) равен этому биту
Пусть n = b020 + b121 + b222 + ... + bi-12i-1,
тогда n2 = (b020 + b121 + b222 + ... + bi-12i-1)2, где i - разрядность числа n.
n2 = b020 + b122 + b224 + ... + bi-122i-2 + b041(b120 + b221 + b322
+ ... + bi-12i-2) + b142(b220 + b321 + b422 + ... + bi-12i-3) + b243(b320 + b421 + b522 + ... + bi-12i-4) + ... + bi-24i-1bi-1,
n2 = b0((1 << 0) | (n >> 1) << 2) + b1((1 << 1) | (n >> 2) << 4) + b2((1 << 2) | (n >> 3) << 6) + ...
+ bi-2((1 << (i-2) | (n >> (i-1)) << (2i-2)) + bi-1((1 << (i-1) | (n >> i) << (2i).
Пусть jk = ((2n) >> k), тогда bk = ((jk+1 & 1) > 0) = ((jk & 2 > 0),
n2 = (j0 & 2 ? (1 << 0) | (j2 << 2) : 0) + (j1 & 2 ? (1 << 1) | (j3 << 4) : 0)
+ (j2 & 1 ? (1 << 2) | (j4 << 6) : 0) + ... + (ji-2 & 2 ? (1 << (i-2)) | (ji << (2i-2) : 0) + (ji-1 & 2 ? (1 << (i-1) | (ji+1 << (2i) : 0).
Однако jk+2 << 2 = (jk >> 2) << 2 = jk & 0b111...100 = jk & (-4) (сброс нулевого и первого битов), поэтому:
n2 = (j0 & 2 ? (j0 & -4 | 1) << 0 : 0) + (j1 & 2 ? (j1 & -4 | 1) << 2 : 0)
+ (j2 & 1 ? (j2 & -4 | 1) << 4 : 0) + ... + (ji-2 & 2 ? (ji-2 & -4 | 1) << (2i-4) : 0) + (ji-1 & 2 ? (ji-1 & -4 | 1) << (2i -2) : 0).
Программная реализация для чисел формата integer
В программе (на языке PHP) алгоритм пирамидального возведения в квадрат реализован в виде функции pow2()
и протестирован для чисел формата integer. Функция pow2($i)
в 2.5
раза проигрывает по времени стандартной процедуре pow($i, 2)
, использующей аппаратное умножение, что можно считать хорошим результатом. В 64-битовой версии PHP7 проигрыш больше (~10 раз).
Текст программы:
define("N", 99); // Количество точек в тестовом массиве.
define("PR", 1); // PR=0 - тест производительности.
// Вывод одномерного массива
function print_1d($text, $v){
$eps = 1e-7;
print "$text"."[";
$cnt = -1;
foreach($v as $key => $item){
if(++$cnt) print ", ";
$flag = true;
if($cnt && !($cnt%5)) print"<br>";
print ""$key"=>";
if(abs((int)$item - $item) > $eps){
printf("%.3f", $item);
}else{
print $item;
}
}
print "], ";
}
// Пирамидальное возведение в квадрат
function pow2($n){
if($n < 0) $n = -$n;
$res = 0;
for($j = $n << 1, $sh = 0; $j > 1; $j >>= 1, $sh += 2){
$res += ($j & 2) ? ($j & -4 | 1) << $sh : 0;
}
return $res;
}
$test = [];
for($i = 0; $i < N; $i++){
$test[] = mt_rand(-9999, 9999);
}
$time_start = microtime(true);
$pow2_ = [];
for($i = 0; $i < N; $i++){
$pow2_[$test[$i]] = pow2($test[$i]);
}
$time1 = microtime(true) - $time_start;
$time_start = microtime(true);
$pow2__ = [];
for($i = 0; $i < N; $i++){
$pow2__[$test[$i]] = pow($test[$i], 2);
}
$time2 = microtime(true) - $time_start;
printf("Возведение в квадрат для N = %d случайных точек", (int)N);
if(PR > 0){
print_1d("<br>*** Пирамидальный алгоритм ***<br>Результаты:<br>", $pow2_);
print_1d("<br>*** Библиотечный алгоритм ***<br>Результаты:<br>", $pow2__);
}else{
print("<br>*** Пирамидальный алгоритм ***<br>Время выполнения: $time1");
print("<br>*** Библиотечный алгоритм ***<br>Время выполнения: $time2");
}
Проверка на случайной выборке:
Возведение в квадрат для N = 99 случайных точек
*** Пирамидальный алгоритм ***
Результаты:
["8154"=>66487716, "993"=>986049, "-3669"=>13461561, "7855"=>61701025, "-7597"=>57714409,
"1568"=>2458624, "-9029"=>81522841, "-6646"=>44169316, "-4934"=>24344356, "-3954"=>15634116,
"-5303"=>28121809, "70"=>4900, "1661"=>2758921, "-3628"=>13162384, "5734"=>32878756,
"4802"=>23059204, "-1886"=>3556996, "4269"=>18224361, "9638"=>92891044, "-7879"=>62078641,
"9218"=>84971524, "-2221"=>4932841, "-7153"=>51165409, "5861"=>34351321, "-3656"=>13366336,
"3058"=>9351364, "-8879"=>78836641, "593"=>351649, "329"=>108241, "-4801"=>23049601,
"-8637"=>74597769, "4480"=>20070400, "9918"=>98366724, "8584"=>73685056, "-8525"=>72675625,
"1920"=>3686400, "5019"=>25190361, "8697"=>75637809, "-5594"=>31292836, "-8198"=>67207204,
"-7119"=>50680161, "-6288"=>39538944, "131"=>17161, "320"=>102400, "-5410"=>29268100,
"9147"=>83667609, "-7411"=>54922921, "-2211"=>4888521, "3990"=>15920100, "-2077"=>4313929,
"-6588"=>43401744, "-9280"=>86118400, "7357"=>54125449, "-1909"=>3644281, "6109"=>37319881,
"911"=>829921, "8898"=>79174404, "-4262"=>18164644, "6537"=>42732369, "-3567"=>12723489,
"1345"=>1809025, "-9229"=>85174441, "6654"=>44275716, "-9949"=>98982601, "9057"=>82029249,
"-8760"=>76737600, "779"=>606841, "-9905"=>98109025, "45"=>2025, "-8160"=>66585600,
"-6796"=>46185616, "-5270"=>27772900, "8019"=>64304361, "-1922"=>3694084, "1542"=>2377764,
"-2489"=>6195121, "-4835"=>23377225, "4752"=>22581504, "-9621"=>92563641, "782"=>611524,
"6617"=>43784689, "4506"=>20304036, "-7679"=>58967041, "-1476"=>2178576, "-6194"=>38365636,
"-4385"=>19228225, "7394"=>54671236, "1323"=>1750329, "870"=>756900, "-6723"=>45198729,
"4162"=>17322244, "-5436"=>29550096, "-7761"=>60233121, "-596"=>355216, "-7463"=>55696369,
"356"=>126736, "429"=>184041, "-1011"=>1022121, "9811"=>96255721],
*** Библиотечный алгоритм ***
Результаты:
["8154"=>66487716, "993"=>986049, "-3669"=>13461561, "7855"=>61701025, "-7597"=>57714409,
"1568"=>2458624, "-9029"=>81522841, "-6646"=>44169316, "-4934"=>24344356, "-3954"=>15634116,
"-5303"=>28121809, "70"=>4900, "1661"=>2758921, "-3628"=>13162384, "5734"=>32878756,
"4802"=>23059204, "-1886"=>3556996, "4269"=>18224361, "9638"=>92891044, "-7879"=>62078641,
"9218"=>84971524, "-2221"=>4932841, "-7153"=>51165409, "5861"=>34351321, "-3656"=>13366336,
"3058"=>9351364, "-8879"=>78836641, "593"=>351649, "329"=>108241, "-4801"=>23049601,
"-8637"=>74597769, "4480"=>20070400, "9918"=>98366724, "8584"=>73685056, "-8525"=>72675625,
"1920"=>3686400, "5019"=>25190361, "8697"=>75637809, "-5594"=>31292836, "-8198"=>67207204,
"-7119"=>50680161, "-6288"=>39538944, "131"=>17161, "320"=>102400, "-5410"=>29268100,
"9147"=>83667609, "-7411"=>54922921, "-2211"=>4888521, "3990"=>15920100, "-2077"=>4313929,
"-6588"=>43401744, "-9280"=>86118400, "7357"=>54125449, "-1909"=>3644281, "6109"=>37319881,
"911"=>829921, "8898"=>79174404, "-4262"=>18164644, "6537"=>42732369, "-3567"=>12723489,
"1345"=>1809025, "-9229"=>85174441, "6654"=>44275716, "-9949"=>98982601, "9057"=>82029249,
"-8760"=>76737600, "779"=>606841, "-9905"=>98109025, "45"=>2025, "-8160"=>66585600,
"-6796"=>46185616, "-5270"=>27772900, "8019"=>64304361, "-1922"=>3694084, "1542"=>2377764,
"-2489"=>6195121, "-4835"=>23377225, "4752"=>22581504, "-9621"=>92563641, "782"=>611524,
"6617"=>43784689, "4506"=>20304036, "-7679"=>58967041, "-1476"=>2178576, "-6194"=>38365636,
"-4385"=>19228225, "7394"=>54671236, "1323"=>1750329, "870"=>756900, "-6723"=>45198729,
"4162"=>17322244, "-5436"=>29550096, "-7761"=>60233121, "-596"=>355216, "-7463"=>55696369,
"356"=>126736, "429"=>184041, "-1011"=>1022121, "9811"=>96255721],
Тестирование быстродействия (PHP5):
Возведение в квадрат для N = 100000 случайных точек
*** Пирамидальный алгоритм ***
Время выполнения: 1.29107379913
*** Библиотечный алгоритм ***
Время выполнения: 0.49302816391
Рекомендации для алгоритма больших чисел
Работа с большими числами ограничивается следующим небольшим фрагментом:
for($j = 2*$n, $sh = 0; $j > 1; $j >>= 1, $sh += 2){
if($j & 2){
$res += ($j ^ 3 | 1) << $sh;
}
}
Рекомендации для больших чисел:
- Формат чисел - двоичный. При использовании массива слов возникнут проблемы со знаковыми битами и хранением реальной длины числа. Более гибким является символьный формат с кодированием и декодированием через функции
chr()
и ord()
, поскольку в нём имеется доступ к каждому байту.
- При сравнении
$j > 1
достаточно длины чисел в битах. При тестировании бита $j & $2
и формировании нулевого и первого битов достаточно использовать только младший фрагмент числа.
- Операции сдвига и суммирования для
$n
должны быть реализованы в полном объёме.
unsigned int
вместоchar
- единицей вычисления является машинное слово, а не один байт, к тому же при использовании знакового типа могут возникнуть проблемы. Для хранения использовать массив. Для арифметических операций представьте как это всё вычисляет машина, придётся это реализовать через "умножение(и не только) в столбик", самому следить за переносами и прочим. – Yuriy Orlov 22 фев '16 в 9:57