Конвертирование из unsigned long long в плавающее с округлением до ближайшего

Question

Мне нужно написать функцию, которая округляет от unsigned long long до плавающего, и округление должно быть ближе к четному. Я не могу просто выполнить приведение типа C ++, поскольку AFAIK стандарт не определяет округление. Я думал об использовании boost :: numeric, но я не смог найти полезного руководства после прочтения документации. Можно ли это сделать с помощью этой библиотеки? Конечно, если есть альтернатива, я был бы рад ее использовать.

Любая помощь будет высоко ценится.

РЕДАКТИРОВАТЬ: Добавление примера, чтобы сделать вещи немного яснее. Предположим, я хочу преобразовать 0xffffff7fffffffff в его представление с плавающей запятой. Стандарт C ++ допускает одно из:

1. 0x5f7fffff ~ 1.9999999 * 2 ^ 63
2. 0x5f800000 = 2 ^ 64

Теперь, если вы добавите ограничение раунда к ближайшему четному, приемлем только первый результат.

Answer 1

Поскольку в источнике так много битов, которые невозможно представить в float, и вы не можете (очевидно) полагаться на преобразование языка, вам придется сделать это самостоятельно.

Я разработал схему, которая может или не может помочь вам.По сути, в float есть 31 бит для представления положительных чисел, поэтому я выбрал 31 наиболее значимый бит в номере источника.Затем я сохраняю и маскирую все нижние биты.Затем, основываясь на значении младших битов, я округляю «новый» LSB вверх или вниз и, наконец, использую static_cast, чтобы создать float.

Я оставил в некоторых койках, которые вы можете удалить по своему усмотрению.

const unsigned long long mask_bit_count = 31;

float ull_to_float2(unsigned long long val)
{
    // How many bits are needed?
    int b = sizeof(unsigned long long) * CHAR_BIT - 1;
    for(; b >= 0; --b)
    {
        if(val & (1ull << b))
        {
            break;
        }
    }

    std::cout << "Need " << (b + 1) << " bits." << std::endl;

    // If there are few enough significant bits, use normal cast and done.
    if(b < mask_bit_count)
    {
        return static_cast<float>(val & ~1ull);
    }

    // Save off the low-order useless bits:
    unsigned long long low_bits = val & ((1ull << (b - mask_bit_count)) - 1);
    std::cout << "Saved low bits=" << low_bits << std::endl;

    std::cout << val << "->mask->";
    // Now mask away those useless low bits:
    val &= ~((1ull << (b - mask_bit_count)) - 1);
    std::cout << val << std::endl;

    // Finally, decide how to round the new LSB:
    if(low_bits > ((1ull << (b - mask_bit_count)) / 2ull))
    {
        std::cout << "Rounding up " << val;
        // Round up.
        val |= (1ull << (b - mask_bit_count));
        std::cout << " to " << val << std::endl;
    }
    else
    {
        // Round down.
        val &= ~(1ull << (b - mask_bit_count));
    }

    return static_cast<float>(val);
}

Answer 2

Я сделал это в Smalltalk для целого числа произвольной точности (LargeInteger), реализовал и протестировал в Squeak / Pharo / Visualworks / Gnu Smalltalk / Dolphin Smalltalk и даже написал в блоге об этом, если вы можете прочитать код Smalltalk http://smallissimo.blogspot.fr/2011/09/clarifying-and-optimizing.html.
Уловка для ускорения алгоритма заключается в следующем: FPU, совместимый с IEEE 754, округляет в точности результат неточной операции.Таким образом, мы можем позволить себе одну неточную операцию и позволить аппаратным средствам правильно работать.Это позволяет нам легко обрабатывать первые 48 бит.Но мы не можем позволить себе две неточные операции, поэтому нам иногда приходится по-разному заботиться о младших битах ...
Надеюсь, код достаточно документирован:

#include <math.h>
#include <float.h>
float ull_to_float3(unsigned long long val)
{
    int prec=FLT_MANT_DIG ;             // 24 bits, the float precision
    unsigned long long high=val>>prec;  // the high bits above float precision
    unsigned long long mask=(1ull<<prec) - 1 ;      // 0xFFFFFFull a mask for extracting significant bits
    unsigned long long tmsk=(1ull<<(prec - 1)) - 1; // 0x7FFFFFull same but tie bit
    // handle trivial cases, 48 bits or less,
    // let FPU apply correct rounding after exactly 1 inexact operation
    if( high <= mask )
        return ldexpf((float) high,prec) + (float) (val & mask);
    // more than 48 bits,
    // what scaling s is needed to isolate highest 48 bits of val?
    int s = 0;
    for( ; high > mask ; high >>= 1) ++s;
    // high now contains highest 24 bits
    float f_high = ldexpf( (float) high , prec + s );
    // store next 24 bits in mid
    unsigned long long mid = (val >> s) & mask;
    // care of rare case when trailing low bits can change the rounding:
    // can mid bits be a case of perfect tie or perfect zero?
    if( (mid & tmsk) == 0ull )
    {
        // if low bits are zero, mid is either an exact tie or an exact zero
        // else just increment mid to distinguish from such case
        unsigned long long low = val & ((1ull << s) - 1);
        if(low > 0ull) mid++;
    }
    return f_high + ldexpf( (float) mid , s );
}

