Какой первый дубль отклоняется от соответствующего длинного на дельту?

Question

Я хочу знать первый дубль от 0d и выше, который отклоняется по длине «того же значения» на некоторую дельту, скажем, 1e-8. Я терплю неудачу здесь, хотя Я пытаюсь сделать это в C, хотя я обычно использую управляемые языки, на всякий случай. Пожалуйста, помогите.

#include <stdio.h>
#include <limits.h>
#define DELTA 1e-8

int main() {
    double d = 0; // checked, the literal is fine
    long i;
    for (i = 0L; i < LONG_MAX; i++) {
         d=i; // gcc does the cast right, i checked
         if (d-i > DELTA || d-i < -DELTA) {
              printf("%f", d);
              break;
         }
    }
}

Я предполагаю, что проблема в том, что я делаю i на удвоение и, следовательно, d == i, а затем разница всегда равна 0. Как еще я могу это правильно определить - я бы предпочел забавное приведение C вместо сравнения строк , что заняло бы вечность.

ОТВЕТ : именно так, как мы ожидали. 2 ^ 53 + 1 = 9007199254740993 - это первая разница в соответствии со стандартными инструментами C / UNIX / POSIX. Большое спасибо pax за его программу. И я думаю, математика снова побеждает.

Answer 1

Двойные числа в IEE754 имеют точность 52 бита, что означает, что они могут хранить числа с точностью до (как минимум) 2 ⁵¹ .

Если ваши длинные 32-битные, они будут иметь только (положительный) диапазон от 0 до 2 ³¹ , поэтому нет 32-битных длинных, которые не могут быть представлены в точности как двойные. Для 64-битной длины это будет (примерно) 2 ⁵² , поэтому я бы начал там, а не с нуля.

Вы можете использовать следующую программу, чтобы определить, где начинаются сбои. В более ранней версии я опирался на тот факт, что последняя цифра в числе, которое непрерывно удваивается, следует за последовательностью {2,4,8,6}. Однако в конечном итоге я решил использовать известный доверенный инструмент (bc) для проверки всего числа, а не только последней цифры.

Имейте в виду, что это может быть затронуто действиями sprintf(), а не реальной точностью двойных чисел (я не думаю, что лично, поскольку у него не было проблем с определенными числами до 2 ¹⁴³ ).

Это программа:

#include <stdio.h>
#include <string.h>

int main() {
    FILE *fin;
    double d = 1.0; // 2^n-1 to avoid exact powers of 2.
    int i = 1;
    char ds[1000];
    char tst[1000];

    // Loop forever, rely on break to finish.
    while (1) {
        // Get C version of the double.
        sprintf (ds, "%.0f", d);

        // Get bc version of the double.
        sprintf (tst, "echo '2^%d - 1' | bc >tmpfile", i);
        system(tst);
        fin = fopen ("tmpfile", "r");
        fgets (tst, sizeof (tst), fin);
        fclose (fin);
        tst[strlen (tst) - 1] = '\0';

        // Check them.
        if (strcmp (ds, tst) != 0) {
            printf( "2^%d - 1 <-- bc failure\n", i);
            printf( "   got       [%s]\n", ds);
            printf( "   expected  [%s]\n", tst);
            break;
        }

        // Output for status then move to next.
        printf( "2^%d - 1 = %s\n", i, ds);
        d = (d + 1) * 2 - 1;  // Again, 2^n - 1.
        i++;
    }
}

Это продолжается до:

2^51 - 1 = 2251799813685247
2^52 - 1 = 4503599627370495
2^53 - 1 = 9007199254740991
2^54 - 1 <-- bc failure
   got       [18014398509481984]
   expected  [18014398509481983]

, где я ожидал, что это не получится.

В качестве отступления я первоначально использовал числа в форме 2 ⁿ , но я получил:

2^136 = 87112285931760246646623899502532662132736
2^137 = 174224571863520493293247799005065324265472
2^138 = 348449143727040986586495598010130648530944
2^139 = 696898287454081973172991196020261297061888
2^140 = 1393796574908163946345982392040522594123776
2^141 = 2787593149816327892691964784081045188247552
2^142 = 5575186299632655785383929568162090376495104
2^143 <-- bc failure
   got       [11150372599265311570767859136324180752990210]
   expected  [11150372599265311570767859136324180752990208]

с двойным размером 8 байт (проверяется с помощью sizeof). Оказалось, что эти числа имеют двоичную форму "1000...", которую можно представить гораздо дольше двойными. Именно тогда я переключился на использование 2 ⁿ -1, чтобы получить лучший битовый шаблон: все один бит.

Answer 2

Первое длинное «неправильно», когда приведение к двойнику не будет выключено на 1e-8, оно будет выключено на 1. До тех пор, пока двойник может соответствовать длинному в его значении, он будет представлять его точно .

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

В Википедии написано 52 бита в значении, не считая неявного начала 1. Это должно означать, что первый длинный тип, который будет приведен к другому значению, равен 2 ^ 53.

Answer 3

Хотя я не решаюсь упомянуть Fortran 95 и его преемников в этом обсуждении, я упомяну, что Fortran, начиная со стандарта 1990 года, предлагает встроенную функцию SPACING, которая сообщает вам, в чем разница между представимыми REAL относительно данного REAL. Вы можете выполнить бинарный поиск, останавливаясь, когда SPACING (X)> DELTA. Для компиляторов, использующих ту же модель с плавающей запятой, что и интересующую вас (вероятно, стандарт IEEE754), вы должны получить те же результаты.

Answer 4

Я подумал, что двойные числа могут точно представлять все целые числа (в пределах их границ).

Если это не так, то вы захотите привести и i, и d к чему-то с большей точностью, чем любой из них. Возможно, длинный двойной подойдет.