复数的有效幅值比较

import array as arr

#================================================
def Main():

#---- Initialize the Counters

        theCount      = 0
        theWrongCount = 0

        thePresumedZeroCount    = 0
        theCloseButNotZeroCount = 0

        theTestExits = arr.array( "i", [ 0, 0, 0, 0, 0, 0, 0 ] )

#---- Test on a Swept Area

        theLimit = 32
        theLimitSquared = theLimit * theLimit

        theWorstHigh = 1.0
        theWorstLow  = 1.0

        for i1 in range( theLimit ):
          ii1 = i1 * i1
          for q1 in range( theLimit ):
            m1 = ii1 + q1 * q1
            for i2 in range( theLimit ):
              ii2 = i2 * i2
              for q2 in range( theLimit ):
                m2 = ii2 + q2 * q2
                D  = m1 - m2

                theCount += 1

                c, t = CompareMags( i1, q1, i2, q2 )

                if t <= 6:
                   theTestExits[t] += 1

                if c == 2:

                   thePresumedZeroCount += 1
                   if D != 0:
                      theCloseButNotZeroCount += 1

                      Q = float( m1 ) / float( m2 )
                      if Q > 1.0:
                         if theWorstHigh < Q:
                            theWorstHigh = Q
                      else:
                         if theWorstLow > Q:
                            theWorstLow = Q

                      print "%2d  %2d   %2d  %2d   %10.6f" % ( i1, q1, i2, q2, Q )

                elif c == 1:
                   if D <= 0:
                      theWrongCount += 1
                      print "Wrong Less ", i1, q1, i2, q2, D, c
                elif c == 0:
                   if D != 0:
                      theWrongCount += 1
                      print "Wrong Equal", i1, q1, i2, q2, D, c
                elif c == -1:
                   if D >= 0:
                      theWrongCount += 1
                      print "Wrong Great", i1, q1, i2, q2, D, c
                else:
                   theWrongCount += 1
                   print "Invalid c value:", i1, q1, i2, q2, D, c


#---- Calculate the Results

        theSureCount   = ( theCount - thePresumedZeroCount )                    
        theSurePercent = 100.0 * theSureCount / theCount

        theCorrectPercent = 100.0 * ( theSureCount - theWrongCount ) \
                          / theSureCount

        if thePresumedZeroCount > 0:
           theCorrectPresumptionPercent = 100.0 * ( thePresumedZeroCount - theCloseButNotZeroCount ) \
                                        / thePresumedZeroCount
        else:
           theCorrectPresumptionPercent = -1

        theFaultyPercent = 100.0 * theCloseButNotZeroCount / theCount

#---- Report the Results

        print
        print "   Count:", theCount
        print
        print "    Sure:", theSurePercent
        print " Correct:", theCorrectPercent
        print
        print "Presumed:", 100 - theSurePercent
        print "Actually:", theCorrectPresumptionPercent
        print
        print "  Faulty:", theFaultyPercent
        print
        print "    High:", theWorstHigh
        print "     Low:", theWorstLow
        print

#---- Report The Cutoff Values

        pct = 0.0
        f   = 100.0 / theCount

        for t in range( 1, 7 ):
          pct += f * theTestExits[t]
          print "%d %8d  %6.2f" % ( t, theTestExits[t], pct )

        print

#================================================
def CompareMags( I1, Q1, I2, Q2 ):

# This function compares the magnitudes of two 
# integer points and returns a comparison result value
#
# Returns  ( c, t )
#
#    c   Comparison
#
#   -1     | (I1,Q1) |  <  | (I2,Q2) |
#    0     | (I1,Q1) |  =  | (I2,Q2) |
#    1     | (I1,Q1) |  >  | (I2,Q2) |
#    2     | (I1,Q1) | ~=~ | (I2,Q2) |
#
#    t   Exit Test
#
#    1     Primary Determination
#    2     D/S Centers are aligned
#    3     Obvious Answers
#    4     Trivial Matching Gaps
#    5     Opposite Gap Sign Cases
#    6     Same Gap Sign Cases
#   10     Small Angle + Count
#   20     CORDIC + Count
#
# It does not matter if the arguments represent vectors 
# or complex numbers.  Nor does it matter if the calling
# routine considers the integers as fixed point values.


