float-float 模拟CPU和GPU上double双精度计算

已于 2024-07-13 18:47:13 修改
阅读量1.4k 收藏 3
点赞数 1
分类专栏： 数值计算/数值优化 文章标签： c++ 算法 混合精度 CUDA
于 2022-06-08 09:31:31 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/weixin_42849849/article/details/125177997
版权
博客介绍了如何通过FLOAT2结构体在CPU和GPU上模拟双精度浮点数运算，包括加减乘除、开方、倒数等操作。FLOAT2使用两个单精度浮点数来近似表示双精度数值，适用于混合精度计算。测试程序展示了FLOAT2在计算根号2和根号3之和等运算中的应用。
摘要由CSDN通过智能技术生成

float-float

基于Lewis书中代码加工。可以通过使用两个float单精度浮点数模拟CPU和GPU上双精度double类型运算。

在需要混合精度计算逻辑中可以使用，保证特定计算的精度需求。

FLOAT2类型源码(float2.h)
#pragma once 
#ifdef __CUDACC__
#define __HOST__DEVICE__  __host__ __device__
#else 
#define __HOST__DEVICE__ 
#endif

#define __INLINE__ __HOST__DEVICE__ inline


struct FLOAT2{

float x32,y32;

__INLINE__ 
FLOAT2():x32(0.0f),y32(0.0f){}

__INLINE__ explicit 
FLOAT2(float x):x32(x),y32(0.0f){}

__INLINE__  explicit 
FLOAT2(float x,float y):x32(x),y32(y) {}

/**************************************************************
 *      类型转换
 * ************************************************************/
//类型转换: FLOAT2->float
__INLINE__  
operator float(){return x32;}

//类型转换: FLOAT2->double
__INLINE__ 
operator double(){return double(x32)+double(y32);}

//double->FLOAT2, 未考虑溢出
__INLINE__  explicit 
FLOAT2(double a)
{
    x32=a;
    y32=a-x32;
}

__INLINE__ static bool isnan(const FLOAT2 a64)
{
    return std::isnan(a64.x32) || std::isnan(a64.y32);
}
__INLINE__ static bool isinf(const FLOAT2 a64)
{
    return !isnan(a64) &&(std::isinf(a64.x32) || std::isinf(a64.y32));
}
__INLINE__ static bool isnormal(const FLOAT2 a64)
{
    return std::isnormal(a64.x32) && std::isnormal(a64.y32);
}
__INLINE__ static bool isfinite(const FLOAT2 a64)
{
    std::isfinite(a64.x32) && std::isfinite(a64.y32);
}
__INLINE__ static bool iszero(const FLOAT2 a64)
{
    return a64.x32==0.0f && a64.y32==0.0f;
}

/********************************************************
 *                   计算函数
 * *******************************************************/
//==
__INLINE__ static bool EQ64(const FLOAT2 a64, const FLOAT2 b64)
{
    return a64.x32==b64.x32 && a64.y32==b64.y32;
}
//> 
__INLINE__ static bool GT64(const FLOAT2 a64, const FLOAT2 b64)
{
    return a64.x32>b64.x32 || (a64.x32==b64.x32 && a64.y32>b64.y32);
}
//>=
__INLINE__ static bool GE64(const FLOAT2 a64, const FLOAT2 b64)
{
    return a64.x32>=b64.x32 || (a64.x32==b64.x32 && a64.y32>=b64.y32);
}
//<
__INLINE__ static bool LT64(const FLOAT2 a64, const FLOAT2 b64)
{
    return !GE64(a64,b64);
}
//<=
__INLINE__ static bool LE64(const FLOAT2 a64, const FLOAT2 b64)
{
    return !GT64(a64,b64);
}

//-
__INLINE__ static FLOAT2 NEG64(const FLOAT2 a64)
{
    return FLOAT2{-a64.x32,-a64.y32};
}

/************************************************************
 *                    Addition
 * ***********************************************************/

//将a32+b32规格化成FLOAT2
__INLINE__ static FLOAT2 QuickTwoSum(const float a32, const float b32)
{
    // assumes a32 > b32
    float s32, e32 ;

    s32 = a32 + b32 ;
    e32 = b32 - (s32 - a32) ;                   // This is NOT always zero!
    return  FLOAT2{s32, e32} ;
}

__INLINE__ static FLOAT2 TwoSum(const float a32, const float b32)
{
    float s32, v32, e32 ;

    s32 = a32 + b32 ;
    v32 = s32 - a32 ;                           // This is NOT the same as b32!
    e32 = (a32 - (s32 - v32)) + (b32 - v32) ;   // This is NOT always zero!
    return FLOAT2{s32, e32} ;
}

__INLINE__ static FLOAT2 Add64(const FLOAT2 a64, const FLOAT2 b64)
{
    FLOAT2 s64, t64 ;

    s64 = TwoSum(a64.x32, b64.x32) ;
    t64 = TwoSum(a64.y32, b64.y32) ;
    s64.y32 += t64.x32 ;
    s64 = QuickTwoSum(s64.x32, s64.y32) ;
    s64.y32 += t64.y32 ;
    return QuickTwoSum(s64.x32, s64.y32) ;
}

__INLINE__ static FLOAT2 Add64(const FLOAT2 a64, const float b32)
{
    FLOAT2 s64;
    s64 = TwoSum(a64.x32, b32) ;
    s64.y32 += a64.y32 ;
    return QuickTwoSum(s64.x32, s64.y32) ;
}

__INLINE__ static FLOAT2 Add64(const float a32, const FLOAT2 b64)
{
    return Add64(b64,a32);
}

/*********************************************************
 *                   Subtraction
 * *******************************************************/
__INLINE__ static FLOAT2 Sub64(const FLOAT2 a64,const FLOAT2 b64)
{
    return Add64(a64,NEG64(b64));
}

__INLINE__ static FLOAT2 Sub64(const FLOAT2 a64,const float b32)
{
    return Add64(a64, -b32);
}

__INLINE__ static FLOAT2 Sub64(const float a32, const FLOAT2 b64)
{
    return NEG64(Sub64(b64,a32));
}

/*****************************************************************
 *                    Multiplication
 * ****************************************************************/

__INLINE__ static FLOAT2 Split(const float a32)
{
    float p32,ahi32,alo32;
    p32=4097.0f*a32;
    ahi32=p32-(p32-a32);
    alo32=a32-ahi32;

    return FLOAT2(ahi32,alo32);
}

__INLINE__ static FLOAT2 TwoProd(const float a32, const float b32)
{
    FLOAT2 a64, b64 ;
    float p32, e32 ;

    p32 = a32 * b32 ;
    a64 = Split(a32) ;
    b64 = Split(b32) ;
    e32 = ((a64.x32*b64.x32 - p32) + a64.x32*b64.y32 + a64.y32*b64.x32) + a64.y32*b64.y32 ;
    return FLOAT2{p32, e32} ;
}

__INLINE__ static FLOAT2 Mul64(const FLOAT2 a64, const FLOAT2 b64)
{
    FLOAT2 p64 ;

    p64 = TwoProd(a64.x32, b64.x32) ;
    p64.y32 += a64.x32 * b64.y32 ;
    p64.y32 += a64.y32 * b64.x32 ;
    return QuickTwoSum(p64.x32, p64.y32) ;
}

__INLINE__ static FLOAT2 Mul64(const FLOAT2 a64, const float b32)
{
    FLOAT2 p64 ;

    p64 = TwoProd(a64.x32, b32) ;
    p64.y32 += a64.y32 * b32 ;
    return QuickTwoSum(p64.x32, p64.y32) ;
}

__INLINE__ static FLOAT2 Mul64(const float a32, const FLOAT2 b64)
{
    return Mul64(b64,a32);
}

/*****************************************************************
 *                    Division
 * ****************************************************************/

//牛顿迭代法，迭代一次就OK
__INLINE__ static FLOAT2 Div64(const FLOAT2 a64, const FLOAT2 b64)
{
    //check divide by zero
    //ASSERT(!iszero(b64))

    FLOAT2 q64, e64 ;
    float i32, d32 ;

    i32 = 1.0f / b64.x32 ;

    float q32=a64.x32*i32;
    d32=Sub64(a64,Mul64(b64,q32)).x32;
    e64=TwoProd(i32,d32);
    return Add64(q32,e64);
}

__INLINE__ static FLOAT2 Div64(const FLOAT2 a64, const float b32)
{
    //check divide by zero
    //ASSERT(b32!=0.0f)

    FLOAT2 q64, e64 ;
    float i32, d32 ;

    i32 = 1.0f / b32 ;

    float q32=a64.x32*i32;
    d32=Sub64(a64,TwoProd(b32,q32)).x32;
    e64=TwoProd(i32,d32);
    return Add64(q32,e64);
}

__INLINE__ static FLOAT2 Div64(const float a32, const FLOAT2 b64)
{
    //check divide by zero
    //ASSERT(!iszero(b64))
    
    FLOAT2 t64(a32);
    return Div64(t64,b64);
}

/*****************************************************************
 *                    special functions
 * ****************************************************************/
 //TODO: 
 //ln,log10
 //sin,cos,sincos
 //exp,exp1m
 //polyval

//计算倒数
__INLINE__ static FLOAT2 Recip64(const FLOAT2 a64)
{
    FLOAT2 e64 ;
    float i32, d32 ;
    
    //近似值
    i32 = 1.0f / a64.x32 ;
    
     //一次牛顿迭代法: 求零点 f(x)=x^{-1}-a
    //迭代: x=x*(2-a*x)
    return Mul64(i32,Sub64(2.0f,Mul64(a64,i32)));
}

//计算1/sqrt(a)
//x=1/sqrt(a)=> f(x)=x^{-2}-a 求零点
//迭代: x=x*(1.5-0.5*a*x*x)
//https://math.stackexchange.com/questions/607794/newton-raphson-for-reciprocal-square-root
 __INLINE__ static FLOAT2 Rsqrt64(const FLOAT2 a64)
 {
     //近似值
     float i32=rsqrt(a64.x32);

     //一次牛顿迭代
     //迭代: x=x*(1.5-0.5*a*x*x)=0.5*x*(3-a*x*x)
     FLOAT2 t64=Mul64(i32,Mul64(i32,a64));
     return Mul64(0.5f*i32, Sub64(3.0f, t64));
 }

