एक सही ढंग से गोलाकार/लगभग सही ढंग से गोलाकार फ़्लोटिंग-पॉइंट क्यूबिक रूट

Question

कंप्यूटिंग का मानना ​​है कि CRlibm में पाए गए मानक लाइब्रेरी फ़ंक्शंस सही ढंग से गोल किए गए हैं। फिर एक डबल-परिशुद्धता इनपुट के सही गोलाकार क्यूबिक रूट की गणना कैसे करेगा?

^{यह प्रश्न अक्सर पूछे जाने वाले प्रश्नों को उद्धृत करने के लिए "वास्तविक समस्या [I] चेहरा" नहीं है। यह इस तरह होमवर्क की तरह थोड़ा सा है। लेकिन क्यूबिक रूट अक्सर एक बार पाया जाता है और कोई कल्पना कर सकता है कि यह प्रश्न वास्तविक समस्या है जिसे किसी का सामना करना पड़ता है।}

के बाद से, "सबसे अच्छा स्टैक ओवरफ़्लो सवाल उन में स्रोत कोड का एक सा है" यहाँ स्रोत कोड का एक सा है:

y = pow(x, 1./3.); 

ऊपर की गणना नहीं करता है एक सही ढंग से घन जड़ गोल क्योंकि 1/3 double के रूप में बिल्कुल प्रतिनिधित्व नहीं किया जा सकता है।

---

अतिरिक्त नोट्स:

एक article का वर्णन करता है एक फ्लोटिंग प्वाइंट घन जड़ गणना करने के लिए कैसे, लेकिन अनुशंसित न्यूटन- Raphson एल्गोरिथ्म के अंतिम यात्रा (रों) के लिए उच्च परिशुद्धता के लिए किया जाना करने के लिए होगा एक सही गोलाकार डबल-परिशुद्धता क्यूबिक रूट की गणना करने के लिए एल्गोरिदम। यह गणना करने का सबसे अच्छा तरीका हो सकता है, लेकिन मैं अभी भी एक शॉर्टकट की तलाश कर रहा हूं जो मौजूदा सही गोलाकार मानकीकृत कार्यों का लाभ उठाएगा।

सी 99 में cbrt() फ़ंक्शन शामिल है, लेकिन सभी कंपाइलरों के लिए or even faithful सही ढंग से गोल होने की उम्मीद नहीं की जा सकती है। सीआरएलआईबीएम के डिजाइनरों ने प्रदान किए गए कार्यों की सूची में cbrt() शामिल करना चुना होगा, लेकिन उन्होंने नहीं किया। सही गोलाकार गणित कार्यों के अन्य पुस्तकालयों में उपलब्ध कार्यान्वयन के संदर्भ स्वागत हैं।

Answer 1

मुझे डर है कि मुझे नहीं पता कि सही ढंग से गोलाकार डबल-परिशुद्धता घन रूट की गारंटी कैसे दी जाए, लेकिन सवाल में उल्लेखित लगभग एकदम सही ढंग से गोल किया जा सकता है। दूसरे शब्दों में, अधिकतम त्रुटि 0.5 ulp के बहुत करीब प्रतीत होती है।

पीटर Markstein, "IA-64 और प्राथमिक कार्य: स्पीड और प्रेसिजन" (प्रेंटिस-हॉल 2000)

सही ढंग से गोलाई पारस्परिक, वर्गमूल, और पारस्परिक वर्गमूल के लिए कुशल FMA आधारित तकनीकों प्रस्तुत करता है, लेकिन यह इस संबंध में घन रूट को कवर नहीं करता है। सामान्य रूप से मार्कस्टीन के दृष्टिकोण के लिए प्रारंभिक परिणाम की आवश्यकता होती है जो अंतिम दौर अनुक्रम से पहले 1 ulp के भीतर सटीक है।मेरे पास घन जड़ों के गोलाकार करने के लिए अपनी तकनीक का विस्तार करने के लिए गणितीय हैथल नहीं है, लेकिन ऐसा लगता है कि यह सिद्धांत रूप में संभव होना चाहिए, कुछ हद तक पारस्परिक वर्ग रूट के समान चुनौती है।

