I am working on an algorithm which requires calculations in large numbers, upto e+30. I am using a 32 bit system with compiler support of 32 bits for long/float/double. So far, by searching online, I've learned that single-precision floating points (FPs) can be used for Double-precision FPs.
From this question asked by someone earlier (Emulate “double” using 2 “float”s) I found this paper which has the algorithm to work with Double-precision FPs in GPUs. It is too confusing for me to implement in C. I just need four basic mathematical operations. Is there any way I could find an example for this which will help me understand it better?
Thanks in advance.
Here is the Code I am working on. It might have errors i can not see, any suggestions would be appreciated to rectify error but that is preety much what I am trying to implement. In the algorithm, POLYNOMIAL_ORDER should be able to go up to forth order (can settle at Third order if the standard deviation is smaller). Few things I am not sure about are 1) Procedures make_float() and make_float() are correct or not, 2) Use of make_float() in the program.
#define POLYNOMIAL_ORDER (3)
#define TC_TABLE_SIZE (14)
typedef struct vector_float2{
float x;
float y;
}float2;
typedef struct
{
float tc0;
float tc1;
float tc2;
float tc3;
}POLYNOMIALS;
typedef struct {
int16_t Temp;
int16_t Comp;
} TempCompPair;
volatile TempCompPair TCtable[TC_TABLE_SIZE] = {{22452,1651},
{25318,1444},
{28268,1133},
{31120,822},
{34027,511},
{36932,185},
{39770,-81},
{42685,-288},
{45531,-407},
{48425,-632},
{51401,-703},
{54460,-1143},
{57202,-1420},
{60027,-1652}};
POLYNOMIALS polynomials;
float matrix[TC_TABLE_SIZE][TC_TABLE_SIZE] = {0};
float average[TC_TABLE_SIZE] = {0};
float make_float(float x, float y)
{
return x+y;
}
float2 make_float2(float a, float b)
{
float2 f2 = {a,b};
return f2;
}
float2 quickTwoSum(float a, float b)
{
float s = a+b;
float e = b - (s - a);
float2 result = {s, e};
return result;
}
float2 twoSum(float a, float b)
{
volatile float s = a + b;
float v = s - a;
float e = (a - (s - v)) + (b - v);
float2 result = {s , e};
return result;
}
float2 df64_add(float2 a, float2 b)
{
float2 s,t;
s = twoSum(a.x, b.x);
t = twoSum(a.y, b.y);
s.y += t.x;
s = quickTwoSum(s.x, s.y);
s.y += t.y;
s = quickTwoSum(s.x, s.y);
return s;
}
float2 split(float a)
{
const float split = 4097; //(1<<12) + 1
float t = a *split;
float a_hi = t - (t - a);
float a_lo = a - a_hi;
float2 result = {a_hi, a_lo};
return result;
}
float2 twoProd(float a, float b)
{
float p = a*b;
float2 aS = split(a);
float2 bS = split(b);
float err = ((aS.x * bS.x - p)
+ aS.x * bS.y + aS.y * bS.x)
+ aS.y * bS.y;
float2 result = {p, err};
return result;
}
float2 df64_mult(float2 a, float2 b)
{
float2 p;
p = twoProd(a.x,b.x);
p.y += a.x * b.y;
p.y += a.y * b.x;
p = quickTwoSum(p.x,p.y);
return p;
}
float2 calculate_power(float base, int pow)
{
int i = 0;
float2 base_f2 = make_float2(base,0);
float2 result_f2 = {1,0};
if(pow == 0)
{
return result_f2;
}
if(pow > 0)
{
if(pow == 1)
{
return base_f2;
}
else
{
for(i = 0; i < pow; i++)
{
result_f2 = df64_mult(result_f2,base_f2);
}
return result_f2;
}
}
else
{
return result_f2;
//Mechanism for negative powers
}
}
void TComp_Polynomial()
{
int i;
int j;
int k;
int size;
float temp;
float2 sum = {0,0};
float2 result0 = {0,0};
float2 result1 = {0,0};
float x[TC_TABLE_SIZE];
float y[TC_TABLE_SIZE];
for(i = 0; i < TC_TABLE_SIZE; i++)
{
x[i] = (float) TCtable[i].Temp;
y[i] = (float) TCtable[i].Comp;
}
size = i;
for(i = 0; i <= POLYNOMIAL_ORDER; i++)
{
for(j = 0; j <= POLYNOMIAL_ORDER; j++)
{
sum.x = 0;
sum.y = 0;
for(k = 0; k < size; k++)
{
// Expression simplified below: **sum += pow(x[k],i+j)**
result0 = calculate_power(x[k], i+j);
sum = df64_add(result0,sum);
}
matrix[i][j] = make_float(sum.x,sum.y);
}
}
for(i = 0; i <= POLYNOMIAL_ORDER; i++)
{
sum.x = 0;
sum.y = 0;
for(j = 0; j < size; j++)
{
// Expression simplified below: **sum += y[j] * pow(x[j],i)**
result0 = calculate_power(x[j], i);
result1 = df64_mult( result0 , make_float2(y[j],0) );
sum = df64_add(result1,sum);
}
average[i] = make_float(sum.x,sum.y);
}
for(i = 0; i <= POLYNOMIAL_ORDER; i++)
{
for(j = 0; j <= POLYNOMIAL_ORDER; j++)
{
if(j != i)
{
if(matrix[i][i]!= 0)
{
temp = matrix[j][i]/matrix[i][i];
}
for(k = i; k < POLYNOMIAL_ORDER; k++)
{
matrix[j][k] -= temp*matrix[i][k];
}
average[j] -= temp*average[i];
}
}
}
if(matrix[0][0] != 0)
{
polynomials.tc0 = average[0]/matrix[0][0];
}
if(matrix[1][1] != 0)
{
polynomials.tc1 = average[1]/matrix[1][1];
}
if(matrix[2][2] != 0)
{
polynomials.tc2 = average[2]/matrix[2][2];
}
if(matrix[3][3] != 0)
{
polynomials.tc3 = average[3]/matrix[3][3];
}
}
and then use the struct polynomials.tc0/1/2/3 in below expression
// Y = T^3 * X3 + T^2 * X2 + T^1 * X1 + X0 ;
double calculate_equation(uint16_t TEMP)
{
double Y;
if(POLYNOMIAL_ORDER == 1)
{
Y = polynomials.tc1*(double)TEMP + polynomials.tc0;
}
else if(POLYNOMIAL_ORDER == 2)
{
Y = (polynomials.tc2 * (double)TEMP + polynomials.tc1)*(double)TEMP + polynomials.tc0;
}
else if(POLYNOMIAL_ORDER == 3)
{
Y = ((polynomials.tc3 * (double)TEMP + polynomials.tc2)*(double)TEMP + polynomials.tc1)*(double)TEMP + polynomials.tc0;
}
else if(POLYNOMIAL_ORDER == 4)
{
Y = (((polynomials.tc4 * (double)TEMP + polynomials.tc3)*(double)TEMP + polynomials.tc2)*(double)TEMP + polynomials.tc1)*(double)TEMP + polynomials.tc0;
}
return Y;
}
And standard Deviation is calculated is as follows:
//sqrt(sigma(error^2))
for(i = 0; i < TC_TABLE_SIZE; i++)
{
actual_comp[i] =(int) calculate_equation(TCtable[i].Temp);
error[i] = TCtable[i].Comp - actual_comp[i] ;
error_sqr += error[i]*error[i];
printf("%u\t%d\t\t%e\n", TCtable[i].Temp, TCtable[i].Comp, actual_comp[i] );
}
error_sqrt = sqrt(error_sqr);
Reference: http://hal.archives-ouvertes.fr/docs/00/06/33/56/PDF/float-float.pdf Guillaume Da Graça, David Defour Implementation of float-float operators on graphics hardware, 7th conference on Real Numbers and Computers, RNC7.
