Напишите свою собственную реализацию функции математического пола, C
Я думал о функции floor, доступной в math.h. Его очень легко использовать:
#include <stdio.h>
#include <math.h>
int main(void)
{
for (double a = 12.5; a < 13.4; a += 0.1)
printf("floor of %.1lf is %.1lf\n", a, floor(a));
return 0;
}
Что, если я хотел бы написать свою собственную реализацию? Будет ли это выглядеть так:
#include <stdio.h>
#include <math.h>
double my_floor(double num)
{
return (int)num;
}
int main(void)
{
double a;
for (a = 12.5; a < 13.4; a += 0.1)
printf("floor of %.1lf is %.1lf\n", a, floor(a));
printf("\n\n");
for (a = 12.5; a < 13.4; a += 0.1)
printf("floor of %.1lf is %.1lf\n", a, my_floor(a));
return 0;
}
?
Кажется, он не работает с отрицательными числами (my_floor), но второй кажется прекрасным (my_floor_2):
#include <stdio.h>
#include <math.h>
double my_floor(double num)
{
return (int)num;
}
double my_floor_2(double num)
{
if(num < 0)
return (int)num - 1;
else
return (int)num;
}
int main(void)
{
double a1 = -12.5;
printf("%lf\n", floor(a1));
printf("%lf\n", my_floor(a1));
printf("%lf\n", my_floor_2(a1));
return 0;
}
выход программы:
![img]()
Является ли один из них в конечном итоге правильным или нет?
Ответы
Ответ 1
Нет, вы не можете справиться с этим так. Лучший способ написать свою собственную реализацию - украсть ту, что из стандартной библиотеки C на вашей платформе. Но обратите внимание, что могут содержать специфические нюансы платформы, поэтому они могут быть не переносимыми.
Функция стандартной библиотеки C floor, как правило, умна тем, что она не работает, перейдя к интегральному типу. Если это произошло, вы рискуете переполнением целого числа signed, поведение которого равно undefined. (Обратите внимание, что наименьший возможный диапазон для int составляет от -32767 до +32767).
Точная реализация также зависит от схемы с плавающей точкой, используемой на вашей платформе.
Для платформы с плавающей точкой IEEE754 и типа long long вы можете использовать эту схему:
- Если величина числа больше 53-й степени 2, верните ее назад (поскольку она уже является интегральной).
- Else, добавьте к 64-битовому типу (длинный длинный) и верните его.
Ответ 2
Обе ваши попытки имеют ограничения:
- Если значение
double находится за пределами диапазона типа int, преобразование в int определяется реализацией.
- Если значение
double отрицательное, но целое, возврат (int)num - 1 неверен.
Вот (почти) переносная версия, которая пытается обрабатывать все случаи:
double my_floor_2(double num) {
if (num >= LLONG_MAX || num <= LLONG_MIN || num != num) {
/* handle large values, infinities and nan */
return num;
}
long long n = (long long)num;
double d = (double)n;
if (d == num || num >= 0)
return d;
else
return d - 1;
}
Это должно быть правильно, если тип long long имеет больше битов значений, чем тип double, что имеет место в большинстве современных систем.
Ответ 3
В С++ и 32-битной арифметике это можно сделать, например, следующим образом:
//---------------------------------------------------------------------------
// IEEE 754 double MSW masks
const DWORD _f64_sig =0x80000000; // sign
const DWORD _f64_exp =0x7FF00000; // exponent
const DWORD _f64_exp_sig=0x40000000; // exponent sign
const DWORD _f64_exp_bia=0x3FF00000; // exponent bias
const DWORD _f64_exp_lsb=0x00100000; // exponent LSB
const DWORD _f64_exp_pos= 20; // exponent LSB bit position
const DWORD _f64_man =0x000FFFFF; // mantisa
const DWORD _f64_man_msb=0x00080000; // mantisa MSB
const DWORD _f64_man_bits= 52; // mantisa bits
// IEEE 754 single masks
const DWORD _f32_sig =0x80000000; // sign
const DWORD _f32_exp =0x7F800000; // exponent
const DWORD _f32_exp_sig=0x40000000; // exponent sign
const DWORD _f32_exp_bia=0x3F800000; // exponent bias
const DWORD _f32_exp_lsb=0x00800000; // exponent LSB
const DWORD _f32_exp_pos= 23; // exponent LSB bit position
const DWORD _f32_man =0x007FFFFF; // mantisa
const DWORD _f32_man_msb=0x00400000; // mantisa MSB
const DWORD _f32_man_bits= 23; // mantisa bits
//---------------------------------------------------------------------------
double f64_floor(double x)
{
const int h=1; // may be platform dependent MSB/LSB order
const int l=0;
union _f64 // semi result
{
double f; // 64bit floating point
DWORD u[2]; // 2x32 bit uint
} y;
DWORD m,a;
int sig,exp,sh;
y.f=x;
// extract sign
sig =y.u[h]&_f64_sig;
// extract exponent
exp =((y.u[h]&_f64_exp)>>_f64_exp_pos)-(_f64_exp_bia>>_f64_exp_pos);
// floor bit shift
sh=_f64_man_bits-exp; a=0;
if (exp<0)
{
a=y.u[l]|(y.u[h]&_f64_man);
if (sig) return -1.0;
return 0.0;
}
// LSW
if (sh>0)
{
if (sh<32) m=(0xFFFFFFFF>>sh)<<sh; else m=0;
a=y.u[l]&(m^0xFFFFFFFF); y.u[l]&=m;
}
// MSW
sh-=32;
if (sh>0)
{
if (sh<_f64_exp_pos) m=(0xFFFFFFFF>>sh)<<sh; else m=_f64_sig|_f64_exp;
a|=y.u[h]&(m^0xFFFFFFFF); y.u[h]&=m;
}
if ((sig)&&(a)) y.f--;
return y.f;
}
//---------------------------------------------------------------------------
float f32_floor(float x)
{
union // semi result
{
float f; // 32bit floating point
DWORD u; // 32 bit uint
} y;
DWORD m,a;
int sig,exp,sh;
y.f=x;
// extract sign
sig =y.u&_f32_sig;
// extract exponent
exp =((y.u&_f32_exp)>>_f32_exp_pos)-(_f32_exp_bia>>_f32_exp_pos);
// floor bit shift
sh=_f32_man_bits-exp; a=0;
if (exp<0)
{
a=y.u&_f32_man;
if (sig) return -1.0;
return 0.0;
}
if (sh>0)
{
if (sh<_f32_exp_pos) m=(0xFFFFFFFF>>sh)<<sh; else m=_f32_sig|_f32_exp;
a|=y.u&(m^0xFFFFFFFF); y.u&=m;
}
if ((sig)&&(a)) y.f--;
return y.f;
}
//---------------------------------------------------------------------------
Цель состоит в том, чтобы сделать маску, которая очистит десятичные биты от мантиссы, а в случае отрицательного ввода и ненулевые очищенные биты уменьшат результат. Чтобы получить доступ к отдельным битам, вы можете преобразовать значение с плавающей запятой в интегральное представление с использованием объединения (например, в примере) или вместо этого использовать указатели.
Я тестировал это в простом VCL приложении:
float f32;
double f64;
AnsiString txt="";
// 64 bit
txt+="[double]\r\n";
for (f64=-10.0;f64<=10.0;f64+=0.1)
if (fabs(floor(f64)-f64_floor(f64))>1e-6)
{
txt+=AnsiString().sprintf("%5.3lf %5.3lf %5.3lf\r\n",f64,floor(f64),f64_floor(f64));
f64_floor(f64);
}
for (f64=1;f64<=1e307;f64*=1.1)
{
if (fabs(floor( f64)-f64_floor( f64))>1e-6) { txt+=AnsiString().sprintf("%lf lf lf\r\n", f64,floor( f64),f64_floor( f64));
f64_floor( f64); }
if (fabs(floor(-f64)-f64_floor(-f64))>1e-6) { txt+=AnsiString().sprintf("%lf lf lf\r\n",-f64,floor(-f64),f64_floor(-f64));
f64_floor(-f64); }
}
// 32 bit
txt+="[float]\r\n";
for (f32=-10.0;f32<=10.0;f32+=0.1)
if (fabs(floor(f32)-f32_floor(f32))>1e-6)
{
txt+=AnsiString().sprintf("%5.3lf %5.3lf %5.3lf\r\n",f32,floor(f32),f32_floor(f32));
f32_floor(f32);
}
for (f32=1;f32<=1e37;f32*=1.1)
{
if (fabs(floor( f32)-f32_floor( f32))>1e-6) { txt+=AnsiString().sprintf("%lf lf lf\r\n", f32,floor( f32),f32_floor( f32));
f32_floor( f32); }
if (fabs(floor(-f32)-f32_floor(-f32))>1e-6) { txt+=AnsiString().sprintf("%lf lf lf\r\n",-f32,floor(-f32),f32_floor(-f32));
f32_floor(-f32); }
}
mm_log->Lines->Add(txt);
без разностного результата (поэтому во всех проверенных случаях он соответствует значениям math.h floor(). Если вы хотите сделать этот снимок вне VCL, просто измените AnsiString на любой тип строки вы получили под рукой и изменили вывод с TMemo::mm_log на все, что у вас было (например, console cout или что-то еще)
Двойной вызов fxx_floor() в случае разницы - для целей отладки (вы можете поместить точку останова и шаг в случае ошибки напрямую).
[Примечание]
Остерегайтесь того, что порядок слов (MSW, LSW) зависит от платформы, поэтому вы должны соответствующим образом скорректировать константы h,l. Этот код не оптимизирован, поэтому его легко понять, поэтому не ожидайте, что он будет быстрым.