#---- Ensure the Points are in the First Quadrant WLOG

        a1 = abs( I1 )
        b1 = abs( Q1 )

        a2 = abs( I2 )
        b2 = abs( Q2 )

#---- Ensure they are in the Lower Half (First Octant) WLOG

        if b1 > a1:
           a1, b1 = b1, a1

        if b2 > a2:
           a2, b2 = b2, a2

#---- Primary Determination

        if a1 > a2:
           if a1 + b1 >= a2 + b2:
              return 1, 1
           else:
              thePresumedResult = 1
              da = a1 - a2
              sa = a1 + a2
              db = b2 - b1
              sb = b2 + b1
        elif a1 < a2:
           if a1 + b1 <= a2 + b2:
              return -1, 1
           else:
              thePresumedResult = -1
              da = a2 - a1
              sa = a2 + a1
              db = b1 - b2
              sb = b1 + b2
        else:
           if b1 > b2:
              return 1, 1
           elif b1 < b2:
              return -1, 1
           else:
              return 0, 1

#---- Bring Factors into 1/2 to 1 Ratio Range

        db, sb = sb, db

        while da < sa:
            da += da
            sb += sb
            if sb > db:
               db, sb = sb, db

#---- Ensure the [b] Factors are Both Even or Odd

        if ( ( sb + db ) & 1 ) > 0:
           da += da
           sa += sa
           db += db
           sb += sb

#---- Calculate Arithmetic Mean and Radius of [b] Factors

        p = ( db + sb ) >> 1
        r = ( db - sb ) >> 1

#---- Calculate the Gaps from the [b] mean and [a] values

        g = da - p
        h = p - sa

#---- If the mean of [b] is centered in (the mean of) [a]

        if g == h:
           if g == r:
              return 0, 2;
           elif g > r:
              return -thePresumedResult, 2
           else:
              return thePresumedResult, 2

#---- Weed Out the Obvious Answers

        if g > h:
           if r > g and r > h:
              return thePresumedResult, 3
        else:             
           if r < g and r < h:
              return -thePresumedResult, 3

#---- Calculate Relative Gaps

        vg = g - r
        vh = h - r

#---- Handle the Trivial Matching Gaps

        if vg == 0:
           if vh > 0:
              return -thePresumedResult, 4
           else:
              return thePresumedResult, 4

        if vh == 0:
           if vg > 0:
              return thePresumedResult, 4
           else:
              return -thePresumedResult, 4


#---- Handle the Gaps with Opposite Sign Cases

        if vg < 0:
           if vh > 0:
              return -thePresumedResult, 5
        else:
           if vh < 0:
              return thePresumedResult, 5

#---- Handle the Gaps with the Same Sign (using numerators)

        theSum = da + sa

        if g < h:
           theBound = ( p << 4 ) - p  
           theMid   = theSum << 3

           if theBound > theMid:
              return -thePresumedResult, 6
        else:
           theBound = ( theSum << 4 ) - theSum
           theMid   = p << 5

           if theBound > theMid:
              return thePresumedResult, 6

#---- Return to IQ Space under XY Names

        x1 = theSum
        x2 = da - sa

        y2 = db + sb
        y1 = db - sb

#---- Ensure Points are in Lower First Quadrant (First Octant)

        if x1 < y1:
           x1, y1 = y1, x1

        if x2 < y2:
           x2, y2 = y2, x2

#---- Variation of Olli's CORDIC to Finish

        for theTryLimit in range( 10 ):
            c, x1, y1, x2, y2 = Iteration( x1, y1, x2, y2, thePresumedResult )
            if c != 2:
               break

            if theTryLimit > 3:   
               print "Many tries needed!", theTryLimit, x1, y1, x2, y2

        return c, 20

#================================================
def Iteration( x1, y1, x2, y2, argPresumedResult ):

