اختيار التقديرات الأولى الجيدة لقسم جولدشميدت

Question

أقوم بحساب المعاملة المتبادلة الثابتة في Q22.10 مع قسم جولدشميت للاستخدام في برمجيات بلدي على الذراع.

يتم ذلك بمجرد تعيين البسط على 1 ، أي أن البسط يصبح العددية في التكرار الأول. أن نكون صادقين ، أنا نوع من اتباع خوارزمية ويكيبيديا عمياء هنا. تقول المقالة إنه إذا تم تحجيم المقام في نطاق نصف المفتوح (0.5 ، 1.0] ، يمكن أن يستند التقدير الأول الجيد إلى المقام وحده: دع F هو العددية المقدرة و D يكون المقام ، ثم F = 2 - د.

ولكن عند القيام بذلك ، أفقد الكثير من الدقة. قل إذا كنت أرغب في العثور على المتبادل 512.00002f. من أجل خفض الرقم ، أفقد 10 بت من الدقة في جزء الكسر ، الذي تم نقله. لذلك ، أسئلتي هي:

• هل هناك طريقة لاختيار تقدير أفضل لا يتطلب التطبيع؟ لماذا ا؟ لما لا؟ دليل رياضي على سبب هذا أو غير ممكن سيكون رائعًا.
• أيضًا ، هل من الممكن تحويلي التقديرات الأولى حتى تتقارب السلسلة بشكل أسرع؟ الآن ، يتقارب بعد التكرار الرابع في المتوسط. في ARM ، يكون هذا حوالي 50 دورة أسوأ حالة ، وهذا لا يأخذ محاكاة CLZ/BSR في الاعتبار ، ولا بحث الذاكرة. إذا كان ذلك ممكنًا ، أود أن أعرف ما إذا كان القيام بذلك يزيد من الخطأ ، ومقدار المبلغ.

ها هو testcase الخاص بي. ملاحظة: تنفيذ البرنامج لـ clz على السطر 13 من رسالتي هنا. يمكنك استبداله بأمر جوهري إذا أردت. clz يجب إرجاع عدد الأصفار الرائدة ، و 32 للقيمة 0.

#include <stdio.h>
#include <stdint.h>

const unsigned int BASE = 22ULL;

static unsigned int divfp(unsigned int val, int* iter)
{
  /* Numerator, denominator, estimate scalar and previous denominator */
  unsigned long long N,D,F, DPREV;
  int bitpos;

  *iter = 1;
  D = val;
  /* Get the shift amount + is right-shift, - is left-shift. */
  bitpos = 31 - clz(val) - BASE;
  /* Normalize into the half-range (0.5, 1.0] */
  if(0 < bitpos)
    D >>= bitpos;
  else
    D <<= (-bitpos);

  /* (FNi / FDi) == (FN(i+1) / FD(i+1)) */
  /* F = 2 - D */
  F = (2ULL<<BASE) - D;
  /* N = F for the first iteration, because the numerator is simply 1.
     So don't waste a 64-bit UMULL on a multiply with 1 */
  N = F;
  D = ((unsigned long long)D*F)>>BASE;

  while(1){
    DPREV = D;
    F = (2<<(BASE)) - D;
    D = ((unsigned long long)D*F)>>BASE;
    /* Bail when we get the same value for two denominators in a row.
      This means that the error is too small to make any further progress. */
    if(D == DPREV)
      break;
    N = ((unsigned long long)N*F)>>BASE;
    *iter = *iter + 1;
  }
  if(0 < bitpos)
    N >>= bitpos;
  else
    N <<= (-bitpos);
  return N;
}


int main(int argc, char* argv[])
{
  double fv, fa;
  int iter;
  unsigned int D, result;

  sscanf(argv[1], "%lf", &fv);

  D = fv*(double)(1<<BASE);
  result = divfp(D, &iter); 

  fa = (double)result / (double)(1UL << BASE);
  printf("Value: %8.8lf 1/value: %8.8lf FP value: 0x%.8X\n", fv, fa, result);
  printf("iteration: %d\n",iter);

  return 0;
}

Answer 1

