Как вы печатаете точное значение числа с плавающей запятой? - PullRequest
31 голосов
/ 09 июля 2010

Прежде всего, это не вопрос новичка с плавающей запятой. Я знаю, что результаты арифметики с плавающей запятой (не говоря уже о трансцендентных функциях) обычно не могут быть представлены точно, и что большинство конечных десятичных чисел не могут быть представлены точно как двоичные числа с плавающей запятой.

Тем не менее, каждое возможное значение с плавающей запятой в точности соответствует диадическому рациональному числу (рациональное число p/q, где q - степень 2), которое, в свою очередь, имеет точное десятичное представление.

Мой вопрос: как вы находите это точное десятичное представление эффективно? sprintf и аналогичные функции обычно указываются с точностью до числа значащих цифр для однозначного определения исходного значения с плавающей запятой; они не обязательно печатают точное десятичное представление. Я знаю один алгоритм, который использовал, но он очень медленный, O(e^2), где e - показатель степени. Вот схема:

  1. Преобразование мантиссы в десятичное целое число. Вы можете сделать это, разделив биты для непосредственного чтения мантиссы, или вы можете написать грязный цикл с плавающей запятой, который сначала умножает значение на степень два, чтобы поместить его в диапазон 1 <= x <10, а затем вытягивает за раз от цифры, приводя к int, вычитая и умножая на 10. </li>
  2. Примените показатель степени путем многократного умножения или деления на 2. Это операция над строкой сгенерированных вами десятичных цифр. Каждые ~ 3 умножения добавят дополнительную цифру слева. Каждое отдельное деление добавит дополнительную цифру справа.

Это действительно лучшее из возможных? Я сомневаюсь в этом, но я не эксперт с плавающей запятой, и я не могу найти способ выполнить вычисления base-10 для представления числа с плавающей запятой без возможности получения неточных результатов (умножение или деление на все, кроме степени 2, является операцией с потерями для чисел с плавающей запятой, если только вы не знаете, что у вас есть свободные биты для работы).

Ответы [ 10 ]

34 голосов
/ 09 июля 2010

Этот вопрос имеет бюрократическую и алгоритмическую части.Число с плавающей запятой хранится внутри как (2 e × m), где e - показатель степени (сам двоичный), а m - мантисса.Бюрократическая часть вопроса заключается в том, как получить доступ к этим данным, но Р., похоже, больше интересует алгоритмическая часть вопроса, а именно преобразование (2 e × m) в дробь (a / b) вдесятичная форма.Ответ на бюрократический вопрос на нескольких языках - «frexp» (это интересная деталь, которую я не знал до сегодняшнего дня).

Это правда, что на первый взгляд требуется O (e 2 ) работает только для записи 2 e в десятичном виде, и еще больше времени для мантиссы.Но, благодаря волшебству алгоритма быстрого умножения Schonhage-Strassen , вы можете сделать это за Õ (e) времени, где тильда означает «до логарифмических факторов».Если вы рассматриваете Schonhage-Strassen как магию, то не так сложно думать о том, что делать.Если e чётно, вы можете рекурсивно вычислить 2 e / 2 , а затем возвести его в квадрат с помощью быстрого умножения.С другой стороны, если e нечетно, вы можете рекурсивно вычислить 2 e-1 и затем удвоить его.Вы должны быть осторожны, чтобы проверить, есть ли версия Schonhage-Strassen в базе 10. Хотя это не широко документировано, это может быть сделано в любой базе.

Преобразование очень длинной мантиссы из двоичной в базовую10 не совсем тот же вопрос, но он имеет аналогичный ответ.Вы можете разделить мантиссу на две половины, m = a 2 k + b.Затем рекурсивно преобразуйте a и b в основание 10, преобразуйте 2 ^ k в основание 10 и выполните другое быстрое умножение для вычисления m в основании 10.

Абстрактный результат, лежащий в основе всего этого, заключается в том, что вы можете преобразовать целые числа изодна база к другой за время Õ (N).

