Frage Schnellster Weg, um den ganzzahligen Teil von sqrt (n) zu erhalten?


Wie wir wissen, ob n ist also kein perfektes Quadrat sqrt(n) wäre keine ganze Zahl. Da ich nur den ganzzahligen Teil brauche, fühle ich das anzurufen sqrt(n) wäre nicht so schnell, da es Zeit braucht, um auch den Bruchteil zu berechnen.

Also meine Frage ist,

Können wir nur den ganzzahligen Teil von erhalten? sqrt (n) ohne den tatsächlichen Wert von zu berechnen sqrt(n)? Der Algorithmus sollte schneller sein als sqrt(n) (definiert in <math.h> oder <cmath>)

Wenn möglich, können Sie den Code schreiben asm blockiere auch.


59
2018-02-08 06:42


Ursprung


Antworten:


Ich würde es versuchen Schnelle Inverse Quadratwurzel Trick.

Es ist ein Weg, um eine sehr gute Annäherung davon zu bekommen 1/sqrt(n) ohne irgendeine Verzweigung, basierend auf etwas Bit-Twiddling, also nicht portierbar (insbesondere zwischen 32-Bit- und 64-Bit-Plattformen).

Sobald Sie es erhalten, müssen Sie nur das Ergebnis umkehren und nimmt den ganzzahligen Teil.

Es könnte natürlich schnellere Tricks geben, da dieser hier ein bisschen herumläuft.

BEARBEITEN: Machen wir das!

Zuerst ein kleiner Helfer:

// benchmark.h
#include <sys/time.h>

template <typename Func>
double benchmark(Func f, size_t iterations)
{
  f();

  timeval a, b;
  gettimeofday(&a, 0);
  for (; iterations --> 0;)
  {
    f();
  }
  gettimeofday(&b, 0);
  return (b.tv_sec * (unsigned int)1e6 + b.tv_usec) -
         (a.tv_sec * (unsigned int)1e6 + a.tv_usec);
}

Dann der Hauptteil:

#include <iostream>

#include <cmath>

#include "benchmark.h"

class Sqrt
{
public:
  Sqrt(int n): _number(n) {}

  int operator()() const
  {
    double d = _number;
    return static_cast<int>(std::sqrt(d) + 0.5);
  }

private:
  int _number;
};

// http://www.codecodex.com/wiki/Calculate_an_integer_square_root
class IntSqrt
{
public:
  IntSqrt(int n): _number(n) {}

  int operator()() const 
  {
    int remainder = _number;
    if (remainder < 0) { return 0; }

    int place = 1 <<(sizeof(int)*8 -2);

    while (place > remainder) { place /= 4; }

    int root = 0;
    while (place)
    {
      if (remainder >= root + place)
      {
        remainder -= root + place;
        root += place*2;
      }
      root /= 2;
      place /= 4;
    }
    return root;
  }

private:
  int _number;
};

// http://en.wikipedia.org/wiki/Fast_inverse_square_root
class FastSqrt
{
public:
  FastSqrt(int n): _number(n) {}

  int operator()() const
  {
    float number = _number;

    float x2 = number * 0.5F;
    float y = number;
    long i = *(long*)&y;
    //i = (long)0x5fe6ec85e7de30da - (i >> 1);
    i = 0x5f3759df - (i >> 1);
    y = *(float*)&i;

    y = y * (1.5F - (x2*y*y));
    y = y * (1.5F - (x2*y*y)); // let's be precise

    return static_cast<int>(1/y + 0.5f);
  }

private:
  int _number;
};


int main(int argc, char* argv[])
{
  if (argc != 3) {
    std::cerr << "Usage: %prog integer iterations\n";
    return 1;
  }

  int n = atoi(argv[1]);
  int it = atoi(argv[2]);

  assert(Sqrt(n)() == IntSqrt(n)() &&
          Sqrt(n)() == FastSqrt(n)() && "Different Roots!");
  std::cout << "sqrt(" << n << ") = " << Sqrt(n)() << "\n";

  double time = benchmark(Sqrt(n), it);
  double intTime = benchmark(IntSqrt(n), it);
  double fastTime = benchmark(FastSqrt(n), it);

  std::cout << "Number iterations: " << it << "\n"
               "Sqrt computation : " << time << "\n"
               "Int computation  : " << intTime << "\n"
               "Fast computation : " << fastTime << "\n";

  return 0;
}

Und die Ergebnisse:

sqrt(82) = 9
Number iterations: 4096
Sqrt computation : 56
Int computation  : 217
Fast computation : 119

// Note had to tweak the program here as Int here returns -1 :/
sqrt(2147483647) = 46341 // real answer sqrt(2 147 483 647) = 46 340.95
Number iterations: 4096
Sqrt computation : 57
Int computation  : 313
Fast computation : 119

