FMA演算

FMA演算（Fused Multiply-Add、融合積和演算)とは、乗算と加算を1命令で行う演算です。

• 乗算と加算を2命令で実行するより速い
• 乗算と加算を2命令で実行する場合より丸め誤差が小さい

数式で表すと x * y + z となります。この数式は「行列の乗算」や「ベクトルの内積」でよく使われます。行列とベクトルといえばSIMD演算の効果が絶大なので昨今のCPUのSIMD命令やGPUでは「FMA演算をSIMD実行できる機能（1命令で複数データにFMA演算を行う機能）」が定番機能となっています。

このFMA演算をCPUやGPUにおいて1命令で行うことで途中の積算を丸めずに積和演算を最終演算結果の誤差を小さくする工夫が考案されており、IEEE 754規格の2008年改訂版で標準化されています。これに準拠したCPUやGPUでの挙動は常に一定となります。

C言語

C言語だと以下のような感じです。標準的な環境ではmath.hにfma関数があります。

#include <stdio.h>
#include <math.h>
int main() {
    double a = 2.0;
    double b = 3.0;
    double c = 4.0;
    double result;

    // FMAを使用して a * b + c を計算
    result = fma(a, b, c);

    printf("Result of FMA: %f\n", result);
    return 0;
}

ベクトルの内積

ベクトルの内積は「FMAの繰り返し」で計算できます。 ベクトルの内積は3DCGにおけるライティングやテクスチャマッピングのUV座標変換で多用されます。

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

#define SIZE 3

int main() {
    double vec1[SIZE] = {1.0, 2.0, 3.0};
    double vec2[SIZE] = {4.0, 5.0, 6.0};
    double dot_product = 0.0;

    for (int i = 0; i < SIZE; i++) {
        dot_product = fma(vec1[i], vec2[i], dot_product);
    }

    printf("Dot product: %f\n", dot_product);

    return 0;
}

行列の乗算

行列の乗算は各要素の積とその総和を計算するためFMAの特性が非常に役立ちます。 行列の乗算は3DCGにおいては主に各種座標変換（ワールド座標変換やビュー座標変換など）で使われます。 また昨今の機械学習においては演算のほとんどが行列の乗算という状況です。

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

#define SIZE 3

void matrix_multiply(double mat1[SIZE][SIZE], double mat2[SIZE][SIZE], double result[SIZE][SIZE]) {
    for (int i = 0; i < SIZE; i++) {
        for (int j = 0; j < SIZE; j++) {
            result[i][j] = 0.0;
            for (int k = 0; k < SIZE; k++) {
                result[i][j] = fma(mat1[i][k], mat2[k][j], result[i][j]);
            }
        }
    }
}

int main() {
    double mat1[SIZE][SIZE] = {
        {1.0, 2.0, 3.0},
        {4.0, 5.0, 6.0},
        {7.0, 8.0, 9.0}
    };
    double mat2[SIZE][SIZE] = {
        {9.0, 8.0, 7.0},
        {6.0, 5.0, 4.0},
        {3.0, 2.0, 1.0}
    };
    double result[SIZE][SIZE];

    matrix_multiply(mat1, mat2, result);

    printf("Result matrix:\n");
    for (int i = 0; i < SIZE; i++) {
        for (int j = 0; j < SIZE; j++) {
            printf("%f ", result[i][j]);
        }
        printf("\n");
    }

    return 0;
}