Вычисление целочисленного корня
От: ravik  
Дата: 10.02.21 13:37
Оценка: :)
Наткнулся на статью на Хабре, в которой автор из-за отсутствия под рукой софта (он программировал 12-разрядный контроллер), пустился во все тяжкие, чтобы заменить не хватающее ему взятие корня неким прикладным маневром. Заинтересовался и набросал решение взятия целочисленного корня из произвольной длины массива. Для иллюстрации здесь решение редуцировано до 32-битных значений (см. в скрытом тексте, описание алгоритма — в комментариях), чтобы сам алгоритм не затерялся в каше обработки пересечений границ блоков.

  Скрытый текст
void TestSQR(void)
{
    unsigned int pv = 14777;    //подкоренное значение
    unsigned int sqr = 0;        //извлеченный корень
    __asm
    {
        mov edx, pv;            //загружаем исходное значение
        xor edi, edi;            //в EDI будет накапливаться вычисляемый корень
        bsr ecx, edx;            //вычисляем позицию старшего значащего разряда (ПСЗР)
        jz QUIT;
        shr ecx, 1;
        mov edi, 1;
        shl edi, cl;            //получаем установленный старший бит корня
        shl edi, cl;            //
        sub edx, edi;            //эти 3 команды сокращают цикл на 2 итерации
        shr edi, cl;            //
        mov ch, cl;
START:  bsr eax, edx;            //вычисляем ПСЗР уменьшающегося с каждым шагом цикла pv
        jz QUIT;                //pv обнулился - значит, корень извлекся нацело
        mov cl, al;
        sub cl, ch;
        jz QUIT;                //без этой строки алгоритм ошибается на корне из 2 и 3
        jc QUIT;
        dec cl;
        jz LAST_BIT;            //если остаток pv >= 2 * sqr + 1, то добавляем к sqr 1
        dec cl;                    //в CL вычислен индекс разряда корня, в который надо добавить 1
        mov eax, 1;                
        shl eax, cl;            //в итерации до возврата к метке START
        mov ebx, edi;            //происходит вычитание из текущего значения pv
        shl ebx, cl;            //вычисленного на лету значения 2 * sqr + 1 
        sub edx, ebx;            //и смещенного в правильную по отношению к pv позицию.
        lea ebx, [edi+eax];        //Вместо умножения на 2 производится 2 вычитания:
        mov edi, ebx;            //сначала без добавленной единицы в правильный разряд,
        shl ebx, cl;            //затем уже вместе с ней
        sub edx, ebx;
        jmp START;
LAST_BIT:                        //проверка, что остаток не превышает 2 * sqr + 1
        lea ebx, [edi*2];
        cmp edx, ebx;
        jbe QUIT;
        inc ebx;
        sub edx, ebx;
        inc edi;
QUIT:
        mov sqr, edi;
        mov pv, edx;
    }
    WCHAR sz[256];
    wsprintf(sz, L"корень: %d\nостаток: %d\n", sqr, pv);
    ::OutputDebugString(sz);
    return;
}


Но возник вопрос разбирающимся в матобеспечении. Взаправдашний алгоритм вычисления квадратного корня программным способом, например, в матбиблиотеках, — примерно похож?! Я до сего момента довольствовался пионерским методом усреднения множителей, знаю еще метод цепных дробей... А вот вышеприведенного велосипеда — не встречал... Удовлетворите любопытство, плж!
Отредактировано 12.02.2021 5:56 ravik . Предыдущая версия . Еще …
Отредактировано 10.02.2021 14:02 ravik . Предыдущая версия .
Отредактировано 10.02.2021 13:53 ravik . Предыдущая версия .
Re: Вычисление целочисленного корня
От: ravik  
Дата: 11.02.21 07:18
Оценка: 5 (2)
Дело в том, что разрабатывая этот очевидный велосипед, я пришел к интересным формулировкам, к собственно извлечению квадратного корня отношения не имеющим. Никогда не пробовали отскочивший на вашу сторону стола целлулоидный шарик прижать ракеткой к поверхности? Амплитуда, частота, всё такое... Когда-то давно программировал умножение задом-наперед — ну возникла потребность оценивать, что больше x^n, или y^m при огромных значениях как X и Y, так и показателей степеней. "Задом-наперед" — это когда ты получаешь сначала старшие биты и их индексы x^n и y^m, и чаще всего сравнение быстро заканчивается. У одного из значений индекс старшего разряда больше?! Ну так все-же понятно! Ах, индексы одинаковы? Так у одного второй по старшинству разряд будет ноль, а у второго 1... И так далее.

