Самый быстрый способ определить, является ли целочисленный квадратный корень целым числом

Question

Я ищу самый быстрый способ определить, является ли значение long идеальным квадратом (то есть его квадратный корень является другим целым числом):

1. Я сделал это простым способом, используя встроенный Math.sqrt() функции, но мне интересно, если есть способ сделать это быстрее, ограничить себя только целочисленным доменом.
2. Ведение справочной таблицы нецелесообразно (поскольку существует около 2 ^31,5 целые числа, площадь которых меньше 2 ⁶³ ).

Вот очень простой и понятный способ, которым я делаю это сейчас:

public final static boolean isPerfectSquare(long n)
{
  if (n < 0)
    return false;

  long tst = (long)(Math.sqrt(n) + 0.5);
  return tst*tst == n;
}

_{Примечание: я использую эту функцию во многих Project Euler задачах. Так что больше никому не придется поддерживать этот код. И этот вид микрооптимизации может реально изменить ситуацию, поскольку одна из задач состоит в том, чтобы выполнить каждый алгоритм менее чем за минуту, и в некоторых задачах эту функцию нужно будет вызывать миллионы раз.}

---

Я пробовал разные варианты решения проблемы:

• После исчерпывающего тестирования я обнаружил, что добавление 0.5 к результату Math.sqrt () необязательно, по крайней мере, на моей машине.
• Быстрый обратный квадратный корень был быстрее, но он дал неверные результаты для n> = 410881. Однако, как подсказывает BobbyShaftoe , мы можем использовать хак FISR для n < 410881.
• Метод Ньютона был немного медленнее, чем Math.sqrt(). Вероятно, это связано с тем, что Math.sqrt() использует что-то похожее на метод Ньютона, но реализовано в аппаратном обеспечении, поэтому оно намного быстрее, чем в Java. Кроме того, метод Ньютона все еще требовал использования двойных чисел.
• Модифицированный метод Ньютона, который использовал несколько трюков, так что использовалась только целочисленная математика, требовал некоторых хаков, чтобы избежать переполнения (я хочу, чтобы эта функция работала со всеми положительными 64-битными целыми числами со знаком), и он все еще был медленнее Math.sqrt().
• Бинарная отбивная была еще медленнее. Это имеет смысл, потому что двоичной отбивке в среднем потребуется 16 проходов, чтобы найти квадратный корень 64-битного числа.
• Согласно тестам Джона, использование or операторов в C ++ быстрее, чем использование switch, но в Java и C #, похоже, нет никакой разницы между or и switch.
• Я также попытался создать таблицу поиска (в качестве частного статического массива из 64 логических значений). Тогда вместо оператора switch или or я бы просто сказал if(lookup[(int)(n&0x3F)]) { test } else return false;. К моему удивлению, это было (немного) медленнее. Это связано с тем, что границы массивов проверяются в Java .

Answer 1

Целочисленная задача заслуживает целочисленного решения. Таким образом

Выполните бинарный поиск по (неотрицательным) целым числам, чтобы найти наибольшее целое число t, такое что t**2 <= n. Затем проверьте, является ли r**2 = n точно. Это занимает время O (log n).

Если вы не знаете, как выполнить двоичный поиск натуральных чисел, потому что множество не ограничено, это легко. Вы начинаете с вычисления вашей возрастающей функции f (выше f(t) = t**2 - n) по степеням два. Когда вы видите, что это становится положительным, вы нашли верхнюю границу. Тогда вы можете сделать стандартный бинарный поиск.

Answer 2

Было отмечено, что последние d цифры идеального квадрата могут принимать только определенные значения. Последние d цифры (в базе b) числа n совпадают с остатком, когда n делится на bd, т.е. в нотации C n % pow(b, d).

