कहन एल्गोरिथम: उत्पादों का सटीक अंतर कैसे प्राप्त करें

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

इस प्रक्रिया में, मैं StackOverflow पर इस पोस्ट पर आया, जिससे मैंने सटीक गणना के लिए एक सुरुचिपूर्ण एल्गोरिथ्म के बारे में सीखा $$$a \ b b-c \ d d$।

लेकिन एल्गोरिथ्म के साथ आगे बढ़ने से पहले, आपको यह समझने की आवश्यकता है कि अभिव्यक्ति में इतनी चालाक क्या है $$$a \ b b-c \ d d$? लेना $$$a = 33962.035$। $$$b = -30438.8$। $$$c = 41563.4$और $$$d = -24871.969$। (ये असली मूल्य हैं जो मुझे pbrt के लॉन्च के दौरान मिले थे ।) 32-बिट फ्लोट वैल्यू के लिए, हमें यह मिलता है: $$$a \ बार b = -1.03376365 \ गुना 10 ^ 9$और $$$c \ टाइम्स d = -1.03376352 \ 10 बार ^ 9$। हम घटाव को पूरा करते हैं, और हम प्राप्त करते हैं $$$-128$। लेकिन अगर आप दोहरी सटीकता के साथ गणना करते हैं, और अंत में उन्हें फ्लोट में परिवर्तित करते हैं, तो आप प्राप्त करते हैं $$$-75.1656$। क्या हुआ?

समस्या यह है कि प्रत्येक कार्य का मूल्य नीचे की रेखा से बहुत आगे जा सकता है $$$-1 \ _ 10 ^ 9$, जहां प्रतिनिधित्व योग्य अस्थायी बिंदु मानों के बीच की दूरी बहुत बड़ी है - 64. जब गोलाई है $$$a \ _ b$और $$$c \ _ d$व्यक्तिगत रूप से, निकटतम प्रतिनिधित्व योग्य फ्लोट में, वे संख्याओं में बदल जाते हैं जो 64 के गुणक हैं। बदले में, उनका अंतर 64 का गुणक होगा, और कोई उम्मीद नहीं होगी कि यह बन जाएगा $$$-75.1656$से अधिक निकट $$$-64$। हमारे मामले में, यह परिणाम और भी अधिक था कि दोनों कार्यों को कैसे पूरा किया गया $$$-1 \ _ 10 ^ 9$। हम सीधे अच्छे पुराने विनाशकारी कमी ^{1 का} सामना करेंगे।

यहाँ ² से बेहतर समाधान है:

inline float DifferenceOfProducts(float a, float b, float c, float d) { float cd = c * d; float err = std::fma(-c, d, cd); float dop = std::fma(a, b, -cd); return dop + err; } 

DifferenceOfProducts() गणना करता है $$$a \ b b-c \ d d$एक तरह से जो भयावह संकुचन से बचा जाता है। इस तकनीक का वर्णन पहली बार प्रसिद्ध विलियम कहन ने लेख ऑन द कॉस्ट ऑफ फ्लोटिंग-पॉइंटकंप्यूटेशन विदाउट एक्सट्रा-प्रीसीज़ अरिथमेटिक में किया था । यह ध्यान देने योग्य है कि कहन का काम समग्र रूप से पढ़ना दिलचस्प है, उनके पास फ्लोटिंग पॉइंट्स की दुनिया की वर्तमान स्थिति के साथ-साथ गणितीय और तकनीकी विचारों पर कई टिप्पणियां हैं। यहाँ उनके निष्कर्षों में से एक है:

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

उनकी बुद्धि को श्रद्धांजलि देने के बाद, हम डिफरेंसऑफ DifferenceOfProducts() लौटते हैं: इस फ़ंक्शन की महारत का आधार फ्यूज्ड मल्टी-ऐड, एफएमए ³ निर्देशों का उपयोग है। गणितीय दृष्टिकोण से, FMA(a,b,c) है $$$a + बार b + c$इसलिए, पहली बार में ऐसा लगता है कि यह ऑपरेशन केवल माइक्रोपॉटिमाइजेशन के रूप में उपयोगी है: दो के बजाय एक निर्देश। हालाँकि एफ.एम.
इसकी एक विशेष संपत्ति है - यह केवल एक बार बंद हो जाती है।

