आईईईई 754 अनुरूप sqrt () दोहरे प्रकार के लिए कार्यान्वयन

Aug 17 2020

मैं double __ieee754_sqrt(double x)फंक्शन को लागू करने की कोशिश कर रहा हूं, जो पहली सन्निकटन प्राप्त करने के लिए हार्डवेयर निर्देश का उपयोग करता है:

double __ieee754_sqrt(double x) {
    double z;
    /* get reciprocal of the square root (6.75 bits accuracy) */
    __asm(" QSEED.DF %0,%1 \n": "=e" (z):"e" (x):);
    z = 1 / z;
    z = ( z + x / z) / 2; /* 1st Newton-Raphson iteration */
    z = ( z + x / z) / 2; /* 2nd Newton-Raphson iteration */
    z = ( z + x / z) / 2; /* 3rd Newton-Raphson iteration */
    z = ( z + x / z) / 2; /* 4th Newton-Raphson iteration */
    return z;
}

हालाँकि, paranoia.c ( लिंक , लिंक ) परीक्षण की शिकायत करता है:

Square root is neither chopped nor correctly rounded.
Observed errors run from -6.0493828e-01 to 5.0000000e-01 ulps.

प्रश्न: कैसे के लिए अतिरिक्त तर्क को लागू करने के लिए chopping and correct rounding?

यूपीडी हार्डवेयर मूल रूप से समर्थन नहीं करता है sqrt()। हार्डवेयर केवल स्क्वायर रूट (6.75 बिट्स सटीकता) के पारस्परिक प्राप्त करने का समर्थन करता है।

UPD2।

खेतों में प्रयुक्त njuffa समाधान (बहुत धन्यवाद!) मामूली परिवर्तन के साथ: उपयोग qseeddf()करने के बजाय qseedf()=> उपयोग fma()के बजाय fmaf()। क्यों? क्योंकि यह double<=>floatरूपांतरणों को छोड़ देता है और इसलिए तेजी से होता है।
हां, फ़्यूज्ड मल्टीप्ले-ऐड निर्देश (FMA) हार्डवेयर द्वारा समर्थित हैं।
चर्चा में भाग लेने और विस्तृत उत्तरों के लिए सभी को धन्यवाद!
इस विषय में रुचि रखने वाले सभी के लिए, यहां sqrt()कार्यान्वयन की सूची दी गई है :
1. Cygwin गणित से। पुस्तकालय ( libm) cygwin-snapshot-20200710-1/newlib/libm/math/e_sqrt.c: कॉपीराइट Copyright (C) 1993 by Sun Microsystems।
2. GNU C लाइब्रेरी से ( glibc):
  1. glibc-2.31/sysdeps/ieee754/dbl-64/e_sqrt.c: हकदार IBM Accurate Mathematical Library।
  2. glibc-2.31/sysdeps/powerpc/fpu/e_sqrt.c: __builtin_fma()कार्यों का उपयोग करना।

जवाब

2 njuffa Aug 18 2020 at 06:38

किसी के स्वयं के कार्यान्वयन के निर्माण से पहले, उपयुक्त और अच्छी तरह से परीक्षण किए गए ओपन-सोर्स कोड उपलब्ध है या नहीं, यह जांचने के लिए इंटरनेट पर खोज करना उचित है।

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

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

नीचे अनुकरणीय आईएसओ-सी 99 कार्यान्वयन में मैं दिखा रहा हूं कि उन लाइनों के साथ एक सही ढंग से गोल डबल-सटीक वर्गमूल कैसे लागू किया जा सकता है। मुझे लगता है कि यह सोचते हैं हूँ doubleआईईईई-754 के लिए नक्शे binary64और कहा कि floatआईईईई-754 के लिए नक्शे binary32। मैं sqrtIEEE-754 राउंड-टू-निकटतम या सम-मोड के साथ कार्यान्वित होने के लिए प्रतिबंधित कर रहा हूं ।

बहुत महत्वपूर्ण बात मैं यह सोचते हैं कि प्रक्रिया हार्डवेयर जुड़े हुए प्रदान करता निर्देश गुणा-जोड़ सकते हैं और इन सही ढंग से मानक गणित पुस्तकालय कार्यों के माध्यम से संपर्क में हैं कि fmafऔर fma। टिप्पणियों में मैंने एफएमए की उपलब्धता के अनुसार ओपी से स्पष्टीकरण मांगा था, लेकिन प्रतिक्रिया उपलब्ध होने से पहले कोड पर शुरू करने का फैसला किया। एफएमए के बिना कार्यान्वयन संभव है लेकिन बहुत अधिक चुनौतीपूर्ण है, और एक पर्याप्त रूप से पूर्ण उपचार एक स्टैकओवरफ्लो उत्तर के स्थान से अधिक होगा।