Бонус: этот код должен округляться в соответствии с вашимРежим округления FPU, какой бы он ни был, поскольку мы косвенно использовали FPU для выполнения округления с помощью операции +.
Однако остерегайтесь агрессивной оптимизации в стандартах Если вы всегда хотите округлить до ближайшего четного, независимо от текущего режима округления, то вам придется увеличивать старшие биты, когда:

• средние биты> tie, где tie = 1ull << (prec-1); </li>
• средние биты == tie и (младшие биты> 0 или старшие биты нечетны).

РЕДАКТИРОВАТЬ:
Если вы придерживаетесь разрыва связи с округлением до ближайшего четного, то другим решением является использование ШУЧУКА РАСШИРЕНИЯ СУММЫ несмежных частей (fhigh, flow) и (fmid), см. http://www -2.cs.cmu.edu/afs/cs/project/quake / public /apers / robust-arithmetic.ps :

#include <math.h>
#include <float.h>
float ull_to_float4(unsigned long long val)
{
    int prec=FLT_MANT_DIG ;             // 24 bits, the float precision
    unsigned long long mask=(1ull<<prec) - 1 ; // 0xFFFFFFull a mask for extracting significant bits
    unsigned long long high=val>>(2*prec);     // the high bits
    unsigned long long mid=(val>>prec) & mask; // the mid bits
    unsigned long long low=val & mask;         // the low bits
    float fhigh = ldexpf((float) high,2*prec);
    float fmid  = ldexpf((float) mid,prec);
    float flow  = (float) low;
    float sum1 = fmid + flow;
    float residue1 = flow - (sum1 - fmid);
    float sum2 = fhigh + sum1;
    float residue2 = sum1 - (sum2 - fhigh);
    return (residue1 + residue2) + sum2;
}

Это делает алгоритм без ответвлений с немного большим количеством операций.Он может работать с другими режимами округления, но я позволю вам проанализировать документ, чтобы убедиться ...

Answer 3

То, что возможно между 8-байтовыми целыми числами и форматом с плавающей запятой, легко объяснить, но не так легко реализовать!

Следующий абзац касается того, что представимо в 8-байтовых целых числах со знаком.

Все положительные целые числа от 1 (2 ^ 0) до 16777215 (2 ^ 24-1) точно представлены в iEEE754 с одинарной точностью (с плавающей точкой).Или, если быть точным, все числа между 2 ^ 0 и 2 ^ 24-2 ^ 0 с шагом 2 ^ 0.Следующий диапазон точно представимых натуральных чисел составляет от 2 ^ 1 до 2 ^ 25-2 ^ 1 с шагом 2 ^ 1 и т. Д. До 2 ^ 39 до 2 ^ 63-2 ^ 39 с шагом 2 ^ 39.

Беззнаковые 8-байтовые целочисленные значения могут быть выражены с точностью до 2 ^ 64-2 ^ 40 с шагом 2 ^ 40.

Здесь формат одиночного пречисона не останавливается, а распространяется на всевплоть до диапазона от 2 ^ 103 до 2 ^ 127-2 ^ 103 с шагом 2 ^ 103.

Для 4-байтовых целых (длинных) самый высокий диапазон с плавающей запятой составляет от 2 ^ 7 до 2 ^ 31-2 ^ 7 с шагом 2 ^ 7.

В архитектуре x86 наибольшим целочисленным типом, поддерживаемым набором команд с плавающей запятой, является 8-байтовое целое число со знаком.2 ^ 64-1 не могут быть загружены обычными средствами.

Это означает, что для данного приращения диапазона, выраженного как "2 ^ i, где i является целым числом> 0", все целые числа, которые заканчиваются битовой комбинацией 0x1 вверхдо 2 ^ i-1 не будет точно отображаться в пределах этого диапазона в плавающем значении. Это означает, что то, что вы называете округлением вверх, на самом деле зависит от того, в каком диапазоне вы работаете. Бесполезно пытаться округлить до 1 (2^ 0) или 16 (2 ^ 4), если гранулярность диапазона, в котором вы находитесь, составляет 2 ^ 19.

Дополнительное следствие того, что вы предлагаете сделать (округление 2 ^ 63-1 до 2 ^ 63) может привести к переполнению (длинный целочисленный формат), если вы попытаетесь выполнить следующее преобразование: longlong_int = (long long) ((float) 2 ^ 63).

Проверьте эту небольшую программу, которую я написал (на C)что должно помочь проиллюстрировать, что возможно, а что нет.

int main (void)
{
  __int64 basel=1,baseh=16777215,src,dst,j;
  float cnvl,cnvh,range;
  int i=0;

  while (i<40)
  {
    src=basel<<i;
    cnvl=(float) src;
    dst=(__int64) cnvl;    /* compare dst with basel */

    src=baseh<<i;
    cnvh=(float) src;
    dst=(__int64) cnvh;    /* compare dst with baseh */

    j=basel;
    while (j<=baseh)
    {
      range=(float) j;
      dst=(__int64) range;

      if (j!=dst) dst/=0;

      j+=basel;
    }

    ++i;
  }
  return i;
}

Эта программа показывает представимые диапазоны целочисленных значений.Их перекрывают: например, 2 ^ 5 представимо во всех диапазонах с нижней границей 2 ^ b, где 1 =