Более быстрый алгоритм 16-битного умножения для 8-битного MCU

Ответ 1

Резюме

• Рассмотрим замену a и b (исходное предложение)
• Попытка избежать условных переходов (Неудачная оптимизация)
• Изменение формы входной формулы (оценка 35%)
• Удаление дублированного сдвига
• Развертывание цикла: "Оптимальная" сборка
• Убедить компилятор предоставить оптимальную сборку

<ч/" > 1. Рассмотрим замену a и b

Одним из улучшений было бы сначала сравнить a и b и обменять их, если a<b: вы должны использовать как b меньший из двух, так что у вас есть минимальное количество циклов. Обратите внимание, что вы можете избежать подкачки путем дублирования кода (if (a<b) затем перейти к разделу зеркального кода), но я сомневаюсь, что это стоит.

<ч/" > 2. Попытка избежать условных переходов (Не успешная оптимизация)

Try:

uint16_t umul16_(uint16_t a, uint16_t b)
{
    ///Here swap if necessary
    uint16_t accum=0;

    while (b) {
        accum += ((b&1) * uint16_t(0xffff)) & a; //Hopefully this multiplication is optimized away
        b>>=1;
        a+=a;
    }

    return accum;
}

От обратной связи Sergio это не приносило улучшений.

<ч/" > 3. Изменение формы входной формулы

Учитывая, что целевая архитектура имеет в основном только 8-битные команды, если вы отделяете верхний и нижний 8 бит входных переменных, вы можете написать:

a = a1 * 0xff + a0;
b = b1 * 0xff + b0;

a * b = a1 * b1 * 0xffff + a0 * b1 * 0xff + a1 * b0 * 0xff + a0 * b0

Теперь здорово, что мы можем выбросить термин a1 * b1 * 0xffff, потому что 0xffff отправит его из вашего регистра.

(16bit) a * b = a0 * b1 * 0xff + a1 * b0 * 0xff + a0 * b0

Кроме того, термин a0*b1 и a1*b0 можно рассматривать как 8-битные умножения, поскольку 0xff: любая часть, превышающая 256, будет отправлена ​​из регистра.

До сих пор интересно!... Но здесь наступает реальная поразительная: a0 * b0 следует рассматривать как умножение на 16 бит, так как вам придется сохранять все результирующие биты. a0 должен быть сохранен на 16 бит, чтобы позволить сдвиг влево. Это умножение имеет половину итераций a * b, оно находится в части 8 бит (из-за b0), но вы все равно должны учитывать упомянутые выше 2 8-битные умножения и окончательную композицию результата. Нам нужна дальнейшая перестройка!

Итак, теперь я собираю b0.

(16bit) a * b = a0 * b1 * 0xff + b0 * (a0 + a1 * 0xff)

Но

(a0 + a1 * 0xff) = a

Итак, получаем:

(16bit) a * b = a0 * b1 * 0xff + b0 * a

Если N были циклами исходного a * b, теперь первое слагаемое представляет собой 8-битное умножение с N/2 циклами, а второе - умножение на 16 бит * 8 бит с N/2 циклами. Рассматривая M количество инструкций на итерацию в оригинальной a*b, 8-битная итерация 8 бит имеет половину инструкций, а 16 бит * 8 бит около 80% от M (одна команда смены меньше для b0 по сравнению с b). Соединяясь, мы имеем сложность N/2*M/2+N/2*M*0.8 = N*M*0.65, поэтому ожидаемая экономия ~ 35% по отношению к оригиналу N*M. Звучит многообещающий.

Это код:

uint16_t umul16_(uint16_t a, uint16_t b)
{
    uint8_t res1 = 0;

    uint8_t a0 = a & 0xff; //This effectively needs to copy the data
    uint8_t b0 = b & 0xff; //This should be optimized away
    uint8_t b1 = b >>8; //This should be optimized away

    //Here a0 and b1 could be swapped (to have b1 < a0)
    while (b1) {///Maximum 8 cycles
        if ( (b1 & 1) )
            res1+=a0;
        b1>>=1;
        a0+=a0;
    }

    uint16_t res = (uint16_t) res1 * 256; //Should be optimized away, it not even a copy!

    //Here swapping wouldn't make much sense
    while (b0) {///Maximum 8 cycles
        if ( (b0 & 1) )
            res+=a;
        b0>>=1;
        a+=a;
    }

    return res;
}

Кроме того, расщепление в 2 цикла должно теоретически удваивать вероятность пропустить некоторые циклы: N/2 может быть небольшим завышенным.

Небольшое дальнейшее улучшение состоит в том, чтобы избежать последнего ненужного сдвига для переменных a. Небольшая сторона примечания: если либо b0, либо b1 равны нулю, она вызывает 2 дополнительных инструкции. Но он также сохраняет первую проверку b0 и b1, что является самым дорогим, потому что он не может проверить статус zero flag из операции сдвига для условного перехода цикла for.