لم أستطع مقاومة قضاء ساعة في مشكلتك ...

تم وصف هذه الخوارزمية في القسم 5.5.2 من "Arithmetique des indinators" لجان ميشيل مولر (باللغة الفرنسية). إنها في الواقع حالة خاصة لتكرارات نيوتن مع 1 كنقطة انطلاق. يعطي الكتاب صياغة بسيطة للخوارزمية لحساب N/D ، مع D تطبيعها في المدى [1/2،1 [:

e = 1 - D
Q = N
repeat K times:
  Q = Q * (1+e)
  e = e*e

عدد البتات الصحيحة يتضاعف في كل تكرار. في حالة 32 بت ، سيكون 4 تكرارات كافية. يمكنك أيضا التكرار حتى e يصبح صغيرًا جدًا بحيث لا يمكن تعديله Q.

يتم استخدام التطبيع لأنه يوفر عدد أقصى البتات المهمة في النتيجة. من الأسهل أيضًا حساب الخطأ وعدد التكرارات اللازمة عندما تكون المدخلات في نطاق معروف.

بمجرد تطبيع قيمة الإدخال الخاصة بك ، لا تحتاج إلى الإزعاج بقيمة القاعدة حتى تحصل على العكسي. لديك ببساطة رقم 32 بت X تطبيع في النطاق 0x80000000 إلى 0xffffffff ، وحساب تقريب y = 2^64/x (y على الأكثر 2^33).

يمكن تنفيذ هذه الخوارزمية المبسطة لتمثيل Q22.10 الخاص بك على النحو التالي:

// Fixed point inversion
// EB Apr 2010

#include <math.h>
#include <stdio.h>

// Number X is represented by integer I: X = I/2^BASE.
// We have (32-BASE) bits in integral part, and BASE bits in fractional part
#define BASE 22
typedef unsigned int uint32;
typedef unsigned long long int uint64;

// Convert FP to/from double (debug)
double toDouble(uint32 fp) { return fp/(double)(1<<BASE); }
uint32 toFP(double x) { return (int)floor(0.5+x*(1<<BASE)); }

// Return inverse of FP
uint32 inverse(uint32 fp)
{
  if (fp == 0) return (uint32)-1; // invalid

  // Shift FP to have the most significant bit set
  int shl = 0; // normalization shift
  uint32 nfp = fp; // normalized FP
  while ( (nfp & 0x80000000) == 0 ) { nfp <<= 1; shl++; } // use "clz" instead

  uint64 q = 0x100000000ULL; // 2^32
  uint64 e = 0x100000000ULL - (uint64)nfp; // 2^32-NFP
  int i;
  for (i=0;i<4;i++) // iterate
    {
      // Both multiplications are actually
      // 32x32 bits truncated to the 32 high bits
      q += (q*e)>>(uint64)32;
      e = (e*e)>>(uint64)32;
      printf("Q=0x%llx E=0x%llx\n",q,e);
    }
  // Here, (Q/2^32) is the inverse of (NFP/2^32).
  // We have 2^31<=NFP<2^32 and 2^32<Q<=2^33
  return (uint32)(q>>(64-2*BASE-shl));
}

int main()
{
  double x = 1.234567;
  uint32 xx = toFP(x);
  uint32 yy = inverse(xx);
  double y = toDouble(yy);

  printf("X=%f Y=%f X*Y=%f\n",x,y,x*y);
  printf("XX=0x%08x YY=0x%08x XX*YY=0x%016llx\n",xx,yy,(uint64)xx*(uint64)yy);
}

كما هو مذكور في الكود ، فإن الضربات ليست كاملة 32x32-> 64 بت. سوف تصبح E أصغر وأصغر وتناسب في البداية على 32 بت. س ستكون دائما على 34 بت. نحن نأخذ فقط 32 بت من المنتجات.

اشتقاق 64-2*BASE-shl يترك كتمرين للقارئ :-). إذا أصبحت 0 أو سلبية ، فإن النتيجة غير ممثلة (قيمة الإدخال صغيرة جدًا).

تعديل. كمتابعة لتعليقي ، إليك نسخة ثانية ذات 32 بت ضمنية على Q.

uint32 inverse2(uint32 fp)
{
  if (fp == 0) return (uint32)-1; // invalid

  // Shift FP to have the most significant bit set
  int shl = 0; // normalization shift for FP
  uint32 nfp = fp; // normalized FP
  while ( (nfp & 0x80000000) == 0 ) { nfp <<= 1; shl++; } // use "clz" instead
  int shr = 64-2*BASE-shl; // normalization shift for Q
  if (shr <= 0) return (uint32)-1; // overflow

  uint64 e = 1 + (0xFFFFFFFF ^ nfp); // 2^32-NFP, max value is 2^31
  uint64 q = e; // 2^32 implicit bit, and implicit first iteration
  int i;
  for (i=0;i<3;i++) // iterate
    {
      e = (e*e)>>(uint64)32;
      q += e + ((q*e)>>(uint64)32);
    }
  return (uint32)(q>>shr) + (1<<(32-shr)); // insert implicit bit
}

Answer 2

بعض الأفكار بالنسبة لك ، على الرغم من أن لا شيء يحل مشكلتك مباشرة كما هو مذكور.

1. لماذا هذا algo للقسمة؟ معظم الانقسامات التي رأيتها في الذراع تستخدم بعض المتغيرات

   
         adcs hi, den, hi, lsl #1
         subcc hi, hi, den
         adcs lo, lo, lo
   

تتكرر N بتات مع بحث ثنائي من CLZ لتحديد من أين تبدأ. هذا سريع جدا.

2. إذا كانت الدقة مشكلة كبيرة ، فأنت لا تقتصر على 32/64 بت لتقديم تمثيل النقطة الثابتة الخاصة بك. سيكون الأمر أبطأ قليلاً ، ولكن يمكنك القيام بإضافة/ADC أو Sub/SBC لنقل القيم عبر السجلات. تم تصميم Mul/MLA أيضًا لهذا النوع من العمل.

مرة أخرى ، لا توجد إجابات مباشرة لك ، ولكن ربما بعض الأفكار للمضي قدمًا في ذلك. من المحتمل أن تساعدني رؤية رمز الذراع الفعلي قليلاً أيضًا.

Answer 3

مادز ، أنت لا تفقد أي دقة على الإطلاق. عندما تقسم 512.00002f على 2^10 ، فأنت فقط تقلل من الأسس من رقم النقطة العائمة بمقدار 10. يبقى Mantissa كما هو. بالطبع ما لم يضرب الأسعار الحد الأدنى لقيمته ولكن لا ينبغي أن يحدث ذلك لأنك تتوسع إلى (0.5 ، 1].

تحرير: حسنًا ، لذا فأنت تستخدم نقطة عشرية ثابتة. في هذه الحالة ، يجب أن تسمح بتمثيل مختلف للمقام في الخوارزمية. قيمة D من (0.5 ، 1] ليس فقط في البداية ولكن طوال الحساب بأكمله (من السهل إثبات أن x * (2-x) <1 لـ x <1). لذلك يجب أن تمثل القاسم مع العشرية العشرية نقطة في القاعدة = 32. بهذه الطريقة سيكون لديك 32 بت من الدقة طوال الوقت.

تحرير: لتنفيذ هذا ، يجب عليك تغيير الأسطر التالية من الكود الخاص بك:

  //bitpos = 31 - clz(val) - BASE;
  bitpos = 31 - clz(val) - 31;
...
  //F = (2ULL<<BASE) - D;
  //N = F;
  //D = ((unsigned long long)D*F)>>BASE;
  F = -D;
  N = F >> (31 - BASE);
  D = ((unsigned long long)D*F)>>31;
...
    //F = (2<<(BASE)) - D;
    //D = ((unsigned long long)D*F)>>BASE;
    F = -D;
    D = ((unsigned long long)D*F)>>31;
...
    //N = ((unsigned long long)N*F)>>BASE;
    N = ((unsigned long long)N*F)>>31;

أيضًا في النهاية ، سيتعين عليك تغيير N ليس عن طريق BitPOs ولكن بعض القيمة المختلفة التي أنا كسول جدًا لأكتشفها الآن :).