И вот там была такая фишка, с которой я до конца не мог разобраться — проект, сроки и т.д.

На определенной итерации становится очевидно, что предшествующий старший бит уже не поменяется, как бы не складывались дальнейшие обстоятельства. Обычно это происходило на относительно небольшом заглублении от него, максимум 5-7 разрядов. Разумеется, существуют исключения, и ты не имеешь подтверждения, пока не вычислишь результат умножения полностью до гребаного миллионного байта — но на то они и исключения. А обычно — 5-7 разрядов. И все, конкретный разряд отливается в гранит.

И сщас при программировании этого велосипеда, отступая от индекса старшего разряда на одну и ту же величину: Ipv — (Isqr + 2), аж похолодел. А если +3? Или даже +8?! Ёпрст, велико произведение, но отступать-то дальше некуда! Сходимость исчезнет. А 5-7 разрядов это максимум 128 вариантов, каждый из которых фиксирует дребезг чеканутого теннисного мячика так, что мама не горюй. Скажем так, исключения, которые не зафиксируют текущий статус-кво, весьма и весьма немногочисленны. Так немногочисленны, что их и в этих 128 вариантах-то не окажется.

Вот такие размышления, пребывая на излете профессиональной карьеры в восхитительной синекуре...
Re: Вычисление целочисленного корня
От: ravik  
Дата: 04.06.21 11:10
Оценка:
Я извиняюсь, что поднимаю, возможно, мало кому интересную тему, но я тут добился некоторого прогресса... Собственно, алгоритм работал с самого начала правильно, просто я тогда, записав по пунктам, что могло быть в нем оптимизировано, добрался до последнего. И перед тем как прыгнуть, решил с кем-то поболтать, для кого вопросы оптимизации не последние.


Надо оптимизировать вот этот цикл:
J0:     sub ebx, esp;
J1:     inc cl;
        adc edi, edi;
        js J3;
J2:     shl eax, 1;
        adc ebx, ebx;
        setc ch;
        add ebx, esp;
        adc ch, 0;
        jz J0;
        btr edi, cl;
        cmp ebx, edi;
        jb CORRECTION_RESTORE;
        bts edi, cl;
        shr ch, 1;
        jmp J1;

Чего здесь происходит? Перед циклом загружается в EBX:EAX два старших 32-разрядных блока числа длинной арифметики (содержащие значимые разряды, я их называю "вершина числа"). Содержимое вершины "прижимается" в направлении границы старших разрядов EBX. Например, было загружено 0x00008000:00000001, после сдвига влево, например, непосредственно перед циклом сдвиг на 16 разрядов EBX:EAX = 0x80000000:00010000. Здесь надо отметить, что старший значащий разряд вершины не обязан попадать вот прямо в старший регистровый. От расчетного сдвига зависит.

В ESP находится так называемый "измерительный блок". В контексте вершины, он изначально прижат к границе старших разрядов, и все его разряды содержательны с одной оговоркой. Его старший разряд всегда единичный, что позволило мне пройти предыдущий пункт оптимизации, используя отрицательное значение измерительного блока, по ходу цикла не вычитая, а складывая с текущей вершиной (add ebx, esp) и получая нужное значение флага переполнения (CF) в случае равенства.