Если вопрос касается стандартных 64-битных чисел с плавающей запятой, то он слишком мал для причудливого алгоритма Шонхаге-Штрассена.В этом диапазоне вы можете вместо этого сохранить работу с различными трюками.Один из подходов состоит в том, чтобы сохранить все 2048 значений 2 e в справочной таблице, а затем работать в мантиссе с асимметричным умножением (между длинным умножением и коротким умножением).Другой трюк заключается в работе на базе 10000 (или более высокой степени 10, в зависимости от архитектуры) вместо базы 10. Но, как указывает Р. в комментариях, 128-битные числа с плавающей запятой уже позволяют достаточно большим показателям вызыватьвопрос и таблицы поиска и стандартное длинное умножение.На практике длинное умножение является самым быстрым, вплоть до нескольких цифр, затем в значительном среднем диапазоне можно использовать умножение Карацубы или умножение Тоом-Кука , а затем после этоговариант Schonhage-Strassen лучше всего подходит не только в теории, но и на практике.

На самом деле, большой целочисленный пакет GMP уже имеет radi (N) преобразование радиуса времени, а также хорошееэвристика, для которой выбор алгоритма умножения.Единственное различие между их решением и моим состоит в том, что вместо выполнения какой-либо большой арифметики в базе 10 они вычисляют большие степени 10 в базе 2. В этом решении им также требуется быстрое деление, но это можно получить из быстрого умножения в любомнескольких способов.