Wo wie erwartet Schnell Berechnung führt viel besser aus als die Int Berechnung.

Oh, und nebenbei, sqrt ist schneller :)


20
2018-02-08 07:29



Edit: diese Antwort ist dumm - verwenden (int) sqrt(i)

Nach dem Profiling mit richtig die Einstellungen (-march=native -m64 -O3) Das obige war a Menge schneller.


Okay, eine etwas alte Frage, aber die "schnellste" Antwort wurde noch nicht gegeben. Der schnellste (ich glaube) Algorithmus ist der Binary Square Root Algorithmus, der vollständig in dieser Embedded.com Artikel.

Es kommt grundsätzlich darauf an:

unsigned short isqrt(unsigned long a) {
    unsigned long rem = 0;
    int root = 0;
    int i;

    for (i = 0; i < 16; i++) {
        root <<= 1;
        rem <<= 2;
        rem += a >> 30;
        a <<= 2;

        if (root < rem) {
            root++;
            rem -= root;
            root++;
        }
    }

    return (unsigned short) (root >> 1);
}

Auf meiner Maschine (Q6600, Ubuntu 10.10) profilierte ich, indem ich die Quadratwurzel der Zahlen 1-100000000 nahm. Verwenden iqsrt(i) nahm 2750 ms. Verwenden (unsigned short) sqrt((float) i) nahm 3600ms. Dies wurde unter Verwendung von g++ -O3. Verwendung der -ffast-math Kompilier-Option waren die Zeiten 2100ms bzw. 3100ms. Beachten Sie, dass dies ohne die Verwendung einer einzigen Assembler-Zeile ist, so dass es wahrscheinlich noch viel schneller sein könnte.

Der obige Code funktioniert sowohl für C und C ++ als auch für kleinere Syntaxänderungen auch für Java.

Was für eine begrenzte Reichweite noch besser funktioniert, ist eine binäre Suche. Auf meiner Maschine bläst das die Version aus dem Wasser um den Faktor 4. Leider ist die Reichweite sehr begrenzt:

#include <stdint.h>

const uint16_t squares[] = {
    0, 1, 4, 9,
    16, 25, 36, 49,
    64, 81, 100, 121,
    144, 169, 196, 225,
    256, 289, 324, 361,
    400, 441, 484, 529,
    576, 625, 676, 729,
    784, 841, 900, 961,
    1024, 1089, 1156, 1225,
    1296, 1369, 1444, 1521,
    1600, 1681, 1764, 1849,
    1936, 2025, 2116, 2209,
    2304, 2401, 2500, 2601,
    2704, 2809, 2916, 3025,
    3136, 3249, 3364, 3481,
    3600, 3721, 3844, 3969,
    4096, 4225, 4356, 4489,
    4624, 4761, 4900, 5041,
    5184, 5329, 5476, 5625,
    5776, 5929, 6084, 6241,
    6400, 6561, 6724, 6889,
    7056, 7225, 7396, 7569,
    7744, 7921, 8100, 8281,
    8464, 8649, 8836, 9025,
    9216, 9409, 9604, 9801,
    10000, 10201, 10404, 10609,
    10816, 11025, 11236, 11449,
    11664, 11881, 12100, 12321,
    12544, 12769, 12996, 13225,
    13456, 13689, 13924, 14161,
    14400, 14641, 14884, 15129,
    15376, 15625, 15876, 16129,
    16384, 16641, 16900, 17161,
    17424, 17689, 17956, 18225,
    18496, 18769, 19044, 19321,
    19600, 19881, 20164, 20449,
    20736, 21025, 21316, 21609,
    21904, 22201, 22500, 22801,
    23104, 23409, 23716, 24025,
    24336, 24649, 24964, 25281,
    25600, 25921, 26244, 26569,
    26896, 27225, 27556, 27889,
    28224, 28561, 28900, 29241,
    29584, 29929, 30276, 30625,
    30976, 31329, 31684, 32041,
    32400, 32761, 33124, 33489,
    33856, 34225, 34596, 34969,
    35344, 35721, 36100, 36481,
    36864, 37249, 37636, 38025,
    38416, 38809, 39204, 39601,
    40000, 40401, 40804, 41209,
    41616, 42025, 42436, 42849,
    43264, 43681, 44100, 44521,
    44944, 45369, 45796, 46225,
    46656, 47089, 47524, 47961,
    48400, 48841, 49284, 49729,
    50176, 50625, 51076, 51529,
    51984, 52441, 52900, 53361,
    53824, 54289, 54756, 55225,
    55696, 56169, 56644, 57121,
    57600, 58081, 58564, 59049,
    59536, 60025, 60516, 61009,
    61504, 62001, 62500, 63001,
    63504, 64009, 64516, 65025
};