हमेशा की तरह $$$a + बार b + c$पहले गणना की $$$a \ _ b$, और फिर यह मान, जिसे सामान्य स्थिति में फ्लोटिंग पॉइंट प्रारूप में नहीं दिखाया जा सकता है, को निकटतम फ्लोट के लिए गोल किया गया है। फिर इस गोल मूल्य में जोड़ दिया जाता है $$$सी$, और यह परिणाम फिर से निकटतम फ्लोट के लिए गोल है। एफएमए को इस तरह से लागू किया जाता है कि केवल अंत में गोलाई का प्रदर्शन किया जाता है - एक मध्यवर्ती मूल्य $$$a \ _ b$पर्याप्त सटीकता बनाए रखता है, इसलिए इसे जोड़ने के बाद $$$सी$तैयार परिणाम सही मूल्य के सबसे करीब होगा $$$a + बार b + c$तैरने का मूल्य।

एफएमए से निपटने के बाद, हम डिफरेंसऑफ प्रॉडक्ट्स DifferenceOfProducts() लौट आएंगे। फिर, मैं इसकी पहली दो पंक्तियाँ दिखाऊंगा:

  float cd = c * d; float err = std::fma(-c, d, cd); 

पहले गोल मान की गणना करता है $$$c \ _ d$और दूसरा ... घटाना $$$c \ _ d$उनके काम से? यदि आप नहीं जानते कि एफएमए कैसे काम करता है, तो आप सोच सकते हैं कि err हमेशा शून्य होगा। लेकिन एफएमए के साथ काम करते समय, दूसरी पंक्ति वास्तव में गणना मूल्य में गोलाई त्रुटि के मूल्य को निकालती है $$$c \ _ d$और इसे err करने के लिए बचाता है। उसके बाद, आउटपुट बहुत सरल है:

  float dop = std::fma(a, b, -cd); return dop + err; 

दूसरा एफएमए एफएमए का उपयोग करके कार्यों के बीच अंतर की गणना करता है, केवल बहुत अंत में गोलाई का प्रदर्शन करता है। नतीजतन, यह भयावह कमी के लिए प्रतिरोधी है, लेकिन इसे एक गोल मूल्य के साथ काम करने की आवश्यकता है। $$$c \ _ d$। दूसरी पंक्ति में हाइलाइट किए गए त्रुटि को जोड़ते हुए return इस समस्या को "पैच" कर देता है। जीननेरॉड एट अल के एक लेख में। यह दिखाया गया है कि परिणाम 1.5 ulps (यूनिट सटीक उपायों) तक सच है, जो बहुत अच्छा है: FMA और सरल फ्लोटिंग पॉइंट ऑपरेशन 0.5 ulps तक सटीक हैं, इसलिए एल्गोरिथ्म लगभग सही है।

हम एक नए हथौड़ा का उपयोग करते हैं

जब आप DifferenceOfProducts() का उपयोग करने के तरीकों की तलाश शुरू करते हैं, तो यह आश्चर्यजनक रूप से उपयोगी साबित होता है। द्विघात समीकरण के विभेदक की गणना? कॉल DifferenceOfProducts(b, b, 4 * a, c) ⁴ । 2x2 मैट्रिक्स के निर्धारक की गणना करना? एल्गोरिथ्म इस समस्या को हल करेगा। Pbrt के अगले संस्करण में , मुझे इसके लिए लगभग 80 उपयोग मिले। उन सभी में से, वेक्टर उत्पाद का कार्य सबसे प्रिय है। यह हमेशा समस्याओं का एक स्रोत रहा है, जिसके कारण आपको अपने हाथों को ऊपर उठाना पड़ा और कार्यान्वयन में दोगुना उपयोग करना पड़ा ताकि आपत्तिजनक कटौती से बच सकें:

 inline Vector3f Cross(const Vector3f &v1, const Vector3f &v2) { double v1x = v1.x, v1y = v1.y, v1z = v1.z; double v2x = v2.x, v2y = v2.y, v2z = v2.z; return Vector3f(v1y * v2z - v1z * v2y, v1z * v2x - v1x * v2z, v1x * v2y - v1y * v2x); } 

और अब हम फ्लोट के साथ काम करना जारी रख सकते हैं और DifferenceOfProducts() उपयोग कर सकते हैं।

 inline Vector3f Cross(const Vector3f &v1, const Vector3f &v2) { return Vector3f(DifferenceOfProducts(v1.y, v2.z, v1.z, v2.y), DifferenceOfProducts(v1.z, v2.x, v1.x, v2.z), DifferenceOfProducts(v1.x, v2.y, v1.y, v2.x)); } 

