C++で学ぶオイラーの四平方恒等式――2数の積を4つの平方数の和で表す方法

問題の概要

本記事では、2つの数値が与えられたとき、オイラーの四平方恒等式（Euler's Four-Square Identity）を利用して両者の積を求める方法を解説します。

オイラーの四平方恒等式とは、「2つの整数がそれぞれ4つの平方数の和として表せるならば、その積もまた4つの平方数の和として表せる」という有名な定理です。つまり、次の関係が成り立ちます。

a = x1² + x2² + x3² + x4²
b = y1² + y2² + y3² + y4²
a × b = z1² + z2² + z3² + z4²

具体例で理解しよう

入力：

a = 54 = 2×2 + 3×3 + 4×4 + 5×5
b = 10 = 1×1 + 2×2 + 1×1 + 2×2

出力： 1×1 + 1×1 + 3×3 + 23×23

解説：

a と b の積は 540 となり、これは実にさまざまな形で表現できます。ここでは一例として、次の表し方を挙げます。

540 = 1×1 + 1×1 + 3×3 + 23×23 = 1 + 1 + 9 + 529

解法のアプローチ

最も単純な解法は、総当たり（トライアル）方式で4つの平方数の組み合わせをすべて試すことです。具体的には、平方数の候補ごとにネストしたループ（合計4重ループ）を回し、計算結果が積と一致した組み合わせをすべて出力します。

ただしこの方法は計算量が膨大になります。ループ変数の範囲を考慮しない素朴な見積もりでは O((a*b)⁴) オーダー、コードのように √(prod) まで探索する場合でも O((a*b)²) 程度となり、非効率です。

解法1：4重ループによる実装例

サンプルコード

#include <bits/stdc++.h>
using namespace std;

void findEulerSquareNumberValue(int a, int b) {
    int prod = a * b;
    int sumSquare = 0;
    for (int x = 0; x <= sqrt(prod); x++) {
        for (int y = x; y <= sqrt(prod); y++) {
            for (int z = y; z <= sqrt(prod); z++) {
                for (int k = z; k <= sqrt(prod); k++) {
                    sumSquare = (x*x) + (y*y) + (z*z) + (k*k);
                    if (sumSquare == prod) {
                        cout<<"The product "<<a<<"*"<<b<<" = "<<prod<<" represented as ";
                        cout<<x<<"*"<<x<<" + ";
                        cout<<y<<"*"<<y<<" + ";
                        cout<<z<<"*"<<z<<" + ";
                        cout<<k<<"*"<<k<<endl;
                        cout<<endl;
                    }
                }
            }
        }
    }
}

int main() {
    int a = (2*2) + (3*3) + (4*4) + (5*5);
    int b = (1*1) + (2*2) + (1*1) + (2*2);
    cout<<"a = (2*2) + (3*3) + (4*4) + (5*5) = "<<a<<endl;
    cout<<"b = (1*1) + (2*2) + (1*1) + (2*2) = "<<b<<endl;
    findEulerSquareNumberValue(a, b);
    return 0;
}

実行結果

