Блог

Мне нужно проверить, лежит ли точка в R ^ 2 в окружности относительно большого радиуса r (до 10 ^ 5). Очевидно, что обычно я бы просто сравнил внутреннее произведение с r ^ 2, но это происходит во встроенной среде, и это не будет работать со значениями int32_t, которые достаточно велики, поскольку квадратуры будут переполнять тип (максимум 32-битные типы).

Возможные решения:

Я мог бы вручную создать 64-битный продукт из двух 32-битных целых чисел (вероятно, то, что я в конечном итоге сделаю).

Я мог бы разделить все на 10 (или любое значение), а затем выполнить обычное сравнение внутреннего продукта, но я теряю точность.

Я мог бы попытаться проверить внутри n-угольника, вписанного в круг, но это требует много вычислений, таблиц и т. Д., И я все еще теряю точность.

Существует ли алгоритм, который обычно используется для подобных вещей?

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

 int check_distance(int x, int y, int r) {
    return (long long)x * x   (long long)y * y <= (long long)r * r;
}
 

Если сгенерированный код кажется слишком медленным, вы можете добавить тест, чтобы проверить, требуется ли 64-разрядная операция. Предполагая x , y и r являются положительными, вот решение, использующее арифметику без знака и точные типы ширины из <stdint.h> :

 int check_distance(uint32_t x, uint32_t y, uint32_t r) {
    if (x <= 46340 amp;amp; y <= 46340 amp;amp; r <= 0xffff) {
        /* 32-bit unsigned expression does not overflow */
        return x * x   y * y <= r * r;
    } else {
        return (uint64_t)x * x   (uint64_t)y * y <= (uint64_t)r * r;
    }
}
 

Обратите внимание на константу 46340, которая равна floor(sqrt(pow(2, 31))) : если оба x и y больше этого значения, x*x y*y будет превышать 2 ³².

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

 int check_distance(uint32_t x, uint32_t y, uint32_t r) {
    if ((x | y | r) <= 0x7fff) {
        /* 32-bit unsigned expression does not overflow */
        return x * x   y * y <= r * r;
    } else {
        return (uint64_t)x * x   (uint64_t)y * y <= (uint64_t)r * r;
    }
}
 

Затем, если вы действительно не хотите использовать 64-разрядную арифметику компилятора, вы можете написать вычисление явно. Учитывая диапазон значений x y и r указанный как <= 100000 , сдвиг значений вправо на 2 бита сохраняет x и y ниже 46340:

 int check_distance(uint32_t x, uint32_t y, uint32_t r) {
    if (x <= 46340 amp;amp; y1 <= 46340 amp;amp; r1 <= 0xffff) {
        /* 32-bit unsigned expression does not overflow */
        return x * x   y * y <= r * r;
    } else {
        /* shift all values right 2 bits to keep them below 46340 */
        uint32_t x0 = x amp; 3;
        uint32_t y0 = y amp; 3;
        uint32_t r0 = r amp; 3;
        uint32_t x1 = x >> 2;
        uint32_t y1 = y >> 2;
        uint32_t r1 = r >> 2;
        uint32_t x2_lo = x0 * (x0   x1 * 8);
        uint32_t y2_lo = y0 * (y0   y1 * 8);
        uint32_t d2_lo = x2_lo   y2_lo;
        uint32_t d2_hi = x1 * x1   y1 * y1   (d2_lo >> 4);
        uint32_t r2_lo = r0 * (r0   r1 * 8);
        uint32_t r2_hi = r1 * r1   (r2_lo >> 4);
        return d2_hi < r2_hi || (d2_hi == r2_hi amp;amp; (d2_lo amp; 15) <= (r2_lo amp; 15));
    }
}
 

Наконец, сдвиг значений на 5 бит позволяет получать числа до 1000000:

 int check_distance(uint32_t x, uint32_t y, uint32_t r) {
    if (x <= 46340 amp;amp; y1 <= 46340 amp;amp; r1 <= 0xffff) {
        /* 32-bit unsigned expression does not overflow */
        return x * x   y * y <= r * r;
    } else {
        /* shift all values right 5 bits to keep them below 46340 */
        uint32_t x0 = x amp; 31;
        uint32_t y0 = y amp; 31;
        uint32_t r0 = r amp; 31;
        uint32_t x1 = x >> 5;
        uint32_t y1 = y >> 5;
        uint32_t r1 = r >> 5;
        uint32_t x2_lo = x0 * (x0   x1 * 64);
        uint32_t y2_lo = y0 * (y0   y1 * 64);
        uint32_t d2_lo = x2_lo   y2_lo;
        uint32_t d2_hi = x1 * x1   y1 * y1   (d2_lo >> 10);
        uint32_t r2_lo = r0 * (r0   r1 * 64);
        uint32_t r2_hi = r1 * r1   (r2_lo >> 10);
        return d2_hi < r2_hi || (d2_hi == r2_hi amp;amp; (d2_lo amp; 1023) <= (r2_lo amp; 1023));
    }
}
 

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

 int check_distance(uint32_t x, uint32_t y, uint32_t r) {
    while (x > 46340 || y > 46340 || r > 0xffff) {
        x >>= 1;
        y >>= 1;
        r >>= 1;
    }
    /* 32-bit unsigned expression no longer overflows */
    return x * x   y * y <= r * r;
}
 