uint16_t umul16_(uint16_t a, uint16_t b)
{
    uint8_t res1 = 0;

    uint8_t a0 = a & 0xff; //This effectively needs to copy the data
    uint8_t b0 = b & 0xff; //This should be optimized away
    uint8_t b1 = b >>8; //This should be optimized away

    //Here a0 and b1 could be swapped (to have b1 < a0)
    if ( (b1 & 1) )
        res1+=a0;
    b1>>=1;
    while (b1) {///Maximum 7 cycles
        a0+=a0;
        if ( (b1 & 1) )
            res1+=a0;
        b1>>=1;
    }

    uint16_t res = (uint16_t) res1 * 256; //Should be optimized away, it not even a copy!

    //Here swapping wouldn't make much sense
    if ( (b0 & 1) )
        res+=a;
    b0>>=1;
    while (b0) {///Maximum 7 cycles
        a+=a;
        if ( (b0 & 1) )
            res+=a;
        b0>>=1;
    }

    return res;
}

<ч/" > 4. Удаление дублированного сдвига

Есть ли еще место для улучшения? Да, так как байты в a0 меняются два раза. Таким образом, должно быть преимущество в объединении двух циклов. Возможно, было бы немного сложно убедить компилятор делать именно то, что мы хотим, особенно с реестром результатов.

Итак, мы обрабатываем один цикл b0 и b1. Первое, что нужно сделать, это условие выхода цикла? Пока использование b0/b1 очищенного статуса было удобным, поскольку оно позволяет избежать использования счетчика. Кроме того, после правого сдвига флаг может быть уже установлен, если результат операции равен нулю, и этот флаг может разрешить условный переход без дальнейших оценок.

Теперь условие выхода цикла может быть сбой (b0 || b1). Однако это может потребовать дорогостоящих вычислений. Одним из решений является сравнение b0 и b1 и переход к двум различным разделам кода: if b1 > b0 Я проверяю условие на b1, иначе я проверяю условие на b0. Я предпочитаю другое решение с 2 циклами, первый выход, когда b0 равен нулю, второй, когда b1 равен нулю. Будут случаи, в которых я буду выполнять нулевые итерации в b1. Дело в том, что во втором цикле я знаю, что b0 равен нулю, поэтому я могу уменьшить количество выполненных операций.

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

uint16_t umul16_(uint16_t a, uint16_t b)
{
    uint16_t res = 0;

    uint8_t b0 = b & 0xff; //This should be optimized away
    uint8_t b1 = b >>8; //This should be optimized away

    //Swapping probably doesn't make much sense anymore
    if ( (b1 & 1) )
        res+=(uint16_t)((uint8_t)(a && 0xff))*256;
    //Hopefully the compiler understands it has simply to add the low 8bit register of a to the high 8bit register of res

    if ( (b0 & 1) )
        res+=a;

    b1>>=1;
    b0>>=1;
    while (b0) {///N cycles, maximum 7
        a+=a;
        if ( (b1 & 1) )
            res+=(uint16_t)((uint8_t)(a & 0xff))*256;
        if ( (b0 & 1) )
            res+=a;
        b1>>=1;
        b0>>=1; //I try to put as last the one that will leave the carry flag in the desired state
    }

    uint8_t a0 = a & 0xff; //Again, not a real copy but a register selection

    while (b1) {///P cycles, maximum 7 - N cycles
        a0+=a0;
        if ( (b1 & 1) )
            res+=(uint16_t) a0 * 256;
        b1>>=1;
    }
    return res;
}

Спасибо Sergio за предоставленную сборку (-Ofast). На первый взгляд, учитывая возмутительное количество mov в коде, кажется, что компилятор не интерпретировал, поскольку мне нужны подсказки, которые я ему дал, чтобы интерпретировать регистры.

Входы: r22, r23 и r24,25.
Набор инструкций AVR: Краткая ссылка, Подробная документация