Алгоритм простой — в каждой итерации измерительный блок сравнивается с вершиной. Если меньше — вычитается и в результат (EDI) записывается в младший разряд "1". Для этого после вычитания мы должны добраться до инструкции jmp J1. Если меньше вершина, то в результат пишется "0" (jz J0). Запись "0" или "1" производится инструкцией adc edi, edi. После записи результата вершина сдвигается слево двумя инструкциями (shl eax, 1; adc ebx, ebx; — Прим. почему инструкция ADC, а не RCL?! А RCL по-разному работает на процессорах AMD и Intel) После записи "1" производится проверка, достаточен ли результат вычитания для достоверности результата? Для этого временно снимается SWAP-бит в EDI (отвечает за своевременную подкачку в EBX очередного блока — инструкции btr edi, cl/bts edi, cl)

Что тут можно оптимизировать? Собственно, весь цикл прогибался под эту финальную оптимизацию. После вычитания может так получиться, что до ближайшего единичного бита вершины может образоваться несколько нулевых старших разрядов. Это означает, что придется проходить несколько полных итераций цикла с записями "0". Вот такие итерации и надо как-то хитро не то чтобы совсем исключить, но пройти по укороченной схеме без вычитания-восстановления (add ebx, esp; и sub ebx, esp; ). Я так представляю вот прямо после восстановления SWAP-бита в EDI надо чего-то выдумать, а не тупо отправлять все на метку J1. Идея с BSR мне не нравится — сразу отмечу.
Re: Вычисление целочисленного корня
От: kov_serg Россия  
Дата: 04.06.21 11:36
Оценка:
Здравствуйте, ravik, Вы писали:

R>А вот вышеприведенного велосипеда — не встречал... Удовлетворите любопытство, плж!


Копайте там мног интересно

https://ru.wikipedia.org/wiki/%D0%A6%D0%B5%D0%BB%D0%BE%D1%87%D0%B8%D1%81%D0%BB%D0%B5%D0%BD%D0%BD%D1%8B%D0%B9_%D0%BA%D0%B2%D0%B0%D0%B4%D1%80%D0%B0%D1%82%D0%BD%D1%8B%D0%B9_%D0%BA%D0%BE%D1%80%D0%B5%D0%BD%D1%8C
https://stackoverflow.com/questions/1100090/looking-for-an-efficient-integer-square-root-algorithm-for-arm-thumb2

еще тут посмотрите

https://en.wikipedia.org/wiki/Fast_inverse_square_root
https://github.com/id-Software/DOOM-3/blob/master/neo/idlib/math/Math.h#L300
Re[2]: Вычисление целочисленного корня
От: ravik  
Дата: 04.06.21 11:59
Оценка:
Здравствуйте, kov_serg, Вы писали:

_>Копайте там мног интересно


Спасибо большое! Но со времени записи стартового топика 4 месяца назад я много чего почитал. Я уже знаю, что реализую алгоритм Неймана, описанный в книге Генри Уоррена и предложенный аж в 1945 году. Я, правда, немного его продвинул, добавив блок оценки достоверности текущего результата... Но принципиально — идея точно его.


Реализация большая — почти 400 ассемблерных инструкций, прошла все возможные тестирования. Скорость устраивает более чем, осталось выполнить жест художника
Автор: ravik
Дата: 04.06.21

Re[2]: Вычисление целочисленного корня
От: Videoman Россия https://hts.tv/
Дата: 10.11.21 22:28
Оценка:
Здравствуйте, ravik, Вы писали:

R>Когда-то давно программировал умножение задом-наперед...


Про корень к сожалению не подскажу, но вот ваши соображения по умножению задом наперед заинтересовало. Вообще я для своих нужд делал библиотеку для работы с длинными целыми (slimcpplib) и в одной из операций мне как раз нужно было реализовать функцию: mulhi(x, y). Функция должна выдавать старшую половину произведения x на y. Я был уверен что оптимизировать её никак нельзя, так как старшие разряды могут зависеть от младших. Т.е. в процессе обратного умножения мы можем получить 0xfffffff..... и потом где-то в самом младшем разряде у нас перенос и все биты пошли поехали меняться с единиц но нули. Ваша логика мне тоже вроде понята, если у нас есть где-то 0 в разряде, то перенос на нем закончится. И вот я подумал, может быть вы подскажите или у вас есть идея, как можно оптимизировать мою функцию, что бы использовать ваш подход, или я не прав и тут ничего не сделать уже?
Re[3]: Вычисление целочисленного корня
От: ravik  
Дата: 19.11.21 08:01
Оценка:
Здравствуйте, Videoman, Вы писали:

