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

Question

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

书中的一个练习（练习 8-2，第 241 页）要求为 C.

中的 double 变量实现模运算符 (%)

所以：

4.6 % 2.1 = 0.4
4.0 % 3.0 = 1.0

因此它基本上是使用 frexp:

实现 dmod

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

我看过这个 post: 它定义了一个算法来实现这个运算符。

但是这本书建议阅读 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 实现它仍然有效。

Answer 1

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

的功能

实现 fmod() 的最简单方法是使用二进制普通除法，它在每次迭代产生一个商位的循环中产生除法的商。重复直到用尽商的所有整数位，同时保留部分余数。在这个过程结束时，最后的余数就是想要的结果。

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

fmod() 的示例性 ISO-C99 实现如下所示。 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 <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 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;
}

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

How do I implement a modulus operator for double variables using frexp?

c

math

floating-point

c89

modulus