原理参考论文：Floating-Point_Inverse_Square_Root_Algorithm_Based_on_Taylor-Series_Expansion

主要代码如下：

clc
close all
clear

%% 
a_init = [1.125 1.375 1.625 1.875]; %% 区间门限
a = a_init(1);    %% 
x = 1.12;       % 1~2
y1 = 1/sqrt(x);

%% 基于泰勒展开式求开根号倒数
y2 = 1/sqrt(a)*(1-(x-a)/(2*a)+3*(x-a)^2/(8*a^2)-5*(x-a)^3/(16*a^3)+35*(x-a)^4/(128*a^4));

%% 误差对比
sub_out = y1-y2

%% 生成FPGA 数组系数
N = 20;                       %%  小数部分量化位数20位，整数部分量化2位 
a_init_bin = zeros(4,N);
out = zeros(5,N);

for kk=1:4
    a = a_init(kk);
    A(1) = 1/sqrt(a);
    A(2) = A(1)*1/(2*a);
    A(3) = A(1)*3/(8*a^2);
    A(4) = A(1)* 5/(16*a^3);
    A(5) = A(1)* 35/(128*a^4);
    a_init_bin(kk,:) = dectobin(a-1,N); %% 小数转定点数
    B = 1/sqrt(2) * A;  
    for k1=1:5
%         out(k1,:) = dectobin(A(k1),N);  %% even number
        out(k1,:) = dectobin(B(k1),N);  %% odd number
    end
    
   switch kk
    case 1
        init0 = out ;
    case 2
        init1 = out ;
    case 3
        init2 = out ;
    case 4
        init3 = out ;
   end
end

%% 测试部分：matlab浮点计算结果和FPGA定点计算结果对比
xin = [0.625 1.12 1.412 1.566 1.851 3.851 7.851]; % 选取7个值，可代表四个区间,以及奇偶
len = length(xin);

for hh=1:len
    x = xin(hh);
    x_out = 1/sqrt(x); 
    sqrt_matlab_out(hh) = x_out;
end

%% FPGA定点计算结果： compare fpga with matlab
load sqrt_out.txt
fpga_out = sqrt_out.';

for hh=1:len
    fpga_bin = DEC_BIN(fpga_out(hh),N+3); %% 十进制转二进制数

    bit_decimal = fpga_bin ;
    bit1 = 3;     %整数位宽
    bit2 = N;     %小数位宽
    decimal = 0;
    integer = 0;
    %% 二进制数转整数
    for j = 1 : bit1
       integer = integer + bit_decimal(j) * (2^(bit1-j));
    end
    
    %% 二进制数转小数
    for j = 1 : bit2
        decimal = decimal + bit_decimal(j+bit1) * (2^(-j));
    end
    
    fpga_out_dec(hh) = integer + decimal;
end

fpga_sub_matlab = fpga_out_dec - sqrt_matlab_out;

figure
plot(fpga_sub_matlab);
title('fpga和matlab计算误码')
xlabel('point')
ylabel('Approximate errors')

从上图可以看出，matlab计算误码和FPGA计算误差10-6量级，如果位宽再扩大，以及划分更多区间，误差会进一步减低。
如需核心代码可私信。