Added scaling to abval2s, sqrt2s macros.

Details:
- Re-defined abval2s and sqrt2s macros to use scaling to avoid underflow
  and overflow from squaring the real and imaginary components. (This is
  the same technique used to fix recent bugs in invscals/invscaljs and
  inverts.)
This commit is contained in:
Field G. Van Zee
2013-11-02 17:19:40 -05:00
parent c5b1ed9409
commit d70f2b089d
2 changed files with 79 additions and 39 deletions

View File

@@ -84,46 +84,70 @@
#define bli_csabval2s( x, a ) \
{ \
bli_csabsq2s( x, a ); \
bli_sssqrt2s( a, a ); \
float s = bli_fmaxabs( bli_creal(x), bli_cimag(x) ); \
float mag = sqrtf( s ) * \
sqrtf( ( bli_creal(x) / s ) * bli_creal(x) + \
( bli_cimag(x) / s ) * bli_cimag(x) ); \
bli_sssetris( mag, 0.0, a ); \
}
#define bli_zsabval2s( x, a ) \
{ \
bli_zsabsq2s( x, a ); \
bli_sssqrt2s( a, a ); \
double s = bli_fmaxabs( bli_zreal(x), bli_zimag(x) ); \
double mag = sqrt( s ) * \
sqrt( ( bli_zreal(x) / s ) * bli_zreal(x) + \
( bli_zimag(x) / s ) * bli_zimag(x) ); \
bli_dssetris( mag, 0.0, a ); \
}
#define bli_cdabval2s( x, a ) \
{ \
bli_cdabsq2s( x, a ); \
bli_ddsqrt2s( a, a ); \
double s = bli_fmaxabs( bli_creal(x), bli_cimag(x) ); \
double mag = sqrt( s ) * \
sqrt( ( bli_creal(x) / s ) * bli_creal(x) + \
( bli_cimag(x) / s ) * bli_cimag(x) ); \
bli_ddsetris( mag, 0.0, a ); \
}
#define bli_zdabval2s( x, a ) \
{ \
bli_zdabsq2s( x, a ); \
bli_ddsqrt2s( a, a ); \
double s = bli_fmaxabs( bli_zreal(x), bli_zimag(x) ); \
double mag = sqrt( s ) * \
sqrt( ( bli_zreal(x) / s ) * bli_zreal(x) + \
( bli_zimag(x) / s ) * bli_zimag(x) ); \
bli_ddsetris( mag, 0.0, a ); \
}
#define bli_ccabval2s( x, a ) \
{ \
bli_ccabsq2s( x, a ); \
bli_ccsqrt2s( a, a ); \
float s = bli_fmaxabs( bli_creal(x), bli_cimag(x) ); \
float mag = sqrtf( s ) * \
sqrtf( ( bli_creal(x) / s ) * bli_creal(x) + \
( bli_cimag(x) / s ) * bli_cimag(x) ); \
bli_scsetris( mag, 0.0, a ); \
}
#define bli_zcabval2s( x, a ) \
{ \
bli_zcabsq2s( x, a ); \
bli_ccsqrt2s( a, a ); \
double s = bli_fmaxabs( bli_zreal(x), bli_zimag(x) ); \
double mag = sqrt( s ) * \
sqrt( ( bli_zreal(x) / s ) * bli_zreal(x) + \
( bli_zimag(x) / s ) * bli_zimag(x) ); \
bli_dcsetris( mag, 0.0, a ); \
}
#define bli_czabval2s( x, a ) \
{ \
bli_czabsq2s( x, a ); \
bli_zzsqrt2s( a, a ); \
double s = bli_fmaxabs( bli_creal(x), bli_cimag(x) ); \
double mag = sqrt( s ) * \
sqrt( ( bli_creal(x) / s ) * bli_creal(x) + \
( bli_cimag(x) / s ) * bli_cimag(x) ); \
bli_dzsetris( mag, 0.0, a ); \
}
#define bli_zzabval2s( x, a ) \
{ \
bli_zzabsq2s( x, a ); \
bli_zzsqrt2s( a, a ); \
double s = bli_fmaxabs( bli_zreal(x), bli_zimag(x) ); \
double mag = sqrt( s ) * \
sqrt( ( bli_zreal(x) / s ) * bli_zreal(x) + \
( bli_zimag(x) / s ) * bli_zimag(x) ); \
bli_dzsetris( mag, 0.0, a ); \
}

View File

