Databases
StackOverflow2013
Postsdbo
PO
 

Note that there are some explanatory texts on larger screens.

plurals
  1. PO
    primarykey
    data
    text
    <p>I could not resist spending an hour on your problem...</p> <p>This algorithm is described in section 5.5.2 of "Arithmetique des ordinateurs" by Jean-Michel Muller (in french). It is actually a special case of Newton iterations with 1 as starting point. The book gives a simple formulation of the algorithm to compute N/D, with D normalized in range [1/2,1[:</p> <pre><code>e = 1 - D Q = N repeat K times: Q = Q * (1+e) e = e*e </code></pre> <p>The number of correct bits doubles at each iteration. In the case of 32 bits, 4 iterations will be enough. You can also iterate until <code>e</code> becomes too small to modify <code>Q</code>.</p> <p>Normalization is used because it provides the max number of significant bits in the result. It is also easier to compute the error and number of iterations needed when the inputs are in a known range.</p> <p>Once your input value is normalized, you don't need to bother with the value of BASE until you have the inverse. You simply have a 32-bit number X normalized in range 0x80000000 to 0xFFFFFFFF, and compute an approximation of Y=2^64/X (Y is at most 2^33).</p> <p>This simplified algorithm may be implemented for your Q22.10 representation as follows:</p> <pre><code>// Fixed point inversion // EB Apr 2010 #include &lt;math.h&gt; #include &lt;stdio.h&gt; // 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&lt;&lt;BASE); } uint32 toFP(double x) { return (int)floor(0.5+x*(1&lt;&lt;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 &amp; 0x80000000) == 0 ) { nfp &lt;&lt;= 1; shl++; } // use "clz" instead uint64 q = 0x100000000ULL; // 2^32 uint64 e = 0x100000000ULL - (uint64)nfp; // 2^32-NFP int i; for (i=0;i&lt;4;i++) // iterate { // Both multiplications are actually // 32x32 bits truncated to the 32 high bits q += (q*e)&gt;&gt;(uint64)32; e = (e*e)&gt;&gt;(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&lt;=NFP&lt;2^32 and 2^32&lt;Q&lt;=2^33 return (uint32)(q&gt;&gt;(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); } </code></pre> <p>As noted in the code, the multiplications are not full 32x32->64 bits. E will become smaller and smaller and fits initially on 32 bits. Q will always be on 34 bits. We take only the high 32 bits of the products.</p> <p>The derivation of <code>64-2*BASE-shl</code> is left as an exercise for the reader :-). If it becomes 0 or negative, the result is not representable (the input value is too small).</p> <p>EDIT. As a follow-up to my comment, here is a second version with an implicit 32-th bit on Q. Both E and Q are now stored on 32 bits:</p> <pre><code>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 &amp; 0x80000000) == 0 ) { nfp &lt;&lt;= 1; shl++; } // use "clz" instead int shr = 64-2*BASE-shl; // normalization shift for Q if (shr &lt;= 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&lt;3;i++) // iterate { e = (e*e)&gt;&gt;(uint64)32; q += e + ((q*e)&gt;&gt;(uint64)32); } return (uint32)(q&gt;&gt;shr) + (1&lt;&lt;(32-shr)); // insert implicit bit } </code></pre>
    singulars
    1. This table or related slice is empty.
    1. POPicking good first estimates for Goldschmidt division
      singulars
      1. PTQuestion
    1. PTAnswer
    1. USEric Bainville
    1. USEric Bainville
    plurals
    1. This table or related slice is empty.
    1. This table or related slice is empty.
    1. POPicking good first estimates for Goldschmidt division
      singulars
      1. PTQuestion
    1. This table or related slice is empty.
    1. VO
      singulars
      1. PO
      1. This table or related slice is empty.
      1. VTUpMod
    2. VO
      singulars
      1. PO
      1. This table or related slice is empty.
      1. VTUpMod
    3. VO
      singulars
      1. PO
      1. This table or related slice is empty.
      1. VTUpMod
    1. COBut in any event, because of the estimate `2 - D` one have to reserve two bits, right? And there is no way to get a smaller estimate than ~1.111111111 ?
      singulars
      1. PO
      1. This table or related slice is empty.
    2. COHere's your well-deserved bounty by the way.
      singulars
      1. PO
      1. This table or related slice is empty.
    3. COThanks. All values of E fit on 32 bits. You don't need extra bits to store E. Q starts at 2^32, and may actually reach 2^33 in the case where NFP is 0x80000000. Since this case can be handled separately, Q will always fit on 33 bits, with bit 32 always set. So you can actually represent Q on 32 bits with an implicit bit. I forgot to add that 64-2*BASE-shl can be tested before iterating. Another point: when the initial FP uses a small number of bits, you probably don't need 32 significant bits (but with no meaning) in Q.
      singulars
      1. PO
      1. USEric Bainville
 

Querying!

 
Guidance
A row detail

Detail views are divided into sections. All the information in the data section comes from columns in the selected row. The other sections display data from other, related rows.

Related data can be related in a to-one or a to-many fashion. Captions of data related in a to-many fashion link to a list view showing a filtered view of the table.

Try moving around until you find a non-empty to-many entry and click on the label to get to one. You can move back to the root by clicking on the database name in the header.

SQuiL has stopped working due to an internal error.

If you are curious you may find further information in the browser console, which is accessible through the devtools (F12).

Reload