حسابات دقيقة وسريعة لأرقام الفاصلة العائمة باستخدام مثال دالة الجيب. الجزء 3: النقطة الثابتة

نواصل دورة من المحاضرات ( جزء 1 و جزء 2 ). في الجزء 2 ، نظرنا إلى ما يوجد داخل مكتبة libm وفي هذا العمل سنحاول تغيير وظيفة do_sin قليلاً لزيادة دقتها وسرعتها. سأذكر هذه الوظيفة مرة أخرى do_sin ):







كما هو موضح في المقالة السابقة ، الجزء 132-145. تم تنفيذه من أجل x ضمن النطاق [0.126، 0.855469]. حسنا. دعنا نحاول كتابة دالة ستكون أكثر دقة ، وربما أسرع ، ضمن الحدود المعطاة.



الطريقة التي سنستخدمها بها واضحة جدًا. من الضروري توسيع دقة الحسابات لتشمل المزيد من المنازل العشرية. سيكون الحل الواضح هو اختيار النوع المزدوج الطويل ، والعد فيه ، ثم التحويل مرة أخرى. من حيث الدقة ، يجب أن يكون الحل جيدًا ، ولكن من حيث الأداء ، قد تكون هناك مشاكل. مع ذلك ، تعد المضاعفة الطويلة نوعًا غريبًا إلى حد ما من البيانات ولا يمثل دعمها في المعالجات الحديثة أولوية. في تعليمات x86_64 SSE / AVX مع نوع البيانات هذا لا تعمل. المعالج الرياضي سوف "ينفجر".



إذن ماذا يجب أن تختار؟ لنلق نظرة فاحصة على حدود السعة والوظيفة.



هم في المنطقة 1.0. أولئك. في الواقع ، لا نحتاج إلى نقطة عائمة. لنستخدم عددًا صحيحًا 64 بت عند حساب الوظيفة. سيعطينا هذا 10-11 بت إضافية للدقة الأصلية. دعنا نتعرف على كيفية التعامل مع هذه الأرقام. يتم تمثيل الرقم في هذا التنسيق كـ a / d ، حيث يمثل a عددًا صحيحًا ، و d هو القاسم الذي نختاره ثابتًا لجميع المتغيرات ونخزنه "في ذاكرتنا" ، وليس في ذاكرة الكمبيوتر. فيما يلي بعض العمليات لهذه الأرقام:

$$

$$\frac{c}{d}=\frac{a}{d}\pm\frac{b}{d}=\frac{a\pm b}{d}\\ \frac{c}{d}=\frac{a}{d}\cdot\frac{b}{d}=\frac{a \cdot b}{d^2}\\ \frac{c}{d}=\frac{a}{d}\cdot x=\frac{a \cdot x}{d}$$

كما ترون ، لا يوجد شيء معقد في ذلك. تُظهر الصيغة الأخيرة الضرب بأي عدد صحيح. لاحظ أيضًا شيئًا واضحًا إلى حد ما وهو أن نتيجة ضرب متغيرين عدد صحيحين غير مشارَّين بالحجم N غالبًا ما يكون عددًا من الحجم يصل إلى 2 * N شاملًا. يمكن أن تتسبب الإضافة في تجاوز يصل إلى 1 بت إضافي.



لنحاول اختيار المقسوم عليه د . من الواضح ، في العالم الثنائي ، أنه من الأفضل اختياره كقوة ثنائية ، حتى لا يتم تقسيمه ، ولكن فقط نقل السجل ، على سبيل المثال. ما هي قوة اثنين يجب أن تختار؟ ابحث عن التلميح في تعليمات آلة الضرب. على سبيل المثال ، تضاعف تعليمات MUL القياسية في نظام x86 سجلين وتكتب النتيجة في سجلين أيضًا ، حيث يكون أحد السجلات هو "الجزء العلوي" من النتيجة ، والثاني هو الجزء السفلي.



على سبيل المثال ، إذا كان لدينا رقمان من 64 بت ، فستكون النتيجة رقم 128 بت مكتوبًا في سجلين 64 بت. دعنا نسمي RH - الحالة "العلوية" ، و RL - الحالة "السفلية" ¹ . ثم ، رياضيًا ، يمكن كتابة النتيجة كـ$$$R=R_H \cdot 2^{64}+R_L$... الآن نستخدم الصيغ أعلاه ونكتب الضرب لـ$$$d=2^{-64}$

$$

$$\frac{c}{d}=\frac{a}{2^{64}}\cdot\frac{b}{2^{64}}=\frac{a \cdot b}{2^{128}}=\frac{R_H \cdot 2^{64} + R_L}{2^{128}}=\frac{R_H + R_L \cdot 2^{-64}}{2^{64}}$$

وقد اتضح أن نتيجة ضرب هذين العددين ذات النقطة الثابتة هي السجل $$$R=R_H$...



إنه أسهل بالنسبة لنظام Aarch64. تعليمات "UMULH" تضاعف سجلين وتكتب الجزء "العلوي" من عملية الضرب إلى السجل الثالث.



حسنا اذن. لقد حددنا رقمًا ذا نقطة ثابتة ، ولكن لا تزال هناك مشكلة واحدة. الأعداد السالبة. في سلسلة تايلور ، يتم التوسع بعلامة متغيرة. للتعامل مع هذه المشكلة ، نقوم بتحويل صيغة حساب كثير الحدود بطريقة Goner ، إلى الشكل التالي:

$$

$$\sin(x)\approx x(1-x^2(1/3!-x^2(1/5!-x^2(1/7!-x^2\cdot1/9!))))$$

