CORDIC怎么开方 cordic计算平方根

CORDIC(坐标旋转数字算法)是一种计算三角、双曲和其他数学函数的数字算法，每次运算均产生一次结果输出。这使我们能够根据应用需求调整算法精度；增加运算迭代次数可以得到更精确的结果。CORDIC 是只使用加法、减法、移位和查找表实现的简单算法，这种算法在FPGA中实现效率高，在硬件算法实现中经常用到。

本文主要在此下文章介绍CORDIC双曲系统的基础上介绍平方根计算。

HLS / Chisel 实现CORDIC算法双曲系统

原理

在CORDIC算法双曲系统的向量模式中，最终的输出结果如下：

$\begin{aligned} &\left\{\begin{array}{l} x_{n}=K_{n} \sqrt{x_{1}^{2}-y_{1}^{2}} \\ y_{n}=0 \\ z_{n}=z_{1}+\tanh ^{-1}\left(y_{1} / x_{1}\right) \\ K_{n}=\prod_{i=1}^{n-1} \sqrt{1-2^{-2 i}} \end{array}\right.\\ &\sigma_{i}=\left\{\begin{array}{l} +1, y_{i} \leqslant 0 \\ -1, y_{i}>0 \end{array}\right. \end{aligned} \\ -0.784<tanh({z_n}) =\frac {y_1} {x_1}<0.784$

具体的数学推导可以如上引用的文章查看。

我们可以注意到，在这个迭代过程中，$x_n$出现了平方根函数，因此我们可以稍作变换：设a为待求的平方根值，设x的初始值为a+1，y的初始值为a-1，z任意值皆可，经过迭代后的x值为：

$x=K^{*} \sqrt{(a+1)^{2}-(a-1)^{2}}=K^{*} \sqrt{4 a}=2 K^{*} \sqrt{a}$

已知

$\begin{gathered} K=1 / P= 0.8297816201389014 \\ P = 1.205136358399695 \end{gathered}$

则将双曲线系统向量模式的输出$x_n$乘P,再右移一位即可得到a的平方根

$\sqrt{a}=(p*x_n)>>1$

输入范围问题

注意到$-0.784<tanh({z_n}) =\frac {y_1} {x_1}<0.784$，则有$-0.784< \frac {a-1} {a+1}<0.784$,所以x的输入范围需要限制在[0,8]范围内。

此外，由于[0,1]范围内数太小，其在迭代的数值计算中误差太大，最终结果也不理想，所以我们希望输入的范围在[1,8].

为了扩充范围到R+,我们可以将输入a使用移位运算移动N位，N为偶数，移位到[1,8]范围后得到$a`$，最后通过HV模块计算后得到结果$x_n$，在移位N/2恢复即可。

HLS实践

注意：这里只实现了基础班的求平方运算,并未通过二进制移位解决输入问题，读者可以参照解决。

#define THETA_TYPE float
#define COS_SIN_TYPE float
#define NUM_ITERATIONS 50

THETA_TYPE cordic_P_mark_n = 1.205136358399695;
void cordic_hv_square_root(COS_SIN_TYPE a, COS_SIN_TYPE &x_out)
{
    x = a+1;
    y = a-1；
    
    for(int j = 1;j <= NUM_ITERATIONS; j++){
        COS_SIN_TYPE temp_x = x;
        
        if(y < 0){
            x = x + y >> j;
            y = y + temp_x >> j;
            theta = theta - cordic_phase_h[j-1];
        } else{
            x = x - y >> j;
            y = y - temp_x >> j;
            theta = theta + cordic_phase_h[j-1];
        }
    }
    
    x_out = (x * cordic_P_mark_n) >> 1;
}

Chisel实践

定点数定义

首先定义定点数，注意在chisel.experimental包里有定点数的类

/* Cordic 的全局常数定义 */
trait Cordic_Const_Data extends HasDataConfig{
  /* 迭代次数配置 */
  val NUM_ITERATIONS = 20
}

trait HasDataConfig {
  /* 定点数的位宽定义 */ 
  val DataWidth = 32
  val BinaryPoint = 20
}

源码实现

这里在引用的文章基础上改了cordic算法HV模块代码，注意这里只需要计算x，y即可，不需要双曲线系统计算中的z

• 首先实现二进制移位，并生成伴生对象

class shift_range_1_8 extends Module with HasDataConfig{
  /*
   * @x : 输入的待移位数 定点数类型 取值范围R
   * @out : 输出移位后的结果 范围[1,8]
   * @cnt : 输出的移位数量，左移为正，右移为负
   * details: 组合逻辑电路将x通过左右二倍数移位至范围[1,8]，并输出移位数目
  **/
  val io = IO(new Bundle {
    val x: FixedPoint = Input(FixedPoint(DataWidth.W, BinaryPoint.BP))
    val out: FixedPoint = Output(FixedPoint(DataWidth.W, BinaryPoint.BP))
    val cnt: SInt = Output(SInt(log2Ceil(DataWidth).W))
  })

  val temp_x: FixedPoint = io.x

  when(temp_x < FixedPoint.fromDouble(1, DataWidth.W, BinaryPoint.BP)){
    /* 比0.5小，需要左移 */
    val index = VecInit(Seq.fill(BinaryPoint>>1)(0.B))
    for(i <- 0 until (BinaryPoint>>1)){
      when((temp_x << (i<<1)) < FixedPoint.fromDouble(1, DataWidth.W, BinaryPoint.BP)){
        index(i) := 0.B
      }.otherwise{
        index(i) := 1.B
      }
    }
    /* 优先编码器获得首先大于等于1的位置，即cnt的值 */
    val temp_cnt = PriorityEncoder(index.asUInt())
    io.cnt := (temp_cnt<<1).asSInt()
    io.out := (temp_x << (temp_cnt<<1))
  }.elsewhen(temp_x > FixedPoint.fromDouble(8, DataWidth.W, BinaryPoint.BP)){
    /* 比8大，需要右移 */
    val index = VecInit(Seq.fill(DataWidth>>1)(0.B))
    for(i <- 0 until (DataWidth>>1)){
      when((temp_x >> (i<<1)) > FixedPoint.fromDouble(8, DataWidth.W, BinaryPoint.BP)){
        index(i) := 0.B
      }.otherwise{
        index(i) := 1.B
      }
    }
    /* 优先编码器获得首先小于等于8的位置，即cnt的值 */
    val temp_cnt = PriorityEncoder(index.asUInt())
    io.cnt := -((temp_cnt<<1).asSInt())
    io.out := (temp_x >> (temp_cnt<<1))
  }.otherwise{
    io.cnt := 0.S
    io.out := temp_x
  }
}


/*
 * 定义这个类的伴生对象，并定义一个工厂方法来简化模块的例化和连线。
**/
object shift_range_1_8 {
  def apply(x: FixedPoint): (FixedPoint, SInt) = {
    /*
     * @x : 输入的待移位数 定点数类型 取值范围R
     * @out : 输出移位后的结果 范围[0.5,1]
     * @cnt : 输出的移位数量，左移为正，右移为负
     * @flag : 输出x的正负符号分别为1，0
     * details: 将x通过左右移位至范围[0.5,1]，并输出移位数目和是否改变符号
    **/
    val unit = Module(new shift_range_1_8)
    unit.io.x := x
    (unit.io.out, unit.io.cnt)
  }
}

• 实现hv模块中计算平方根源码

class cordic_hv_square_root_origin extends Module with HasDataConfig with Cordic_Const_Data{
  /*
   * @in : 输入的向量横轴坐标 定点数类型 in∈[0,8]
   * @out : 输出的x^0.5 定点数类型
   * details: 利用cordic双曲坐标系的向量模式迭代的标得到(in+1,in-1)坐标的输出x,
              即2K(in^0.5),再乘1/K=P,右移一位即可。注意HV模式限制了y/x的范围,故
              in∈[0,8]
   */
  val io = IO(new Bundle {
    val in: FixedPoint = Input(FixedPoint(DataWidth.W, BinaryPoint.BP))
    val out: FixedPoint = Output(FixedPoint(DataWidth.W, BinaryPoint.BP))
  })

  val x: FixedPoint = io.in + (1.F(DataWidth.W, BinaryPoint.BP))
  val y: FixedPoint = io.in - (1.F(DataWidth.W, BinaryPoint.BP))

  /* 在趋近于无穷大的时候，K值变为恒定 */
  val cordic_P_h: FixedPoint = FixedPoint.fromDouble(1.205136358399695, DataWidth.W, BinaryPoint.BP)

  /* 初始化计算的寄存器数组，形成NUM_ITERATIONS级流水 */
  val current_x: Vec[FixedPoint] = RegInit(VecInit(Seq.fill(NUM_ITERATIONS)(FixedPoint.fromDouble(0.0, DataWidth.W, BinaryPoint.BP)))) // cos
  val current_y: Vec[FixedPoint] = RegInit(VecInit(Seq.fill(NUM_ITERATIONS)(FixedPoint.fromDouble(0.0, DataWidth.W, BinaryPoint.BP)))) // sin

  for (i <- 1 to NUM_ITERATIONS) {
    /*
     * x[i] = K`(x[i-1] + sigma * 2^(-i) * y[i-1])
     * y[i] = K`(y[i-1] + sigma * 2^(-i) * x[i-1])
     * z[i] = z[i-1] - sigma * arctanh(2^(-i))
    **/
    if (i == 1) {
      /* 流水线第一级直接对(io.x,io.y)点做运算 */
      when(y < FixedPoint.fromDouble(0.0, DataWidth.W, BinaryPoint.BP)) {
        current_x(i - 1) := x + (y >> i)
        current_y(i - 1) := y + (x >> i)
      }.otherwise {
        current_x(i - 1) := x - (y >> i)
        current_y(i - 1) := y - (x >> i)
      }
    } else {
      when(current_y(i - 2) < FixedPoint.fromDouble(0.0, DataWidth.W, BinaryPoint.BP)) {
        current_x(i - 1) := current_x(i - 2) + (current_y(i - 2) >> i) // 移位替代乘法
        current_y(i - 1) := current_y(i - 2) + (current_x(i - 2) >> i)
      }.otherwise {
        current_x(i - 1) := current_x(i - 2) - (current_y(i - 2) >> i)
        current_y(i - 1) := current_y(i - 2) - (current_x(i - 2) >> i)
      }
    }
  }

  io.out := (current_x(NUM_ITERATIONS - 1) * cordic_P_h ) >> 1
}

测试

import org.scalatest._
import chisel3._
import chiseltest._
import chisel3.experimental._
import scala.math._

class Cordic_Square_Root_Tester extends FlatSpec with ChiselScalatestTester with Matchers {
  behavior of "mytest2"
  it should "do something" in {
    test(new cordic_hv_square_root) { c =>
      c.io.in.poke(2.F(32.W, 20.BP))
      c.clock.step(20)
      println(s"out ${c.io.out.peek}\n")

      c.io.in.poke(4.F(32.W, 20.BP))
      c.clock.step(20)
      println(s"out ${c.io.out.peek}\n")

      c.io.in.poke(8.F(32.W, 20.BP))
      c.clock.step(20)
      println(s"out ${c.io.out.peek}\n")

      c.io.in.poke(16.F(32.W, 20.BP))
      c.clock.step(20)
      println(s"out ${c.io.out.peek}\n")

      c.io.in.poke(64.F(32.W, 20.BP))
      c.clock.step(20)
      println(s"out ${c.io.out.peek}\n")

      c.io.in.poke(128.F(32.W, 20.BP))
      c.clock.step(20)
      println(s"out ${c.io.out.peek}\n")

      // 注意这里32位不够
//      c.io.in.poke(FixedPoint.fromDouble(65532,32.W, 20.BP))
//      c.clock.step(20)
//      println(s"out ${c.io.out.peek}\n")

      c.io.in.poke(0.00132.F(32.W, 20.BP))
      c.clock.step(20)
      println(s"out ${c.io.out.peek}\n")

      c.io.in.poke(0.6315.F(32.W, 20.BP))
      c.clock.step(20)
      println(s"out ${c.io.out.peek}\n")

      c.io.in.poke(0.000632.F(32.W, 20.BP))
      c.clock.step(20)
      println(s"out ${c.io.out.peek}\n")
    }
  }
}

// 执行 sbt "testOnly Cordic_Square_Root_Tester"