如何使用frexp为双变量实现模运算符? [英] How do I implement a modulus operator for double variables using frexp?
问题描述
我关注的是 Kernighan& Pike的"UNIX编程环境" .
这本书的一个练习(练习8-2，第241页)要求为C中的double变量实现模运算符(％).
An exercise from the book (Exercise 8-2, page 241) asks implementing the modulo operator (%) for double variables in C.
所以:
4.6 % 2.1 = 0.4
4.0 % 3.0 = 1.0
因此，它基本上是使用frexp来实现dmod的:
Therefore it is basically implementing dmod using frexp:
dmod(4.6, 2.1) would return 0.4
dmod(4,0, 3.0) would return 1.0
我看过这篇文章:在定点类型上实现模量定义了实现此运算符的算法.
I have seen this post: implementing a modulus on fixed point type which defines an algorithm to implement this operator.
但是这本书建议读者阅读frexp(3)，因此我想可以使用该功能来做到这一点.
But the book suggests as a hint to read frexp(3), so I guess it is possible to do it using that function.
现在，如果我正确理解了手册页，则该功能将执行类似(伪代码)的操作:
Now if I understood the man page correctly, that function does things like (pseudocode):
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
但是我仍然不知道如何混合这些值以获得模运算符.
But I still can not figure out how to mix those values to get the modulus operator.
这本书很旧，也许有更好的方法可以使用，但是这个问题更具学术性，对于理解如何使用frexp实现它仍然有效.
The book is old and there are probably better ways of doing it, but the question is more academic and still valid to understand how to implement this with frexp.
推荐答案
我不知道作者对浮点数的模数采用什么规范.我在这里假设他们指的是标准C库函数fmod()的功能.
I do not know what specification the authors assumed for modulo of floating-point numbers. I assume here that they are referring to the functionality of the standard C library functionfmod().
实现fmod()的最简单方法是使用二进制长除法，该除法在循环中产生除法的商，该循环在每次迭代中产生一个商位.重复直到所有商的整数位都用完，同时保留部分余数.在过程结束时，最后的余数代表了预期的结果.
The simplest way of implementing fmod() is to use binary longhand division which produces the quotient of the division in a loop that produces one quotient bit per iteration. Repeat until all integer bits of the quotient have been exhausted, while retaining the partial remainder. At the end of the process, the final remainder represents desired result.
要开始进行长期除法，我们必须在开始时将除数与除数正确对齐.这可以通过缩放使得分红> =除数>分红/2. frexp()与ldexp()结合使用可根据指数进行粗略缩放，可能需要根据有效位数(尾数)进行精炼.
To start the longhand division, we have to properly align the divisor with the dividend at the start. This is achieved by scaling such that dividend >= divisor > dividend/2. The use of frexp() in conjunction with ldexp() provides a rough scaling based on the exponents which may have to be refined based on the significands (mantissas).
fmod()的示例性ISO-C99实现如下所示. remainder()的实现看起来很相似，但由于要求将商四舍五入到最接近或偶数，而不是将其截断而变得更复杂.
An exemplary ISO-C99 implementation of fmod()is shown below. Implementation of remainder() would look similar but a bit more complicated due to the requirement to round the quotient to nearest-or-even, rather than truncating it.
#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为双变量实现模运算符?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