sbrs //Tests a single bit in a register and skips the next instruction if the bit is set. Skip takes 2 clocks. 
ldi // Load immediate, 1 clock
sbiw // Subtracts immediate to *word*, 2 clocks

    00000010 <umul16_Antonio5>:
      10:    70 ff           sbrs    r23, 0
      12:    39 c0           rjmp    .+114        ; 0x86 <__SREG__+0x47>
      14:    41 e0           ldi    r20, 0x01    ; 1
      16:    00 97           sbiw    r24, 0x00    ; 0
      18:    c9 f1           breq    .+114        ; 0x8c <__SREG__+0x4d>
      1a:    34 2f           mov    r19, r20
      1c:    20 e0           ldi    r18, 0x00    ; 0
      1e:    60 ff           sbrs    r22, 0
      20:    07 c0           rjmp    .+14         ; 0x30 <umul16_Antonio5+0x20>
      22:    28 0f           add    r18, r24
      24:    39 1f           adc    r19, r25
      26:    04 c0           rjmp    .+8          ; 0x30 <umul16_Antonio5+0x20>
      28:    e4 2f           mov    r30, r20
      2a:    45 2f           mov    r20, r21
      2c:    2e 2f           mov    r18, r30
      2e:    34 2f           mov    r19, r20
      30:    76 95           lsr    r23
      32:    66 95           lsr    r22
      34:    b9 f0           breq    .+46         ; 0x64 <__SREG__+0x25>
      36:    88 0f           add    r24, r24
      38:    99 1f           adc    r25, r25
      3a:    58 2f           mov    r21, r24
      3c:    44 27           eor    r20, r20
      3e:    42 0f           add    r20, r18
      40:    53 1f           adc    r21, r19
      42:    70 ff           sbrs    r23, 0
      44:    02 c0           rjmp    .+4          ; 0x4a <__SREG__+0xb>
      46:    24 2f           mov    r18, r20
      48:    35 2f           mov    r19, r21
      4a:    42 2f           mov    r20, r18
      4c:    53 2f           mov    r21, r19
      4e:    48 0f           add    r20, r24
      50:    59 1f           adc    r21, r25
      52:    60 fd           sbrc    r22, 0
      54:    e9 cf           rjmp    .-46         ; 0x28 <umul16_Antonio5+0x18>
      56:    e2 2f           mov    r30, r18
      58:    43 2f           mov    r20, r19
      5a:    e8 cf           rjmp    .-48         ; 0x2c <umul16_Antonio5+0x1c>
      5c:    95 2f           mov    r25, r21
      5e:    24 2f           mov    r18, r20
      60:    39 2f           mov    r19, r25
      62:    76 95           lsr    r23
      64:    77 23           and    r23, r23
      66:    61 f0           breq    .+24         ; 0x80 <__SREG__+0x41>
      68:    88 0f           add    r24, r24
      6a:    48 2f           mov    r20, r24
      6c:    50 e0           ldi    r21, 0x00    ; 0
      6e:    54 2f           mov    r21, r20
      70:    44 27           eor    r20, r20
      72:    42 0f           add    r20, r18
      74:    53 1f           adc    r21, r19
      76:    70 fd           sbrc    r23, 0
      78:    f1 cf           rjmp    .-30         ; 0x5c <__SREG__+0x1d>
      7a:    42 2f           mov    r20, r18
      7c:    93 2f           mov    r25, r19
      7e:    ef cf           rjmp    .-34         ; 0x5e <__SREG__+0x1f>
      80:    82 2f           mov    r24, r18
      82:    93 2f           mov    r25, r19
      84:    08 95           ret
      86:    20 e0           ldi    r18, 0x00    ; 0
      88:    30 e0           ldi    r19, 0x00    ; 0
      8a:    c9 cf           rjmp    .-110        ; 0x1e <umul16_Antonio5+0xe>
      8c:    40 e0           ldi    r20, 0x00    ; 0
      8e:    c5 cf           rjmp    .-118        ; 0x1a <umul16_Antonio5+0xa>

<ч/" > 5. Развертывание цикла: "Оптимальная" сборка

Со всей этой информацией попробуйте понять, что было бы "оптимальным" решением, учитывая ограничения архитектуры. "Оптимальное" цитируется, потому что "оптимальный" во многом зависит от входных данных и того, что мы хотим оптимизировать. Предположим, что мы хотим оптимизировать количество циклов в худшем случае. Если мы идем в худшем случае, разворот цикла - разумный выбор: мы знаем, что у нас есть 8 циклов, и мы удаляем все тесты, чтобы понять, закончили ли мы (если b0 и b1 равны нулю). До сих пор мы использовали трюк "мы сдвигаемся, и мы проверяем флаг нуля", чтобы проверить, нужно ли нам выйти из цикла. Убрано это требование, мы можем использовать другой трюк: мы сдвигаемся, а проверяем бит переноса (бит, который мы отправили из регистра при переключении) , чтобы понять, должен ли я обновлять результат. Учитывая набор инструкций, в сборке "описательный" код инструкции становятся следующими.

//Input: a = a1 * 256 + a0, b = b1 * 256 + b0
//Output: r = r1 * 256 + r0

Preliminary:
P0 r0 = 0 (CLR)
P1 r1 = 0 (CLR)

Main block:
0 Shift right b0 (LSR)
1 If carry is not set skip 2 instructions = jump to 4 (BRCC)
2 r0 = r0 + a0 (ADD)
3 r1 = r1 + a1 + carry from prev. (ADC)
4 Shift right b1 (LSR)
5 If carry is not set skip 1 instruction = jump to 7 (BRCC)
6 r1 = r1 + a0 (ADD)
7 a0 = a0 + a0 (ADD)  
8 a1 = a1 + a1 + carry from prev. (ADC)

