如何使用frexp实现双变量的模运算符？
如何使用frexp实现双变量的模运算符？
提问于 2020-07-08 00:29:26
我正在跟踪Kernighan&Pike“UNIX编程环境”。
本书中的一个练习(练习8-2，第241页)要求在double中为C变量实现模块化操作符(%)。
所以：
4.6 % 2.1 = 0.4
4.0 % 3.0 = 1.0因此，它基本上是使用dmod实现frexp。
dmod(4.6, 2.1) would return 0.4
dmod(4,0, 3.0) would return 1.0我看过这篇文章：implementing a modulus on fixed point type，它定义了一个实现这个操作符的算法。
但是这本书暗示我们要阅读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实现这个问题仍然是有效的。
回答 1
Stack Overflow用户
回答已采纳
发布于 2020-07-12 19:31:18
我不知道作者为浮点数的模假定了什么规格.我在这里假设它们是指标准C库函数fmod()的功能。
实现fmod()的最简单方法是使用二进制长和除法，它在一个循环中生成除法的商，该循环每次迭代生成一个商位。重复，直到商的所有整数位都用完为止，同时保留部分余数。在过程结束时，最后的余数表示期望的结果。
要开始长线除法，我们必须在开始时适当地将除数与股息对齐。这是通过扩展分红>=除数>股利/2来实现的。将frexp()与ldexp()结合使用提供了一个基于指数的粗略标度，而指数可能需要根据意义(mantissas)进行细化。
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;
}页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/62785780
复制相关文章
相似问题
腾讯云开发者