चूंकि ओपी ने लक्ष्य वास्तुकला को निर्दिष्ट नहीं किया है या शुरुआती सन्निकटन का विवरण प्रदान नहीं किया है, इसलिए मैं नीचे दिए गए अपने स्वयं के शुरुआती सन्निकटन का उपयोग कर रहा हूं जो कि अंतराल [0.25, 1] पर एक बहुपद न्यूनतम अनुमान पर आधारित है, जिसमें सभी गैर-असाधारण तर्क कम हो गए हैं। qseedf()परिणाम लगभग 7 बिट के लिए सटीक हैं, इसलिए ओपी के अंतर्निहित हार्डवेयर की तुलना में थोड़ा बेहतर है। क्या यह अंतर महत्वपूर्ण है, मैं आकलन नहीं कर सकता।

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

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <string.h>
#include <math.h>

/* Approximate 1/sqrt(a) on [0.25, 1] with an accuracy of about 7 bits */
float qseedf (float a)
{
    float r;

    r =             -2.43845296f;
    r = fmaf (r, a,  6.22994471f);
    r = fmaf (r, a, -5.91090727f);
    r = fmaf (r, a,  3.11237526f);
    return r;
}

double my_sqrt (double a)
{    
    const double QNAN_INDEFINITE = 0.0 / 0.0;
    const double half = 0.5;
    const double three_eighth = 0.375;
    double refined_rsqrt_approx, sqrt_approx, sqrt_residual, result, b;
    double rsqrt_approx, rsqrt_approx_err, rsqrt_approx_squared, reduced_arg;
    float argf, approxf, approxf_err;
    int e, t, f;

    /* handle normal cases */
    if ((a >= 0) && (a < INFINITY)) {
        /* compute exponent adjustments */
        b = frexp (a, &e);
        t = e - 2*512;
        f = t / 2;
        t = t - 2 * f;
        f = f + 512;

        /* map argument into the primary approximation interval [0.25,1) */
        reduced_arg = ldexp (b, t);
        
        /* Compute initial low-precision approximation */
        argf = (float)reduced_arg;
        approxf = qseedf (argf);
        
        /* Apply two Newton-Raphson iterations with quadratic convergence */
        approxf_err = fmaf (-argf, approxf * approxf, 1.0f);
        approxf = fmaf (0.5f * approxf, approxf_err, approxf);
        approxf_err = fmaf (-argf, approxf * approxf, 1.0f);
        approxf = fmaf (0.5f * approxf, approxf_err, approxf);
        
        /* rsqrt approximation is now accurate to 1 single-precision ulp */
        rsqrt_approx = (double)approxf;

        /* Perform a Halley iteration wih cubic convergence. Based on the work
           of Peter Markstein. See: Peter Markstein, "IA-64 and Elementary 
           Functions", Prentice Hall 2000
        */
        rsqrt_approx_squared = rsqrt_approx * rsqrt_approx;
        rsqrt_approx_err = fma (-reduced_arg, rsqrt_approx_squared, 1.0);
        refined_rsqrt_approx = fma (fma (rsqrt_approx_err, three_eighth, half), 
                                rsqrt_approx * rsqrt_approx_err, rsqrt_approx);
        sqrt_approx = reduced_arg * refined_rsqrt_approx;
        sqrt_residual = fma (-sqrt_approx, sqrt_approx, reduced_arg);
        result = fma (sqrt_residual, half * refined_rsqrt_approx, sqrt_approx);

        /* map back from primary approximation interval by jamming exponent */
        result = ldexp (result, f);
    } else {
        /* handle special cases */
        result = (a < 0) ? QNAN_INDEFINITE : (a + a);
    }
    return result;
}

/*
  https://groups.google.com/forum/#!original/comp.lang.c/qFv18ql_WlU/IK8KGZZFJx4J
  From: geo <[email protected]>
  Newsgroups: sci.math,comp.lang.c,comp.lang.fortran
  Subject: 64-bit KISS RNGs
  Date: Sat, 28 Feb 2009 04:30:48 -0800 (PST)

  This 64-bit KISS RNG has three components, each nearly
  good enough to serve alone.    The components are:
  Multiply-With-Carry (MWC), period (2^121+2^63-1)
  Xorshift (XSH), period 2^64-1
  Congruential (CNG), period 2^64
*/
static uint64_t kiss64_x = 1234567890987654321ULL;
static uint64_t kiss64_c = 123456123456123456ULL;
static uint64_t kiss64_y = 362436362436362436ULL;
static uint64_t kiss64_z = 1066149217761810ULL;
static uint64_t kiss64_t;
#define MWC64  (kiss64_t = (kiss64_x << 58) + kiss64_c, \
                kiss64_c = (kiss64_x >> 6), kiss64_x += kiss64_t, \
                kiss64_c += (kiss64_x < kiss64_t), kiss64_x)