[Repeat same instructions for another 7 times]

Ветвление принимает 1 команду, если не вызван прыжок, 2 в противном случае. Все остальные инструкции - 1 цикл. Таким образом, состояние b1 не влияет на количество циклов, тогда как у нас есть 9 циклов, если b0 = 1 и 8 циклов, если b0 = 0. Подсчитая инициализацию, 8 итераций и пропустив последнее обновление a0 и a1, в худшем случае (b0 = 11111111b), мы имеем в общей сложности 8 * 9 + 2 - 2 = 72 цикла. Я бы не знал, какая реализация на С++ убедит компилятор в ее создании. Может быть:

 void iterate(uint8_t& b0,uint8_t& b1,uint16_t& a, uint16_t& r) {
     const uint8_t temp0 = b0;
     b0 >>=1;
     if (temp0 & 0x01) {//Will this convince him to use the carry flag?
         r += a;
     }
     const uint8_t temp1 = b1;
     b1 >>=1;
     if (temp1 & 0x01) {
         r+=(uint16_t)((uint8_t)(a & 0xff))*256;
     }
     a += a;
 }

 uint16_t umul16_(uint16_t a, uint16_t b) {
     uint16_t r = 0;
     uint8_t b0 = b & 0xff;
     uint8_t b1 = b >>8;

     iterate(b0,b1,a,r);
     iterate(b0,b1,a,r);
     iterate(b0,b1,a,r);
     iterate(b0,b1,a,r);
     iterate(b0,b1,a,r);
     iterate(b0,b1,a,r);
     iterate(b0,b1,a,r);
     iterate(b0,b1,a,r); //Hopefully he understands he doesn't need the last update for variable a
     return r;
 }

Но, учитывая предыдущий результат, чтобы действительно получить желаемый код, нужно действительно переключиться на сборку!

Наконец, можно было бы также рассмотреть более экстремальную интерпретацию разворота цикла: инструкции sbrc/sbrs позволяют протестировать определенный бит регистра. Поэтому мы можем избежать смещения b0 и b1, и на каждом цикле проверить бит. Единственная проблема заключается в том, что эти инструкции позволяют пропускать следующую команду, а не для пользовательского перехода. Итак, в "повествовательном коде" это будет выглядеть так:

Main block:
0 Test Nth bit of b0 (SBRS). If set jump to 2 (+ 1cycle) otherwise continue with 1
1 Jump to 4 (RJMP)
2 r0 = r0 + a0 (ADD)
3 r1 = r1 + a1 + carry from prev. (ADC)
4 Test Nth bit of (SBRC). If cleared jump to 6 (+ 1cycle) otherwise continue with 5
5 r1 = r1 + a0 (ADD)
6 a0 = a0 + a0 (ADD)  
7 a1 = a1 + a1 + carry from prev. (ADC)

В то время как вторая подстановка позволяет сохранить 1 цикл, нет явного преимущества во второй замене. Тем не менее, я считаю, что код С++ может быть проще интерпретировать для компилятора. Учитывая 8 циклов, инициализацию и пропустить последнее обновление a0 и a1, мы имеем теперь 64 цикла.

Код С++:

 template<uint8_t mask>
 void iterateWithMask(const uint8_t& b0,const uint8_t& b1, uint16_t& a, uint16_t& r) {
     if (b0 & mask)
         r += a;
     if (b1 & mask)
         r+=(uint16_t)((uint8_t)(a & 0xff))*256;
     a += a;
 }

 uint16_t umul16_(uint16_t a, const uint16_t b) {
     uint16_t r = 0;
     const uint8_t b0 = b & 0xff;
     const uint8_t b1 = b >>8;

     iterateWithMask<0x01>(b0,b1,a,r);
     iterateWithMask<0x02>(b0,b1,a,r);
     iterateWithMask<0x04>(b0,b1,a,r);
     iterateWithMask<0x08>(b0,b1,a,r);
     iterateWithMask<0x10>(b0,b1,a,r);
     iterateWithMask<0x20>(b0,b1,a,r);
     iterateWithMask<0x40>(b0,b1,a,r);
     iterateWithMask<0x80>(b0,b1,a,r);

     //Hopefully he understands he doesn't need the last update for a
     return r;
 }

Обратите внимание, что в этой реализации 0x01, 0x02 не являются реальным значением, а всего лишь подсказкой для компилятора, чтобы узнать, какой бит тестировать. Поэтому маска не может быть получена сдвигом вправо: иначе, чем все другие функции, которые были замечены до сих пор, на самом деле это не эквивалентная версия цикла.

Одна большая проблема заключается в том, что

r+=(uint16_t)((uint8_t)(a & 0xff))*256;

Он должен быть просто суммой верхнего регистра r с нижним регистром a. Не интерпретируется, как хотелось бы. Другой вариант:

