From 6876ae0c3b77b3fb08185e8619801343e4388868 Mon Sep 17 00:00:00 2001
From: Angelika Schwarz <17718454+angsch@users.noreply.github.com>
Date: Wed, 20 Sep 2023 19:10:08 +0200
Subject: [PATCH 1/2] Fix division by zero in zrotg

The cases
[ c  s ] * [ 0      ] = [ |db_i| ]
[-s  c ]   [ i*db_i ]   [  0     ]
and
[ c  s ] * [ 0      ] = [ |db_r| ]
[-s  c ]   [ db_r   ]   [  0     ]
computed s incorrectly. To flip the entries of vector,
s should be conjg(db)/|db| and not conjg(db) / da,
where da == 0.0.
---
 interface/zrotg.c | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/interface/zrotg.c b/interface/zrotg.c
index af6f85c1ca..4d2a9d5106 100644
--- 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));

From db3a43c8edeb36ecc9e7cde10b1c06be3f2147fc Mon Sep 17 00:00:00 2001
From: Angelika Schwarz <17718454+angsch@users.noreply.github.com>
Date: Wed, 20 Sep 2023 19:42:13 +0200
Subject: [PATCH 2/2] Simplify rotg

* The check da != ZERO is no longer necessary since there
  is a special case ada == ZERO, where ada = |da|.
* Add the missing check c != ZERO before the division.

Note that with these two changes the long double code
follows the float/double version of the code.
---
 interface/rotg.c | 9 ++-------
 1 file changed, 2 insertions(+), 7 deletions(-)

diff --git a/interface/rotg.c b/interface/rotg.c
index 8d40d9c53c..423ebda21d 100644
--- 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;