 __INLINE__ static FLOAT2 Sqrt64(const FLOAT2 a64)
 {
     return Recip64(Rsqrt64(a64));
 }

}; //struct FLOAT2
测试程序(test_float2.cu)
#include "float2.h"
#include <cstdio>
__INLINE__  static void test_float2()
{
   const double sqrt_2=1.4142135623730950488016887242097;
   const double sqrt_3=1.7320508075688772935274463415059;
   const double sqrt_2_plus_sqrt_3=3.146264369941972560695830907207;
   const float z=1.5f; //3/2

   printf("sqrt_2=%18.16g, sqrt_3=%18.16g\n",sqrt_2,sqrt_3);

   FLOAT2 a64(sqrt_2),b64(sqrt_3);

   auto x64=FLOAT2::Add64(a64,b64);

   double y=x64;

   printf("y=%18.16g, %18.16g\n",y,sqrt_2+sqrt_3);

   printf("%18.16g\n",double(FLOAT2::Sub64(z,FLOAT2(sqrt_2))));

   x64=FLOAT2::Div64(a64,b64);
   y=x64;
   printf("y=%18.16g\n",y);

   x64=FLOAT2::Recip64(FLOAT2(sqrt_2));
   y=x64;
   printf("y=%18.16g\n",y);
   
   x64=FLOAT2::Rsqrt64(FLOAT2(2.0f));
   y=x64;
   printf("y=%18.16g\n",y);

}

