इंटीजर क्यूब रूट

Question

मैं 64-बिट (हस्ताक्षरित) घन जड़ों के लिए तेज़ कोड ढूंढ रहा हूं। (मैं सी का उपयोग कर रहा हूं और जीसीसी के साथ संकलन कर रहा हूं, लेकिन मुझे लगता है कि अधिकांश काम आवश्यक भाषा-और कंपाइलर-अज्ञेयवादी होंगे।) मैं एक 64-बिट असीमित पूर्णांक के साथ इंगित करता हूं।

एक इनपुट को देखते हुए n मैं की आवश्यकता होती है ऐसी है कि

r * r * r <= n && n < (r + 1) * (r + 1) * (r + 1) 

है कि, मैं n का घनमूल चाहते हो (अभिन्न) वापसी मान आर, पूर्णांक। मूल कोड जैसे

return (ulong)pow(n, 1.0/3); 

सीमा के अंत की ओर घूमने के कारण गलत है।

ulong 
cuberoot(ulong n) 
{ 
    ulong ret = pow(n + 0.5, 1.0/3); 
    if (n < 100000000000001ULL) 
     return ret; 
    if (n >= 18446724184312856125ULL) 
     return 2642245ULL; 
    if (ret * ret * ret > n) { 
     ret--; 
     while (ret * ret * ret > n) 
      ret--; 
     return ret; 
    } 
    while ((ret + 1) * (ret + 1) * (ret + 1) <= n) 
     ret++; 
    return ret; 
} 

सही परिणाम देता है, लेकिन इसकी तुलना में धीमी गति से असंगत कोड है।

यह कोड गणित पुस्तकालय के लिए है और इसे विभिन्न कार्यों से कई बार बुलाया जाएगा। गति महत्वपूर्ण है, लेकिन आप गर्म कैश पर भरोसा नहीं कर सकते हैं (इसलिए 2,642,245-एंट्री बाइनरी खोज जैसे सुझाव सही हैं)।

तुलना के लिए, यहां कोड है सही पूर्णांक वर्ग रूट की गणना करता है।

ulong squareroot(ulong a) { 
    ulong x = (ulong)sqrt((double)a); 
    if (x > 0xFFFFFFFF || x*x > a) 
     x--; 
    return x; 
} 

Answer 1

"हैकर डिलाइट" पुस्तक में इस और कई अन्य समस्याओं के लिए एल्गोरिदम है। कोड ऑनलाइन here है। EDIT: वह कोड 64-बिट इनट्स के साथ ठीक से काम नहीं करता है, और 64-बिट के लिए इसे ठीक करने के तरीके पर पुस्तक में निर्देश कुछ उलझन में हैं। उचित 64-बिट कार्यान्वयन (परीक्षण केस सहित) ऑनलाइन here है।

मुझे शक है अपने squareroot समारोह "सही ढंग से" काम करता है कि - यह तर्क के लिए ulong a होना चाहिए, नहीं n :) (लेकिन एक ही दृष्टिकोण sqrt के बजाय cbrt का उपयोग कर काम करेगा, हालांकि सभी वेब सी गणित पुस्तकालयों घनमूल कार्य)।

Answer 2

मैं करूंगा research how to do it by hand, और उसके बाद का अनुवाद है कि एक ऐसे कंप्यूटर एल्गोरिदम में, बल्कि आधार से 10

हम (स्यूडोकोड) की तरह एक एल्गोरिथ्म कुछ के साथ अंत से आधार 2 में काम कर रहे:

Find the largest n such that (1 << 3n) < input. 
result = 1 << n. 
For i in (n-1)..0: 
    if ((result | 1 << i)**3) < input: 
     result |= 1 << i. 

हम (result | 1 << i)**3 की गणना को अनुकूलित करके, बिटटाइव-या समकक्ष के बराबर है, result**3 + 3 * i * result ** 2 + 3 * i ** 2 * result + i ** 3 पर रीफैक्टरिंग, result**3 और result**2 के मानों को कैश करके, और गुणा के बजाय शिफ्ट का उपयोग करके।

Answer 3

आप अपने राउंडिंग त्रुटियों को ठीक करने के लिए एक न्यूटन के कदम की कोशिश कर सकते:

ulong r = (ulong)pow(n, 1.0/3); 
if(r==0) return r; /* avoid divide by 0 later on */ 
ulong r3 = r*r*r; 
ulong slope = 3*r*r; 

ulong r1 = r+1; 
ulong r13 = r1*r1*r1; 

/* making sure to handle unsigned arithmetic correctly */ 
if(n >= r13) r+= (n - r3)/slope; 
if(n < r3) r-= (r3 - n)/slope; 

एक एकल न्यूटन कदम पर्याप्त होना चाहिए, लेकिन आप बंद-एक करके (या शायद अधिक?) त्रुटियों हो सकता है। आप देख सकते हैं/एक अंतिम जाँच & वेतन वृद्धि कदम का उपयोग करने वालों, अपने OQ में के रूप में ठीक:

while(r*r*r > n) --r; 
while((r+1)*(r+1)*(r+1) <= n) ++r; 

या कुछ इस तरह।

(मैं मानता मैं आलसी हूँ, यह ध्यान से निर्धारित करने के लिए जो (यदि जांच & वेतन वृद्धि में से कोई भी) वास्तव में आवश्यक है ... जाँच करने के लिए है करने के लिए सही रास्ता)

Answer 4

तो pow भी है महंगा, आप परिणाम के अनुमान लगाने के लिए एक गिनती-अग्रणी-शून्य निर्देश का उपयोग कर सकते हैं, फिर एक लुकअप टेबल का उपयोग कर सकते हैं, फिर इसे समाप्त करने के लिए कुछ न्यूटन चरण।

int k = __builtin_clz(n); // counts # of leading zeros (often a single assembly insn) 
int b = 64 - k;   // # of bits in n 
int top8 = n >> (b - 8); // top 8 bits of n (top bit is always 1) 
int approx = table[b][top8 & 0x7f]; 

को देखते हुए b और top8, आप (मेरे कोड, 8K प्रविष्टियों में) एक लुकअप तालिका का उपयोग करने के लिए cuberoot(n) एक अच्छा सन्निकटन खोजने के लिए कर सकते हैं। इसे समाप्त करने के लिए कुछ न्यूटन चरणों का उपयोग करें (आने वाले तूफान का जवाब देखें)।

Answer 5

// On my pc: Math.Sqrt 35 ns, cbrt64 <70ns, cbrt32 <25 ns, (cbrt12 < 10ns) 

// cbrt64(ulong x) is a C# version of: 
// http://www.hackersdelight.org/hdcodetxt/acbrt.c.txt  (acbrt1) 

// cbrt32(uint x) is a C# version of: 
// http://www.hackersdelight.org/hdcodetxt/icbrt.c.txt  (icbrt1) 

// Union in C#: 
// http://www.hanselman.com/blog/UnionsOrAnEquivalentInCSairamasTipOfTheDay.aspx 

using System.Runtime.InteropServices; 
[StructLayout(LayoutKind.Explicit)] 
public struct fu_32 // float <==> uint 
{ 
[FieldOffset(0)] 
public float f; 
[FieldOffset(0)] 
public uint u; 
} 

private static uint cbrt64(ulong x) 
{ 
    if (x >= 18446724184312856125) return 2642245; 
    float fx = (float)x; 
    fu_32 fu32 = new fu_32(); 
    fu32.f = fx; 
    uint uy = fu32.u/4; 
    uy += uy/4; 
    uy += uy/16; 
    uy += uy/256; 
    uy += 0x2a5137a0; 
    fu32.u = uy; 
    float fy = fu32.f; 
    fy = 0.33333333f * (fx/(fy * fy) + 2.0f * fy); 
    int y0 = (int)          
     (0.33333333f * (fx/(fy * fy) + 2.0f * fy));  
    uint y1 = (uint)y0;         

    ulong y2, y3; 
    if (y1 >= 2642245) 
    { 
     y1 = 2642245; 
     y2 = 6981458640025; 
     y3 = 18446724184312856125; 
    } 
    else 
    { 
     y2 = (ulong)y1 * y1; 
     y3 = y2 * y1; 
    } 
    if (y3 > x) 
    { 
     y1 -= 1; 
     y2 -= 2 * y1 + 1; 
     y3 -= 3 * y2 + 3 * y1 + 1; 
     while (y3 > x) 
     { 
      y1 -= 1; 
      y2 -= 2 * y1 + 1; 
      y3 -= 3 * y2 + 3 * y1 + 1; 
     } 
     return y1; 
    } 
    do 
    { 
     y3 += 3 * y2 + 3 * y1 + 1; 
     y2 += 2 * y1 + 1; 
     y1 += 1; 
    } 
    while (y3 <= x); 
    return y1 - 1; 
} 

private static uint cbrt32(uint x) 
{ 
    uint y = 0, z = 0, b = 0; 
    int s = x < 1u << 24 ? x < 1u << 12 ? x < 1u << 06 ? x < 1u << 03 ? 00 : 03 : 
                 x < 1u << 09 ? 06 : 09 : 
              x < 1u << 18 ? x < 1u << 15 ? 12 : 15 : 
                 x < 1u << 21 ? 18 : 21 : 
          x >= 1u << 30 ? 30 : x < 1u << 27 ? 24 : 27; 
    do 
    { 
     y *= 2; 
     z *= 4; 
     b = 3 * y + 3 * z + 1 << s; 
     if (x >= b) 
     { 
      x -= b; 
      z += 2 * y + 1; 
      y += 1; 
     } 
     s -= 3; 
    } 
    while (s >= 0); 
    return y; 
} 

private static uint cbrt12(uint x) // x < ~255 
{ 
    uint y = 0, a = 0, b = 1, c = 0; 
    while (a < x) 
    { 
     y++; 
     b += c; 
     a += b; 
     c += 6; 
    } 
    if (a != x) y--; 
    return y; 
}