整数の自乗アルゴリズム（sqr32）

2025-12-03

Takayuki HOSODA

概要

これは、16×16 乗算器しか持たないプロセッサ向けに設計した 32-bit 整数の自乗（square）アルゴリズムと C での実装例です。

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

アルゴリズム

a² = a_h² B² + 2 a_h a_l B + a_l²
B = 2¹⁶

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

C 実装

sqr32 -- 与えられた 32-bit 整数の自乗を 64-bit 符号なし整数として返します。

/* 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 するものとして扱っています。
※ 64-bit 乗算器がある環境では、通常の a * a を使用してください。