r+=(uint16_t) 256 *((uint8_t)(a & 0xff));

<ч/" > 6. Убедить компилятор предоставить оптимальную сборку

Мы также можем сохранить константу a, а вместо нее - результат r. В этом случае мы обрабатываем b, начиная с самого значащего бита. Сложность эквивалентна, но компилятор может быть проще переваривать. Кроме того, на этот раз мы должны быть осторожны, чтобы написать явно последний цикл, который не должен делать дальнейший сдвиг вправо для r.

 template<uint8_t mask>
 void inverseIterateWithMask(const uint8_t& b0,const uint8_t& b1,const uint16_t& a, const uint8_t& a0, uint16_t& r) {
     if (b0 & mask)
         r += a;
     if (b1 & mask)
         r+=(uint16_t)256*a0; //Hopefully easier to understand for the compiler?
     r += r;
 }

 uint16_t umul16_(const uint16_t a, const uint16_t b) {
     uint16_t r = 0;
     const uint8_t b0 = b & 0xff;
     const uint8_t b1 = b >>8;
     const uint8_t a0 = a & 0xff;

     inverseIterateWithMask<0x80>(b0,b1,a,r);
     inverseIterateWithMask<0x40>(b0,b1,a,r);
     inverseIterateWithMask<0x20>(b0,b1,a,r);
     inverseIterateWithMask<0x10>(b0,b1,a,r);
     inverseIterateWithMask<0x08>(b0,b1,a,r);
     inverseIterateWithMask<0x04>(b0,b1,a,r);
     inverseIterateWithMask<0x02>(b0,b1,a,r);

     //Last iteration:
     if (b0 & 0x01)
         r += a;
     if (b1 & 0x01)
         r+=(uint16_t)256*a0;

     return r;
 }

Ответ 2

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

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

uint16_t umul16_(uint16_t op1, uint16_t op2)
{
    uint16_t accum=0;
    uint16_t a, b;
    a=op1; b=op2;
    if( op1<op2 ) { a=op2; b=0p1; } //swap operands to loop on smaller
    while (b) {
        accum += (b&1)?a:0;
        b>>=1;
        a+=a;
    }

    return accum;
}

Много лет назад я написал "четвертое", что способствовало вычислению, а не отраслевому подходу, и это предполагает выбор того, какую ценность использовать,

Вы можете использовать массив, чтобы полностью избежать теста, который ваш компилятор может использовать для генерации в качестве нагрузки со смещения. Определите массив, содержащий два значения, 0 и a, и обновите значение для a в конце цикла,

uint16_t umul16_(uint16_t op1, uint16_t op2)
{
    uint16_t accum=0;
    uint16_t pick[2];
    uint16_t a, b;
    a=op1; b=op2;
    if( op1<op2 ) { a=op2; b=0p1; } //swap operands to loop on smaller
    pick[0]=0;
    pick[1]=a;

    while (b) {
        accum += pick[(b&1)]; //avoid test completely
        b>>=1;
        pick[1] += pick[1]; //(a+=a);
    }

    return accum;
}

Да, зло. Но я обычно не пишу кода.

Изменить - изменено, чтобы добавить swap в цикл на меньшем op1 или op2 (меньшее количество проходов). Это исключило бы полезность тестирования для аргумента = 0.

Ответ 3

Ну, смесь LUT и shift обычно работает

Что-то вдоль линии, умножая 8-битные объекты. Предположим, что они состоят из двух квадов

uint4_t u1, l1, u2, l2;
uint8_t a = 16*u1 + l1;
uint8_t b = 16*u2 + l2;

product = 256*u1*u2 + 16*u1*l2 + 16*u2*l1 + l1*l1;

inline uint4_t hi( uint8_t v ) { return v >> 4; }
inline uint4_t lo( uint8_t v ) { return v & 15; }

inline uint8_t LUT( uint4_t x, uint4_t y ) {
    static uint8_t lut[256] = ...;
    return lut[x | y << 4]
}

uint16_t multiply(uint8_t a, uint8_t b) {
    return (uint16_t)LUT(hi(a), hi(b)) << 8 +
           ((uint16_t)LUT(hi(a), lo(b)) + (uint16_t)LUT(lo(a), hi(b)) << 4 +
           (uint16_t)LUT(lo(a), lo(b));
}

просто заполните lut [] с результатами умножения. В вашем случае в зависимости от памяти вы можете пойти с квадрациклом (размером 256 точек) или с байтами (65536 размер LUT) или что-нибудь среднее между

Ответ 4

Один из подходов - развернуть цикл. У меня нет компилятора для используемой платформы, поэтому я не могу смотреть на сгенерированный код, но такой подход мог бы помочь.

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

(Код примечания не тестировался, с моей головы. Мне интересно, как выглядит сгенерированный код!)

#define UMUL16_STEP(a, b, shift) \
    if ((b) & (1U << (shift))) result += ((a) << (shift)));

