هل هناك أي طريقة للحصول على التقريب الصحيح باستخدام تعليمات i387 fsqrt؟

Question

هل هناك أي طريقة للحصول على التقريب الصحيح باستخدام تعليمات i387 fsqrt؟...

...بصرف النظر عن تغيير وضع الدقة في كلمة التحكم x87 - أعلم أن هذا ممكن، لكنه ليس حلاً معقولاً لأنه يحتوي على مشكلات سيئة من نوع إعادة الدخول حيث سيكون وضع الدقة خاطئًا في حالة مقاطعة عملية sqrt.

المشكلة التي أتعامل معها هي كما يلي:x87 fsqrt ينفذ كود التشغيل عملية جذر تربيعي مستديرة بشكل صحيح (حسب IEEE 754) بدقة سجلات fpu، والتي سأفترض أنها دقة ممتدة (80 بت).ومع ذلك، أريد استخدامه لتنفيذ وظائف جذر تربيعي فعالة مفردة ومزدوجة الدقة مع تقريب النتائج بشكل صحيح (حسب وضع التقريب الحالي).نظرًا لأن النتيجة لها دقة زائدة، فإن الخطوة الثانية هي تحويل النتيجة إلى جولات دقة مفردة أو مزدوجة مرة أخرى، مما قد يؤدي إلى ترك نتيجة غير مقربة بشكل صحيح.

مع بعض العمليات، من الممكن التغلب على هذه المشكلة مع التحيزات.على سبيل المثال، يمكنني تجنب الدقة الزائدة في نتائج الجمع عن طريق إضافة انحياز في شكل قوة اثنين والتي تجبر 52 بتًا كبيرًا من قيمة الدقة المزدوجة على آخر 52 بت من الجزء العشري ذو الدقة الموسعة 63 بت .لكنني لا أرى أي طريقة واضحة للقيام بهذه الخدعة باستخدام الجذر التربيعي.

أي أفكار ذكية؟

(تم وضع علامة C أيضًا لأن التطبيق المقصود هو تنفيذ C sqrt و sqrtf المهام.)

Answer 1

أولاً، دعونا نبتعد عن ما هو واضح:يجب أن تستخدم SSE بدلاً من x87.إس إس إي sqrtss و sqrtsd التعليمات تفعل بالضبط ما تريد، وهي مدعومة على جميع أنظمة x86 الحديثة، كما أنها أسرع بشكل ملحوظ أيضًا.

الآن، إذا كنت مصرًا على استخدام x87، فسأبدأ بالأخبار الجيدة:لا تحتاج إلى القيام بأي شيء من أجل تعويم.انت تحتاج 2p + 2 بتات لحساب جذر تربيعي مستدير بشكل صحيح بتنسيق النقطة العائمة p-bit.لأن 80 > 2*24 + 2, ، سيتم دائمًا تقريب التقريب الإضافي إلى الدقة الفردية بشكل صحيح، وسيكون لديك جذر تربيعي تم تقريبه بشكل صحيح.

الآن الأخبار السيئة: 80 < 2*53 + 2, ، لذلك لا يوجد مثل هذا الحظ للدقة المزدوجة.يمكنني أن أقترح عدة حلول؛إليكم فكرة سهلة ولطيفة من أعلى رأسي.

1. يترك y = round_to_double(x87_square_root(x));
2. استخدم منتج Dekker (الرأس والذيل) للحساب a و b مثل ذلك y*y = a + b بالضبط.
3. حساب المتبقية r = x - a - b.
4. if (r == 0) return y
5. if (r > 0), ، يترك y1 = y + 1 ulp, ، وحساب a1, b1 شارع. y1*y1 = a1 + b1.يقارن r1 = x - a1 - b1 ل r, ، والعودة سواء y أو y1, ، اعتمادًا على أي منها يحتوي على الجزء المتبقي الأصغر (أو الذي يحتوي على بتة ذات ترتيب منخفض صفر، إذا كانت البقايا متساوية في الحجم).
6. if (r < 0), ، افعل نفس الشيء من أجل y1 = y - 1 ulp.

يعالج هذا الإجراء وضع التقريب الافتراضي فقط؛ومع ذلك، في أوضاع التقريب الموجهة، فإن التقريب إلى تنسيق الوجهة يؤدي إلى الإجراء الصحيح.

Answer 2

حسنًا ، أعتقد أن لدي حلًا أفضل:

1. إحصاء - عد y=sqrt(x) بدقة ممتدة (fsqrt).
2. إذا كانت آخر 11 بت ليست كذلك 0x400, ، ما عليك سوى التحويل إلى دقة مضاعفة والعودة.
3. يضيف 0x100-(fpu_status_word&0x200) إلى الكلمة المنخفضة لتمثيل الدقة الموسعة.
4. تحويل إلى دقة مزدوجة والعودة.

تعتمد الخطوة 3 على حقيقة أن البت C1 (0x200) لكلمة الحالة هو 1 إذا وفقط إذا fsqrtتم تقريب نتيجة .وهذا صحيح لأنه، بسبب الاختبار في الخطوة 2، x لم يكن مربعا كاملا.لو كان مربعا كاملا، y لن يكون لها أي أجزاء تتجاوز الدقة المزدوجة.

قد يكون من الأسرع تنفيذ الخطوة 3 مع تشغيل النقطة العائمة المشروطة بدلاً من العمل على تمثيل البت وإعادة التحميل.

إليك الكود (يبدو أنه يعمل في جميع الحالات):

sqrt:
    fldl 4(%esp)
    fsqrt
    fstsw %ax
    sub $12,%esp
    fld %st(0)
    fstpt (%esp)
    mov (%esp),%ecx
    and $0x7ff,%ecx
    cmp $0x400,%ecx
    jnz 1f
    and $0x200,%eax
    sub $0x100,%eax
    sub %eax,(%esp)
    fstp %st(0)
    fldt (%esp)
1:  add $12,%esp
    fstpl 4(%esp)
    fldl 4(%esp)
    ret

Answer 3

قد لا يكون ما تريده، لأنه لا يستفيد من 387 fsqrt التعليمات، ولكن هناك كفاءة بشكل مدهش sqrtf(float) في glibc تم تنفيذها باستخدام حساب عدد صحيح 32 بت.حتى أنه يتعامل مع NaNs وInfs والمستويات الفرعية بشكل صحيح - قد يكون من الممكن إزالة بعض عمليات التحقق هذه باستخدام تعليمات x87 الحقيقية/أعلام كلمات التحكم FP.يرى: glibc-2.14/sysdeps/ieee754/flt-32/e_sqrtf.c

ال dbl-64/e_sqrt.c الكود ليس ودودًا جدًا.من الصعب معرفة الافتراضات التي يتم وضعها في لمحة واحدة.ومن الغريب أن المكتبة i386 sqrt[f|l] التنفيذ فقط اتصل fsqrt, ، ولكن قم بتحميل القيمة بشكل مختلف. flds ل س، fldl لموانئ دبي.