math - 从一对浮子到单个浮子的近似可逆函数？

Question

前段时间，我在一些 CAD 代码中遇到了一对函数，用于将一组坐标(一对 float)编码为单个 float(用作散列键)，然后将单个 float 解压缩回原始对。

前向和后向函数只使用标准的数学运算——没有魔法摆弄 float 的位级表示，没有提取和交错单个数字或类似的东西。显然，这种反转在实践中并不完美，因为从两个 float 到一个 float 会损失相当大的精度，但根据函数的维基百科页面，它应该是完全可逆的，给定无限精度算法。

不幸的是，我不再处理该代码，而且我忘记了该函数的名称，所以我无法再次在维基百科上查找它。有人知道满足该描述的命名数学函数吗？

Answer 1

在评论中，OP 阐明了所需的函数应该将两个 IEEE-754 binary64 操作数映射到一个这样的操作数。实现这一点的一种方法是将每个 binary64( double )数字映射到 [0, 2²⁶-2] 中的正整数，然后使用一个很好的已知配对函数将两个这样的整数映射到区间 [0,2⁵²) 中的单个正整数，这在 binary64 中可以精确表示它有 52 个存储的有效数字(“尾数”)位。

根据 OP 的评论，binary64 操作数的范围不受限制，因此所有二进制数都应该以相同的相对精度表示，我们还需要处理零、无穷大和 NaN。为此，我选择了 log2() 来压缩数据。零被视为最小的 binary64 次正规，0x1.0p-1074，其结果是 0x1.0p-1074 和零都将解压为零。 log2 的结果落在 [-1074, 1024) 范围内。因为我们需要存储操作数的符号，所以我们将对数值偏置 1074，得到 [0, 2098)，然后将其缩放到几乎 [0, 2²⁵), 四舍五入到最接近的整数，最后加上原始binary64操作数的符号。几乎利用整个范围的动机是在范围的顶部留出一点空间用于无穷大和 NaN 的特殊编码(因此使用单一的规范 NaN 编码)。

由于文献中已知的配对函数仅对自然数起作用，因此我们需要从整数到自然数的映射。这很容易通过将负整数映射到奇自然数，而将正整数映射到偶自然数来轻松实现。因此，我们的操作数从 (-2²⁵, +2²⁵) 映射到 [0, 2²⁶-2]。然后配对函数将 [0, 2²⁶-2] 中的两个整数组合成 [0, 2⁵²] 中的单个整数。

文献中已知的不同配对函数的加扰行为不同，这可能会影响问题中提到的散列功能。它们的性能也可能不同。因此，我在下面的代码中为 pair()/unpair() 实现提供了四种不同的配对函数。相应的文献引用请见代码中的注释。

打包操作数的解包涉及以相反顺序应用每个打包步骤的逆。取消配对函数给了我们两个自然整数。这些被映射到两个整数，这两个整数被映射到两个对数值，然后用 exp2() 取幂以恢复原始数字，增加一些工作以获得特殊值和符号正确。

虽然对数以 10^-8 量级的相对精度表示，但最终结果的预期最大相对误差为 10^-5 量级，因为众所周知的取幂误差放大特性。在大量测试中，pack()/unpack() 往返观察到的最大相对误差为 2.167e-5。

下面是我的算法的 ISO C99 实现以及我的测试框架的一部分。这应该可以通过一些努力移植到其他编程语言。

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <math.h>
#include <float.h>

#define SZUDZIK_ELEGANT_PAIRING  (1)
#define ROZSA_PETER_PAIRING      (2)
#define ROSENBERG_STRONG_PAIRING (3)
#define CHRISTOPH_MICHEL_PAIRING (4)

#define PAIRING_FUNCTION  (ROZSA_PETER_PAIRING)

/*
  https://groups.google.com/forum/#!original/comp.lang.c/qFv18ql_WlU/IK8KGZZFJx4J
  From: geo <gmars...@gmail.com>
  Newsgroups: sci.math,comp.lang.c,comp.lang.fortran
  Subject: 64-bit KISS RNGs
  Date: Sat, 28 Feb 2009 04:30:48 -0800 (PST)

  This 64-bit KISS RNG has three components, each nearly
  good enough to serve alone.    The components are:
  Multiply-With-Carry (MWC), period (2^121+2^63-1)
  Xorshift (XSH), period 2^64-1
  Congruential (CNG), period 2^64
*/
static uint64_t kiss64_x = 1234567890987654321ULL;
static uint64_t kiss64_c = 123456123456123456ULL;
static uint64_t kiss64_y = 362436362436362436ULL;
static uint64_t kiss64_z = 1066149217761810ULL;
static uint64_t kiss64_t;
#define MWC64  (kiss64_t = (kiss64_x << 58) + kiss64_c, \
                kiss64_c = (kiss64_x >> 6), kiss64_x += kiss64_t, \
                kiss64_c += (kiss64_x < kiss64_t), kiss64_x)