V>Здравствуйте, ravik, Вы писали:


R>>Когда-то давно программировал умножение задом-наперед...


V>Про корень к сожалению не подскажу, но вот ваши соображения по умножению задом наперед заинтересовало. Вообще я для своих нужд делал библиотеку для работы с длинными целыми (slimcpplib) и в одной из операций мне как раз нужно было реализовать функцию: mulhi(x, y). Функция должна выдавать старшую половину произведения x на y. Я был уверен что оптимизировать её никак нельзя, так как старшие разряды могут зависеть от младших. Т.е. в процессе обратного умножения мы можем получить 0xfffffff..... и потом где-то в самом младшем разряде у нас перенос и все биты пошли поехали меняться с единиц но нули. Ваша логика мне тоже вроде понята, если у нас есть где-то 0 в разряде, то перенос на нем закончится. И вот я подумал, может быть вы подскажите или у вас есть идея, как можно оптимизировать мою функцию, что бы использовать ваш подход, или я не прав и тут ничего не сделать уже?


Привет!

Ну да, оптимизировать невозможно. Но!

bool KMultiple::operator== (IKInt* pi)
{
    if (this->HighBit() != pi->HighBit())
        return false;
    for (int i = m_iOld; i > 0; i--)
        if (this->BitValue(i) != pi->BitValue(i))
            return false;
    return true;
}


... выпендреж с умножением задом наперед делался, чтобы реализовать этот оператор сравнения. Класс KMultiple как раз и реализовывал умножение задом наперед, и делал это не до финала 0xFFFFFFF....F, а только до позиции i, передаваемой в метод KMultiple::BitValue(). По цитируемому тексту метода видно, как алгоритм продвигается от старшего бита к младшему до тех пор, пока биты совпадают с аналогичными по позиции сравниваемого числа. Т.е., умножение задом наперед реализовывалось не абстрактно, а с конкретной прикладной целью.

P.S.: И да KMultiple::BitValue() вычислит-таки самый младший бит, если до него дело дойдет. И все это реализовано по-честному и прекрасно работает до чисел с миллиардом десятичных цифр. Помнится, это было > 300Mb в бинарном эквиваленте.
Re[4]: Вычисление целочисленного корня
От: Videoman Россия https://hts.tv/
Дата: 19.11.21 21:35
Оценка:
Здравствуйте, ravik, Вы писали:

R>...


Ну понятно, т.е. можно оптимизировать если количество бит очень велико. Жаль
Re: Вычисление целочисленного корня
От: fk0 Россия https://fk0.name
Дата: 02.04.22 10:05
Оценка:
Здравствуйте, ravik, Вы писали:

R>Наткнулся на статью на Хабре, в которой автор из-за отсутствия под рукой софта (он программировал 12-разрядный контроллер), пустился во все тяжкие, чтобы заменить не хватающее ему взятие корня неким прикладным маневром. Заинтересовался и набросал решение взятия целочисленного корня из произвольной длины массива. Для иллюстрации здесь решение редуцировано до 32-битных значений (см. в скрытом тексте, описание алгоритма — в комментариях), чтобы сам алгоритм не затерялся в каше обработки пересечений границ блоков.


Нихера не понял, зачем кому-то читать твой непонятный ассемблер? Квадратный корень вычисляется,
оптимальным образом, в зависимости от платформы:

1) при наличии только целочисленного АЛУ -- методом CORDIC ("цифра за цифрой");

2) при наличии FPU -- методом Ньютона например.
 
Подождите ...
Wait...
Пока на собственное сообщение не было ответов, его можно удалить.