inline int isqrt(uint16_t x) {
    const uint16_t *p = squares;

    if (p[128] <= x) p += 128;
    if (p[ 64] <= x) p +=  64;
    if (p[ 32] <= x) p +=  32;
    if (p[ 16] <= x) p +=  16;
    if (p[  8] <= x) p +=   8;
    if (p[  4] <= x) p +=   4;
    if (p[  2] <= x) p +=   2;
    if (p[  1] <= x) p +=   1;

    return p - squares;
}

Eine 32-Bit-Version kann hier heruntergeladen werden: https://gist.github.com/3481770


16
2018-03-14 09:17



Während ich vermute, dass Sie eine Menge von Optionen finden können, indem Sie nach "Fast Integer Square Root" suchen, hier sind einige potentiell neue Ideen, die gut funktionieren könnten (jeder unabhängig, oder vielleicht können Sie sie kombinieren):

  1. Mach ein static const Array aller perfekten Quadrate in der Domäne, die Sie unterstützen möchten, und führen Sie eine schnelle zweiglose binäre Suche darauf durch. Der resultierende Index im Array ist die Quadratwurzel.
  2. Konvertiere die Zahl in Fließkomma und zerlege sie in Mantisse und Exponent. Halbiere den Exponenten und multipliziere die Mantisse mit einem magischen Faktor (deine Aufgabe ist es). Dies sollte Ihnen eine sehr gute Annäherung geben können. Fügen Sie einen letzten Schritt ein, um es anzupassen, wenn es nicht genau ist (oder verwenden Sie es als Ausgangspunkt für die obige binäre Suche).

6
2018-02-08 06:49



Ich denke Google search bietet gute Artikel wie Calculate an integer square root über zu viele mögliche Wege der schnellen Berechnung diskutiert und es gibt gute Nachschlagewerke, ich denke, niemand kann hier besser bieten als sie (und wenn jemand kann zuerst Papier darüber produzieren), aber wenn Sie sie lesen und es gibt Unklarheiten mit Ihnen, dann können wir Ihnen gut helfen.


6
2018-02-08 07:09



Wenn Ihnen eine Annäherung nichts ausmacht, wie wäre es mit dieser ganzzahligen sqrt-Funktion, die ich zusammengeschustert habe?

int sqrti(int x)
{
    union { float f; int x; } v; 

    // convert to float
    v.f = (float)x;

    // fast aprox sqrt
    //  assumes float is in IEEE 754 single precision format 
    //  assumes int is 32 bits
    //  b = exponent bias
    //  m = number of mantissa bits
    v.x  -= 1 << 23; // subtract 2^m 
    v.x >>= 1;       // divide by 2
    v.x  += 1 << 29; // add ((b + 1) / 2) * 2^m

    // convert to int
    return (int)v.f;
}

Es verwendet den hier beschriebenen Algorithmus Wikipedia Artikel. Auf meinem Rechner ist es fast doppelt so schnell wie sqrt :)


5
2018-03-11 00:42



Um integer sqrt auszuführen, können Sie diese Spezialisierung der Newton-Methode verwenden:

Def isqrt(N):

    a = 1
    b = N

    while |a-b| > 1
        b = N / a
        a = (a + b) / 2

    return a

Grundsätzlich gilt für jedes x, dass das sqrt im Bereich (x ... N / x) liegt, also teilen wir dieses Intervall bei jeder Schleife für den neuen Schätzwert. So ähnlich wie binäre Suche, aber es konvergiert schneller.

Dies konvergiert in O (loglog (N)), was sehr schnell ist. Es verwendet auch überhaupt keinen Fließkommawert und es funktioniert auch gut für beliebige Ganzzahlen.


3
2017-08-26 17:55



Warum schlägt niemand die schnellste Methode vor?

Ob:

  1. Der Zahlenbereich ist begrenzt
  2. Speicherverbrauch ist nicht entscheidend
  3. Anwendungsstartzeit ist nicht kritisch

dann erstellen int[MAX_X] gefüllt (beim Start) mit sqrt(x) (Sie müssen die Funktion nicht verwenden sqrt() dafür).

All diese Bedingungen passen gut zu meinem Programm. Insbesondere int[10000000] Array wird konsumieren 40MB.

Was denkst du darüber?


3
2017-08-28 18:30



