डबल-सटीक तर्कों से शुरू होने वाले 80-बिट विस्तारित परिशुद्धता कंप्यूटेशंस की गुण

Question

इंटरपोलेशन फ़ंक्शंस के दो कार्यान्वयन यहां दिए गए हैं। तर्क u1 हमेशा 0. और 1. के बीच है।

#include <stdio.h> 

double interpol_64(double u1, double u2, double u3) 
{ 
    return u2 * (1.0 - u1) + u1 * u3; 
} 

double interpol_80(double u1, double u2, double u3) 
{ 
    return u2 * (1.0 - (long double)u1) + u1 * (long double)u3; 
} 

int main() 
{ 
    double y64,y80,u1,u2,u3; 
    u1 = 0.025; 
    u2 = 0.195; 
    u3 = 0.195; 
    y64 = interpol_64(u1, u2, u3); 
    y80 = interpol_80(u1, u2, u3); 
    printf("u2: %a\ny64:%a\ny80:%a\n", u2, y64, y80); 
} 

साथ 80-बिट long double रों एक सख्त आईईईई 754 मंच पर, interpol_64() में सभी संगणना आईईईई 754 डबल परिशुद्धता के अनुसार किया, और 80-बिट बढ़ाया परिशुद्धता में interpol_80() में कर रहे हैं। कार्यक्रम प्रिंट:

u2: 0x1.8f5c28f5c28f6p-3 
y64:0x1.8f5c28f5c28f5p-3 
y80:0x1.8f5c28f5c28f6p-3 

मैं संपत्ति में दिलचस्पी है "परिणाम समारोह से वापस लौटे हमेशा होता है के बीच u2 और u3"। यह गुण interpol_64() का झूठा है, जैसा ऊपर main() में मानों द्वारा दिखाया गया है।

क्या संपत्ति के पास interpol_80() का सत्य होने का मौका है? यदि यह नहीं है, तो काउंटर उदाहरण क्या है? अगर हम जानते हैं कि u2 != u3 या उनके बीच न्यूनतम दूरी है तो क्या इससे मदद मिलती है? क्या इंटरमीडिएट कंप्यूटेशंस के लिए एक महत्वपूर्ण चौड़ाई निर्धारित करने के लिए कोई तरीका है जिस पर संपत्ति को सत्य होने की गारंटी दी जाएगी?

संपादित करें: मैंने कोशिश किए गए सभी यादृच्छिक मूल्यों पर, संपत्ति को आयोजित किया जब इंटरमीडिएट कंप्यूटेशंस आंतरिक रूप से विस्तारित परिशुद्धता में किया गया था। यदि interpol_80() ने long double तर्क लिया, तो काउंटर-उदाहरण बनाने के लिए अपेक्षाकृत आसान होगा, लेकिन यहां प्रश्न विशेष रूप से ऐसे फ़ंक्शन के बारे में है जो double तर्क लेता है। यह एक काउंटर उदाहरण बनाने के लिए बहुत कठिन बनाता है, अगर कोई है।

---

नोट: एक संकलक पैदा x87 निर्देश interpol_64() और interpol_80() के लिए एक ही कोड उत्पन्न कर सकता है, लेकिन यह मेरे सवाल का स्पर्शरेखा है।

Answer 1

हाँ, interpol_80() सुरक्षित है, चलिए इसे प्रदर्शित करते हैं।

समस्या कहा गया है कि आदानों 64bits फ्लोट कर रहे हैं

rnd64(ui) = ui 

परिणाम वास्तव में है (यह मानते हुए * और + गणितीय क्रियाओं कर रहे हैं)

r = u2*(1-u1)+(u1*u3) 

इष्टतम वापसी मान 64 बिट फ्लोट करने के लिए गोल

है

r64 = rnd64(r) 

जैसा कि हमारे पास ये गुण हैं

u2 <= r <= u3 

यह गारंटी है कि U1, U2, U3 के 80bits को

rnd64(u2) <= rnd64(r) <= rnd64(u3) 
u2 <= r64 <= u3 

रूपांतरण सटीक भी हैं। यहां तक ​​कि सबसे पास करने के लिए

rf = rnd(rnd(u2*rnd(1-u1)) + rnd(u1*u3)) 

मानते हुए दौर, यह सबसे कम 2 ULP बंद सटीक हो जाएगा:

rnd80(ui)=ui 

अब, चलो 0 <= u2 <= u3 मान, तो अयथार्थ नाव संचालन के साथ प्रदर्शन करते हैं ज्यादा से ज्यादा 4 राउंडिंग त्रुटियों की ओर जाता है मूल्य। तो राउंडिंग 64 बिट्स नाव या 80 बिट के साथ किया जाता तैरता:

r - 2 ulp64(r) <= rf64 <= r + 2 ulp64(r) 
r - 2 ulp80(r) <= rf80 <= r + 2 ulp80(r) 

rf64 2 ULP अंतर हो सकता है तो इंटरपोल-64() rnd64(rf80) के बारे में असुरक्षित है, लेकिन क्या?
हम चाहते हैं कि बता सकते हैं:

rnd64(r - 2 ulp80(r)) <= rnd64(rf80) <= rnd64(r + 2 ulp80(r)) 

के बाद से 0 <= u2 <= u3, तो

ulp80(u2) <= ulp80(r) <= ulp80(r3) 
rnd64(u2 - 2 ulp80(u2)) <= rnd64(r - 2 ulp80(r)) <= rnd64(rf80) 
rnd64(u3 + 2 ulp80(u3)) >= rnd64(r + 2 ulp80(r)) >= rnd64(rf80) 