#---- Try to reduce the Magnitudes

        while ( x1 & 1 ) == 0 and \
              ( y1 & 1 ) == 0 and \
              ( x2 & 1 ) == 0 and \
              ( y2 & 1 ) == 0:
            x1 >>= 1
            y1 >>= 1
            x2 >>= 1
            y2 >>= 1

#---- Set the Perpendicular Values (clockwise to downward)

        dx1 =  y1
        dy1 = -x1

        dx2 =  y2
        dy2 = -x2

        sdy = dy1 + dy2

#---- Allocate the Arrays for Length Storage

        wx1 = arr.array( "i" )
        wy1 = arr.array( "i" )
        wx2 = arr.array( "i" )
        wy2 = arr.array( "i" )

#---- Locate the Search Range

        thePreviousValue = x1 + x2  # Guaranteed Big Enough

        for theTries in range( 10 ):
            wx1.append( x1 )
            wy1.append( y1 )
            wx2.append( x2 )
            wy2.append( y2 )

            if x1 > 0x10000000 or x2 > 0x10000000:
               print "Danger, Will Robinson!"
               break

            theValue = abs( y1 + y2 + sdy )

            if theValue > thePreviousValue:
               break

            thePreviousValue = theValue

            x1 += x1
            y1 += y1
            x2 += x2
            y2 += y2

#---- Prepare for the Search

        theTop = len( wx1 ) - 1

        thePivot = theTop - 1

        x1 = wx1[thePivot]
        y1 = wy1[thePivot]
        x2 = wx2[thePivot]
        y2 = wy2[thePivot]

        theValue = abs( y1 + y2 + sdy )

#---- Binary Search

        while thePivot > 0:
            thePivot -= 1

            uy1 = y1 + wy1[thePivot]
            uy2 = y2 + wy2[thePivot]

            theUpperValue = abs( uy1 + uy2 + sdy )

            ly1 = y1 - wy1[thePivot]
            ly2 = y2 - wy2[thePivot]

            theLowerValue = abs( ly1 + ly2 + sdy )

            if theUpperValue < theLowerValue:
               if theUpperValue < theValue:
                  x1 += wx1[thePivot]
                  x2 += wx2[thePivot]
                  y1  = uy1
                  y2  = uy2

                  theValue = theUpperValue

            else:
               if theLowerValue < theValue:
                  x1 -= wx1[thePivot]
                  x2 -= wx2[thePivot]
                  y1  = ly1
                  y2  = ly2

                  theValue = theLowerValue

#---- Apply the Rotation

        x1 += dx1
        y1 += dy1

        x2 += dx2
        y2 += dy2

#---- Bounce Points Below the Axis to Above

        if y1 < 0:
           y1 = -y1

        if y2 < 0:
           y2 = -y2

#---- Comparison Determination

        c = 2

        if x1 > x2:
           if x1 + y1 >= x2 + y2:
              c = argPresumedResult
        elif x1 < x2:
           if x1 + y1 <= x2 + y2:
              c = -argPresumedResult
        else:
           if y1 > y2:
              c = argPresumedResult
           elif y1 < y2:
              c = -argPresumedResult
           else:
              c =  0

#---- Exit

        return c, x1, y1, x2, y2

#================================================
def MyVersionOfOllis( I1, Q1, I2, Q2 ):

# Returns  ( c, t )
#
#    c   Comparison
#
#   -1     | (I1,Q1) |  <  | (I2,Q2) |
#    1     | (I1,Q1) |  >  | (I2,Q2) |
#
#    t   Exit Test
#
#    1     (Partial) Primary Determination
#    2     CORDIC Loop + 1
#    6     Terminating Guess

#---- Set Extent Parameter

        maxIterations = 4

#---- Ensure the Points are in the First Quadrant WLOG

        I1 = abs( I1 )
        Q1 = abs( Q1 )

        I2 = abs( I2 )
        Q2 = abs( Q2 )

#---- Ensure they are in the Lower Half (First Octant) WLOG

        if Q1 > I1:
           I1, Q1 = Q1, I1

        if Q2 > I2:
           I2, Q2 = Q2, I2

