#include <stdio.h>
#include <stdint.h>
#include <math.h>
#include <assert.h>
typedef int64_t rtnum_t;
typedef int64_t rtden_t;
static bool get_rational(double x, rtnum_t& num, rtden_t& den)
{
// 出力パラメータ初期化
num = (rtnum_t)0;
den = (rtden_t)0;
// 負号保存
const bool bNeg = (x < 0.0);
x = fabs(x);
// 指数部取得
int exp;
(void)frexp(x, &exp);
const int IEEE754_exp = exp - 1;
// 誤差の指数取得
// これは仮数部最下位bitが2の何乗かを表す。
// 倍精度浮動小数点数の仮数部(小数部のみ表現)のビット数は52なので、52を引く。
int err_exp = IEEE754_exp - 52;
if (err_exp < -1022) {
// 倍精度浮動小数点数で表現可能な指数の下限が精度の下限。
// (2^(-1022)未満は非正規化数)
err_exp = -1022;
}
// 誤差取得
// これは仮数部最下位bitが表す数にあたる。取得にあたり、
// pow()は中で多項式近似しているかもしれないので信用しないことにする。
double err = 1.0;
if (err_exp >= 0) {
for (int i = 0; i < err_exp; i++) {
err *= 2.0; // 2のべきなので計算誤差0。
}
}
else {
for (int i = 0; i < -err_exp; i++) {
err /= 2.0; // 2のべきなので計算誤差0。
}
}
// (とは言いつつ、実は一致するらしい。)
assert(err == pow(2.0, err_exp));
// 誤差が1の位より大きいケースには対応しない。
if (err > 1.0) {
assert(0);
return false;
}
// 分子取得
// max(0, x)∈[ 10^exp10, 10^(exp10 + 1) )となるexp10を求める。
double base10 = 1.0;
int exp10 = 0;
while (base10 <= x) {
base10 *= 10.0; // 整数倍なので計算誤差0。
exp10++;
}
while (base10 > x && exp10 > 0) {
assert(base10 > 0.0 && fmod(base10, 10.0) == 0.0);
base10 /= 10.0; // 10の倍数を10で割るので計算誤差0。
exp10--;
}
// 整数部出力
// 上の桁から求めていく。
while (base10 > err && exp10 > 0) {
const double tmp = x / base10;
const double intpart = floor(tmp);
const rtnum_t digit = (rtnum_t)intpart % 10;
num *= 10;
num += digit;
assert(base10 > 0.0 && fmod(base10, 10.0) == 0.0);
base10 /= 10.0; // 10の倍数を10で割るので計算誤差0。
exp10--;
}
// 1の桁
{
const double tmp = x / base10;
const double intpart = floor(tmp);
const rtnum_t digit = (rtnum_t)intpart % 10;
num *= 10;
num += digit;
exp10--;
}
// 小数部出力
// 上の桁から求めていく。
double x_frac = x - floor(x);
err *= 10.0; // 整数倍なので計算誤差0。
while (err < 1.0) {
x_frac *= 10.0; // 整数倍なので計算誤差0。
const double intpart = floor(x_frac);
const rtnum_t digit = (rtnum_t)intpart % 10;
num *= 10;
num += digit;
exp10--;
err *= 10.0; // 整数倍なので計算誤差0。
}
// この時点でexp10は10進数表記したときの誤差の桁を指しており、
// 誤差±10^(exp10)の10進数表記が得られたことになる。
// 余分な0または9を削除
exp10++;
if (num % 10 == 9) {
// 余分な9を削除。
// 2進数のIEEE 754形式xから10進数表現numに直したとき
// 運悪く9の循環小数になったときここに来る。
// 削除といっても+1するだけで、後続の「0」削除処理が仕事をしてくれる。
num++;
}
// 余分な0を削除
while (num % 10 == 0) {
num /= 10;
exp10++;
}
// exp10 > 0のケースには対応しない。
if (exp10 > 0) {
assert(0);
return false;
}
// 負号復元
if (bNeg) {
num = -num;
}
// 分母確定
den = (rtden_t)1;
for (int i = 0; i < -exp10; i++) {
den *= 10;
}
return true;
}
void RationalTest()
{
// 有理数化テスト
{
const char* x_str = "0.1";
const double x = 0.1;
rtnum_t num;
rtden_t den;
get_rational(x, num, den);
printf("x=%s: rt=%lld/%lld\n", x_str, num, den);
}
{
const char* x_str = "0.11";
const double x = 0.11;
rtnum_t num;
rtden_t den;
get_rational(x, num, den);
printf("x=%s: rt=%lld/%lld\n", x_str, num, den);
}
{
const char* x_str = "0.111";
const double x = 0.111;
rtnum_t num;
rtden_t den;
get_rational(x, num, den);
printf("x=%s: rt=%lld/%lld\n", x_str, num, den);
}
{
const char* x_str = "0.1111";
const double x = 0.1111;
rtnum_t num;
rtden_t den;
get_rational(x, num, den);
printf("x=%s: rt=%lld/%lld\n", x_str, num, den);
}
{
const char* x_str = "0.9";
const double x = 0.9;
rtnum_t num;
rtden_t den;
get_rational(x, num, den);
printf("x=%s: rt=%lld/%lld\n", x_str, num, den);
}
{
const char* x_str = "0.99";
const double x = 0.99;
rtnum_t num;
rtden_t den;
get_rational(x, num, den);
printf("x=%s: rt=%lld/%lld\n", x_str, num, den);
}
{
const char* x_str = "0.999";
const double x = 0.999;
rtnum_t num;
rtden_t den;
get_rational(x, num, den);
printf("x=%s: rt=%lld/%lld\n", x_str, num, den);
}
{
const char* x_str = "0.9999";
const double x = 0.9999;
rtnum_t num;
rtden_t den;
get_rational(x, num, den);
printf("x=%s: rt=%lld/%lld\n", x_str, num, den);
}
{
const char* x_str = "0.33333";
const double x = 0.33333;
rtnum_t num;
rtden_t den;
get_rational(x, num, den);
printf("x=%s: rt=%lld/%lld\n", x_str, num, den);
}
{
const char* x_str = "123.34567";
const double x = 123.34567;
rtnum_t num;
rtden_t den;
get_rational(x, num, den);
printf("x=%s: rt=%lld/%lld\n", x_str, num, den);
}
{
const char* x_str = "123.345678";
const double x = 123.345678;
rtnum_t num;
rtden_t den;
get_rational(x, num, den);
printf("x=%s: rt=%lld/%lld\n", x_str, num, den);
}
{
const char* x_str = "123.3456789";
const double x = 123.3456789;
rtnum_t num;
rtden_t den;
get_rational(x, num, den);
printf("x=%s: rt=%lld/%lld\n", x_str, num, den);
}
{
const char* x_str = "123.345699999999";
const double x = 123.345699999999;
rtnum_t num;
rtden_t den;
get_rational(x, num, den);
printf("x=%s: rt=%lld/%lld\n", x_str, num, den);
}
{
const char* x_str = "12.34567";
const double x = 12.34567;
rtnum_t num;
rtden_t den;
get_rational(x, num, den);
printf("x=%s: rt=%lld/%lld\n", x_str, num, den);
}
{
const char* x_str = "12";
const double x = 12;
rtnum_t num;
rtden_t den;
get_rational(x, num, den);
printf("x=%s: rt=%lld/%lld\n", x_str, num, den);
}
{
const char* x_str = "123E-5";
const double x = 123E-5;
rtnum_t num;
rtden_t den;
get_rational(x, num, den);
printf("x=%s: rt=%lld/%lld\n", x_str, num, den);
}
}
int main()
{
RationalTest();
}