Выполнение 32-битной арифметики на встроенном устройстве

Вы написали:

16 наиболее значимых должны быть точными.

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

Недавно я работал над реализацией с фиксированной точкой квадратный корень. Он предназначен для работы с числами в диапазоне [0, 1), представленными в 0.16 формат, т.е. без знака, 0 бит целая часть, 16 бит дробная часть. Она основана на полиномиальном приближении

√x ≈ 0,1908901 + 1,5392302 x − 1,4475870 x² + 1,0217778 x³ − 0,3043110 x⁴

что весьма неплохо для x ∈ [0.25, 1]. Числа в (0, 0.25) могут быть переместился в этот интервал путем последовательных умножений на 4, а затем разделив результат на 2. Другими словами, используя тот факт, что

√x = √(4 x) / 2

Максимальная ошибка составляет 2,56e-4, т.е. 16,8 ulps (единиц в последнем место). В худшем случае для выполнения потребуется 298 циклов ЦП. Примечание что для достижения этой производительности мне пришлось реализовать умножение в ассемблере. Иначе gcc дал бы мне неэффективный код.

Вот код:

/*
 * Convert x to a 16-bit fixed point number
 * with n bits after the binary point.
 */
#define FIXED(x, n) ((uint16_t)((float)(x) * (1UL << (n)) + .5))

/*
 * Fixed point multiplication.
 *
 * Multiply two fixed point numbers in unsigned 0.16 format.
 * Returns result in the same format.
 * Rounds to nearest, ties rounded up.
 */
static uint16_t mul_fix_u16(uint16_t x, uint16_t y)
{
    uint16_t result;
#if defined(__AVR_HAVE_MUL__) && defined(__AVR_HAVE_MOVW__)
    /* Optimized ASM version. */
    asm volatile(
        "mul  %B1, %B2\n\t"
        "movw %A0, r0\n\t"
        "ldi  r19, 0x80\n\t"
        "clr  r18\n\t"
        "mul  %A1, %A2\n\t"
        "add  r19, r1\n\t"
        "adc  %A0, r18\n\t"
        "adc  %B0, r18\n\t"
        "mul  %B1, %A2\n\t"
        "add  r19, r0\n\t"
        "adc  %A0, r1\n\t"
        "adc  %B0, r18\n\t"
        "mul  %A1, %B2\n\t"
        "add  r19, r0\n\t"
        "adc  %A0, r1\n\t"
        "adc  %B0, r18\n\t"
        "clr  r1"
        : "=&r" (result)
        : "r" (x), "r" (y)
        : "r18", "r19"
    );
#else
    /* Generic C version. Compiles to inefficient 32 bit code. */
    result = ((uint32_t) x * y + 0x8000) >> 16;
#endif
    return result;
}

/*
 * Fixed point square root.
 * Argument and result in unsigned 0.16 format.
 *
 * If defined(REALLY_ACCURATE):
 *   max error = 7.67e-6
 *   execution time:
 *     937   cycles (worst case)
 *     828.1 cycles (average)
 * else
 *   max error = 2.56e-4
 *   execution time:
 *     298   cycles (worst case)
 *     184.7 cycles (average)
 */
uint16_t sqrt_fix(uint16_t x)
{
    if (x == 0) return 0;
    uint8_t n;
#ifdef REALLY_ACCURATE
    if (x == 0xffff) return 0xffff;
    uint16_t x0 = x;
#endif
    for (n = 0; x < 0x4000; n++) x <<= 2;  // scale the operand
    uint16_t y;
    y = FIXED(0.30435, 15);
    y = FIXED(1.02173, 15) - mul_fix_u16(x, y);
    y = FIXED(1.44754, 15) - mul_fix_u16(x, y);
    y = FIXED(1.53925, 15) - mul_fix_u16(x, y);
    y = FIXED(0.19089, 16) + (mul_fix_u16(x, y) << 1);
    while (n--) y >>= 1;                   // scale the result
#ifdef REALLY_ACCURATE
    y = (y + ((uint32_t)x0 << 16) / y + 1) / 2;  // Babylonian iteration
#endif
    return y;
}

Если вам действительно нужен очень точный результат, вы можете #define REALLY_ACCURATE. Это добавит одну итерацию Вавилонский метод, который по сути является методом Ньютона, примененным к квадратный корень:

y ← (y + x/y) / 2

где y — это наше приближение √x. Это выглядит как легкая победа: просто деление, сложение и сдвиг бита. И это делает результат практически идеально, с максимальной ошибкой 7,67e-6 (0,502 ulps). Но время выполнения увеличивается более чем в три раза, при 937 циклов в худшем случае. Это потому, что чип AVR не имеет аппаратная поддержка для разделения и программная альтернатива, предоставляемая компилятор ужасно медленный.