In vielen Fällen wird sogar ein exakter ganzzahliger sqrt-Wert nicht benötigt, genug, um eine gute Annäherung davon zu erhalten. (Zum Beispiel passiert es oft bei der DSP-Optimierung, wenn das 32-Bit-Signal auf 16 Bit oder 16 Bit auf 8 Bit komprimiert werden soll, ohne dass die Genauigkeit um den Nullpunkt herum verloren geht).

Ich habe diese nützliche Gleichung gefunden:

k = ceil(MSB(n)/2); - MSB(n) is the most significant bit of "n"


sqrt(n) ~= 2^(k-2)+(2^(k-1))*n/(2^(2*k))); - all multiplications and divisions here are very DSP-friendly, as they are only 2^k.

Diese Gleichung erzeugt eine glatte Kurve (n, sqrt (n)), ihre Werte unterscheiden sich nicht sehr von reellem sqrt (n) und können daher nützlich sein, wenn die ungefähre Genauigkeit ausreicht.


2
2017-08-08 12:36



Wenn Sie Leistung für die Berechnung der Quadratwurzel benötigen, werden Sie wahrscheinlich viele berechnen. Warum also nicht die Antwort zwischenspeichern? Ich kenne den Bereich für N in Ihrem Fall nicht und auch nicht, wenn Sie die Quadratwurzel derselben Ganzzahl mehrfach berechnen, aber wenn ja, können Sie das Ergebnis bei jedem Aufruf der Methode zwischenspeichern (in einem Array wäre das) der effizienteste, wenn nicht zu groß).


1
2018-02-08 08:24



Das ist so kurz, dass es 99% Inline gibt.

static inline int sqrtn(int num) {
    int i;
    __asm__ (
        "pxor %%xmm0, %%xmm0\n\t"   // clean xmm0 for cvtsi2ss
        "cvtsi2ss %1, %%xmm0\n\t"   // convert num to float, put it to xmm0
        "sqrtss %%xmm0, %%xmm0\n\t" // square root xmm0
        "cvttss2si %%xmm0, %0"      // float to int
        :"=r"(i):"r"(num):"%xmm0"); // i: result, num: input, xmm0: scratch register
    return i;
}

Warum sauber? xmm0? Dokumentation von cvtsi2ss

Der Zieloperand ist ein XMM-Register. Das Ergebnis wird im niedrigen Doppelwort des Zieloperanden gespeichert, und die oberen drei Doppelwörter bleiben unverändert.

GCC Intrinsic-Version (läuft nur auf GCC):

#include <xmmintrin.h>
int sqrtn2(int num) {
    register __v4sf xmm0 = {0, 0, 0, 0};
    xmm0 = __builtin_ia32_cvtsi2ss(xmm0, num);
    xmm0 = __builtin_ia32_sqrtss(xmm0);
    return __builtin_ia32_cvttss2si(xmm0);
}

Intel Intrinsic-Version (getestet auf GCC, Clang, ICC):

#include <xmmintrin.h>
int sqrtn2(int num) {
    register __m128 xmm0 = _mm_setzero_ps();
    xmm0 = _mm_cvt_si2ss(xmm0, num);
    xmm0 = _mm_sqrt_ss(xmm0);
    return _mm_cvt_ss2si(xmm0);
}

^^^^ Alle von ihnen erfordert SSE 1. (nicht einmal SSE 2)


1
2018-06-29 16:51



Auf meinem Computer mit gcc, mit -ffast-math, dauert das Konvertieren einer 32-Bit-Ganzzahl in float und das Verwenden von sqrtf 1,2 s pro 10 ^ 9 Ops (ohne -ffast-math dauert es 3,54 s).

Der folgende Algorithmus verwendet 0,87 s pro 10 ^ 9 auf Kosten einer gewissen Genauigkeit: Fehler können bis zu -7 oder +1 betragen, obwohl der RMS-Fehler nur 0,79 beträgt:

uint16_t SQRTTAB[65536];

inline uint16_t approxsqrt(uint32_t x) { 
  const uint32_t m1 = 0xff000000;
  const uint32_t m2 = 0x00ff0000;
  if (x&m1) {
    return SQRTTAB[x>>16];
  } else if (x&m2) {
    return SQRTTAB[x>>8]>>4;
  } else {
    return SQRTTAB[x]>>8;
  }
}

Die Tabelle besteht aus:

void maketable() {
  for (int x=0; x<65536; x++) {
    double v = x/65535.0;
    v = sqrt(v);
    int y = int(v*65535.0+0.999);
    SQRTTAB[x] = y;
  }
}

Ich fand, dass die Verfeinerung der Bisektion mit weiteren if-Anweisungen die Genauigkeit verbessert, aber es verlangsamt auch die Dinge bis zu dem Punkt, dass sqrtf schneller ist, zumindest mit -frast-math.


0
2018-02-10 15:38