@@ -55,15 +55,19 @@
}
#define bli_cssqrt2s( x, a ) \
{ \
float mag = sqrtf( bli_creal(x) * bli_creal(x) + \
bli_cimag(x) * bli_cimag(x) ); \
float s = bli_fmaxabs( bli_creal(x), bli_cimag(x) ); \
float mag = sqrtf( s ) * \
sqrtf( ( bli_creal(x) / s ) * bli_creal(x) + \
( bli_cimag(x) / s ) * bli_cimag(x) ); \
\
(a) = sqrtf( ( mag + bli_creal(x) ) / 2.0F ); \
(a) = sqrtf( ( mag + bli_creal(x) ) / 2.0 ); \
}
#define bli_zssqrt2s( x, a ) \
{ \
double mag = sqrt( bli_zreal(x) * bli_zreal(x) + \
bli_zimag(x) * bli_zimag(x) ); \
double s = bli_fmaxabs( bli_zreal(x), bli_zimag(x) ); \
double mag = sqrt( s ) * \
sqrt( ( bli_zreal(x) / s ) * bli_zreal(x) + \
( bli_zimag(x) / s ) * bli_zimag(x) ); \
\
(a) = sqrt( ( mag + bli_zreal(x) ) / 2.0 ); \
}
@@ -79,15 +83,19 @@
}
#define bli_cdsqrt2s( x, a ) \
{ \
float mag = sqrtf( bli_creal(x) * bli_creal(x) + \
bli_cimag(x) * bli_cimag(x) ); \
double s = bli_fmaxabs( bli_creal(x), bli_cimag(x) ); \
double mag = sqrt( s ) * \
sqrt( ( bli_creal(x) / s ) * bli_creal(x) + \
( bli_cimag(x) / s ) * bli_cimag(x) ); \
\
(a) = sqrtf( ( mag + bli_creal(x) ) / 2.0F ); \
(a) = sqrt( ( mag + bli_creal(x) ) / 2.0 ); \
}
#define bli_zdsqrt2s( x, a ) \
{ \
double mag = sqrt( bli_zreal(x) * bli_zreal(x) + \
bli_zimag(x) * bli_zimag(x) ); \
double s = bli_fmaxabs( bli_zreal(x), bli_zimag(x) ); \
double mag = sqrt( s ) * \
sqrt( ( bli_zreal(x) / s ) * bli_zreal(x) + \
( bli_zimag(x) / s ) * bli_zimag(x) ); \
\
(a) = sqrt( ( mag + bli_zreal(x) ) / 2.0 ); \
}
@@ -105,16 +113,20 @@
}
#define bli_ccsqrt2s( x, a ) \
{ \
float mag = sqrtf( bli_creal(x) * bli_creal(x) + \
bli_cimag(x) * bli_cimag(x) ); \
float s = bli_fmaxabs( bli_creal(x), bli_cimag(x) ); \
float mag = sqrtf( s ) * \
sqrtf( ( bli_creal(x) / s ) * bli_creal(x) + \
( bli_cimag(x) / s ) * bli_cimag(x) ); \
\
bli_creal(a) = sqrtf( ( mag + bli_creal(x) ) / 2.0F ); \
bli_cimag(a) = sqrtf( ( mag - bli_cimag(x) ) / 2.0F ); \
bli_creal(a) = sqrtf( ( mag + bli_creal(x) ) / 2.0 ); \
bli_cimag(a) = sqrtf( ( mag - bli_cimag(x) ) / 2.0 ); \
}
#define bli_zcsqrt2s( x, a ) \
{ \
double mag = sqrt( bli_zreal(x) * bli_zreal(x) + \
bli_zimag(x) * bli_zimag(x) ); \
double s = bli_fmaxabs( bli_zreal(x), bli_zimag(x) ); \
double mag = sqrt( s ) * \
sqrt( ( bli_zreal(x) / s ) * bli_zreal(x) + \
( bli_zimag(x) / s ) * bli_zimag(x) ); \
\
bli_creal(a) = sqrt( ( mag + bli_zreal(x) ) / 2.0 ); \
bli_cimag(a) = sqrt( ( mag - bli_zimag(x) ) / 2.0 ); \
@@ -129,20 +141,24 @@
#define bli_dzsqrt2s( x, a ) \
{ \
bli_zreal(a) = sqrt( (x) ); \
bli_zimag(a) = 0.0F; \
bli_zimag(a) = 0.0; \
}
#define bli_czsqrt2s( x, a ) \
{ \
float mag = sqrtf( bli_creal(x) * bli_creal(x) + \
bli_cimag(x) * bli_cimag(x) ); \
double s = bli_fmaxabs( bli_creal(x), bli_cimag(x) ); \
double mag = sqrt( s ) * \
sqrt( ( bli_creal(x) / s ) * bli_creal(x) + \
( bli_cimag(x) / s ) * bli_cimag(x) ); \
\
bli_zreal(a) = sqrtf( ( mag + bli_creal(x) ) / 2.0F ); \
bli_zimag(a) = sqrtf( ( mag - bli_cimag(x) ) / 2.0F ); \
bli_zreal(a) = sqrt( ( mag + bli_creal(x) ) / 2.0 ); \
bli_zimag(a) = sqrt( ( mag - bli_cimag(x) ) / 2.0 ); \
}
#define bli_zzsqrt2s( x, a ) \
{ \
double mag = sqrt( bli_zreal(x) * bli_zreal(x) + \
bli_zimag(x) * bli_zimag(x) ); \
double s = bli_fmaxabs( bli_zreal(x), bli_zimag(x) ); \
double mag = sqrt( s ) * \
sqrt( ( bli_zreal(x) / s ) * bli_zreal(x) + \
( bli_zimag(x) / s ) * bli_zimag(x) ); \
\
bli_zreal(a) = sqrt( ( mag + bli_zreal(x) ) / 2.0 ); \
bli_zimag(a) = sqrt( ( mag - bli_zimag(x) ) / 2.0 ); \