uint16_t umul16(uint16_t a, uint16_t b)
{
    uint16_t result = 0;

    UMUL16_STEP(a, b, 0);
    UMUL16_STEP(a, b, 1);
    UMUL16_STEP(a, b, 2);
    UMUL16_STEP(a, b, 3);
    UMUL16_STEP(a, b, 4);
    UMUL16_STEP(a, b, 5);
    UMUL16_STEP(a, b, 6);
    UMUL16_STEP(a, b, 7);
    UMUL16_STEP(a, b, 8);
    UMUL16_STEP(a, b, 9);
    UMUL16_STEP(a, b, 10);
    UMUL16_STEP(a, b, 11);
    UMUL16_STEP(a, b, 12);
    UMUL16_STEP(a, b, 13);
    UMUL16_STEP(a, b, 14);
    UMUL16_STEP(a, b, 15);

    return result;
}

Update:

В зависимости от того, что делает ваш компилятор, макрос UMUL16_STEP может измениться. Альтернативой может быть:

#define UMUL16_STEP(a, b, shift) \
    if ((b) & (1U << (shift))) result += (a); (a) << 1;

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

Мое предположение о том, как должен выглядеть ассемблер на бит, r0: r1 - результат, r2: r3 - a и r4: r5 - b:

sbrc r4, 0
add r0, r2
sbrc r4, 0
addc r1, r3
lsl r2
rol r3

Это должно выполняться в постоянное время без ветки. Проверьте биты в r4, а затем проверьте бит в r5 для более высоких восьми бит. Это должно выполнить умножение в 96 циклов на основе моего чтения руководства по набору инструкций.

Ответ 5

Неответ, ассемблер tinyARM (веб-документ) вместо С++ или C. Я изменил довольно общий multiply-by-squares-lookup для скорости (< 50 циклов, за исключением вызовов и возвратных издержек) за счет установки только в AVR с размером не менее 1 Кбайт ОЗУ с использованием 512 выровненных байтов для таблицу нижней половины квадратов. На частоте 20 МГц это прекрасно соответствовало бы временному пределу 2 max 3 usec, который еще не появился в правильном вопросе, но Серхио Формиггини хотел 16 МГц. По состоянию на 2015/04 год есть только один ATtiny от Atmel с такой большой оперативной памятью, и это указано до 8 МГц & hellip; (Роллинг вашего "собственного" (например, от OpenCores), ваш FPGA, вероятно, имеет кучу быстрых множителей (18 × 18 бит кажутся популярными), если не процессорные ядра.)
Посмотрите на смену и добавьте сдвиг фактора влево, разверните 16 × 16 → 16 и/или улучшите его (например, wiki post). (Возможно, вы можете создать ответ на этот вопрос в вики-сообществе, заданный в вопросе.)

.def    a0  = r16   ; factor low byte
.def    a1  = r17
#warning two warnings about preceding definitions of
#warning  r16 and r17 are due and may as well be ignored
.def    a   = r16   ; 8-bit factor
.def    b   = r17   ; 8-bit factor ; or r18, rather?
.def    b0  = r18   ; factor low byte
.def    b1  = r19
.def    p0  = r20   ; product low byte
.def    p1  = r21

; "squares table" SqTab shall be two 512 Byte tables of
;  squares of 9-bit natural numbers, divided by 4

; Idea: exploit p = a * b = Squares[a+b] - Squares[a-b]

init:
    ldi     r16, 0x73
    ldi     r17, 0xab
    ldi     r18, 23
    ldi     r19, 1
    ldi     r20, HIGH(SRAM_SIZE)
    cpi     r20, 2
    brsh    fillSqTable ; ATtiny 1634?
    rjmp    mpy16T16
fillSqTable:
    ldi     r20, SqTabH
    subi    r20, -2
    ldi     zh, SqTabH
    clr     zl
; generate sqares by adding up odd numbers starting at 1 += -1
    ldi     r22, 1
    clr     r23
    ser     r26
    ser     r27
fillLoop:
    add     r22, r26
    adc     r23, r27
    adiw    r26, 2
    mov     r21, r23
    lsr     r21         ; get bits 9:2
    mov     r21, r22
    ror     r21
    lsr     r21
    bst     r23, 1
    bld     r21, 7
    st      z+, r21
    cp      zh, r20
    brne    fillLoop
    rjmp    mpy16F16

; assembly lines are marked up with cycle count
;  and (latest) start cycle in block.
; If first line in code block, the (latest) block start cycle
;  follows; else if last line, the (max) block cycle total

