Напишите свою собственную реализацию функции пола математики, C

Question

Я думал о функции 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; 
} 

выхода программы:

Является ли один из них в конечном итоге правильно или нет?

Answer 1

Нет, вы не можете справиться с этим таким образом. Лучший способ написать свою собственную реализацию - украсть ту, что из стандартной библиотеки C на вашей платформе. Но обратите внимание, что могут содержать специфические нюансы платформы, поэтому они могут быть не переносимыми.

Функция стандартной библиотеки C floor, как правило, умна тем, что она не работает путем преобразования в интегральный тип. Если бы это произошло, вы столкнулись бы с риском переполнения целых чисел signed, поведение которого не определено. (Обратите внимание, что наименьший возможный диапазон для int составляет от -32767 до +32767).

Точная реализация также зависит от схемы с плавающей точкой, используемой на вашей платформе.

Для платформы с использованием IEEE754 с плавающей точкой, а long long типа вы могли принять эту схему:

1. Если величина числа больше, чем 53-й степени 2, вернуть его обратно (как он уже интегрален).
2. Иначе, бросьте на 64-битный тип (длинный длинный) и верните его обратно.

Answer 2

В C++ и 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.hfloor() значения.Если вы хотите, чтобы дать этому выстрел вне VCL затем просто изменить AnsiString к любому типу строки вы получили в руки и изменить выход из TMemo::mm_log все, что вы получили (например, консоли cout или любой другой)

Двойной призвании fxx_floor() в случае разницы для целей отладки (вы можете поместить точку останова и шаг в случае ошибки непосредственно).

[Примечания]

Берегитесь порядком слов (MSW, LSW) является платформо-зависимым, так что вы должны настроить h,l константы соответственно. Этот код не оптимизирован, поэтому его легко понять, поэтому не ожидайте, что он будет быстрым.

Answer 3

Обе ваши попытки имеют ограничения:

• Если значение 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, что имеет место на самом современном системы.