पोस्ट की शुरुआत से वह चालाक उदाहरण वास्तव में एक वेक्टर काम का हिस्सा है। एक निश्चित चरण में, pbrt कोड को वैक्टर के वेक्टर उत्पाद की गणना करने की आवश्यकता होती है $$$(33962.035, 41563.4, 7706.415)$और $$$(- 24871.969, -30438.8, -5643.727)$। जब फ्लोट का उपयोग करके गणना की जाती है, तो हमें एक वेक्टर मिलेगा $$$(1552, -1248, -128)$। (सामान्य नियम: यदि फ्लोटिंग पॉइंट गणनाओं में जहां बड़ी संख्याएं शामिल हैं, तो आपको बहुत बड़े पूर्णांक मान नहीं मिलते हैं, तो यह लगभग निश्चित संकेत है कि एक भयावह कमी हुई है।)

दोहरी परिशुद्धता के साथ, वेक्टर उत्पाद है $$$(1556.0276, -1257.5151, -75.1656)$। हम देखते हैं कि फ्लोट के साथ $$$x$सामान्य लग रहा है $$$य$पहले से ही बहुत बुरा है $$$z$एक आपदा है जो एक समाधान खोजने के लिए प्रेरणा बन गई है। और हम क्या DifferenceOfProducts() साथ परिणाम प्राप्त करेंगे। उत्पाद DifferenceOfProducts() और फ्लोट मान? $$$(1556.0276, -1257.5153, -75.1656)$। अर्थ $$$x$और $$$z$डबल परिशुद्धता के अनुरूप है, और $$$य$थोड़ा ऑफसेट - इसलिए अतिरिक्त ulp से आया।

गति के बारे में क्या? डिफरऑफप्रोडक्ट्स DifferenceOfProducts() दो एफएमए करता है, साथ ही गुणा और जोड़ भी। एक भोले एल्गोरिथ्म को एक एफएमए और एक गुणा के साथ लागू किया जा सकता है, जो, ऐसा प्रतीत होता है, आधे से अधिक समय लेना चाहिए। लेकिन व्यवहार में, यह पता चला है कि रजिस्टरों से मूल्यों को प्राप्त करने के बाद, यह मायने नहीं रखता है: मेरे लैपटॉप पर आयोजित सिंथेटिक बेंचमार्क में, DifferenceOfProducts() भोली एल्गोरिथ्म की तुलना में केवल 1.09 गुना अधिक महंगा है। डबल प्रिसिजन ऑपरेशन 2.98 गुना धीमा था।

जैसे ही आप एक भयावह कमी की संभावना के बारे में सीखते हैं, कोड में सभी प्रकार के निर्दोष दिखने वाले भाव संदिग्ध लगने लगते हैं। DifferenceOfProducts() उनमें से अधिकांश के लिए एक अच्छा इलाज लगता है। इसका उपयोग करना आसान है, और हमारे पास इसका उपयोग न करने का कोई विशेष कारण नहीं है।

नोट

1. विभिन्न संकेतों के साथ मात्रा घटाना या एक ही संकेत के साथ मूल्यों को जोड़ने पर भयावह कमी कोई समस्या नहीं है। हालांकि, विभिन्न संकेतों के साथ मूल्यों को जोड़ने पर यह एक समस्या बन सकती है। यही है, राशियों को मतभेद के समान संदेह के साथ माना जाना चाहिए।
2. पाठक के लिए एक अभ्यास के रूप में, मैं SumOfProducts() फ़ंक्शन लिखने का सुझाव देता हूं, जो भयावह संकुचन से सुरक्षा प्रदान करता है। यदि आप कार्य को जटिल करना चाहते हैं, तो बताएं कि dop + err == dop DifferenceOfProducts() , dop + err == dop में क्यों, क्योंकि संकेत a*b और c*d अलग हैं।
3. एफएमए निर्देश एक दशक से अधिक समय तक जीपीयू पर उपलब्ध रहा है, और कम से कम पांच वर्षों के लिए अधिकांश सीपीयू पर। सीपीयू का उपयोग करते समय सीधे उन्हें उत्पन्न करने के लिए सीपीयू में संकलक झंडे जोड़ना आवश्यक हो सकता है std::fma() ; gcc और -march=native , -march=native काम करता है।
4. IEEE फ़्लोटिंग पॉइंट फॉर्मेट में, दो की शक्तियों द्वारा गुणन ठीक-ठीक किया जाता है, इसलिए 4 * a किसी भी गोल त्रुटि का कारण नहीं बनता है, जब तक कि अतिप्रवाह नहीं होता है।