#define XSH64  (kiss64_y ^= (kiss64_y << 13), kiss64_y ^= (kiss64_y >> 17), \
                kiss64_y ^= (kiss64_y << 43))
#define CNG64  (kiss64_z = 6906969069ULL * kiss64_z + 1234567ULL)
#define KISS64 (MWC64 + XSH64 + CNG64)

double uint64_as_double (uint64_t a)
{
    double r;
    memcpy (&r, &a, sizeof r);
    return r;
}

#define LOG2_BIAS    (1074.0)
#define CLAMP_LOW    (exp2(-LOG2_BIAS))
#define SCALE        (15993.5193125)
#define NAN_ENCODING (33554430)
#define INF_ENCODING (33554420)

/* check whether argument is an odd integer */
int is_odd_int (double a)
{
    return (-2.0 * floor (0.5 * a) + a) == 1.0;
}

/* compress double-precision number into an integer in (-2**25, +2**25) */
double compress (double x)
{
    double t;

    t = fabs (x);
    t = (t < CLAMP_LOW) ? CLAMP_LOW : t;
    t = rint ((log2 (t) + LOG2_BIAS) * SCALE);
    if (isnan (x)) t = NAN_ENCODING;
    if (isinf (x)) t = INF_ENCODING;
    return copysign (t, x);
}

/* expand integer in (-2**25, +2**25) into double-precision number */
double decompress (double x)
{
    double s, t;

    s = fabs (x);
    t = s / SCALE;
    if (s == (NAN_ENCODING)) t = NAN;
    if (s == (INF_ENCODING)) t = INFINITY;
    return copysign ((t == 0) ? 0 : exp2 (t - LOG2_BIAS), x);
}

/* map whole numbers to natural numbers. Here: (-2^25, +2^25) to [0, 2^26-2] */
double map_Z_to_N (double x)
{
    return (x < 0) ? (-2 * x - 1) : (2 * x);
}

/* Map natural numbers to whole numbers. Here: [0, 2^26-2] to (-2^25, +2^25) */
double map_N_to_Z (double x)
{
    return is_odd_int (x) ? (-0.5 * (x + 1)) : (0.5 * x);
}

#if PAIRING_FUNCTION == SZUDZIK_ELEGANT_PAIRING
/* Matthew Szudzik, "An elegant pairing function." In Wolfram Research (ed.) 
   Special NKS 2006 Wolfram Science Conference, pp. 1-12. 

   Here: map two natural numbers in [0, 2^26-2] to natural number in [0, 2^52),
   and vice versa
*/
double pair (double x, double y)
{
    return (x != fmax (x, y)) ? (y * y + x) : (x * x + x + y);
}

void unpair (double z, double *x, double *y)
{
    double sqrt_z = trunc (sqrt (z));
    double sqrt_z_diff = z - sqrt_z * sqrt_z;
    *x = (sqrt_z_diff < sqrt_z) ? sqrt_z_diff : sqrt_z;
    *y = (sqrt_z_diff < sqrt_z) ? sqrt_z : (sqrt_z_diff - sqrt_z);
}

#elif PAIRING_FUNCTION == ROZSA_PETER_PAIRING
/*
  Rozsa Peter, "Rekursive Funktionen" (1951), p. 44. Via David R. Hagen's blog,
  https://drhagen.com/blog/superior-pairing-function/

  Here: map two natural numbers in [0, 2^26-2] to natural number in [0, 2^52),
  and vice versa
*/
double pair (double x, double y)
{
    double mx = fmax (x, y);
    double mn = fmin (x, y);
    double sel = (mx == x) ? 0 : 1;
    return mx * mx + mn * 2 + sel;
}

void unpair (double z, double *x, double *y)
{
    double sqrt_z = trunc (sqrt (z));
    double sqrt_z_diff = z - sqrt_z * sqrt_z;
    double half_diff = trunc (sqrt_z_diff * 0.5);
    *x = is_odd_int (sqrt_z_diff) ? half_diff : sqrt_z;
    *y = is_odd_int (sqrt_z_diff) ? sqrt_z : half_diff;
}

