银河系.site

这几天需要用到求一个数组的中值，普遍的做法是先排序，然后取中间那个值。
但是这样做，冗余度太高，我并不需要其他部分的顺序。

于是找到一个猜测法，基本思想是，先猜测一个值a，求大于a和小于a的数的个数。
如果小于a的数比较多，就减少猜测值。反之则增加猜测值。

上代码

    Real , Intent(IN) :: x(n)
    Real, Parameter :: AFac = 1.5, amp = 1.5
    Integer :: np, nm, j, n
    Real :: a, eps, ap, am, sum, sumx, xp, xm, dum
    a   = (x(1)+x(n))/2.0 !//猜测值
    eps = abs(x(n)-x(1))
    ap  = Huge(ap)
    am  = -Huge(ap)
    Do
      sum  = 0.
      sumx = 0.
      np   = 0
      nm   = 0
      xp   = Huge(xp)
      xm   = -Huge(xm)
      Do j = 1, n !//循环求大于猜测值和小于猜测值的个数
        If (x(j)/=a) Then
          If (x(j)>a) Then
            np = np + 1
            If (x(j)<xp) xp = x(j)
          Else If (x(j)<a) Then
            nm = nm + 1
            If (x(j)>xm) xm = x(j)
          End If
          dum = 1./(eps+abs(x(j)-a))
          sum = sum + dum
          sumx = sumx + x(j)*dum
        End If
      End Do
      If (np-nm>=2) Then !猜测值太低了
        am = a
        xmed = xp + max(0., sumx/sum-a)*amp
        If (xmed>ap) xmed = (a+ap)/2
        eps = AFac*abs(xmed-a)
        if(eps<1.0e-6) exit !//精度不需要那么高，不猜了
        a = xmed
      Else If (nm-np>=2) Then !猜测值太高了
        ap = a
        xmed = xm + min(0., sumx/sum-a)*amp
        If (xmed<am) xmed = (a+am)/2
        eps = AFac*abs(xmed-a)
        if(eps<1.0e-6) exit !//精度不需要那么高，不猜了
        a = xmed
      Else !猜测值合适
        If (mod(n,2)==0) Then !//偶数个数求平均
          If (np==nm) Then
            xmed = (xp+xm)/2
          Else If (np>nm) Then
            xmed = (a+xp)/2
          Else
            xmed = (xm+a)/2
          End If
        Else !//奇数个数
          If (np==nm) Then
            xmed = a
          Else If (np>nm) Then
            xmed = xp
          Else
            xmed = xm
          End If
        End If
        exit
      End If      
    End Do
  End Function median

采用中二猜测法，我的测试中可以提升2-3倍的效率。数组越大，越显著。