Skip to content

Fix scale factor calculations in CROTG/ZROTG and reformat the file #6123

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
martin-frbg merged 2 commits into from
Oct 14, 2026
Merged

Fix scale factor calculations in CROTG/ZROTG and reformat the file #6123

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
226cee4
Fix scale factor calculations and use same the SAFMAX value as Refere…
martin-frbg Oct 13, 2026
6a13afa
Merge branch 'OpenMathLib:develop' into issue6120
martin-frbg Oct 13, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
191 changes: 98 additions & 93 deletions interface/zrotg.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -43,65 +43,69 @@ void CNAME(void *VDA, void *VDB, FLOAT *C, void *VS) {
IDEBUG_START;

FUNCTION_PROFILE_START();

long double safmax = 1./safmin;
safmax *= 2.;
do {
if (db_r == ZERO && db_i == ZERO) {
*C = ONE;
*(S + 0) = ZERO;
*(S + 1) = ZERO;
break;
}
if (db_r == ZERO && db_i == ZERO) {
*C = ONE;
*(S + 0) = ZERO;
*(S + 1) = ZERO;
break;
}

#if defined DOUBLE
long double rtmax = safmax /DBL_EPSILON;
#else
long double rtmax = safmax /FLT_EPSILON;
#endif

long double safmax = 1./safmin;
#if defined DOUBLE
long double rtmax = safmax /DBL_EPSILON;
#else
long double rtmax = safmax /FLT_EPSILON;
#endif
*(S1 + 0) = *(DB + 0);
*(S1 + 1) = *(DB + 1) *-1;
if (da_r == ZERO && da_i == ZERO) {
*C = ZERO;
if (db_r == ZERO) {
(*DA) = fabsl(db_i);
*(S1 + 0) = *(DB + 0);
*(S1 + 1) = *(DB + 1) *-1;

if (da_r == ZERO && da_i == ZERO) {
*C = ZERO;
if (db_r == ZERO) {
(*DA) = fabsl(db_i);
*S = *S1 /(*DA);
*(S+1) = *(S1+1) /(*DA);
break;
} else if ( db_i == ZERO) {
*DA = fabsl(db_r);
} else if ( db_i == ZERO) {
*DA = fabsl(db_r);
*S = *S1 /(*DA);
*(S+1) = *(S1+1) /(*DA);
break;
} else {
long double g1 = MAX( fabsl(db_r), fabsl(db_i));
rtmax =sqrt(safmax/2.);
if (g1 > rtmin && g1 < rtmax) { // unscaled
d = sqrt(adb);
*S = *S1 /d;
*(S+1) = *(S1+1) /d;
*DA = d ;
*(DA+1) = ZERO;
break;
} else { // scaled algorithm
long double u = MIN ( safmax, MAX ( safmin, g1));
FLOAT gs_r = db_r/u;
FLOAT gs_i = db_i/u;
d = sqrt ( gs_r*gs_r + gs_i*gs_i);
*S = gs_r / d;
*(S + 1) = (gs_i * -1) / d;
*DA = d * u;
*(DA+1) = ZERO;
break;
}
}
} else {
FLOAT f1 = MAX ( fabsl(da_r), fabsl(da_i));
FLOAT g1 = MAX ( fabsl(db_r), fabsl(db_i));
rtmax = sqrt(safmax / 4.);
if ( f1 > rtmin && f1 < rtmax && g1 > rtmin && g1 < rtmax) { //unscaled
long double h = ada + adb;
double adahsq = sqrt(ada * h);
if (ada >= h *safmin) {
long double g1 = MAX( fabsl(db_r), fabsl(db_i));
rtmax =sqrt(safmax/2.);
if (g1 > rtmin && g1 < rtmax) { // unscaled
d = sqrt(adb);
*S = *S1 /d;
*(S+1) = *(S1+1) /d;
*DA = d ;
*(DA+1) = ZERO;
break;
} else { // scaled algorithm
long double u = MIN ( safmax, MAX ( safmin, g1));
FLOAT gs_r = db_r/u;
FLOAT gs_i = db_i/u;
d = sqrt ( gs_r*gs_r + gs_i*gs_i);
*S = gs_r / d;
*(S + 1) = (gs_i * -1) / d;
*DA = d * u;
*(DA+1) = ZERO;
break;
}
}

} else {

FLOAT f1 = MAX ( fabsl(da_r), fabsl(da_i));
FLOAT g1 = MAX ( fabsl(db_r), fabsl(db_i));
rtmax = sqrt(safmax / 4.);
if ( f1 > rtmin && f1 < rtmax && g1 > rtmin && g1 < rtmax) { //unscaled
long double h = ada + adb;
double adahsq = sqrt(ada * h);
if (ada >= h *safmin) {
*C = sqrt(ada/h);
*R = *DA / *C;
*(R+1) = *(DA+1) / *C;
Expand All @@ -113,55 +117,56 @@ void CNAME(void *VDA, void *VDB, FLOAT *C, void *VS) {
*S = *S1 * (*R/h) - *(S1+1) * (*(R+1)/h);
*(S+1) = *S1 * (*(R+1)/h) + *(S1+1) * (*(R)/h);
}
} else {
*C = ada / adahsq;
if (*C >= safmin) {
*R = *DA / *C;
*(R+1) = *(DA+1) / *C;
} else {
*C = ada / adahsq;
if (*C >= safmin) {
*R = *DA / *C;
*(R+1) = *(DA+1) / *C;
} else {
*R = *DA * (h / adahsq);
*(R+1) = *(DA+1) * (h / adahsq);
*R = *DA * (h / adahsq);
*(R+1) = *(DA+1) * (h / adahsq);
}
*S = *S1 * ada / adahsq;
*(S+1) = *(S1+1) * ada / adahsq;
}
*DA=*R;
*(DA+1)=*(R+1);
break;
} else { // scaled
FLOAT fs_r, fs_i, gs_r, gs_i;
*S = *S1 * (*DA) / adahsq;
*(S+1) = *(S1+1) * (*DA) / adahsq;
}
*DA=*R;
*(DA+1)=*(R+1);
break;
} else { // scaled
FLOAT fs_r, fs_i, gs_r, gs_i;
long double v,w,f2,g2,h;
long double u = MIN ( safmax, MAX ( safmin, MAX(f1,g1)));
gs_r = db_r/u;
gs_i = db_i/u;
g2 = sqrt ( gs_r*gs_r + gs_i*gs_i);
if (f1 /u < rtmin) {
gs_r = db_r/u;
gs_i = db_i/u;
g2 = gs_r*gs_r + gs_i*gs_i;
if (f1 /u < rtmin) {
v = MIN (safmax, MAX (safmin, f1));
w = v / u;
fs_r = *DA/ v;
fs_i = *(DA+1) / v;
f2 = sqrt ( fs_r*fs_r + fs_i*fs_i);
h = f2 * w * w + g2;
} else { // use same scaling for both
f2 = fs_r*fs_r + fs_i*fs_i;
h = f2 * w * w + g2;
} else { // use same scaling for both
w = 1.;
fs_r = *DA/ u;
fs_i = *(DA+1) / u;
f2 = sqrt ( fs_r*fs_r + fs_i*fs_i);
f2 = fs_r*fs_r + fs_i*fs_i;
h = f2 + g2;
}
if ( f2 >= h * safmin) {
*C = sqrt ( f2 / h );
*DA = fs_r / *C;
*(DA+1) = fs_i / *C;
rtmax *= 2;
if ( f2 > rtmin && h < rtmax) {
*S = gs_r * (fs_r /sqrt(f2*h)) - gs_i * (fs_i / sqrt(f2*h));
*(S+1) = gs_r * (fs_i /sqrt(f2*h)) + gs_i * -1. * (fs_r / sqrt(f2*h));
} else {
*S = gs_r * (*DA/h) - gs_i * (*(DA+1) / h);
*(S+1) = gs_r * (*(DA+1) /h) + gs_i * -1. * (*DA / h);
}
} else { // intermediates might overflow
}
if ( f2 >= h * safmin) {
*C = sqrt ( f2 / h );
*DA = fs_r / *C;
*(DA+1) = fs_i / *C;
rtmax *= 2;
if ( f2 > rtmin && h < rtmax) {
FLOAT f2h=sqrt(f2*h);
*S = gs_r * (fs_r /f2h) + gs_i * (fs_i / f2h);
*(S+1) = gs_r * (fs_i /f2h) - gs_i * -1. * (fs_r / f2h);
} else {
*S = gs_r * (*DA/h) - gs_i * (*(DA+1) / h);
*(S+1) = gs_r * (*(DA+1) /h) + gs_i * -1. * (*DA / h);
}
} else { // intermediates might overflow
d = sqrt ( f2 * h);
*C = f2 /d;
if (*C >= safmin) {
Expand All @@ -173,14 +178,14 @@ void CNAME(void *VDA, void *VDB, FLOAT *C, void *VS) {
}
*S = gs_r * (fs_r /d) - gs_i * (fs_i / d);
*(S+1) = gs_r * (fs_i /d) + gs_i * -1. * (fs_r / d);
}
*C *= w;
*DA *= u;
*(DA+1) *= u;
break;
}
*C *= w;
*DA *= u;
*(DA+1) *= u;
break;
}
} while (0);
}
} while (0);

FUNCTION_PROFILE_END(4, 4, 4);
IDEBUG_END;
Expand Down
Loading