#elif PAIRING_FUNCTION == ROSENBERG_STRONG_PAIRING
/*
  A. L. Rosenberg and H. R. Strong, "Addressing arrays by shells", 
  IBM Technical Disclosure Bulletin, Vol. 14, No. 10, March 1972,
  pp. 3026-3028.

  Arnold L. Rosenberg, "Allocating storage for extendible arrays," 
  Journal of the ACM, Vol. 21, No. 4, October 1974, pp. 652-670.
  Corrigendum, Journal of the ACM, Vol. 22, No. 2, April 1975, p. 308.

  Matthew P. Szudzik, "The Rosenberg-Strong Pairing Function", 2019
  https://arxiv.org/abs/1706.04129

  Here: map two natural numbers in [0, 2^26-2] to natural number in [0, 2^52),
  and vice versa
 */
double pair (double x, double y)
{
    double mx = fmax (x, y);
    return mx * mx + mx + x - y;
}

void unpair (double z, double *x, double *y)
{
    double sqrt_z = trunc (sqrt (z));
    double sqrt_z_diff = z - sqrt_z * sqrt_z;
    *x = (sqrt_z_diff < sqrt_z) ? sqrt_z_diff : sqrt_z;
    *y = (sqrt_z_diff < sqrt_z) ? sqrt_z : (2 * sqrt_z - sqrt_z_diff);
}
#elif PAIRING_FUNCTION == CHRISTOPH_MICHEL_PAIRING
/*
  Christoph Michel, "Enumerating a Grid in Spiral Order", September 7, 2016,
  https://cmichel.io/enumerating-grid-in-spiral-order. Via German Wikipedia,
  https://de.wikipedia.org/wiki/Cantorsche_Paarungsfunktion

  Here: map two natural numbers in [0, 2^26-2] to natural number in [0, 2^52),
  and vice versa
*/
double pair (double x, double y)
{
    double mx = fmax (x, y);
    return mx * mx + mx + (is_odd_int (mx) ? (x - y) : (y - x));
}

void unpair (double z, double *x, double *y)
{
    double sqrt_z = trunc (sqrt (z));
    double sqrt_z_diff = z - sqrt_z * (sqrt_z + 1);
    double min_clamp = fmin (sqrt_z_diff, 0);
    double max_clamp = fmax (sqrt_z_diff, 0);
    *x = is_odd_int (sqrt_z) ? (sqrt_z + min_clamp) : (sqrt_z - max_clamp);
    *y = is_odd_int (sqrt_z) ? (sqrt_z - max_clamp) : (sqrt_z + min_clamp);
}
#else
#error unknown PAIRING_FUNCTION
#endif

/* Lossy pairing function for double precision numbers. The maximum round-trip
   relative error is about 2.167e-5
*/
double pack (double a, double b)
{
    double c, p, q, s, t;

    p = compress (a);
    q = compress (b);
    s = map_Z_to_N (p);
    t = map_Z_to_N (q);
    c = pair (s, t);
    return c;
}

/* Unpairing function for double precision numbers. The maximum round-trip
   relative error is about 2.167e-5 */
void unpack (double c, double *a, double *b)
{
    double s, t, u, v;

    unpair (c, &s, &t);
    u = map_N_to_Z (s);
    v = map_N_to_Z (t);
    *a = decompress (u);
    *b = decompress (v);
}

int main (void)
{
    double a, b, c, ua, ub, relerr_a, relerr_b;
    double max_relerr_a = 0, max_relerr_b = 0; 

#if PAIRING_FUNCTION == SZUDZIK_ELEGANT_PAIRING
    printf ("Testing with Szudzik's elegant pairing function\n");
#elif PAIRING_FUNCTION == ROZSA_PETER_PAIRING
    printf ("Testing with Rozsa Peter's pairing function\n");
#elif PAIRING_FUNCTION == ROSENBERG_STRONG_PAIRING
    printf ("Testing with Rosenberg-Strong pairing function\n");
#elif PAIRING_FUNCTION == CHRISTOPH_MICHEL_PAIRING
    printf ("Testing with C. Michel's spiral pairing function\n");
#else
#error unkown PAIRING_FUNCTION
#endif

    do {
        a = uint64_as_double (KISS64);
        b = uint64_as_double (KISS64);

        c = pack (a, b);
        unpack (c, &ua, &ub);

        if (!isnan(ua) && !isinf(ua) && (ua != 0)) {
            relerr_a = fabs ((ua - a) / a);
            if (relerr_a > max_relerr_a) {
                printf ("relerr_a= %15.8e  a=% 23.16e  ua=% 23.16e\n", 
                        relerr_a, a, ua);
                max_relerr_a = relerr_a;
            }
        }
        if (!isnan(ub) && !isinf(ub) && (ub != 0)) {
            relerr_b = fabs ((ub - b) / b);
            if (relerr_b > max_relerr_b) {
                printf ("relerr_b= %15.8e  b=% 23.16e  ub=% 23.16e\n", 
                        relerr_b, b, ub);
                max_relerr_b = relerr_b;
            }
        }
    } while (1);
    return EXIT_SUCCESS;
}