16 голосов
/ 10 июля 2010

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

  1. Функция Дэвида Гея dtoa() в dtoa.c (http://www.netlib.org/fp/dtoa.c).

  2. Функция ___printf_fp() в файле /stdio-common/printf_fp.c в glibc (например, http://ftp.gnu.org/gnu/glibc/glibc-2.11.2.tar.gz,).

Обе будут печатать столько цифр, сколько вы запрашиваете в %f типе printf (как я писал здесь: http://www.exploringbinary.com/print-precision-of-dyadic-fractions-varies-by-language/ и http://www.exploringbinary.com/print-precision-of-floating-point-integers-varies-too/).

5 голосов
/ 09 июля 2010

Была проделана большая работа по печати чисел с плавающей точкой.Золотой стандарт - распечатывать десятичный эквивалент минимальной длины, чтобы при считывании десятичного эквивалента вы получали то же число с плавающей запятой, с которого вы начали, независимо от режима округления во время обратного чтения.Вы можете прочитать об алгоритме в превосходной статье от Burger и Dybvig .

3 голосов
/ 09 июля 2010

Хотя это C # и ваш вопрос помечен C, у Джона Скита есть код для преобразования double в его точное представление в виде строки: http://www.yoda.arachsys.com/csharp/DoubleConverter.cs

С первого взгляда не слишком сложно портировать на C, и даже легче писать на C ++.

При дальнейшем размышлении выясняется, что алгоритм Джона также является O (e ^ 2), поскольку он также проходит по показателю степени. Однако это означает, что алгоритм равен O (log (n) ^ 2) (где n - число с плавающей запятой), и я не уверен, что вы можете преобразовать из базы 2 в базу 10 быстрее, чем время в квадрате логарифма.

2 голосов
/ 09 июля 2010

Будучи не экспертом по плавающей запятой, я бы предпочел использовать хорошо протестированную библиотеку с открытым исходным кодом.

GNU MPFR - хороший вариант.

Библиотека MPFR - это библиотека C для вычислений с плавающей точкой с множественной точностью и правильным округлением.Основная цель MPFR - предоставить библиотеку для вычислений с плавающей точкой с множественной точностью, которая является одновременно эффективной и имеет четко определенную семантику.

1 голос
/ 13 июля 2010

sprintf и аналогичные функции обычно указываются с точностью до числа значащих цифр для однозначного определения исходного значения с плавающей запятой;они не обязательно печатают точное десятичное представление.

Вы можете запросить более значимые цифры, чем по умолчанию:

printf("%.100g\n", 0.1);

печатает 0.1000000000000000055511151231257827021181583404541015625.

0 голосов
/ 23 августа 2013

Есть 3 способа

  1. печать чисел в bin или hex

    Это самый точный способ. Я предпочитаю hex, потому что это больше похоже на базу 10 для чтения / ощущения, например F.8h = 15.5 без потери точности.

  2. печать в dec, но только соответствующие цифры

    Под этим я подразумеваю только те цифры, которые могут содержать 1 в вашем номере как можно ближе.

    num из целых цифр просты и точны (без потери точности):

    // n10 - base 10 integer digits
    // n2  - base 2 integer digits
    n10=log10(2^n2)
    n10=log2(2^n2)/log2(10)
    n10=n2/log2(10)
    n10=ceil(n2*0.30102999566398119521373889472449)
    // if fist digit is 0 and n10 > 1 then n10--
    

    num из дробные цифры более хитро и эмпирически я нашел это:

    // n10 - base 10 fract. digits
    // n2  - base 2 fract. digits >= 8
    n10=0; if (n02==8) n10=1;
    else if (n02==9) n10=2;
    else if (n02> 9)
        {
        n10=((n02-9)%10);
             if (n10>=6) n10=2;
        else if (n10>=1) n10=1;
        n10+=2+(((n02-9)/10)*3);
        }
    

    если вы создадите таблицу зависимостей n02 <-> n10, то увидите, что константа 0.30102999566398119521373889472449 все еще присутствует, но при старте с 8 битов, поскольку less не может представлять 0.1 с удовлетворительной точностью (0.85 - 1.15). из-за отрицательных показателей базы 2 зависимость не линейная, а паттерны. Этот код работает для небольшого количества битов (<=52), но при большем количестве битов может быть ошибка, поскольку используемый шаблон не соответствует log10(2) или 1/log2(10) точно.

    для большего количества битов я использую это:

    n10=7.810+(9.6366363636363636363636*((n02>>5)-1.0));
    

    но эта формула выровнена на 32 бита !!! а также большая ошибка объявления количества битов к нему

    P.S. дальнейший анализ двоичного представления десятичных чисел

    0.1
    0.01
    0.001
    0.0001
    ...
    

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

    для ясности:

    8 bin digits -> 1 dec digits
    9 bin digits -> 2 dec digits
    10 bin digits -> 3 dec digits
    11 bin digits -> 3 dec digits
    12 bin digits -> 3 dec digits
    13 bin digits -> 3 dec digits
    14 bin digits -> 3 dec digits
    15 bin digits -> 4 dec digits
    16 bin digits -> 4 dec digits
    17 bin digits -> 4 dec digits
    18 bin digits -> 4 dec digits
    19 bin digits -> 5 dec digits
    20 bin digits -> 6 dec digits
    21 bin digits -> 6 dec digits
    22 bin digits -> 6 dec digits
    23 bin digits -> 6 dec digits
    24 bin digits -> 6 dec digits
    25 bin digits -> 7 dec digits
    26 bin digits -> 7 dec digits
    27 bin digits -> 7 dec digits
    28 bin digits -> 7 dec digits
    29 bin digits -> 8 dec digits
    30 bin digits -> 9 dec digits
    31 bin digits -> 9 dec digits
    32 bin digits -> 9 dec digits
    33 bin digits -> 9 dec digits
    34 bin digits -> 9 dec digits
    35 bin digits -> 10 dec digits
    36 bin digits -> 10 dec digits
    37 bin digits -> 10 dec digits
    38 bin digits -> 10 dec digits
    39 bin digits -> 11 dec digits
    40 bin digits -> 12 dec digits
    41 bin digits -> 12 dec digits
    42 bin digits -> 12 dec digits
    43 bin digits -> 12 dec digits
    44 bin digits -> 12 dec digits
    45 bin digits -> 13 dec digits
    46 bin digits -> 13 dec digits
    47 bin digits -> 13 dec digits
    48 bin digits -> 13 dec digits
    49 bin digits -> 14 dec digits
    50 bin digits -> 15 dec digits
    51 bin digits -> 15 dec digits
    52 bin digits -> 15 dec digits
    53 bin digits -> 15 dec digits
    54 bin digits -> 15 dec digits
    55 bin digits -> 16 dec digits
    56 bin digits -> 16 dec digits
    57 bin digits -> 16 dec digits
    58 bin digits -> 16 dec digits
    59 bin digits -> 17 dec digits
    60 bin digits -> 18 dec digits
    61 bin digits -> 18 dec digits
    62 bin digits -> 18 dec digits
    63 bin digits -> 18 dec digits
    64 bin digits -> 18 dec digits
    

    И наконец не забудьте округлить срезанные цифры !!! Это означает, что если цифра после последней соответствующей цифры равна >=5 в декабре, то последняя последняя соответствующая цифра должна быть +1 ... и если она уже равна 9, то вы должны перейти к предыдущей цифре и т. Д. *

  3. печать точного значения

    Чтобы вывести точное значение дробного двоичного числа , просто выведите дробные n цифр, где n - это количество дробных битов, поскольку представленное значение является суммой этих значений, поэтому число дробные десятичные дроби не могут быть больше, чем num дробных цифр LSB . Материал выше (bullet # 2 ) имеет отношение к сохранению десятичного числа в float (или печати только соответствующих десятичных знаков).

    отрицательные степени двух точных значений ...

    2^- 1 = 0.5
    2^- 2 = 0.25
    2^- 3 = 0.125
    2^- 4 = 0.0625
    2^- 5 = 0.03125
    2^- 6 = 0.015625
    2^- 7 = 0.0078125
    2^- 8 = 0.00390625
    2^- 9 = 0.001953125
    2^-10 = 0.0009765625
    2^-11 = 0.00048828125
    2^-12 = 0.000244140625
    2^-13 = 0.0001220703125
    2^-14 = 0.00006103515625
    2^-15 = 0.000030517578125
    2^-16 = 0.0000152587890625
    2^-17 = 0.00000762939453125
    2^-18 = 0.000003814697265625
    2^-19 = 0.0000019073486328125
    2^-20 = 0.00000095367431640625
    

    теперь отрицательные степени 10 печатаются с использованием стиля точных значений для 64-битных doubles:

    10^+ -1 = 0.1000000000000000055511151231257827021181583404541015625
            = 0.0001100110011001100110011001100110011001100110011001101b
    10^+ -2 = 0.01000000000000000020816681711721685132943093776702880859375
            = 0.00000010100011110101110000101000111101011100001010001111011b
    10^+ -3 = 0.001000000000000000020816681711721685132943093776702880859375
            = 0.000000000100000110001001001101110100101111000110101001111111b
    10^+ -4 = 0.000100000000000000004792173602385929598312941379845142364501953125
            = 0.000000000000011010001101101110001011101011000111000100001100101101b
    10^+ -5 = 0.000010000000000000000818030539140313095458623138256371021270751953125
            = 0.000000000000000010100111110001011010110001000111000110110100011110001b
    10^+ -6 = 0.000000999999999999999954748111825886258685613938723690807819366455078125
            = 0.000000000000000000010000110001101111011110100000101101011110110110001101b
    10^+ -7 = 0.0000000999999999999999954748111825886258685613938723690807819366455078125
            = 0.0000000000000000000000011010110101111111001010011010101111001010111101001b
    10^+ -8 = 0.000000010000000000000000209225608301284726753266340892878361046314239501953125
            = 0.000000000000000000000000001010101111001100011101110001000110000100011000011101b
    10^+ -9 = 0.0000000010000000000000000622815914577798564188970686927859787829220294952392578125
            = 0.0000000000000000000000000000010001001011100000101111101000001001101101011010010101b
    10^+-10 = 0.00000000010000000000000000364321973154977415791655470655996396089904010295867919921875
            = 0.00000000000000000000000000000000011011011111001101111111011001110101111011110110111011b
    10^+-11 = 0.00000000000999999999999999939496969281939810930172340963650867706746794283390045166015625
            = 0.00000000000000000000000000000000000010101111111010111111111100001011110010110010010010101b
    10^+-12 = 0.00000000000099999999999999997988664762925561536725284350612952266601496376097202301025390625
            = 0.00000000000000000000000000000000000000010001100101111001100110000001001011011110101000010001b
    10^+-13 = 0.00000000000010000000000000000303737455634003709136034716842278413651001756079494953155517578125
            = 0.00000000000000000000000000000000000000000001110000100101110000100110100001001001011101101000001b
    10^+-14 = 0.000000000000009999999999999999988193093545598986971343290729163921781719182035885751247406005859375
            = 0.000000000000000000000000000000000000000000000010110100001001001101110000110101000010010101110011011b
    10^+-15 = 0.00000000000000100000000000000007770539987666107923830718560119501514549256171449087560176849365234375
            = 0.00000000000000000000000000000000000000000000000001001000000011101011111001111011100111010101100001011b
    10^+-16 = 0.00000000000000009999999999999999790977867240346035618411149408467364363417573258630000054836273193359375
            = 0.00000000000000000000000000000000000000000000000000000111001101001010110010100101111101100010001001101111b
    10^+-17 = 0.0000000000000000100000000000000007154242405462192450852805618492324772617063644020163337700068950653076171875
            = 0.0000000000000000000000000000000000000000000000000000000010111000011101111010101000110010001101101010010010111b
    10^+-18 = 0.00000000000000000100000000000000007154242405462192450852805618492324772617063644020163337700068950653076171875
            = 0.00000000000000000000000000000000000000000000000000000000000100100111001001011101110100011101001001000011101011b
    10^+-19 = 0.000000000000000000099999999999999997524592683526013185572915905567688179926555402943222361500374972820281982421875
            = 0.000000000000000000000000000000000000000000000000000000000000000111011000001111001001010011111011011011010010101011b
    10^+-20 = 0.00000000000000000000999999999999999945153271454209571651729503702787392447107715776066783064379706047475337982177734375
            = 0.00000000000000000000000000000000000000000000000000000000000000000010111100111001010000100001100100100100100001000100011b
    

    теперь отрицательные степени 10 напечатаны в стиле только соответствующих десятичных цифр (я к этому привык) для 64-битных doubles:

    10^+ -1 = 0.1
    10^+ -2 = 0.01
    10^+ -3 = 0.001
    10^+ -4 = 0.0001
    10^+ -5 = 0.00001
    10^+ -6 = 0.000001
    10^+ -7 = 0.0000001
    10^+ -8 = 0.00000001
    10^+ -9 = 0.000000001
    10^+-10 = 0.0000000001
    10^+-11 = 0.00000000001
    10^+-12 = 0.000000000001
    10^+-13 = 0.0000000000001
    10^+-14 = 0.00000000000001
    10^+-15 = 0.000000000000001
    10^+-16 = 0.0000000000000001
    10^+-17 = 0.00000000000000001
    10^+-18 = 0.000000000000000001
    10^+-19 = 0.0000000000000000001
    10^+-20 = 0.00000000000000000001
    

    надеюсь, это поможет:)

0 голосов
/ 09 июля 2010

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

Т.е.

10^2 = 2^6 + 2^5 + 2^2

Тогдасумма:

mantissa<<6 + mantissa<<5 + mantissa<<2

Я думаю, что разбить это будет на O (n) на количество цифр, сдвиг O (1), а суммирование O (n)цифры ...

Вам понадобится целочисленный класс, достаточно большой, чтобы хранить результаты, конечно ...

Дайте мне знать - мне любопытно, действительноя думаю.: -)

0 голосов
/ 09 июля 2010

Если вы хотите получить более точные результаты, почему бы не использовать математику с фиксированной точкой?Конверсии быстрые.Ошибка известна и может быть исправлена.Не точный ответ на ваш вопрос, но другая идея для вас.

0 голосов
/ 09 июля 2010

Ты не. Самое близкое, что вы можете найти, это сбросить байты.

...