Максимальный радиус, который вы можете иметь, используя описанные вами предпосылки 100000 , равен , что приводит к квадрату модуля 10 000 000 000 , для которого требуется (поскольку все числа положительные) 35 целое число бит ( unsigned квадратного радиуса), которое должно быть представлено.

Основываясь на этих предпосылках и на том факте, что у вас нет простого способа получить 64-битные целые числа, и имея довольно низкие дополнительные биты, мы можем масштабировать результаты на четыре бита в норме (два бита в исходных координатах), чтобы достичь полной емкости для обработки 100000 координат до 32-битного целого числа без знака.

В моем первом издании этого ответа я предположил, что для обработки полного набора значений было достаточно только одного сдвига в координатах (два бита в вычисленной норме), и учитывалась потеря 1 бита точности, но я ошибался, и потребовался один дополнительный бит. Необходимо сдвинуть результаты как минимум на три бита, чтобы вместить полный набор входных данных, поэтому я решил масштабировать координаты на два бита, и поэтому результаты будут масштабироваться на четыре. Поскольку я решил также всегда масштабировать и возвращать долю квадратной единицы в виде целого числа в диапазоне от 0 до 15 (в шестнадцатых долях квадратной единицы). Таким образом, вы добьетесь точных результатов, сравнив сначала целые части двух точек и используя дробные части, в случае совпадения целых частей. Это делает вычисления и значение возвращаемых результатов более согласованными, чем ранее, и дает вам полную точность с целочисленными координатами.

Вы запросили рабочую реализацию, поэтому я разместил ее для вас ниже:

 #include <stdio.h>
#include <stdint.h>

/* calculate the square of a divided by four number and
 * accumulate the fraction (in sixteenths of a square unit)
 * into the reference pointed by frac_p.  */
uint32_t
square_of_div16(uint32_t x, int *frac_p)
{
    /* we use (IP   FP)^2 = IP^2   2*IP*FP   FP^2 */

    uint32_t int_part    = x >>  2;                /* divide by four */
    uint32_t frac_part   = x amp; 0x3;                /* mod 4 */
    uint32_t int_result  = int_part  * int_part;   /* square of IP */
    int      frac_result = frac_part * frac_part;  /* square of FP */
    uint32_t mixed_prod  = int_part  * frac_part;  /* IP*FP */

    int_result   =  mixed_prod >> 1;
    frac_result  = (mixed_prod amp; 1) << 3;
    if (frac_result >= 0x10) { /* carry process */
        int_result   = frac_result >> 4;
        frac_result amp;= 0x0f;
    }
    if (frac_p) *frac_p  = frac_resu< /* accumulate */
    return int_resu<
}

/* this calculates the squared norm scaled to one sixteenth
 * of the original coordinates (scaled by one fourth).
 * The ref_fraction pointer is a reference of a variable to
 * accumulate the fraction sixteenths of a square unit.  If
 * you are not interested in the fraction value, you can just
 * pass NULL as parameter. */
uint32_t
norm_scaled(uint32_t x, uint32_t y, int *ref_fraction)
{
    int fraction = 0;
    uint32_t result = 0;

    result  = square_of_div16(x, amp;fraction);
    result  = square_of_div16(y, amp;fraction);

    if (ref_fraction)
        *ref_fraction  = fraction; /* the excess */

    return resu<
}

/* TEST MAIN PROGRAM.  Just input pairs of coordinates in the
 * same line (separated by spaces) and calculate the squared
 * norm of the vector, scaled by 1/16 (accumulating the
 * fraction of the value in 1/16s of a square unit in the
 * location referenced.  This is done using double floating
 * point numbers and uint32_t integers. */
int main()
{
    char line[256];
    while (fgets(line, sizeof line, stdin) != NULL) {
        int x = 0, y = 0, fraction = 0;

        sscanf(line, "%u%u", amp;x, amp;y);

        uint32_t norm_16th = norm_scaled(x, y, amp;fraction);

        printf("Trying (%u, %u) => %u (fraction = %d/16)n",
                x, y, norm_16th, fraction);

        double norm_sq_16th
            = (double) x/4.0 * (double)x/4.0
              (double) y/4.0 * (double)y/4.0;

        printf("squared norm scaled: %.8fn", norm_sq_16th);
    }
    printf("Program finishedn");
}
 

Функция square_of_div16 вычисляет масштабированный модуль, деленный на 16 числа, поэтому мы можем использовать его для вычисления квадратов x и y координат. Функция принимает указатель frac_p на целочисленную переменную для хранения дробной части (в шестнадцатых долях квадратной единицы)

Затем функция norm_scaled вычисляет норму, используя square_of_div16 функцию и добавляя оба результата. Дробная часть накапливается для обоих вызовов, а результат накапливается в указанной переменной по указателю ref_fraction . Здесь выполняется обработка переноса, чтобы дать правильные результаты.

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