如何使用 frexp 实现双变量的模运算符？

我正在关注Kernighan&Pike“UNIX 编程环境”.

书中的一个练习（练习 8-2，第 241 页）要求实现模运算符 (%)double变量在C.

So:

4.6 % 2.1 = 0.4
4.0 % 3.0 = 1.0

因此基本上是在实施dmod using frexp:

dmod(4.6, 2.1) would return 0.4
dmod(4,0, 3.0) would return 1.0

我看过这个帖子：在定点类型上实现模数它定义了一个算法来实现这个运算符。

但这本书建议作为阅读的提示frexp(3)，所以我想可以使用该函数来做到这一点。

现在，如果我正确理解了手册页，该函数会执行类似（伪代码）的操作：

a,b -- double variables
a_exp,b_exp -- integer exponents for frexp
a_x = frexp(a,&a_exp) --> a = a_x * 2^a_exp
b_x = frexp(b,&b_exp) --> b = b_x * 2^b_exp
c=a/b
c_exp -- integer exponent for frexp
c_x = frexp(c,&c_exp) --> c = c_x * 2^c_exp

但我仍然不知道如何混合这些值来获得模运算符。

这本书很旧，可能有更好的方法来做到这一点，但问题更具学术性，并且对于理解如何实现这一点仍然有效frexp.

我不知道作者对浮点数模假设的规范是什么。我在这里假设它们指的是标准 C 库函数的功能fmod().

最简单的实现方式fmod()是使用二进制普通除法，它在循环中产生除法的商，每次迭代产生一位商位。重复此操作，直到商的所有整数位都用完，同时保留部分余数。在该过程结束时，最终的余数代表期望的结果。

要开始普通除法，我们必须在开始时将除数与被除数正确对齐。这是通过缩放使得被除数 >= 除数 > 被除数/2 来实现的。指某东西的用途frexp()和这个结合ldexp()提供基于指数的粗略缩放，可能需要根据有效数（尾数）进行细化。

示例性 ISO-C99 实施fmod()如下所示。实施remainder()看起来类似，但有点复杂，因为需要将商舍入到最接近的或偶数，而不是截断它。

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

/* returns the floating-point remainder of a/b (rounded towards zero) */
double my_fmod (double a, double b)
{
    const double NAN_INDEFINITE = 0.0 / 0.0;
    double r;
    if (isnan (a) || isnan (b)) {
        r = a + b;
    } else if (isinf (a) || (b == 0.0)) {
        r = NAN_INDEFINITE;
    } else {
        double fa, fb, dividend, divisor;
        int expo_a, expo_b;
        fa = fabs (a);
        fb = fabs (b);
        if (fa >= fb) {
            dividend = fa;
            /* normalize divisor */
            (void)frexp (fa, &expo_a);
            (void)frexp (fb, &expo_b);
            divisor = ldexp (fb, expo_a - expo_b);
            if (divisor <= 0.5 * dividend) {
                divisor += divisor;
            }
            /* compute quotient one bit at a time */
            while (divisor >= fb) {
                if (dividend >= divisor) {
                    dividend -= divisor;
                }
                divisor *= 0.5;
            }
            /* dividend now represents remainder */
            r = copysign (dividend, a);
        } else {
            r = a;
        }
    }
    return r;
}

/*
  From: geo <[email protected]>
  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 int64_as_double (int64_t a)
{
    double r;
    memcpy (&r, &a, sizeof r);
    return r;
}

int32_t double_as_int64 (double a)
{
    int64_t r;
    memcpy (&r, &a, sizeof r);
    return r;
}

int main (void)
{
    double a, b, res, ref;
    uint64_t i = 0;
    do {
        a = int64_as_double (KISS64);
        b = int64_as_double (KISS64);
        ref = fmod (a, b);
        res = my_fmod (a, b);
        if (double_as_int64 (res) != double_as_int64 (ref)) {
            printf ("error: a=% 23.16e b=% 23.16e res=% 23.16e ref=% 23.16e\n", a, b, res, ref);
            return EXIT_FAILURE;
        }
        i++;
        if (!(i & 0xfffff)) printf ("\r%llu", i);
    } while (i);
    return EXIT_SUCCESS;
}