🌐English

概要

これは、16×16 乗算器しか持たない プロセッサ向けに設計した 整数の自乗アルゴリズムと、その C による実装例 (sqr32) です。

自乗の対称性を利用し、Karatsuba 法 的な分解によって、 1 回の 32×32-bit 乗算を 3 回の 16×16-bit 乗算に置き換えています。

また、付録として同じ分解に基づく 18-bit 自乗器（sqr18u）の Verilog 実装例も示します。

※ 64-bit 乗算器がある環境では、通常の a * a を使用してください。
FPGA においても 18×18 の DSP 乗算器が利用可能な場合は、単純な乗算の方が一般に効率的です。 ここで示した分解は、DSP ブロックが使えない場合や節約したい場合に有効です。

アルゴリズム

ここで a_H, a_L は入力値を 上位・下位 16-bit に分割したものです。

C 実装

/* sqr32 - calculate squre of given integer value
 * Rev.1.0 (2025-12-22) (c) 2009-2025 Takayuki HOSODA
 * SPDX-License-Identifier: BSD-3-Clause
 */
#include <inttypes.h>
uint64_t sqr32(int32_t a);
#ifdef HAVE_NO_MUL64
uint64_t sqr32(int32_t a) {
    uint32_t        ax, ah, al;
    uint64_t        bh, bm, bl;
    // Karatsuba method: a * a = a_h^2 * B^2 + 2 * a_h * a_l * B + a_l^2
    ax = (uint32_t)(a < 0 ? -a : a);       // assume two's complement; INT32_MIN wraps
    al = (ax      ) & 0xffff;              // lower 16-bit of input
    ah = (ax >> 16) & 0xffff;              // upper 16-bit of input
    bh = (uint64_t)(ah * ah) << (2 * 16);  // 16x16 multiplier and 64-bit shifter
    bm = (uint64_t)(ah * al) << (1 + 16);  // 16x16 multiplier and 64-bit shifter
    bl = (uint64_t)(al * al);              // 16x16 multiplier
    return bl + bm + bh;                   // 64-bit adder
}
#else
uint64_t sqr32(int32_t x) {
    int64_t eax = (int64_t)x;
    return (uint64_t)(eax * eax);
}
#endif

※ two's complement を前提とし、INT32_MIN は wrap するものとして扱っています。

追補 — Verilog での実装

この実装は 18×18 乗算器を回避し、9×9 乗算へ分解することで、LUT ベースのロジックに効率よくマップされる構成としています。

sqr18u — 与えられた 18-bit の符号なし整数の自乗を 36-bit 符号なし整数として返します。

/* sqr18u -- unsigned 18-bit squarer without hardware 18x18 multiplier
 * Rev.1.0 (2025-12-22) (c) 2009-2025 Takayuki HOSODA
 * SPDX-License-Identifier: BSD-3-Clause
 */
module sqr18u (input  wire [17:0] a, output wire [35:0] y);
    // split input into high / low 9-bit parts
    wire [8:0] ah = a[17:9];
    wire [8:0] al = a[ 8:0];
    // 3x (9x9, LUT-based) multiplications instead of 18x18 multiplication. 
    wire [17:0] p_hh = ah * ah;
    wire [17:0] p_ll = al * al;
    wire [17:0] p_hl = ah * al;
    // Square-specialized Karatsuba-like decomposition:
    // a^2 = al^2 + 2*ah*al*2^9 + ah^2*2^18
    assign y = (p_ll)
             + (p_hl << (1 + 9))
             + (p_hh << (9 + 9)); // 36-bit adder
endmodule

SEE ALSO

• Karatsuba 法による乗算