Линейная интерполяция с плавающей запятой

Если не принимать во внимание разницу в точности, это выражение эквивалентно

float lerp(float a, float b, float f)
{
    return a + f * (b - a);
}

Это 2 сложения/вычитания и 1 умножение вместо 2 сложений/вычитаний и 2 умножения.

Предполагая, что математика с плавающей запятой доступна, алгоритм OP является хорошим и всегда превосходит альтернативный a + f * (b - a) из-за потери точности, когда a и b значительно различаются по величине.

Например:

// OP's algorithm
float lint1 (float a, float b, float f) {
    return (a * (1.0f - f)) + (b * f);
}

// Algebraically simplified algorithm
float lint2 (float a, float b, float f) {
    return a + f * (b - a);
}

В этом примере, предполагая, что 32-битные числа с плавающей точкой lint1(1.0e20, 1.0, 1.0) правильно вернут 1,0, тогда как lint2 неправильно вернет 0,0.

Большая часть потерь точности приходится на операторы сложения и вычитания, когда операнды значительно различаются по величине. В приведенном выше случае виновниками являются вычитание в b - a и сложение в a + f * (b - a). Алгоритм OP не страдает от этого, поскольку компоненты полностью перемножаются перед добавлением.

Для случая a=1e20, b=1 приведен пример разных результатов. Программа испытаний:

#include <stdio.h>
#include <math.h>

float lint1 (float a, float b, float f) {
    return (a * (1.0f - f)) + (b * f);
}

float lint2 (float a, float b, float f) {
    return a + f * (b - a);
}

int main () {
    const float a = 1.0e20;
    const float b = 1.0;
    int n;
    for (n = 0; n <= 1024; ++ n) {
        float f = (float)n / 1024.0f;
        float p1 = lint1(a, b, f);
        float p2 = lint2(a, b, f);
        if (p1 != p2) {
            printf("%i %.6f %f %f %.6e\n", n, f, p1, p2, p2 - p1);
        }
    }
    return 0;
}

Вывод, немного скорректированный под форматирование:

    f            lint1               lint2             lint2-lint1
0.828125  17187500894208393216  17187499794696765440  -1.099512e+12
0.890625  10937500768952909824  10937499669441282048  -1.099512e+12
0.914062   8593750447104196608   8593749897348382720  -5.497558e+11
0.945312   5468750384476454912   5468749834720641024  -5.497558e+11
0.957031   4296875223552098304   4296874948674191360  -2.748779e+11
0.972656   2734375192238227456   2734374917360320512  -2.748779e+11
0.978516   2148437611776049152   2148437474337095680  -1.374390e+11
0.986328   1367187596119113728   1367187458680160256  -1.374390e+11
0.989258   1074218805888024576   1074218737168547840  -6.871948e+10
0.993164    683593798059556864    683593729340080128  -6.871948e+10
1.000000                     1                     0  -1.000000e+00

Если вы используете микроконтроллер без FPU, то операции с плавающей запятой будут очень дорогими. Легко может быть в двадцать раз медленнее для операции с плавающей запятой. Самое быстрое решение — просто выполнить всю математику, используя целые числа.

Количество мест после фиксированной двоичной точки (http://blog.credland.net/2013/09/binary-fixed-point-explanation.html?q=fixed+binary+point): XY_TABLE_FRAC_BITS.

Вот функция, которую я использую:

inline uint16_t unsignedInterpolate(uint16_t a, uint16_t b, uint16_t position) {
    uint32_t r1;
    uint16_t r2;

    /* 
     * Only one multiply, and one divide/shift right.  Shame about having to
     * cast to long int and back again.
     */

    r1 = (uint32_t) position * (b-a);
    r2 = (r1 >> XY_TABLE_FRAC_BITS) + a;
    return r2;    
}

Со встроенной функцией это должно быть прибл. 10-20 циклов.

Если у вас есть 32-битный микроконтроллер, вы сможете использовать большие целые числа и получать большие числа или большую точность без ущерба для производительности. Эта функция использовалась в 16-битной системе.

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

Стоит отметить, что стандартные формулы линейной интерполяции f1(t)=a+t(ba), f2(t)=b-(ba)(1-t) и f3(t)=a(1- t)+bt не гарантирует монотонности при использовании арифметики с плавающей запятой. В частности, если a != b, не гарантируется, что f1(1.0) == b или f2(0.0) == a, а для a == b не гарантируется, что f3(t) будет равно a , когда 0 ‹ t ‹ 1.

Эта функция работала для меня на процессорах, поддерживающих IEEE754 с плавающей запятой, когда мне нужно, чтобы результаты были монотонными (я использую ее с двойной точностью, но с плавающей запятой тоже должно работать):

double lerp(double a, double b, double t) 
{
    if (t <= 0.5)
        return a+(b-a)*t;
    else
        return b-(b-a)*(1.0-t);
}

Начиная с C++20 вы можете использовать std::lerp(), который, вероятно, будет наилучшая возможная реализация для вашей цели.

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

int lerp_int(int a, int b, float f)
{
    //float diff = (float)(b-a);
    //float frac = f*diff;
    //return a + (int)frac;
    return a + (int)(f * (float)(b-a));
}

Это умножает два броска и один поплавок. Если на вашей платформе приведение выполняется быстрее, чем сложение/вычитание с плавающей запятой, и если вам полезен целочисленный ответ, это может быть разумной альтернативой.