__global__ void kernel_float2()
{
    test_float2();
}

void 
test_driver()
{
   //==========================================================
   printf("----------->On CPU:\n");
   test_float2();
   
   //============================================================
   printf("------------>On GPU: \n");
   kernel_float2<<<1,1>>>();
}

int main()
{
    test_driver();
    return(0);
}

//编译: nvcc -O3 -use_fast_math test_float2.cu

/***********************************************************
matlab:

vpa(sqrt(2)) #1.4142135623730950488016887242097
vpa(sqrt(3)) #1.7320508075688772935274463415059
vpa(sqrt(2)+sqrt(3)) #3.146264369941972560695830907207
vpa(3/2-sqrt(2)) #0.085786437626904854525378141261172
vpa(sqrt(2/3)) #0.81649658092772603273242802490196 

**************************************************************/
网络资源
float-float模拟double运算
ieee float, fp16,bfp16,tfp32

float16

bfloat16

混合精度计算
weixin_42849849
关注
专栏目录
如何通过CPU多核并行和GPU多核并行两种方式提升GPU的运算性能 Accelerating GPUs with Multithreading Programming
程序员光剑
08-05 1328
GPU(Graphics Processing Unit)多核并行处理单元是当前科技领域里最热门的硬件之一。如今GPU已经集成在笔记本电脑、平板电脑、服务器和手机等各类设备中。由于其高性能的计算能力及其独特的编程模型（比如CUDA），使得GPU在科技界占有重要地位。同时，GPU的多线程编程技术也越来越受到关注。本文将从多线程编程的基本概念出发，介绍基于CUDA编程模型的多线程编程技术。然后，介绍如何通过CPU多核并行和GPU多核并行两种方式提升GPU的运算性能。
CUDA：使用Warp的双精度GEM计算的实例
希望我的博客，能帮上你解决学习中工作中所遇到的问题
12-21 166
CUDA：使用Warp的双精度GEM计算的实例
用2个float模拟double
weixin_34342992的博客
01-19 393
在有些设备上只有float没有double，比如前几代GPU、部分移动设备等。当非得用到double精度的时候该怎么办？ 我记得去年在某个地方见到过用2个float模拟double的作法，经过一番玩命地搜索，得来全不费功夫，就在CUDA SDK的Mandelbrot例子里找到了2个float模拟double乘法的函数。甚至，GTX280上的double也是类似的方法模拟出来的，所 以慢的惊人...
CpuFloat.apk
09-24
显示CPU GPU 频率，浮窗显示 适用于Android版本 是一个工具apk 不知道出处，如有侵权，联系删。
用单精度浮点模拟双精度的运算，实现扩展浮点精度的算法
qq_25675517的博客
11-04 1837
实际上，任何精度的浮点运算也可以像BigInteger那样计算，可以用多个浮点数据表示更高精度的浮点数据。而这就要求实现带余数的浮点四则运算。 带余数的浮点四则运算的定义： 加法和减法，余数和结果的和能精确表示运算值。 乘法和除法，余数和结果必须同号，对于乘法，余数和结果的和能精确表示运算值，而对于除法，余数就是原来的意思。 不多解释，代码如下：（可能存在错误） #include <iostream> #include <xmmintrin.h> using na.
基准测试移动 GPU 中的浮点精度 - 第 1 部分
weixin_30950607的博客
03-29 231
原文:Benchmarking floating-point precision in mobile GPUs 作者：Tom Olson/2013年5月29日下午2:39 / 谈到GPU性能时（在Mali官方网站有很多讨论），我们通常会谈论速度。在之前的博文中，我们谈论了每秒可在屏幕中显示的像素数，每秒可以宣称绘制的三角形数（不要问），最近我们还谈论了每秒可以运行的...
2023.6.13-【double】双精度浮点型：求两个整数的平均值，以及double型数据的加减乘除
最新发布
cheng19960315的博客
06-13 1441
利用双精度浮点型数据【double】来实现两个整数的平均值计算。
c++单精度和双精度的运算例子（计算机占位存储）
Briwisdom的博客
02-15 8880
计算机中float, double类型数据分别占据4，8个字节，其中float类型和double可以表示的小数位数不同，导致了精度不同。double的精度更高。 计算机中数据的表示由：符号位， 指数位，尾数位组成。比如一个float类型数字的二进制由左到右依次是符号位，指数位，尾数位。 类型 符号位 指数位 尾数位 总位数（bit) float 1 8 23 32 double 1 11 52 64 数字1.4在计算机中
gpu浮点计算能力floaps_GPU中双精度（double)浮点运算的实现！
weixin_39823299的博客
12-22 2689
当前NVIDA的GPU芯片仅支持单精度(float)浮点运算，对一些应用来说精度可能不够用，一些关键的步骤可能需要双精度运算，才能保证程序的正常执行。对此，本人尝试用两个单精度浮点数数来代表一个双精度浮点数：//类型定义, 对于GPU编程，我们可以直接用内置类型float2, 这里的定义只是为了清晰起见。typedef struct { float x, y;} float2双精度和float2之...
gpu浮点计算能力floaps_认识GPU浮点计算精度
weixin_39626586的博客
12-22 1076
● 认识GPU浮点计算精度DirectX 9.0时代提出的Shader Model 2.0最重要的一点改进是增加对浮点数据的处理功能，以前GPU只能对整数进行处理，改进后提高渲染精度，使最终处理的色彩格式达到电影级别。Shader Model 2.0时代突破了以前限制PC图形图象质量在数学上的精度障碍，它的每条渲染流水线都升级为128位浮点颜色，让游戏程序设计师们更容易更轻松的创造出更漂亮的效果，...
nVIDIA显卡CUDA性能测试工具，可测试nVIDIA显卡的浮点运算性能
12-20
nVIDIA显卡CUDA性能测试工具，可测试nVIDIA显卡的参数，nVIDIA显卡的CUDA浮点运算性能
GPU双精度测试：Unity GPU双精度运算测试
02-20
GPU双精度测试 Unity GPU双精度运算测试 预览
C++矩阵乘法计算 || GPU && CPU 实现
春江水暖鸭先知
09-03 5152
废话 使用英伟达加速矩阵运算一般使用cublas库进行操作，在cublas库中，使用cublasSgemv()和cublasSgemm()分别进行矩阵向量间的乘法运算与矩阵矩阵间的乘法运算。本文将具体的解释上述两个函数的参数以及具体的应用例子。参照官方解释 矩阵乘法函数解释 仅仅使用单精度(float)进行解释，若想尝试双精度(double),可以进入官网查看相应函数. ...
gpu浮点计算能力floaps_基准测试移动 GPU 中的浮点精度 - 第 2 部分
weixin_29689845的博客
01-13 671
投稿人：，2013年6月11日这是有关GPU中浮点质量的一系列博文中的第二篇，我的灵感源自发表于的文章。在中，我宣称许多程序员其实并不真正了解浮点数字，也指出如果您准备将它用于比较棘手的东西，那么最好先准备好更加深入地了解其运作原理，甚至要超过您所希望达到的深度。我介绍了Stuart的测试，也说明它披露了GPU片段着色器中所用浮点精度的位数。这很有趣，但这个测试告诉我们的不止这些。在本文中，我...
单精度和双精度的运算问题
qq_44792712的博客
08-02 2711
你好，我的朋友 学习目标： 单精度和双精度的运算问题 问题引出： 当我们在运行System.out.println(1.0-0.66); 结果会是什么呢？ 问题解决： 可能许多小伙伴会觉得输出0.34吧，但是事实果真如此吗？ 首先，大家都知道，小数类型的数据默认是Double类型，也就是说1.0是双精度型。而计算机底层是以二进制的原理存储数据，包括运算数据。减法在计算机中也是加法的过程。所以当我们将0.66转换成二进制时是除不尽的，此时就是问题所在，所以我们输出的结果是：0.3399999999999
GPU编程和流式多处理器（四）
吴建明wujianming_110117
01-05 406
GPU编程和流式多处理器（四） 3.2. 单精度（32位） 单精度浮点支持是GPU计算的主力军。GPU已经过优化，可以在此数据类型上原生提供高性能，不仅适用于核心标准IEEE操作（例如加法和乘法），还适用于非标准操作（例如对先验的近似（例如sin（）和log（）））。32位值与整数保存在同一寄存器文件中，因此单精度浮点值和32位整数（使用__float_as_int（）和__int_as_float（））之间的强制转换是免费的。 加法，乘法和乘加 编译器自动将浮点值的+，–和运算符转换为加，乘和乘加指令。所
getUsed是干嘛的？float64是在计算cpu的什么信息？
阿里云开发者
05-05 411
是一种数据类型，通常用于存储浮点数。在计算 CPU 信息时，可能会使用一些浮点数，例如 CPU 的占用率等信息，指获取某个资源（例如 CPU、内存等）的使用量，常常用于计算资源利用率或监控资源使用情况。的精度较高，可以使计算结果更加准确。
C语言不同数据类型间的混合运算+常见数据类型
热门推荐
Fengjingdisan的博客
07-29 2万+
1zai在 在程序中经常会遇到不同类型的数据进行运算，若一个运算符两侧的数据类型不同，则先自动进行类型转换，使两者具有同一类型，然后进行运算，现将规律总结如下： 1. +. -. *. /运算的两个数中有一个数为float或double型，结果是double型，因为系统将所有float型数据都先转换为double型，然后进行运算； 2. 若int型与float或double型数据
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值