Потерялся 1 бит в long double

224
25 октября 2017, 05:57

По следам вопроса о битовом представлении вещественных чисел и моего ответа на него.

Хочу программно для любого вещественного типа определить, сколько бит в нём отводится под мантиссу, а сколько под экспоненту. Для этого написал следующий код (в нём бит под знак считается отдельно от мантиссы, поэтому числа на 1 меньше):

https://ideone.com/YuIWNc - код на Си (float, double, long double)
https://ideone.com/342B4S - код на Си++ (float, double, long double)
https://ideone.com/VURQnw - код на Си++ (float, double, long double, __float128)

#include <cstdio>
template <typename typed> void count(unsigned *result_m, unsigned *result_e)
{
  typed x = 1, exp;
  unsigned res, e;
  for (res=0; x!=0; ++res) x/=2;
  for (exp=1,e=0; exp*2<res; ++e) exp*=2;
  *result_e = e+1;
  *result_m = res-exp+1;
}
int main(void)
{
  unsigned f_m, f_e, d_m, d_e, ld_m, ld_e, f128_m, f128_e;
  count<float>(&f_m, &f_e);
  count<double>(&d_m, &d_e);
  count<long double>(&ld_m, &ld_e);
  count<__float128>(&f128_m, &f128_e);
  printf("              S    M   E   SZ\n");
  printf("float:        1  %3u  %2u  %3u\n",    f_m,    f_e, 8 * sizeof(float));
  printf("double:       1  %3u  %2u  %3u\n",    d_m,    d_e, 8 * sizeof(double));
  printf("long double:  1  %3u  %2u  %3u\n",   ld_m,   ld_e, 8 * sizeof(long double));
  printf("__float128:   1  %3u  %2u  %3u\n", f128_m, f128_e, 8 * sizeof(__float128));
}

Получается так:

              S    M   E   SZ
float:        1   23   8   32
double:       1   52  11   64
long double:  1   63  15  128
__float128:   1  112  15  128

Для float, double и даже __float128 всё работает (Википедия, IEEE 754-2008).
А вот с long double возникают проблемы:

  1. 1+63+15 = 79 - 79 бит. Вместо 80. Где ещё один бит?
  2. long double представляет собой 10-байтные числа, но sizeof вернул 16.
    Как-то можно получить 10?
Answer 1

Один бит потерялся из-за того, что на платформе x86 80-битное плавающее значение имеет одно принципиальное отличие в представлении от 32- и 64-битных IEEE754 плавающих значений (float и double).

float и double используют представление с неявной ведущей единицей в мантиссе. То есть в нормализованном представлении старшая единица в мантиссе не хранится явно, а лишь подразумевается. А вот в расширенном 80-битном плавающем типе long double эта ведущая единица в мантиссе всегда хранится явно.

Из-за этого и возникает разница.

Для типов float и double ваш первый цикл сначала итерирует через нормализованные представления числа, в которых явная мантисса всегда равна нулю, а экспонента уменьшается от половины своего максимального значения (127 для float) до значения 1:

// Для `float`
// Нормализованные представления: мантисса равна 0, а экспонента убывает от 127 до 1
0x3F800000
...
0x00800000  <- после 126 делений

После этого ваш цикл продолжает итерировать через денормализованные представления числа, в которых экспонента равна 0, а по мантиссе вправо движется одинокая единица. Когда эта одинокая единица вылетит за правый край мантиссы, x станет равен нулю и цикл завершится

// Денормализованные представления: экспонента равна 0, а мантисса состоит
// из движущейся вправо единицы
0x00400000
0x00200000
...
0x00000001
0x00000000  <- после 150 делений

Заметьте, в float и double единица в мантиссе возникает только в самом первом денормализованном значении и проходит через все биты мантиссы. Получается, что количество денормализованных ненулевых значений в этом случае равно количеству битов в мантиссе.

Однако при использовании long double единица в старшем бите мантиссы явно присутствовала всегда, с самого начала. Когда в вашем цикле экспонента long double достигает нуля и цикл начинает подсчитывать денормализованные значения long double, единица в мантиссе не возникает "из никуда" в старшую позицию мантиссы (как это было в в float и double), а уже присутствует в старшей позиции изначально и "стартует" именно оттуда. Из-за этого та часть цикла, которая считает денормализованные значения, делает на одну итерацию меньше.

Кстати, странная манера складывать в одну сумму - res - половину диапазона экспоненты и ширину мантиссы чревата проблемами. Вы потом вычисляете значение log2 res и ожидаете, что эта величина будет правильно описывать количество битов в экспоненте. Однако если в некотором гипотетическом плавающем типе мантисса окажется очень широкой, то величина log2 res может оказаться ошибочной.