;**************************************************************
;*
;* "mpy16F16" - 16x16->16 Bit Unsigned Multiplication
;*                        using table lookup
;* Sergio Formiggini special edition
;* Multiplies  two 16-bit register values a1:a0 and b1:b0.
;* The result is placed in p1:p0.
;*
;* Number of flash words: 318 + return = 
;*                       (40 + 256(flash table) + 22(RAM init))
;* Number of cycles     : 49 + return
;* Low  registers used  : None
;* High registers used  : 7+2 (a1:a0, b1:b0, p1:p0, sq;
;*                             + Z(r31:r30))
;* RAM bytes used       : 512 (squares table)
;*
;**************************************************************
mpy16F16:
    ldi     ZH, SqTabH>>1;1 0   0   squares table>>1
    mov     ZL, a0      ; 1 1
    add     ZL, b0      ; 1 2       a0+b0
    rol     ZH          ; 1 3       9 bit offset
    ld      p0, Z       ; 2 4       a0+b0l          1
    lpm     p1, Z       ; 3 6   9   a0+b0h          2

    ldi     ZH, SqTabH  ; 1 0   9   squares table

    mov     ZL, a1      ; 1 0   10
    sub     ZL, b0      ; 1 1       a1-b0
    brcc    noNegF10    ; 1 2
    neg     ZL          ; 1 3
noNegF10:
    ld      sq, Z       ; 2 4       a1-b0l          3
    sub     p1, sq      ; 1 6   7

    mov     ZL, a0      ; 1 0   17
    sub     ZL, b1      ; 1 1       a0-b1
    brcc    noNegF01    ; 1 2
    neg     ZL          ; 1 3
noNegF01:
    ld      sq, Z       ; 2 4       a0-b1l          4
    sub     p1, sq      ; 1 6   7

    mov     ZL, a0      ; 1 0   24
    sub     ZL, b0      ; 1 1       a0-b0
    brcc    noNegF00    ; 1 2
    neg     ZL          ; 1 3
noNegF00:
    ld      sq, Z       ; 2 4       a0-b0l          5
    sub     p0, sq      ; 1 6
    lpm     sq, Z       ; 3 7       a0-b0h          6*
    sbc     p1, sq      ; 1 10  11

    ldi     ZH, SqTabH>>1;1 0   35
    mov     ZL, a1      ; 1 1
    add     ZL, b0      ; 1 2       a1+b0
    rol     ZH          ; 1 3
    ld      sq, Z       ; 2 4       a1+b0l          7
    add     p1, sq      ; 1 6   7

    ldi     ZH, SqTabH>>1;1 0   42
    mov     ZL, a0      ; 1 1
    add     ZL, b1      ; 1 2       a0+b1
    rol     ZH          ; 1 3
    ld      sq, Z       ; 2 4       a0+b1l          8
    add     p1, sq      ; 1 6   7

    ret                 ;       49

.CSEG
.org 256; words?!
SqTableH:
.db   0,   0,   0,   0,   0,   0,   0,   0,   0,   0
.db   0,   0,   0,   0,   0,   0,   0,   0,   0,   0
.db   0,   0,   0,   0,   0,   0,   0,   0,   0,   0
.db   0,   0,   1,   1,   1,   1,   1,   1,   1,   1
.db   1,   1,   1,   1,   1,   1,   2,   2,   2,   2
.db   2,   2,   2,   2,   2,   2,   3,   3,   3,   3
.db   3,   3,   3,   3,   4,   4,   4,   4,   4,   4
.db   4,   4,   5,   5,   5,   5,   5,   5,   5,   6
.db   6,   6,   6,   6,   6,   7,   7,   7,   7,   7
.db   7,   8,   8,   8,   8,   8,   9,   9,   9,   9
.db   9,   9,  10,  10,  10,  10,  10,  11,  11,  11
.db  11,  12,  12,  12,  12,  12,  13,  13,  13,  13
.db  14,  14,  14,  14,  15,  15,  15,  15,  16,  16
.db  16,  16,  17,  17,  17,  17,  18,  18,  18,  18
.db  19,  19,  19,  19,  20,  20,  20,  21,  21,  21
.db  21,  22,  22,  22,  23,  23,  23,  24,  24,  24
.db  25,  25,  25,  25,  26,  26,  26,  27,  27,  27
.db  28,  28,  28,  29,  29,  29,  30,  30,  30,  31
.db  31,  31,  32,  32,  33,  33,  33,  34,  34,  34
.db  35,  35,  36,  36,  36,  37,  37,  37,  38,  38
.db  39,  39,  39,  40,  40,  41,  41,  41,  42,  42
.db  43,  43,  43,  44,  44,  45,  45,  45,  46,  46
.db  47,  47,  48,  48,  49,  49,  49,  50,  50,  51
.db  51,  52,  52,  53,  53,  53,  54,  54,  55,  55
.db  56,  56,  57,  57,  58,  58,  59,  59,  60,  60
.db  61,  61,  62,  62,  63,  63,  64,  64,  65,  65
.db  66,  66,  67,  67,  68,  68,  69,  69,  70,  70
.db  71,  71,  72,  72,  73,  73,  74,  74,  75,  76
.db  76,  77,  77,  78,  78,  79,  79,  80,  81,  81
.db  82,  82,  83,  83,  84,  84,  85,  86,  86,  87
.db  87,  88,  89,  89,  90,  90,  91,  92,  92,  93
.db  93,  94,  95,  95,  96,  96,  97,  98,  98,  99
.db 100, 100, 101, 101, 102, 103, 103, 104, 105, 105
.db 106, 106, 107, 108, 108, 109, 110, 110, 111, 112
.db 112, 113, 114, 114, 115, 116, 116, 117, 118, 118
.db 119, 120, 121, 121, 122, 123, 123, 124, 125, 125
.db 126, 127, 127, 128, 129, 130, 130, 131, 132, 132
.db 133, 134, 135, 135, 136, 137, 138, 138, 139, 140
.db 141, 141, 142, 143, 144, 144, 145, 146, 147, 147
.db 148, 149, 150, 150, 151, 152, 153, 153, 154, 155
.db 156, 157, 157, 158, 159, 160, 160, 161, 162, 163
.db 164, 164, 165, 166, 167, 168, 169, 169, 170, 171
.db 172, 173, 173, 174, 175, 176, 177, 178, 178, 179
.db 180, 181, 182, 183, 183, 184, 185, 186, 187, 188
.db 189, 189, 190, 191, 192, 193, 194, 195, 196, 196
.db 197, 198, 199, 200, 201, 202, 203, 203, 204, 205
.db 206, 207, 208, 209, 210, 211, 212, 212, 213, 214
.db 215, 216, 217, 218, 219, 220, 221, 222, 223, 224
.db 225, 225, 226, 227, 228, 229, 230, 231, 232, 233
.db 234, 235, 236, 237, 238, 239, 240, 241, 242, 243
.db 244, 245, 246, 247, 248, 249, 250, 251, 252, 253
.db 254, 255
; word addresses, again?!
.equ SqTabH = (high(SqTableH) << 1)