#---- (Partial) Primary Determination

        if I1 < I2 and  I1 + Q1 <= I2 + Q2:
           return -1, 1

        if I1 > I2 and  I1 + Q1 >= I2 + Q2:
           return 1, 1

#---- CORDIC Loop

        Q1pow1 = Q1 >> 1
        I1pow1 = I1 >> 1
        Q2pow1 = Q2 >> 1
        I2pow1 = I2 >> 1

        Q1pow2 = Q1 >> 3
        I1pow2 = I1 >> 3
        Q2pow2 = Q2 >> 3
        I2pow2 = I2 >> 3

        for n in range ( 1, maxIterations+1 ):
            newI1 = I1 + Q1pow1
            newQ1 = Q1 - I1pow1
            newI2 = I2 + Q2pow1
            newQ2 = Q2 - I2pow1

            I1 = newI1
            Q1 = abs( newQ1 )
            I2 = newI2
            Q2 = abs( newQ2 )

            if I1 <= I2 - I2pow2:
               return -1, 1 + n

            if I2 <= I1 - I1pow2:
               return 1, 1 + n

            Q1pow1 >>= 1
            I1pow1 >>= 1
            Q2pow1 >>= 1
            I2pow1 >>= 1

            Q1pow2 >>= 2
            I1pow2 >>= 2
            Q2pow2 >>= 2
            I2pow2 >>= 2

#---- Terminating Guess

        Q1pow1 <<= 1
        Q2pow1 <<= 1

        if I1 + Q1pow1 < I2 + Q2pow1:
           return -1, 6
        else:   
           return 1, 6

#================================================
Main()

Swap db, sb   # d is now the larger quantity

While da < sa
  da =<< 1
  sb =<< 1
  If sb > db Then Swap db, sb # s is the smaller quantity
EndWhile

$ p如果$ d_b = 2 s_b $，则为/ 3 $。

让
$$
\开始{aligned}
g＆= d_a-p \\
h＆= p-s_a \\
\ end {aligned}
$$

将它们插入可能已经翻倍的$ D $公式。

$$
\ begin {aligned}
D 2 ^ k＆=（p + g）（ph）-（p + r）（pr）\\
＆= [p ^ 2 +（gh ）p-gh]-[p ^ 2-r ^ 2] \\
＆=（gh）p + [r ^ 2- gh] \\
\ end {aligned}
$$

如果$ g = h $则很简单：如果$ r = g $相等。如果$ r> g $则1较大，否则2较大。

让
$$
\ begin {aligned}
v_g＆= g-r \\
v_h＆= h-r \\
\ end {aligned}

求值$ D2 ^ k的RHS的两个项$ equation。

$$
\ begin {aligned}
r ^ 2-gh＆= r ^ 2-（r + v_g）（r + v_h）\\
＆= -v_g v_h-（v_g + v_h）r \\
\ end {aligned}
$$

and

$$ g- h = v_g-v_h $$

插入它们。

$$
\ begin {aligned}
D 2 ^ k＆=（gh ）p + [r ^ 2- gh] \\
＆=（v_g-v_h）p-v_g v_h-（v_g + v_h）r \\
＆= v_g（pr）-v_h（p + r）-v_g v_h \\
＆= v_g s_b-v_h d_b-\ left（\ frac {v_h v_g} {2} + \ frac {v_h v_g} {2} \ right）\\
＆= v_g（s_b- \ frac {v_h} {2}）-v_h（d_b + \ frac {v_g} {2}）\\
\ end {aligned}
$$

将两边都乘以2以除去分数。

$$
\ begin {aligned}
D 2 ^ {k + 1}＆= v_g （2s_b-v_h）-v_h（2d_b + v_g）\\
\ end {aligned}
$$

如果$ v_g $或$ v_h $为零，则

同样，如果$ v_g $和$ v_h $具有相反的符号，则D的符号确定也很小。

仍然在最后一个工作sliver ......


非常接近。

theHandledPercent 98.6582716049

theCorrectPercent 100.0

稍​​后再继续.......欢迎任何人为相同的标志盒找到正确的处理逻辑。


前一天，又迈出了一大步。