Answer 2

Справка по C/C++ сообщает удивительную подробность:

long double — тип с плавающей точкой повышенной точности. Не обязательно отображается на типы IEEE-754. Обычно 80-битный тип с плавающей точкой формата x87 на архитектурах x86 и x86-64.

Получается, ответственность за реализацию лежит целиком на компиляторе.

Чтобы проверить, куда и какой бит делся, можно добавить распечатку в процесс "обнаружения" всех бит:

#include <cstdio>
#include <cstring>
#include <cinttypes>
template <typename typed> void count(unsigned *result_m, unsigned *result_e)
{
  typed x = 1, exp;
  unsigned res, e;
  std::uint32_t bytes[4] {0,0,0,0};
  for (res=0; x!=0; ++res) {
    x/=2;
    std::memcpy(bytes, &x, sizeof(typed));
    printf("x(%3d): %08x %08x %08x %08x\n", res, bytes[0], bytes[1], bytes[2], bytes[3]);
  }
  for (exp=1,e=0; exp*2<res; ++e) exp*=2;
  *result_e = e+1;
  *result_m = res-exp+1;
}
int main(void)
{
  unsigned f_m, f_e, d_m, d_e, ld_m, ld_e;
  count<float>(&f_m, &f_e);
  count<double>(&d_m, &d_e);
  count<long double>(&ld_m, &ld_e);
  printf("              S   M   E   SZ\n");
  printf("float:        1  %2u  %2u  %3u\n",  f_m,  f_e, 8 * sizeof(float));
  printf("double:       1  %2u  %2u  %3u\n",  d_m,  d_e, 8 * sizeof(double));
  printf("long double:  1  %2u  %2u  %3u\n", ld_m, ld_e, 8 * sizeof(long double));
}

Так, например, для 32-битного числа:

x(  0): 3f000000
x(  1): 3e800000
x(  2): 3e000000
x(  3): 3d800000
...с номера 126 наступает денормализация:
x(124): 01000000
x(125): 00800000
x(126): 00400000
x(127): 00200000
x(128): 00100000
x(129): 00080000
x(130): 00040000
...и далее перебор заканчивается:
x(147): 00000002
x(148): 00000001
x(149): 00000000

Где бит пропал - там и заковыка. Аналогично будут видны все незадействованные биты, и так можно оценить эффективный размер типа.

Продолжение для long double:

x(    0): 00000000 80000000 00003ffe 00000000
x(    1): 00000000 80000000 00003ffd 00000000
x(    2): 00000000 80000000 00003ffc 00000000
x(    3): 00000000 80000000 00003ffb 00000000
...денормализация:
x(16379): 00000000 80000000 00000003 00000000
x(16380): 00000000 80000000 00000002 00000000
x(16381): 00000000 80000000 00000001 00000000
x(16382): 00000000 40000000 00000000 00000000
x(16383): 00000000 20000000 00000000 00000000
...
x(16410): 00000000 00000004 00000000 00000000
x(16411): 00000000 00000002 00000000 00000000
x(16412): 00000000 00000001 00000000 00000000
x(16413): 80000000 00000000 00000000 00000000
x(16414): 40000000 00000000 00000000 00000000
...финал:
x(16443): 00000002 00000000 00000000 00000000
x(16444): 00000001 00000000 00000000 00000000
x(16445): 00000000 00000000 00000000 00000000

Итак, один бит действительно всегда выставлен в 1 для нормальных чисел, в соответствии с английской википедией на тему Extended Precision - и становится 0 после денормализации. В "стандартных" 32-битных и 64-битных представлениях этот бит опускается

READ ALSO
set и unordered_set

set и unordered_set

Все мы знаем, что std::set - это красно-чёрное дерево (соответственно со сложностью поиска O(log n)), а std::unordered_set - хэш-таблица с поиском за константуКакие...

254
Увеличение времени исполнения кода

Увеличение времени исполнения кода

Столкнулся со следующей проблемой: Есть код, написанный не одним человеком, который работает с использованием MPIВ части последовательного...

202
Высвобождение памяти

Высвобождение памяти

Не могу найти где забываю высвободить памятьПодскажите,пожалуйста

254
Выполнение 2х функций в одном файле javascript

Выполнение 2х функций в одном файле javascript

Изобретаю велосипед, и столкнулась с проблемойЕсть код, в котором должны обрабатываться 2 клика

207