.DSEG
RAMTab .BYTE 512

Ответ 6

Наконец-то ответ, если назойливый: я не смог (пока) получить AVR-C-компилятор из GCC, поместив его в код 8K. (Для получения ассемблера см. Умножение AVR: Нет Holds Barred).
Подход - это то, что все, кто использовал устройство Duff, попытались сделать вторую попытку:
используйте переключатель. Используя макросы, исходный код выглядит совершенно безвредным, если массируется:

#define low(mp)     case mp: p = a0 * (uint8_t)(mp) << 8; break
#define low4(mp)    low(mp); low(mp + 1); low(mp + 2); low(mp + 3)
#define low16(mp)   low4(mp); low4(mp + 4); low4(mp + 8); low4(mp + 12)
#define low64(mp)   low16(mp); low16(mp + 16); low16(mp + 32); low16(mp + 48)
#if preShift
# define CASE(mp)   case mp: return p + a * (mp)
#else
# define CASE(mp)   case mp: return (p0<<8) + a * (mp)
#endif
#define case4(mp)   CASE(mp); CASE(mp + 1); CASE(mp + 2); CASE(mp + 3)
#define case16(mp)  case4(mp); case4(mp + 4); case4(mp + 8); case4(mp + 12)
#define case64(mp)  case16(mp); case16(mp + 16); case16(mp + 32); case16(mp + 48)

extern "C" __attribute__ ((noinline))
 uint16_t mpy16NHB16(uint16_t a, uint16_t b)
{
    uint16_t p = 0;
    uint8_t b0 = (uint8_t)b, b1 = (uint8_t)(b>>8);
    uint8_t a0 = (uint8_t)a, p0;

    switch (b1) {
        case64(0);
        case64(64);
        case64(128);
        case64(192);
    }
#if preShift
    p = p0 << 8;
#endif
#if preliminaries
    if (0 == b0) {
        p = -a;
        if (b & 0x8000)
            p += a << 9;
        if (b & 0x4000)
            p += a << 8;
        return p;
    }
    while (b0 & 1) {
        a <<= 1;
        b0 >>= 1;
    }
#endif
    switch (b0) {
        low64(0);
        low64(64);
        low64(128);
        low64(192);
    }
    return ~0;
}
int main(int ac, char const *const av[])
{
    char buf[22];
    for (uint16_t a = 0 ; a < a+1 ; a++)
      for (uint16_t m = 0 ; m <= a ; m++)
        puts(itoa(//shift4(ac)+shift3MaskAdd((uint16_t)av[0], ac)
    //      +shift4Add(ac, (uint16_t)av[0])
    //           + mpy16NHB16(ac, (ac + 105))
                 mpy16NHB16(a, m), buf, 10));
}