تأكد من أنه رياضيًا مطابقًا تمامًا للصيغة الأصلية. لكن في كل قوس يوجد رقم مثل$$$1/(2n + 1)! - x^2\cdot(\cdots)$دائما إيجابية. أولئك. يسمح هذا التحويل بتقييم التعبير كأعداد صحيحة بدون إشارة.

constexpr mynumber toint    = {{0x00000000, 0x43F00000}};  /*  18446744073709551616 = 2^64     */
constexpr mynumber todouble = {{0x00000000, 0x3BF00000}};  /*  ~5.42101086242752217003726400434E-20 = 2^-64     */

double sin_e7(double xd) {
  uint64_t x = xd * toint.x;
  uint64_t xx = mul2(x, x);
  uint64_t res = tsx[19]; 
  for(int i = 17; i >= 3; i -= 2) {
    res = tsx[i] - mul2(res, xx);
  }
  res = mul2(res, xx);
  res = x - mul2(x, res);
  return res * todouble.x;
}

أين $$$tsx[i]=1/i!$بتنسيق نقطة ثابتة. هذه المرة ، للراحة ، قمت بنشر كل التعليمات البرمجية على fast_sine github ، وتخلصت من quadmath للتوافق مع clang and arm. وقمت بتغيير طريقة حساب الخطأ قليلاً.



يتم عرض مقارنة بين دالة الجيب القياسية ودالة النقطة الثابتة في الجدولين أدناه. يوضح الجدول الأول دقة الحساب (وهي نفسها تمامًا بالنسبة إلى x86_64 و ARM). الجدول الثاني هو مقارنة الأداء.







أثناء الاختبار ، تم حساب قيمة الجيب "الحقيقية" باستخدام مكتبة MPFR... تم اعتبار الحد الأقصى لقيمة ULP على أنه أقصى انحراف عن القيمة "الحقيقية". النسبة المئوية للأخطاء - عدد الحالات التي لم تتطابق فيها القيمة المحسوبة لدالة الجيب بواسطتنا أو بواسطة libm مع قيمة الجيب المقرب إلى الضعف. يُظهر متوسط ​​قيمة الانحراف "اتجاه" خطأ الحساب: المبالغة في تقدير القيمة أو التقليل من شأنها. كما ترى من الجدول ، تميل وظيفتنا إلى المبالغة في تقدير قيم الجيب. يمكن إصلاح هذا! من قال إن قيم tsx يجب أن تكون مساوية تمامًا لمعاملات سلسلة تايلور. تقترح فكرة واضحة إلى حد ما نفسها لتغيير قيم المعاملات من أجل تحسين دقة التقريب وإزالة المكون الثابت للخطأ. من الصعب إجراء مثل هذا الاختلاف بشكل صحيح. ولكن يمكننا محاولة. لنأخذ على سبيل المثالالقيمة الرابعة من صفيف معاملات tsx (tsx [3]) وتغيير الرقم الأخير a إلى f. دعونا نعيد تشغيل البرنامج ونرى الدقة (sin_e7a). انظروا انها قليلا لكنها زادت! نضيف هذه الطريقة إلى حصالةنا.



الآن دعونا نرى ما هو الأداء. للاختبار ، أخذت ما كان موجودًا في متناول اليد i5 mobile و raspberry الرابع (Raspberry PI 4 8GB) ، GCC10 من توزيع Ubuntu 20.04 x64 لكلا النظامين.







أنا لا أدعي أن أكون أكثر دقة في هذه القياسات. يمكن إجراء اختلافات تصل إلى عدة عشرات بالمائة اعتمادًا على حمل المعالج. يمكن إجراء الاستنتاج الرئيسي على هذا النحو. لا يؤدي التبديل إلى الحساب الصحيح إلى زيادة الأداء في المعالجات الحديثة ² . إن العدد الذي لا يمكن تصوره من الترانزستورات في المعالجات الحديثة يجعل من الممكن إجراء حسابات معقدة بسرعة. لكن ، أعتقد أنه في معالجات مثل Intel Atom ، وكذلك على وحدات التحكم الضعيفة ، يمكن أن يوفر هذا الأسلوب مكاسب كبيرة في الأداء. ما رأيك؟



على الرغم من أن هذا الأسلوب أدى إلى زيادة الدقة ، فإن هذا الاكتساب في الدقة يبدو أكثر إثارة للاهتمام من كونه مفيدًا. من حيث الأداء ، يمكن لهذا النهج أن يجد نفسه ، على سبيل المثال ، في إنترنت الأشياء. ولكن بالنسبة للحوسبة عالية الأداء ، لم تعد سائدة. في عالم اليوم ، تفضل SSE / AVX / CUDA استخدام حساب الوظيفة المتوازية. وفي حساب الفاصلة العائمة. لا توجد نظائر متوازية لوظيفة MUL. الوظيفة نفسها هي بالأحرى تكريم للتقاليد.



في الفصل التالي ، سأصف كيف يمكنك استخدام AVX بشكل فعال لإجراء العمليات الحسابية. مرة أخرى ، دعنا ندخل في كود libm ونحاول تحسينه.



¹ لا توجد سجلات بهذه الأسماء في معالجات معروفة لي. تم اختيار الأسماء على سبيل المثال.

²وتجدر الإشارة هنا إلى أن جهاز ARM الخاص بي مجهز بأحدث إصدار من المعالج الرياضي. إذا قام المعالج بمحاكاة حسابات الفاصلة العائمة ، فقد تكون النتائج مختلفة بشكل كبير.