बिट-वार एल्गोरिदम सही गोल के साथ जड़ों की गणना के लिए आसानी से उधार देते हैं। चूंकि आईईईई -754 राउंडिंग मोड के लिए निकटतम-या-यहां तक ​​कि नहीं हो सकता है, इसलिए किसी को बस गणना करने की आवश्यकता होती है जब तक कि उसने सभी मंटिसा बिट्स और एक राउंड बिट का उत्पादन नहीं किया हो। बिनोमियल प्रमेय के आधार पर स्क्वायर रूट के लिए बिट-वार एल्गोरिदम, दोनों गैर-बहाली और पुनर्विक्रय दोनों प्रकारों में अच्छी तरह से जाना जाता है और हार्डवेयर कार्यान्वयन का आधार रहा है। द्विपद प्रमेय के माध्यम से एक ही दृष्टिकोण घनमूल के लिए काम करता है, और वहाँ एक अल्पज्ञात कागज कि एक गैर बहाल कार्यान्वयन का विवरण देता है है: वर्ग मूल और घन मूल निकाला जा रहा है के लिए

एच पेंग, "एल्गोरिदम , "कार्यवाही 5 वें आईईईई अंतर्राष्ट्रीय संगोष्ठी पर कंप्यूटर अंकगणित, पीपी 121-126, 1 9 81.

सर्वश्रेष्ठ मैं इसके साथ प्रयोग करने से कह सकता हूं यह पूर्णांक से घन जड़ों के निष्कर्षण के लिए पर्याप्त काम करता है। चूंकि प्रत्येक पुनरावृत्ति केवल एक ही परिणाम उत्पन्न करता है लेकिन यह बिल्कुल तेज़ नहीं है। फ्लोटिंग-पॉइंट अंकगणित में अनुप्रयोगों के लिए इसमें कुछ बहीखाता चर का उपयोग करने की कमी है जिसके लिए अंतिम परिणाम की बिट्स की संख्या लगभग दोगुनी होती है। इसका मतलब है कि एक डबल-परिशुद्धता घन रूट को लागू करने के लिए 128-बिट पूर्णांक अंकगणितीय का उपयोग करने की आवश्यकता होगी।

नीचे मेरा सी 99 कोड Halley's rational method for the cube root पर आधारित है जिसमें घन अभिसरण का अर्थ है कि शुरुआती अनुमान को प्रत्येक पुनरावृत्ति में मान्य अंकों के ट्रिपल की संख्या के रूप में बहुत सटीक नहीं होना चाहिए। गणना विभिन्न तरीकों से व्यवस्थित की जा सकती है। आम तौर पर यह संख्यानुसार फायदेमंद है के रूप में

new_guess पुनरावृत्ति योजनाओं की व्यवस्था करने के: = old_guess + सुधार

एक पर्याप्त करीब प्रारंभिक अनुमान के लिए के बाद से

, correctionold_guess की तुलना में काफी छोटा होता है। एक्स * (एक्स - - क)/(2 * x + y)

इस विशेष व्यवस्था है = एक्स:

एक्स: यह घनमूल के लिए निम्नलिखित यात्रा योजना की ओर जाता है Kahan's notes on cube root में भी सूचीबद्ध है। FMA (fused-multiply add) के उपयोग के लिए स्वाभाविक रूप से उधार देने का इसका और लाभ है। एक दोष यह है कि 2 * x की गणना अतिप्रवाह हो सकती है, इसलिए डबल-परिशुद्धता इनपुट डोमेन के कम से कम हिस्से के लिए एक तर्क कमी योजना आवश्यक है। मेरे कोड में मैं केवल पर गैर-असाधारण इनपुट पर तर्क कमी को लागू करता हूं, आईईईई -754 डबल-परिशुद्धता संचालन के एक्सपोनेंट्स के सीधे हेरफेर के आधार पर।

