квадратный корень на асме AVR

[avreal]

Дык и так можно

Код: Выделить всё

        subi    mask0, 0
        sbci    mask1, 0
        sbci    mask2, 0
        sbci    mask3, 0
        brne loop

А если маска в R24...R27, то и так тоже

Код: Выделить всё

        sbiw    mask0, 0
        sbci    mask2, 0
        sbci    mask3, 0
        brne loop

p.s. В любом случае -- это прямой перенос на asm алгоритма, сделанного под C.
Который делалася так, чтобы продемонстрировать идею РПП и получить разумно-быструю функцию для любой архитектуры без необходимости что-то дотачивать (это всё поначалу вообще на AT89C55 гонялось).
На асме можно и по-другому сделать, посмотрите тему в RU.ALGORITHMS -- с чего тема началась. Там «цифра за цифрой», что не очень уобно реализуется на C, но на асме, мне кажется, может позволить частично объединить регистры (как при делении, когда две переменные сдвигаются вместе и частное вдвигается на место делимого).
Кажется, именно такой алгоритм на асме для AVR где-то вте же годы написал Александр Труш, только там из 3-байтового числа корень извлекался.

[avreal]

Сейчас скомпилил в avr-gcc указанный С-шный исходник, так тело цикла 1:1 тому, что я дал выше (начиная от метки loop и заканичая сдвигом mask на два), а вот с организацией цикла он прави не прав одновремённо. У AVR-то регистров немеряно.
Конец цикла он проверяет не по 0 в mask, а тупо по счётчику, так как при данном начальном значении mask и сдвигах на 2 тело цикла будт пройдено 16 раз (кстати, в дів раза меньше, чем при делении 32/32, так что способ извлечения корня y = (x + y/x) / 2, требующий нескольких делений, явно отстаёт). Но он почемуто решил, что счётчик нужно сделать двухбайтовым.
Так что можно ещё так изменить

Код: Выделить всё

	... инициализация
	ldi	counter, 16 	; <-----------

loop:
	movw ...
	...
skip:
	...
	rol	mask3		; конец сдвига маски

	dec	counter 	; <-----------
	brne loop

Ещё тройка тактов в цикле экономится, на всём корне при 16МГц тактовой ещё микросекунды три выиграется. С учётом всех изменений микросекунд в 40 в наихудшем случае влезет. В среднем, думаю, 35.
И именно из этого уже ничего не выжать. Надо смтреть в сторону «цифра за цифрой» (бит за битом), читайте Винокурова в той теме.

[orinoko]

Жалко столько регистров использовать )). У меня ж они не только для расчёта используются. Ща проверим быстродействие.
И всё таки какой древний алгоритм, а раз выжил - значит хороший.

[orinoko]

Программа была переделана по вашим рекомендациям. Скорость расчёта улучшилась с 65 до 35 мкс. Вот конечная версия. На этом можно остановиться, я думаю. Всем большое спасибо.

[avreal]

Ну так если регистры очень нужны где-то ещё, то всё равно будет по скорости лучше и по объёму кода не хуже в начале четыре регистра push-нуть, задействовать их для вычислений, потом назад pop-нуть.

Вот глянул первоначальный вариант - там sqr инициализируетя нулями (4*sts = 16байт, но половина из них уйдёт и на инициализацию регистров) потом в цикле дважды считывается (8*lds = 32 байта) и один раз сохраняется (4*sts = 16 байт). По об́ъёму кода это во много раз больше, чем 4*push+4*pop (16байт). Ну а по скорости так само собой. Доп расхода памяи на стек нет, на это уйдёт та память, которая была нужна для временной переменной sqr, причём в другое время этот стек будет для чего-то другого использован.

[avreal]

И всё таки какой древний алгоритм

Если алгоритм древний, то я тогда какой?

[orinoko]

И всё таки какой древний алгоритм

Если алгоритм древний, то я тогда какой?

А если это вы автор этого алгоритма, то вы замечательный человек . Надеюсь, вы мне простите тогда, что я заюзал ваши авторские права.

[Psych]

народ мне интересно а если число 33 битное, то какое число нужно записать в mask по первому примеру? Я так догадываюсь что 0x400000000??

[orinoko]

Не факт. И, мне почему-то кажется, что исходное число должно иметь чётное количество бит

[avreal]

Число с нечётным числом бит всегда можно превратить в число с чётным их числом, добавив слева нулевой бит

В маске поднят младший бит из самой старшей пары бит.
Для 10-битного числа маска будет 0x100
33-битное число рассматривайте как 34-битное, маска будет 0x100000000ULL

p.s. Маска представляет собой квадрат числа с единственной единичкой. Начальная маска -- квадрат числа с единичкой в старшем разряде (максимального) ожидаемого корня.
Так и получается, что маска прыгает через два бита и всегда в младшем бите из пары.
0x8000 (старшая единичка в 16-битном числе результата) в квадрате давало маску 0x40000000.
Если описания по ссылкам почитать, это должно быть ясно.

