Merge pull request OpenMathLib#4235 from angsch/develop

Fix division by zero in [z]rotg
martin-frbg · Oct 9, 2023 · 4a0f863 · 4a0f863
2 parents 617294b + db3a43c
commit 4a0f863
Show file tree

Hide file tree

Showing 2 changed files with 7 additions and 12 deletions.
diff --git a/interface/rotg.c b/interface/rotg.c
@@ -66,13 +66,8 @@ void CNAME(FLOAT *DA, FLOAT *DB, FLOAT *C, FLOAT *S){
     c = da / r;
     s = db / r;
     z = ONE;
-    if (da != ZERO) {
-      if (ada > adb){
-	z = s;
-      } else {
-	z = ONE / c;
-      }
-    }
+    if (ada > adb) z = s;
+    if ((ada <= adb) && (c != ZERO)) z = ONE / c;
 
     *C = c;
     *S = s;

diff --git a/interface/zrotg.c b/interface/zrotg.c
@@ -61,16 +61,16 @@ void CNAME(void *VDA, void *VDB, FLOAT *C, void *VS) {
 		*(S1 + 0) = *(DB + 0);
 		*(S1 + 1) = *(DB + 1) *-1;
 	if (da_r == ZERO && da_i == ZERO) {
-	    *C = ZERO;	  
+	    *C = ZERO;
 	    if (db_r == ZERO) {
 		    (*DA) = fabsl(db_i);
-		*S = *S1 /da_r;
-		*(S+1) = *(S1+1) /da_r;
+		*S = *S1 /(*DA);
+		*(S+1) = *(S1+1) /(*DA);
 		return;
 	    } else if ( db_i == ZERO) {
 		    *DA = fabsl(db_r);
-		*S = *S1 /da_r;
-		*(S+1) = *(S1+1) /da_r;
+		*S = *S1 /(*DA);
+		*(S+1) = *(S1+1) /(*DA);
 		return;
 	    } else {
 	        long double g1 = MAX( fabsl(db_r), fabsl(db_i));