原始符号确定方程可以这样分解：

$$
\开始{对齐}
D＆= d_a s_a-d_b s_b \\
＆= \ left（\ sqrt {d_a s_a}-\ sqrt {d_b s_b} \ right）\ left（\ sqrt {d_a s_a} + \ sqrt {d_b s_b} \ right）\\
\ end {aligned}
$$

总和始终为正，因此不会影响D的符号。差异因子是两个几何平均值的差。

几何平均值可以用算术平均值来近似。这是“ alpha max加beta min算法”背后的工作原理。算术平均值也是几何平均值的上限。

由于$ s $值的下限是$ d / 2 $，因此可以为几何平均值建立一个大致的下限，而无需进行大量计算。
＆\ ge \ frac {d} {2} \\
\ sqrt {ds}＆\ ge \ frac {d} {\ sqrt {2}} \\
＆= \ frac {\ frac {d} {\ sqrt {2}}} {{d + s）/ 2} \ cdot \ frac { d + s} {2} \\
＆= \ frac {\ sqrt {2}} {1 + s / d} \ cdot \ frac {d + s} {2} \\
＆ \ ge \ frac {\ sqrt {2}} {1 + 1/2} \ cdot \ frac {d + s} {2} \\
＆= \ frac {2} {3} \ sqrt {2 } \ cdot \ frac {d + s} {2} \\
＆\约0.9428 \ cdot \ frac {d + s} {2} \\
＆> \ frac {15} {16 } \ cdot \ frac {d + s} {2} \\
\ end {aligned}
$$
如果a的算术平均值大于b的算术平均值，则上限b的几何平均值的a小于a的几何平均值的下限，这意味着b必须小于a。反之亦然。

这可以处理许多以前未处理的案例。结果现在是：

theHandledPercent 99.52

theCorrectPercent 100.0

上面的源代码已更新。

那些仍未处理的是“太接近通话”。他们可能需要查找表来解决。请继续关注.....


嘿，丹，

好吧，我会缩短它，但没有一个是多余的。除了第一部分外，这就是使球滚动的原因。因此，最重要的摘要将差不多长。我确实打算写一篇博客文章。这是一个有趣的练习，比我最初想象的要深入得多。

我确实整理了有关Olli坡度边界的旁注。实际上需要执行多少操作。叙述中的数学只是操作的理由。

真正的“赢家”应该是最有效的算法。真正的测试是在同一平台上编程并在此处进行测试的两种方法。现在，我可以告诉你，我的代码（用C语言编写）将仅由于我正在使用整数原型制作而他正在使用它，而他正在使用具有大量昂贵操作的浮点数。

我现在的想法是，我未处理的其余0.5％案例最好使用CORDIC迭代方法进行处理。我将尝试以整数形式实现Olli方法的一种变体，包括他的创造力变化的斜率。 “太接近通话”类别的确应该非常接近。

我同意您的观点，Olli做得很好。我从他那里学到了很多东西。

最后，难道我们都不都是赢家吗？


Dan，

您对CORDIC的信念令人鼓舞。我实现了与Olli有所不同的无损CORDIC，但可能等效。到目前为止，我还没有找到您断言这是最终解决方案的说法。我不会发布代码，因为还应该再进行一次测试来限制它。我将区域限制为四分之一圆（相当于一个完整的圆），并且对它们的评分有所不同。

我建议奥利（Olli）也是这样。

这是我目前的结果：



在0.08％的情况下，假设为零，实际为87.6％。这意味着只有0.01％的答案是错误的，即错误地假定为零。对于那些商（I1 ^ 2 + Q1 ^ 2）/（I2 ^ 2 + Q2 ^ 2）。显示了高和低值。取平方根即可得到实际的幅度比。

大约83％的案件是由主要裁定所捕获的，不需要任何进一步的处理。这样可以节省很多工作。在大约0.7％的案例中需要使用CORDIC。 （在以前的测试中为0.5％。）

If ( (db+sb) & 1 ) > 0 Then
   da =<< 1
   sa =<< 1
   db =<< 1
   sb =<< 1
EndIf