a = (2*2) + (3*3) + (4*4) + (5*5) = 54
b = (1*1) + (2*2) + (1*1) + (2*2) = 10
The product 54*10 = 540 represented as 1*1 + 1*1 + 3*3 + 23*23
The product 54*10 = 540 represented as 1*1 + 3*3 + 13*13 + 19*19
The product 54*10 = 540 represented as 1*1 + 5*5 + 15*15 + 17*17
The product 54*10 = 540 represented as 1*1 + 7*7 + 7*7 + 21*21
The product 54*10 = 540 represented as 1*1 + 9*9 + 13*13 + 17*17
The product 54*10 = 540 represented as 2*2 + 4*4 + 6*6 + 22*22
The product 54*10 = 540 represented as 2*2 + 4*4 + 14*14 + 18*18
The product 54*10 = 540 represented as 2*2 + 6*6 + 10*10 + 20*20
The product 54*10 = 540 represented as 2*2 + 12*12 + 14*14 + 14*14
The product 54*10 = 540 represented as 3*3 + 3*3 + 9*9 + 21*21
The product 54*10 = 540 represented as 3*3 + 7*7 + 11*11 + 19*19
The product 54*10 = 540 represented as 3*3 + 9*9 + 15*15 + 15*15
The product 54*10 = 540 represented as 3*3 + 11*11 + 11*11 + 17*17
The product 54*10 = 540 represented as 4*4 + 10*10 + 10*10 + 18*18
The product 54*10 = 540 represented as 5*5 + 5*5 + 7*7 + 21*21
The product 54*10 = 540 represented as 5*5 + 11*11 + 13*13 + 15*15
The product 54*10 = 540 represented as 6*6 + 6*6 + 12*12 + 18*18
The product 54*10 = 540 represented as 7*7 + 7*7 + 9*9 + 19*19
The product 54*10 = 540 represented as 7*7 + 9*9 + 11*11 + 17*17
The product 54*10 = 540 represented as 9*9 + 11*11 + 13*13 + 13*13
The product 54*10 = 540 represented as 10*10 + 10*10 + 12*12 + 14*14

解法2：3重ループへの最適化

時間計算量を抑えるもう一つの方法が、3重のネストループに削減する手法です。3つの平方数を固定した時点で、残りの値が完全平方数かどうかを判定します。完全平方数であればその平方根が4つ目の値となり、そうでなければその組み合わせには解が存在しないと判断できます。これにより第4のループが不要になり、計算量は O((a*b)^3/2) まで削減され、解法の効率が大幅に向上します。

サンプルコード

#include <bits/stdc++.h>
using namespace std;

void findEulerSquareNumberValue(int a, int b) {
    int prod = a * b;
    int sumSquare = 0;
    for (int x = 0; x <= sqrt(prod); x++) {
        for (int y = x; y <= sqrt(prod); y++) {
            for (int z = y; z <= sqrt(prod); z++) {
                sumSquare = (x*x) + (y*y) + (z*z);
                float k = sqrt(prod - sumSquare);
                if (floor(k) == ceil(k)) {
                    cout<<"The product "<<a<<"*"<<b<<" = "<<prod<<" represented as ";
                    cout<<x<<"*"<<x<<" + ";
                    cout<<y<<"*"<<y<<" + ";
                    cout<<z<<"*"<<z<<" + ";
                    cout<<k<<"*"<<k<<endl;
                    cout<<endl;
                }
            }
        }
    }
}

int main() {
    int a = (2*2) + (3*3) + (4*4) + (5*5);
    int b = (1*1) + (2*2) + (1*1) + (2*2);
    cout<<"a = (2*2) + (3*3) + (4*4) + (5*5) = "<<a<<endl;
    cout<<"b = (1*1) + (2*2) + (1*1) + (2*2) = "<<b<<endl;
    findEulerSquareNumberValue(a, b);
    return 0;
}

実行結果（抜粋）

a = (2*2) + (3*3) + (4*4) + (5*5) = 54
b = (1*1) + (2*2) + (1*1) + (2*2) = 10
The product 54*10 = 540 represented as 1*1 + 1*1 + 3*3 + 23*23
The product 54*10 = 540 represented as 1*1 + 1*1 + 23*23 + 3*3
The product 54*10 = 540 represented as 1*1 + 3*3 + 13*13 + 19*19
The product 54*10 = 540 represented as 1*1 + 3*3 + 19*19 + 13*13
The product 54*10 = 540 represented as 1*1 + 5*5 + 15*15 + 17*17
The product 54*10 = 540 represented as 1*1 + 5*5 + 17*17 + 15*15
...
（以下、要素の並び順が異なる組み合わせも含めて多数の出力が続きます）

まとめ

オイラーの四平方恒等式を用いれば、それぞれ4つの平方数の和で表せる2つの整数の積も、必ず4つの平方数の和として表せることが保証されます。C++では4重ループの総当たりでも解けますが、残り1つの値を「完全平方数かどうかの判定」で求めることで、ループを1段減らして計算量を大幅に削減できるのがポイントです。数値が大きくなるほど、この最適化の効果は顕著に現れます。