Это можно обобщить на любой модуль m, т.е. n % m можно использовать для исключения некоторого процента чисел из идеальных квадратов. Модуль, который вы используете в настоящее время, равен 64, что позволяет 12, т.е. 19% остатков, как возможные квадраты. С небольшим кодированием я нашел модуль 110880, который позволяет только 2016, т.е. 1,8% остатков в качестве возможных квадратов. Таким образом, в зависимости от стоимости операции модуля (т. Е. Деления) и поиска в таблице по сравнению с квадратным корнем на вашем компьютере, использование этого модуля может быть быстрее.

Кстати, если у Java есть способ хранить упакованный массив битов для таблицы поиска, не используйте его. В наши дни 110880 32-разрядных слов - это не много ОЗУ, и выбор машинного слова будет быстрее, чем выборка одного бита.

Answer 3

Для производительности вам очень часто приходится идти на некоторые компромиссы. Другие выражали различные методы, однако вы заметили, что хак Кармака был быстрее до определенных значений N. Затем вы должны проверить «n», и если оно меньше, чем число N, используйте хак Кармака, иначе используйте какой-то другой описанный метод. в ответах здесь.

Answer 4

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

• Мод-256 тест
• Неточный тест mod-3465 (исключает целочисленное деление за счет некоторых ложных срабатываний)
• Квадратный корень с плавающей точкой, округлить и сравнить с входным значением

Я также экспериментировал с этими модификациями, но они не помогли производительности:

• Дополнительный тест мод-255
• Деление входного значения на степени 4
• Быстрый обратный квадратный корень (для работы при больших значениях N требуется 3 итерации, что достаточно, чтобы сделать его медленнее, чем аппаратная функция квадратного корня.)

public class SquareTester {

    public static boolean isPerfectSquare(long n) {
        if (n < 0) {
            return false;
        } else {
            switch ((byte) n) {
            case -128: case -127: case -124: case -119: case -112:
            case -111: case -103: case  -95: case  -92: case  -87:
            case  -79: case  -71: case  -64: case  -63: case  -60:
            case  -55: case  -47: case  -39: case  -31: case  -28:
            case  -23: case  -15: case   -7: case    0: case    1:
            case    4: case    9: case   16: case   17: case   25:
            case   33: case   36: case   41: case   49: case   57:
            case   64: case   65: case   68: case   73: case   81:
            case   89: case   97: case  100: case  105: case  113:
            case  121:
                long i = (n * INV3465) >>> 52;
                if (! good3465[(int) i]) {
                    return false;
                } else {
                    long r = round(Math.sqrt(n));
                    return r*r == n; 
                }
            default:
                return false;
            }
        }
    }

    private static int round(double x) {
        return (int) Double.doubleToRawLongBits(x + (double) (1L << 52));
    }

    /** 3465<sup>-1</sup> modulo 2<sup>64</sup> */
    private static final long INV3465 = 0x8ffed161732e78b9L;

    private static final boolean[] good3465 =
        new boolean[0x1000];

    static {
        for (int r = 0; r < 3465; ++ r) {
            int i = (int) ((r * r * INV3465) >>> 52);
            good3465[i] = good3465[i+1] = true;
        }
    }

}

Answer 5

Следующее упрощение решения maaartinus, по-видимому, позволяет сократить время выполнения на несколько процентов, но я недостаточно хорош для сравнительного анализа, чтобы создать эталон, которому я могу доверять:

long goodMask; // 0xC840C04048404040 computed below
{
    for (int i=0; i<64; ++i) goodMask |= Long.MIN_VALUE >>> (i*i);
}