#define XSH64  (kiss64_y ^= (kiss64_y << 13), kiss64_y ^= (kiss64_y >> 17), \
                kiss64_y ^= (kiss64_y << 43))
#define CNG64  (kiss64_z = 6906969069ULL * kiss64_z + 1234567ULL)
#define KISS64 (MWC64 + XSH64 + CNG64)

int main (void)
{
    const uint64_t N = 10000000000ULL; /* desired number of test cases */
    double arg, ref, res;
    uint64_t argi, refi, resi, count = 0;
    double spec[] = {0, 1, INFINITY, NAN};

    printf ("test a few special cases:\n");
    for (int i = 0; i < sizeof (spec)/sizeof(spec[0]); i++) {
        printf ("my_sqrt(%22.13a) = %22.13a\n", spec[i], my_sqrt(spec[i]));
        printf ("my_sqrt(%22.13a) = %22.13a\n", -spec[i], my_sqrt(-spec[i]));
    }
    
    printf ("test %llu random cases:\n", N);
    do {
        count++;
        argi = KISS64;
        memcpy (&arg, &argi, sizeof arg);
        res = my_sqrt (arg);
        ref = sqrt (arg);
        memcpy (&resi, &res, sizeof resi);
        memcpy (&refi, &ref, sizeof refi);
        if (resi != refi) {
            printf ("\rerror @ arg=%22.13a  res=%22.13a  ref=%22.13a\n",
                    arg, res, ref);
            return EXIT_FAILURE;
        }
        if ((count & 0xfffff) == 0) printf ("\r[%llu]", count);
    } while (count < N);
    printf ("\r[%llu]", count);
    printf ("\ntests PASSED\n");
    return EXIT_SUCCESS;
}

उपरोक्त कार्यक्रम का आउटपुट इस तरह दिखना चाहिए:

test a few special cases:
my_sqrt(  0x0.0000000000000p+0) =   0x0.0000000000000p+0
my_sqrt( -0x0.0000000000000p+0) =  -0x0.0000000000000p+0
my_sqrt(  0x1.0000000000000p+0) =   0x1.0000000000000p+0
my_sqrt( -0x1.0000000000000p+0) =  -0x1.#IND000000000p+0
my_sqrt(  0x1.#INF000000000p+0) =   0x1.#INF000000000p+0
my_sqrt( -0x1.#INF000000000p+0) =  -0x1.#IND000000000p+0
my_sqrt(  0x1.#QNAN00000000p+0) =   0x1.#QNAN00000000p+0
my_sqrt( -0x1.#QNAN00000000p+0) =  -0x1.#QNAN00000000p+0
test 10000000000 random cases:
[10000000000]
tests PASSED

1 RickJames Aug 18 2020 at 01:08

z = 1 / z;
z = ( z + x / z) / 2; /* 1st Newton-Raphson iteration */
...

z = 1 / z;
z += ( x / z - z) * 0.5; /* 1st Newton-Raphson iteration */
...

यह तेज हो सकता है।

और एक पुनरावृति को जल्दी से रोकें (मुझे लगता है।)

जब आप रुकें, तुलना करें z*zऔर x। z*zकोई तुलना में छोटा होगा (मुझे लगता है कि) x। 1ulp से घटाएँ zऔर z*zबनाम जाँचें x। यह "सही गोलाई" की सही जाँच नहीं है, लेकिन इसके बीच zऔर तय करने के लिए "काफी अच्छा" हो सकता है z - 1ulp।

चूँकि आपको त्रुटियों की इतनी बड़ी रेंज मिली है, इसलिए मुझे चिंता है कि राउंडिंग, या यहाँ तक कि सटीक होने पर फ़्लोटिंग पॉइंट 'हार्डवेयर' का बाकी हिस्सा खराब है।

उफ़, मैं भूल गया। आपको एक अनुमान देने का एक कारण था 1/z- लगभग 1 / z जारी रखें; आप इसे विभाजित करने के बजाय कई गुना कर सकते हैं, जिससे (अधिकांश हार्डवेयर में) काफी तेजी से और संभवतः कम राउंडऑफ के साथ।

z = ( z + x * z) * 0.5; /* 1st Newton-Raphson iteration */
...
z = 1 / z;

इसके अलावा, देखें कि क्या किसी के लिए गुणा करने के बजाय घातांक को घटाने का कोई तरीका है / 2।