[Psych]

Хм я так понимаю у этого алгоритма есть небольшя погрешность?

[Psych]

Народ хелп!! При числе 2С52С0000 и больше выдаёт какую-то хрень((

[Jack_A]

Хм я так понимаю у этого алгоритма есть небольшя погрешность?

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

Народ хелп!! При числе 2С52С0000 и больше выдаёт какую-то хрень((

Стартертопика интересовал 32-битный аргумент. А число 2С52С0000 уже 34-битное.

[avreal]

Народ хелп!! При числе 2С52С0000 и больше выдаёт какую-то хрень((

Алгоритм работает нормально, если какая-то программа выдаёт хрень -- надо искать в ней ошибку.
Раз число вылезло за 32 бита, нужно подправить типы и маску.

Код: Выделить всё

#include <inttypes.h>
#include <stdio.h>

uint32_t isqrt64(uint64_t from) {
     uint64_t mask = 0x4000000000000000ULL;
     uint64_t sqr = 0, temp;
     do {
         temp = sqr | mask;
         sqr >>= 1;
         if( temp <= from ) {
             sqr |= mask;
             from -= temp;
         }
     } while( mask >>= 2 );
     if( sqr < from && sqr < UINT32_MAX) ++sqr;
     return (uint32_t)sqr;
}

#define TST(a) printf("sqrt(0x%" PRIX64 ") = 0x%" PRIX32 "\n", (uint64_t)a, isqrt64(a) );

int main()
{
        TST( 4 );
        TST( 0x10000UL );
        TST( 0x100000000ULL );
        TST( 0x2C52C0000ULL );
        TST( 0x3FFFFFFFFULL );
        TST( 0x4F0000000ULL );
}

И выдача программы:

Код: Выделить всё

sqrt(0x4) = 0x2
sqrt(0x10000) = 0x100
sqrt(0x100000000) = 0x10000
sqrt(0x2C52C0000) = 0x1AA16
sqrt(0x3FFFFFFFF) = 0x20000
sqrt(0x4F0000000) = 0x238D8

Всё честно.

Вот делаем маску изначально меньше и ограничиваем маской диапазон входных чисел в 34 бита, куда влезет 2C52C0000 (макс. значение 34-битного входа 3 FFFF FFFF, а с учётом округления алгоритм будет выдавать правильные ответы до 4 0002 0000, но это ему просто повезёт по причине реально 64-битных вычислений):

Код: Выделить всё

#include <inttypes.h>
#include <stdio.h>

uint32_t isqrt64(uint64_t from) {
     uint64_t mask = 0x100000000ULL;
     uint64_t sqr = 0, temp;
     do {
         temp = sqr | mask;
         sqr >>= 1;
         if( temp <= from ) {
             sqr |= mask;
             from -= temp;
         }
     } while( mask >>= 2 );
     if( sqr < from && sqr < UINT32_MAX) ++sqr;
     return (uint32_t)sqr;
}

#define TST(a) printf("sqrt(0x%" PRIX64 ") = 0x%" PRIX32 "\n", (uint64_t)a, isqrt64(a) );

int main()
{
        TST( 4 );
        TST( 0x10000UL );
        TST( 0x100000000ULL );
        TST( 0x2C52C0000ULL );
        TST( 0x3FFFFFFFFULL );
        TST( 0x4F0000000ULL );
}

Код: Выделить всё

sqrt(0x4) = 0x2
sqrt(0x10000) = 0x100
sqrt(0x100000000) = 0x10000
sqrt(0x2C52C0000) = 0x1AA16
sqrt(0x3FFFFFFFF) = 0x20000
sqrt(0x4F0000000) = 0x20000

В последней строке ответ неправильный, но тут уже не выдержан диапазон аргумента, а алгоритм ни при чём.

[Psych]

Блин уже голова трещит....и все равно нифга не канает((

Собственно делаю вычиление: x=sqrt(a^2+b^2+c^2)

Теперь предположим что a=b=c=#F500 --- все нормально, результат правильный(и до него тоже)

А если a=b=c=#F600 () и далее , то все уже фигня выходит((

Маску менял и так и эдак...одна фигня((

[Psych]

Народ поясните подробнее while(mask>>=2)?

[Psych]

Все вопрос решен. Сам даже незнаю как, но все заработало

[avreal]

Это плохо. То, что непонятно как само заработало, с тем же успехо само и перестанет. Кто в доме хозяин?

[Psych]

Это плохо. То, что непонятно как само заработало, с тем же успехо само и перестанет. Кто в доме хозяин?

Че эт перестанет то?! . Уже проверил на всяких разных числах. Просто я чето поменял... и забыл че поменял))

[avreal]

Это и плохо, что забыл Грабли нужно помнить.
В следующий раз при редактировании может нечаянно назад поменяться. Впрочем, тогда и память может освежиться.