public boolean isSquare(long x) {
    // This tests if the 6 least significant bits are right.
    // Moving the to be tested bit to the highest position saves us masking.
    if (goodMask << x >= 0) return false;
    // Remove an even number of trailing zeros, leaving at most one.
    x >>= (Long.numberOfTrailingZeros(x) & (-2);
    // Repeat the test on the 6 least significant remaining bits.
    if (goodMask << x >= 0 | x <= 0) return x == 0;
    // Do it in the classical way.
    // The correctness is not trivial as the conversion from long to double is lossy!
    final long tst = (long) Math.sqrt(x);
    return tst * tst == x;
}

Стоит проверить, как пропустить первый тест,

if (goodMask << x >= 0) return false;

повлияет на производительность.

Answer 6

Вы должны избавиться от 2-степенной части N с самого начала.

2-е редактирование Волшебное выражение для м ниже должно быть

m = N - (N & (N-1));

а не как написано

Конец второго редактирования

m = N & (N-1); // the lawest bit of N
N /= m;
byte = N & 0x0F;
if ((m % 2) || (byte !=1 && byte !=9))
  return false;

1-е редактирование:

Незначительное улучшение:

m = N & (N-1); // the lawest bit of N
N /= m;
if ((m % 2) || (N & 0x07 != 1))
  return false;

Конец первого редактирования

Теперь продолжайте как обычно. Таким образом, к тому времени, когда вы доберетесь до части с плавающей запятой, вы уже избавились от всех чисел, чья 2-степенная часть нечетна (примерно половина), и тогда вы будете считать только 1/8 того, что осталось. То есть вы запускаете часть с плавающей запятой на 6% чисел.

Answer 7

Мне нравится идея использовать почти правильный метод для некоторых входных данных. Вот версия с более высоким «смещением». Код, кажется, работает и проходит мой простой тестовый пример.

Просто замените ваш:

if(n < 410881L){...}

код с этим:

if (n < 11043908100L) {
    //John Carmack hack, converted to Java.
    // See: http://www.codemaestro.com/reviews/9
    int i;
    float x2, y;

    x2 = n * 0.5F;
    y = n;
    i = Float.floatToRawIntBits(y);
    //using the magic number from 
    //http://www.lomont.org/Math/Papers/2003/InvSqrt.pdf
    //since it more accurate
    i = 0x5f375a86 - (i >> 1);
    y = Float.intBitsToFloat(i);
    y = y * (1.5F - (x2 * y * y));
    y = y * (1.5F - (x2 * y * y)); //Newton iteration, more accurate

    sqrt = Math.round(1.0F / y);
} else {
    //Carmack hack gives incorrect answer for n >= 11043908100.
    sqrt = (long) Math.sqrt(n);
}

Answer 8

Project Euler упоминается в тегах, и многие из проблем в нем требуют проверки номера >> 2 ^ 64. Большинство из упомянутых выше оптимизаций не работают легко, когда вы работаете с 80-байтовым буфером.

Я использовал java BigInteger и слегка модифицированную версию метода Ньютона, которая лучше работает с целыми числами. Проблема заключалась в том, что точные квадраты n ^ 2 сходились к (n-1) вместо n, потому что n ^ 2-1 = (n-1) (n + 1), и окончательная ошибка была всего на один шаг ниже конечного делителя алгоритм прекращен. Это было легко исправить, добавив один к исходному аргументу перед вычислением ошибки. (Добавьте два для корней куба и т. Д.)

Одним приятным атрибутом этого алгоритма является то, что вы можете сразу сказать, является ли число идеальным квадратом - конечная ошибка (не исправление) в методе Ньютона будет равна нулю. Простая модификация также позволяет вам быстро вычислить floor (sqrt (x)) вместо ближайшего целого числа. Это удобно при нескольких проблемах Эйлера.

Answer 9

Это доработка от десятичного к двоичному алгоритму старого калькулятора Марчанта (извините, у меня нет ссылки) в Ruby, адаптированном специально для этого вопроса:

def isexactsqrt(v)
    value = v.abs
    residue = value
    root = 0
    onebit = 1
    onebit <<= 8 while (onebit < residue)
    onebit >>= 2 while (onebit > residue)
    while (onebit > 0)
        x = root + onebit
        if (residue >= x) then
            residue -= x
            root = x + onebit
        end
        root >>= 1
        onebit >>= 2
    end
    return (residue == 0)
end

Вот пример чего-то подобного (пожалуйста, не голосуйте за стиль кодирования / запахи или неуклюжий ввод-вывод - это алгоритм, который имеет значение, а C ++ не мой родной язык) В этом случае мы ищем остаток == 0:

#include <iostream>  

using namespace std;  
typedef unsigned long long int llint;

class ISqrt {           // Integer Square Root
    llint value;        // Integer whose square root is required
    llint root;         // Result: floor(sqrt(value))
    llint residue;      // Result: value-root*root
    llint onebit, x;    // Working bit, working value

public:

    ISqrt(llint v = 2) {    // Constructor
        Root(v);            // Take the root 
    };

    llint Root(llint r) {   // Resets and calculates new square root
        value = r;          // Store input
        residue = value;    // Initialise for subtracting down
        root = 0;           // Clear root accumulator

        onebit = 1;                 // Calculate start value of counter
        onebit <<= (8*sizeof(llint)-2);         // Set up counter bit as greatest odd power of 2 
        while (onebit > residue) {onebit >>= 2; };  // Shift down until just < value

        while (onebit > 0) {
            x = root ^ onebit;          // Will check root+1bit (root bit corresponding to onebit is always zero)
            if (residue >= x) {         // Room to subtract?
                residue -= x;           // Yes - deduct from residue
                root = x + onebit;      // and step root
            };
            root >>= 1;
            onebit >>= 2;
        };
        return root;                    
    };
    llint Residue() {           // Returns residue from last calculation
        return residue;                 
    };
};

int main() {
    llint big, i, q, r, v, delta;
    big = 0; big = (big-1);         // Kludge for "big number"
    ISqrt b;                            // Make q sqrt generator
    for ( i = big; i > 0 ; i /= 7 ) {   // for several numbers
        q = b.Root(i);                  // Get the square root
        r = b.Residue();                // Get the residue
        v = q*q+r;                      // Recalc original value
        delta = v-i;                    // And diff, hopefully 0
        cout << i << ": " << q << " ++ " << r << " V: " << v << " Delta: " << delta << "\n";
    };
    return 0;
};

Answer 10

Вызов sqrt, как уже упоминалось, не совсем точен, но интересно и поучительно, что он не уносит другие ответы с точки зрения скорости. В конце концов, последовательность инструкций языка ассемблера для sqrt крошечная. У Intel есть инструкция по аппаратному обеспечению, которая не используется Java, я полагаю, потому что она не соответствует IEEE.

Так почему же это медленно? Потому что Java на самом деле вызывает подпрограмму C через JNI, и это на самом деле медленнее, чем вызов подпрограммы Java, которая сама по себе медленнее, чем встроенная. Это очень раздражает, и Java должна была придумать лучшее решение, то есть, при необходимости, создание вызовов библиотеки с плавающей запятой. Ну хорошо.

Я подозреваю, что в C ++ все сложные альтернативы будут терять скорость, но я не проверял их все. То, что я сделал, и что люди Java найдут полезными, - это простой взлом, расширение тестирования специального случая, предложенного А. Рексом. Используйте одно длинное значение в качестве битового массива, который не проверяется по границам. Таким образом, у вас есть 64-битный логический поиск.

typedef unsigned long long UVLONG
UVLONG pp1,pp2;

void init2() {
  for (int i = 0; i < 64; i++) {
    for (int j = 0; j < 64; j++)
      if (isPerfectSquare(i * 64 + j)) {
    pp1 |= (1 << j);
    pp2 |= (1 << i);
    break;
      }
   }
   cout << "pp1=" << pp1 << "," << pp2 << "\n";  
}


inline bool isPerfectSquare5(UVLONG x) {
  return pp1 & (1 << (x & 0x3F)) ? isPerfectSquare(x) : false;
}

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

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

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

Процедура init2 вызывается один раз для инициализации статических значений pp1 и pp2. Обратите внимание, что в моей реализации на C ++ я использую unsigned long long, поэтому, поскольку вы подписаны, вам придется использовать оператор >>>.

Нет внутренней необходимости проверять границы массива, но оптимизатор Java должен довольно быстро разобраться с этим, поэтому я не виню их за это.