सौभाग्य से, सीमा (u2-ulp64(u2)/2 , u2+ulp64(u2)/2) में हर संख्या की तरह हम पाते हैं

rnd64(u2 - 2 ulp80(u2)) = u2 
rnd64(u3 + 2 ulp80(u3)) = u3 

ulp80(x)=ulp62(x)/2^(64-53) के बाद से

हम इस प्रकार U2 < = U3 < = 0 के लिए सबूत

u2 <= rnd64(rf80) <= u3 

मिलता है, हम आसानी से एक ही सबूत आवेदन कर सकते हैं।

अध्ययन का अंतिम मामला यू 2 < = 0 < = u3 है। यदि हम 2 बड़े मान घटाते हैं, तो परिणाम उल (बड़े-बड़े)/2 के बजाय उल (बड़ा)/2 बंद हो सकता है ...
इस प्रकार इस दावे को हम बना अब और नहीं रखता है:

r - 2 ulp64(r) <= rf64 <= r + 2 ulp64(r) 

सौभाग्य से, u2 <= u2*(1-u1) <= 0 <= u1*u3 <= u3 और इस

u2 <= rnd(u2*rnd(1-u1)) <= 0 <= rnd(u1*u3) <= u3 

गोलाई इस प्रकार के बाद संरक्षित है के बाद से जोड़ा मात्रा विपरीत संकेत के होते हैं:

u2 <= rnd(u2*rnd(1-u1)) + rnd(u1*u3) <= u3 

राउंडिंग के बाद भी जाता है, इसलिए हम एक बार फिर

की गारंटी दे सकते हैं

u2 <= rnd64(rf80) <= u3 

QED

पूरा हम denormal आदानों (क्रमिक अधःप्रवाह) का ध्यान देना चाहिए होने के लिए, लेकिन मुझे आशा है कि आपको लगता है कि तनाव परीक्षणों के साथ शातिर नहीं होगा। मैं प्रदर्शित नहीं होगा कि उन के साथ होता है ...

संपादित:

यहाँ एक अनुवर्ती के रूप में निम्नलिखित धारणा थोड़ा समीपवर्ती था और कुछ टिप्पणियों पर जेनरेट किया गया 0 < = u2 < = U3

r - 2 ulp80(r) <= rf80 <= r + 2 ulp80(r) 

हम निम्नलिखित असमानताओं लिख सकते हैं:

rnd(1-u1) <= 1 
rnd(1-u1) <= 1-u1+ulp(1)/4 
u2*rnd(1-u1) <= u2 <= r 
u2*rnd(1-u1) <= u2*(1-u1)+u2*ulp(1)/4 
u2*ulp(1) < 2*ulp(u2) <= 2*ulp(r) 
u2*rnd(1-u1) < u2*(1-u1)+ulp(r)/2 

के लिए अगले राउंडिंग आपरेशन, हम योग के दूसरे भाग के लिए

ulp(u2*rnd(1-u1)) <= ulp(r) 
rnd(u2*rnd(1-u1)) < u2*(1-u1)+ulp(r)/2 + ulp(u2*rnd(1-u1))/2 
rnd(u2*rnd(1-u1)) < u2*(1-u1)+ulp(r)/2 + ulp(r)/2 
rnd(u2*rnd(1-u1)) < u2*(1-u1)+ulp(r) 

उपयोग करते हैं, हमने:

u1*u3 <= r 
rnd(u1*u3) <= u1*u3 + ulp(u1*u3)/2 
rnd(u1*u3) <= u1*u3 + ulp(r)/2 

rnd(u2*rnd(1-u1))+rnd(u1*u3) < u2*(1-u1)+u1*u3 + 3*ulp(r)/2 
rnd(rnd(u2*rnd(1-u1))+rnd(u1*u3)) < r + 3*ulp(r)/2 + ulp(r+3*ulp(r)/2)/2 
ulp(r+3*ulp(r)/2) <= 2*ulp(r) 
rnd(rnd(u2*rnd(1-u1))+rnd(u1*u3)) < r + 5*ulp(r)/2 

मैं मूल दावा साबित नहीं हुआ है, लेकिन वह दूर नहीं ...

Answer 2

interpol_64 में हानि-परिशुद्धता का मुख्य स्रोत गुणा है। दो 53-बिट मंटिसस गुणा करने से 105- या 106-बिट उत्पन्न होता है (उच्च बिट कैर्री पर निर्भर करता है) मंटिसा। 80-बिट विस्तारित सटीक मान में फ़िट होने के लिए यह बहुत बड़ा है, इसलिए सामान्य रूप से, आपके पास 80-बिट संस्करण में भी हानि-परिशुद्धता होगी। जब ऐसा होता है तो सटीक मात्रा में बहुत मुश्किल होती है; सबसे ज्यादा कहना आसान है कि यह तब होता है जब गोल त्रुटियां जमा होती हैं। ध्यान दें कि दो शर्तों को जोड़ते समय एक छोटा राउंडिंग चरण भी है।

अधिकांश लोगों को शायद वैसे ही जैसे एक समारोह के साथ इस समस्या का समाधान होगा:

double interpol_64(double u1, double u2, double u3) 
{ 
    return u2 + u1 * (u3 - u2); 
} 

लेकिन ऐसा लगता है कि आप राउंडिंग मुद्दों, नहीं बेहतर कार्यान्वयन में अंतर्दृष्टि की तलाश कर रहे हैं।