अंतराल [0.125,1) प्राथमिक अनुमान अंतराल के रूप में कार्य करता है। एक बहुपद मिनिमैक्स सन्निकटन का उपयोग किया जाता है जो प्रारंभिक अनुमान [0.5,1] में देता है। संकीर्ण सीमा गणना के निम्न-सटीकता भागों के लिए सिंगल-प्रेसिजन अंकगणितीय के उपयोग की सुविधा प्रदान करती है।

मैं अपने कार्यान्वयन की त्रुटि सीमाओं के बारे में कुछ भी साबित नहीं कर सकता, हालांकि, एक संदर्भ कार्यान्वयन (लगभग 200 बिट्स के लिए सटीक) के खिलाफ 200 मिलियन से अधिक यादृच्छिक परीक्षण वैक्टरों के साथ परीक्षण करने से पता चलता है कि आधा उलझन में कोई त्रुटि नहीं है, अधिकतम त्रुटि 0.5 उल के करीब होना चाहिए।

double my_cbrt (double a) 
{ 
    double b, u, v, r; 
    float bb, uu, vv; 
    int e, f, s; 

    if ((a == 0.0) || isinf(a) || isnan(a)) { 
     /* handle special cases */ 
     r = a + a; 
    } else { 
     /* strip off sign-bit */ 
     b = fabs (a); 
     /* compute exponent adjustments */ 
     b = frexp (b, &e); 
     s = e - 3*342; 
     f = s/3; 
     s = s - 3 * f; 
     f = f + 342; 
     /* map argument into the primary approximation interval [0.125,1) */ 
     b = ldexp (b, s); 
     bb = (float)b; 
     /* approximate cube root in [0.125,1) with relative error 5.22e-3 */ 
     uu =    0x1.2f32c0p-1f; 
     uu = fmaf (uu, bb, -0x1.62cc2ap+0f); 
     uu = fmaf (uu, bb, 0x1.7546e0p+0f); 
     uu = fmaf (uu, bb, 0x1.5d0590p-2f); 
     /* refine cube root using two Halley iterations w/ cubic convergence */ 
     vv = uu * uu; 
     uu = fmaf (fmaf (vv, uu, -bb)/fmaf (vv, 2.0f*uu, bb), -uu, uu); 
     u = (double)uu; 
     v = u * u; // this product is exact 
     r = fma (fma (v, u, -b)/fma (v, 2.0*u, b), -u, u); 
     /* map back from primary approximation interval by jamming exponent */ 
     r = ldexp (r, f); 
     /* restore sign bit */ 
     r = copysign (r, a); 
    } 
    return r; 
} 

Answer 2

यह देखते हुए कि वक्र x = y^3 पर बहुत आसानी से गणना करने योग्य तर्कसंगत बिंदु हैं, मैं एस^3 ~ एक्स के बारे में कम करने के लिए प्रेरित हूं, जिसमें तर्कसंगत और कुछ बिट्स चौड़े हैं। तो फिर तुम है:

cbrt(x) = s * cbrt(1 + (x - s^3)/s) 

स्पष्ट बात तो अपने पसंदीदा श्रृंखला सन्निकटन का उपयोग कर संशोधन पद का मूल्यांकन, और सिर पूंछ FMA गणित के माध्यम से एक अवशिष्ट की गणना एक ULP द्वारा परिणाम ऊपर या नीचे टक्कर करने के लिए यदि आवश्यक हो तो (यदि आप करने के लिए है जाहिर है, ज्यादातर समय पूर्ण गणना की आवश्यकता नहीं होगी)।

यह पूरी तरह से सवाल की भावना में नहीं है, लेकिन निश्चित रूप से काम करने के लिए बनाया जा सकता है, और इस तरह आवश्यक सीमाओं को साबित करना बहुत आसान है। उम्मीद है कि कोई और कुछ और चालाक सुझा सकता है (मैंने पहले से ही महीने के लिए अपनी चतुरता का उपयोग किया)।