Fix catan, catanh spurious overflows (bug 15409).

7 files changed, 338 insertions, 90 deletions

diff --git a/math/libm-test.inc b/math/libm-test.inc
index a62ae9a11f..b91035b489 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -4364,6 +4364,70 @@ catan_test (void)
   TEST_c_c (catan, 0.5L, -0x1p8191L, 1.570796326794896619231321691639751442099L, -1.833603867554847165621412392048483165956e-2466L);
   TEST_c_c (catan, -0.5L, -0x1p8191L, -1.570796326794896619231321691639751442099L, -1.833603867554847165621412392048483165956e-2466L);
 #endif
+  TEST_c_c (catan, 0x1p100L, 0.0L, 1.570796326794896619231321691638962581193L, 0.0L);
+  TEST_c_c (catan, 0x1p100L, -0.0L, 1.570796326794896619231321691638962581193L, -0.0L);
+  TEST_c_c (catan, -0x1p100L, 0.0L, -1.570796326794896619231321691638962581193L, 0.0L);
+  TEST_c_c (catan, -0x1p100L, -0.0L, -1.570796326794896619231321691638962581193L, -0.0L);
+  TEST_c_c (catan, 0.0L, 0x1p100L, 1.570796326794896619231321691639751442099L, 7.888609052210118054117285652827862296732e-31L);
+  TEST_c_c (catan, -0.0L, 0x1p100L, -1.570796326794896619231321691639751442099L, 7.888609052210118054117285652827862296732e-31L);
+  TEST_c_c (catan, 0.0L, -0x1p100L, 1.570796326794896619231321691639751442099L, -7.888609052210118054117285652827862296732e-31L);
+  TEST_c_c (catan, -0.0L, -0x1p100L, -1.570796326794896619231321691639751442099L, -7.888609052210118054117285652827862296732e-31L);
+  TEST_c_c (catan, 0x1.fp127L, 0.0L, 1.570796326794896619231321691639751442096L, 0.0L);
+  TEST_c_c (catan, 0x1.fp127L, -0.0L, 1.570796326794896619231321691639751442096L, -0.0L);
+  TEST_c_c (catan, -0x1.fp127L, 0.0L, -1.570796326794896619231321691639751442096L, 0.0L);
+  TEST_c_c (catan, -0x1.fp127L, -0.0L, -1.570796326794896619231321691639751442096L, -0.0L);
+  TEST_c_c (catan, 0.0L, 0x1.fp127L, 1.570796326794896619231321691639751442099L, 3.033533808573645181854803967025150136306e-39L, UNDERFLOW_EXCEPTION_FLOAT);
+  TEST_c_c (catan, -0.0L, 0x1.fp127L, -1.570796326794896619231321691639751442099L, 3.033533808573645181854803967025150136306e-39L, UNDERFLOW_EXCEPTION_FLOAT);
+  TEST_c_c (catan, 0.0L, -0x1.fp127L, 1.570796326794896619231321691639751442099L, -3.033533808573645181854803967025150136306e-39L, UNDERFLOW_EXCEPTION_FLOAT);
+  TEST_c_c (catan, -0.0L, -0x1.fp127L, -1.570796326794896619231321691639751442099L, -3.033533808573645181854803967025150136306e-39L, UNDERFLOW_EXCEPTION_FLOAT);
+  TEST_c_c (catan, 0x1.fp127L, 0x1.fp127L, 1.570796326794896619231321691639751442097L, 1.516766904286822590927401983512575068153e-39L, UNDERFLOW_EXCEPTION_FLOAT);
+  TEST_c_c (catan, 0x1.fp127L, -0x1.fp127L, 1.570796326794896619231321691639751442097L, -1.516766904286822590927401983512575068153e-39L, UNDERFLOW_EXCEPTION_FLOAT);
+  TEST_c_c (catan, -0x1.fp127L, 0x1.fp127L, -1.570796326794896619231321691639751442097L, 1.516766904286822590927401983512575068153e-39L, UNDERFLOW_EXCEPTION_FLOAT);
+  TEST_c_c (catan, -0x1.fp127L, -0x1.fp127L, -1.570796326794896619231321691639751442097L, -1.516766904286822590927401983512575068153e-39L, UNDERFLOW_EXCEPTION_FLOAT);
+#ifndef TEST_FLOAT
+  TEST_c_c (catan, 0x1p900L, 0.0L, 1.570796326794896619231321691639751442099L, 0.0L);
+  TEST_c_c (catan, 0x1p900L, -0.0L, 1.570796326794896619231321691639751442099L, -0.0L);
+  TEST_c_c (catan, -0x1p900L, 0.0L, -1.570796326794896619231321691639751442099L, 0.0L);
+  TEST_c_c (catan, -0x1p900L, -0.0L, -1.570796326794896619231321691639751442099L, -0.0L);
+  TEST_c_c (catan, 0.0L, 0x1p900L, 1.570796326794896619231321691639751442099L, 1.183052186166774710972751597518026531652e-271L);
+  TEST_c_c (catan, -0.0L, 0x1p900L, -1.570796326794896619231321691639751442099L, 1.183052186166774710972751597518026531652e-271L);
+  TEST_c_c (catan, 0.0L, -0x1p900L, 1.570796326794896619231321691639751442099L, -1.183052186166774710972751597518026531652e-271L);
+  TEST_c_c (catan, -0.0L, -0x1p900L, -1.570796326794896619231321691639751442099L, -1.183052186166774710972751597518026531652e-271L);
+  TEST_c_c (catan, 0x1.fp1023L, 0.0L, 1.570796326794896619231321691639751442099L, 0.0L);
+  TEST_c_c (catan, 0x1.fp1023L, -0.0L, 1.570796326794896619231321691639751442099L, -0.0L);
+  TEST_c_c (catan, -0x1.fp1023L, 0.0L, -1.570796326794896619231321691639751442099L, 0.0L);
+  TEST_c_c (catan, -0x1.fp1023L, -0.0L, -1.570796326794896619231321691639751442099L, -0.0L);
+  TEST_c_c (catan, 0.0L, 0x1.fp1023L, 1.570796326794896619231321691639751442099L, 5.742126086470197117652213464083623391533e-309L, UNDERFLOW_EXCEPTION_DOUBLE);
+  TEST_c_c (catan, -0.0L, 0x1.fp1023L, -1.570796326794896619231321691639751442099L, 5.742126086470197117652213464083623391533e-309L, UNDERFLOW_EXCEPTION_DOUBLE);
+  TEST_c_c (catan, 0.0L, -0x1.fp1023L, 1.570796326794896619231321691639751442099L, -5.742126086470197117652213464083623391533e-309L, UNDERFLOW_EXCEPTION_DOUBLE);
+  TEST_c_c (catan, -0.0L, -0x1.fp1023L, -1.570796326794896619231321691639751442099L, -5.742126086470197117652213464083623391533e-309L, UNDERFLOW_EXCEPTION_DOUBLE);
+  TEST_c_c (catan, 0x1.fp1023L, 0x1.fp1023L, 1.570796326794896619231321691639751442099L, 2.871063043235098558826106732041811695767e-309L, UNDERFLOW_EXCEPTION_DOUBLE);
+  TEST_c_c (catan, 0x1.fp1023L, -0x1.fp1023L, 1.570796326794896619231321691639751442099L, -2.871063043235098558826106732041811695767e-309L, UNDERFLOW_EXCEPTION_DOUBLE);
+  TEST_c_c (catan, -0x1.fp1023L, 0x1.fp1023L, -1.570796326794896619231321691639751442099L, 2.871063043235098558826106732041811695767e-309L, UNDERFLOW_EXCEPTION_DOUBLE);
+  TEST_c_c (catan, -0x1.fp1023L, -0x1.fp1023L, -1.570796326794896619231321691639751442099L, -2.871063043235098558826106732041811695767e-309L, UNDERFLOW_EXCEPTION_DOUBLE);
+#endif
+#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384
+  TEST_c_c (catan, 0x1p10000L, 0.0L, 1.570796326794896619231321691639751442099L, 0.0L);
+  TEST_c_c (catan, 0x1p10000L, -0.0L, 1.570796326794896619231321691639751442099L, -0.0L);
+  TEST_c_c (catan, -0x1p10000L, 0.0L, -1.570796326794896619231321691639751442099L, 0.0L);
+  TEST_c_c (catan, -0x1p10000L, -0.0L, -1.570796326794896619231321691639751442099L, -0.0L);
+  TEST_c_c (catan, 0.0L, 0x1p10000L, 1.570796326794896619231321691639751442099L, 5.012372749206452009297555933742977749322e-3011L);
+  TEST_c_c (catan, -0.0L, 0x1p10000L, -1.570796326794896619231321691639751442099L, 5.012372749206452009297555933742977749322e-3011L);
+  TEST_c_c (catan, 0.0L, -0x1p10000L, 1.570796326794896619231321691639751442099L, -5.012372749206452009297555933742977749322e-3011L);
+  TEST_c_c (catan, -0.0L, -0x1p10000L, -1.570796326794896619231321691639751442099L, -5.012372749206452009297555933742977749322e-3011L);
+  TEST_c_c (catan, 0x1.fp16383L, 0.0L, 1.570796326794896619231321691639751442099L, 0.0L);
+  TEST_c_c (catan, 0x1.fp16383L, -0.0L, 1.570796326794896619231321691639751442099L, -0.0L);
+  TEST_c_c (catan, -0x1.fp16383L, 0.0L, -1.570796326794896619231321691639751442099L, 0.0L);
+  TEST_c_c (catan, -0x1.fp16383L, -0.0L, -1.570796326794896619231321691639751442099L, -0.0L);
+  TEST_c_c (catan, 0.0L, 0x1.fp16383L, 1.570796326794896619231321691639751442099L, 8.676395208031209048419813722120651877672e-4933L, UNDERFLOW_EXCEPTION);
+  TEST_c_c (catan, -0.0L, 0x1.fp16383L, -1.570796326794896619231321691639751442099L, 8.676395208031209048419813722120651877672e-4933L, UNDERFLOW_EXCEPTION);
+  TEST_c_c (catan, 0.0L, -0x1.fp16383L, 1.570796326794896619231321691639751442099L, -8.676395208031209048419813722120651877672e-4933L, UNDERFLOW_EXCEPTION);
+  TEST_c_c (catan, -0.0L, -0x1.fp16383L, -1.570796326794896619231321691639751442099L, -8.676395208031209048419813722120651877672e-4933L, UNDERFLOW_EXCEPTION);
+  TEST_c_c (catan, 0x1.fp16383L, 0x1.fp16383L, 1.570796326794896619231321691639751442099L, 4.338197604015604524209906861060325938836e-4933L, UNDERFLOW_EXCEPTION);
+  TEST_c_c (catan, 0x1.fp16383L, -0x1.fp16383L, 1.570796326794896619231321691639751442099L, -4.338197604015604524209906861060325938836e-4933L, UNDERFLOW_EXCEPTION);
+  TEST_c_c (catan, -0x1.fp16383L, 0x1.fp16383L, -1.570796326794896619231321691639751442099L, 4.338197604015604524209906861060325938836e-4933L, UNDERFLOW_EXCEPTION);
+  TEST_c_c (catan, -0x1.fp16383L, -0x1.fp16383L, -1.570796326794896619231321691639751442099L, -4.338197604015604524209906861060325938836e-4933L, UNDERFLOW_EXCEPTION);
+#endif
 
   TEST_c_c (catan, 0.75L, 1.25L, 1.10714871779409050301706546017853704L, 0.549306144334054845697622618461262852L);
   TEST_c_c (catan, -2, -3, -1.4099210495965755225306193844604208L, -0.22907268296853876629588180294200276L);
@@ -4496,6 +4560,70 @@ catanh_test (void)
   TEST_c_c (catanh, 0.5L, -0x1p8191L, 1.681051571556046753131338908660876301299e-4932L, -1.570796326794896619231321691639751442099L, UNDERFLOW_EXCEPTION);
   TEST_c_c (catanh, -0.5L, -0x1p8191L, -1.681051571556046753131338908660876301299e-4932L, -1.570796326794896619231321691639751442099L, UNDERFLOW_EXCEPTION);
 #endif
+  TEST_c_c (catanh, 0x1p100L, 0.0L, 7.888609052210118054117285652827862296732e-31L, 1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, 0x1p100L, -0.0L, 7.888609052210118054117285652827862296732e-31L, -1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, -0x1p100L, 0.0L, -7.888609052210118054117285652827862296732e-31L, 1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, -0x1p100L, -0.0L, -7.888609052210118054117285652827862296732e-31L, -1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, 0.0L, 0x1p100L, 0.0L, 1.570796326794896619231321691638962581193L);
+  TEST_c_c (catanh, -0.0L, 0x1p100L, -0.0L, 1.570796326794896619231321691638962581193L);
+  TEST_c_c (catanh, 0.0L, -0x1p100L, 0.0L, -1.570796326794896619231321691638962581193L);
+  TEST_c_c (catanh, -0.0L, -0x1p100L, -0.0L, -1.570796326794896619231321691638962581193L);
+  TEST_c_c (catanh, 0x1.fp127L, 0.0L, 3.033533808573645181854803967025150136306e-39L, 1.570796326794896619231321691639751442099L, UNDERFLOW_EXCEPTION_FLOAT);
+  TEST_c_c (catanh, 0x1.fp127L, -0.0L, 3.033533808573645181854803967025150136306e-39L, -1.570796326794896619231321691639751442099L, UNDERFLOW_EXCEPTION_FLOAT);
+  TEST_c_c (catanh, -0x1.fp127L, 0.0L, -3.033533808573645181854803967025150136306e-39L, 1.570796326794896619231321691639751442099L, UNDERFLOW_EXCEPTION_FLOAT);
+  TEST_c_c (catanh, -0x1.fp127L, -0.0L, -3.033533808573645181854803967025150136306e-39L, -1.570796326794896619231321691639751442099L, UNDERFLOW_EXCEPTION_FLOAT);
+  TEST_c_c (catanh, 0.0L, 0x1.fp127L, 0.0L, 1.570796326794896619231321691639751442096L);
+  TEST_c_c (catanh, -0.0L, 0x1.fp127L, -0.0L, 1.570796326794896619231321691639751442096L);
+  TEST_c_c (catanh, 0.0L, -0x1.fp127L, 0.0L, -1.570796326794896619231321691639751442096L);
+  TEST_c_c (catanh, -0.0L, -0x1.fp127L, -0.0L, -1.570796326794896619231321691639751442096L);
+  TEST_c_c (catanh, 0x1.fp127L, 0x1.fp127L, 1.516766904286822590927401983512575068153e-39L, 1.570796326794896619231321691639751442097L, UNDERFLOW_EXCEPTION_FLOAT);
+  TEST_c_c (catanh, 0x1.fp127L, -0x1.fp127L, 1.516766904286822590927401983512575068153e-39L, -1.570796326794896619231321691639751442097L, UNDERFLOW_EXCEPTION_FLOAT);
+  TEST_c_c (catanh, -0x1.fp127L, 0x1.fp127L, -1.516766904286822590927401983512575068153e-39L, 1.570796326794896619231321691639751442097L, UNDERFLOW_EXCEPTION_FLOAT);
+  TEST_c_c (catanh, -0x1.fp127L, -0x1.fp127L, -1.516766904286822590927401983512575068153e-39L, -1.570796326794896619231321691639751442097L, UNDERFLOW_EXCEPTION_FLOAT);
+#ifndef TEST_FLOAT
+  TEST_c_c (catanh, 0x1p900L, 0.0L, 1.183052186166774710972751597518026531652e-271L, 1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, 0x1p900L, -0.0L, 1.183052186166774710972751597518026531652e-271L, -1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, -0x1p900L, 0.0L, -1.183052186166774710972751597518026531652e-271L, 1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, -0x1p900L, -0.0L, -1.183052186166774710972751597518026531652e-271L, -1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, 0.0L, 0x1p900L, 0.0L, 1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, -0.0L, 0x1p900L, -0.0L, 1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, 0.0L, -0x1p900L, 0.0L, -1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, -0.0L, -0x1p900L, -0.0L, -1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, 0x1.fp1023L, 0.0L, 5.742126086470197117652213464083623391533e-309L, 1.570796326794896619231321691639751442099L, UNDERFLOW_EXCEPTION_DOUBLE);
+  TEST_c_c (catanh, 0x1.fp1023L, -0.0L, 5.742126086470197117652213464083623391533e-309L, -1.570796326794896619231321691639751442099L, UNDERFLOW_EXCEPTION_DOUBLE);
+  TEST_c_c (catanh, -0x1.fp1023L, 0.0L, -5.742126086470197117652213464083623391533e-309L, 1.570796326794896619231321691639751442099L, UNDERFLOW_EXCEPTION_DOUBLE);
+  TEST_c_c (catanh, -0x1.fp1023L, -0.0L, -5.742126086470197117652213464083623391533e-309L, -1.570796326794896619231321691639751442099L, UNDERFLOW_EXCEPTION_DOUBLE);
+  TEST_c_c (catanh, 0.0L, 0x1.fp1023L, 0.0L, 1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, -0.0L, 0x1.fp1023L, -0.0L, 1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, 0.0L, -0x1.fp1023L, 0.0L, -1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, -0.0L, -0x1.fp1023L, -0.0L, -1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, 0x1.fp1023L, 0x1.fp1023L, 2.871063043235098558826106732041811695767e-309L, 1.570796326794896619231321691639751442099L, UNDERFLOW_EXCEPTION_DOUBLE);
+  TEST_c_c (catanh, 0x1.fp1023L, -0x1.fp1023L, 2.871063043235098558826106732041811695767e-309L, -1.570796326794896619231321691639751442099L, UNDERFLOW_EXCEPTION_DOUBLE);
+  TEST_c_c (catanh, -0x1.fp1023L, 0x1.fp1023L, -2.871063043235098558826106732041811695767e-309L, 1.570796326794896619231321691639751442099L, UNDERFLOW_EXCEPTION_DOUBLE);
+  TEST_c_c (catanh, -0x1.fp1023L, -0x1.fp1023L, -2.871063043235098558826106732041811695767e-309L, -1.570796326794896619231321691639751442099L, UNDERFLOW_EXCEPTION_DOUBLE);
+#endif
+#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384
+  TEST_c_c (catanh, 0x1p10000L, 0.0L, 5.012372749206452009297555933742977749322e-3011L, 1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, 0x1p10000L, -0.0L, 5.012372749206452009297555933742977749322e-3011L, -1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, -0x1p10000L, 0.0L, -5.012372749206452009297555933742977749322e-3011L, 1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, -0x1p10000L, -0.0L, -5.012372749206452009297555933742977749322e-3011L, -1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, 0.0L, 0x1p10000L, 0.0L, 1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, -0.0L, 0x1p10000L, -0.0L, 1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, 0.0L, -0x1p10000L, 0.0L, -1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, -0.0L, -0x1p10000L, -0.0L, -1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, 0x1.fp16383L, 0.0L, 8.676395208031209048419813722120651877672e-4933L, 1.570796326794896619231321691639751442099L, UNDERFLOW_EXCEPTION);
+  TEST_c_c (catanh, 0x1.fp16383L, -0.0L, 8.676395208031209048419813722120651877672e-4933L, -1.570796326794896619231321691639751442099L, UNDERFLOW_EXCEPTION);
+  TEST_c_c (catanh, -0x1.fp16383L, 0.0L, -8.676395208031209048419813722120651877672e-4933L, 1.570796326794896619231321691639751442099L, UNDERFLOW_EXCEPTION);
+  TEST_c_c (catanh, -0x1.fp16383L, -0.0L, -8.676395208031209048419813722120651877672e-4933L, -1.570796326794896619231321691639751442099L, UNDERFLOW_EXCEPTION);
+  TEST_c_c (catanh, 0.0L, 0x1.fp16383L, 0.0L, 1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, -0.0L, 0x1.fp16383L, -0.0L, 1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, 0.0L, -0x1.fp16383L, 0.0L, -1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, -0.0L, -0x1.fp16383L, -0.0L, -1.570796326794896619231321691639751442099L);
+  TEST_c_c (catanh, 0x1.fp16383L, 0x1.fp16383L, 4.338197604015604524209906861060325938836e-4933L, 1.570796326794896619231321691639751442099L, UNDERFLOW_EXCEPTION);
+  TEST_c_c (catanh, 0x1.fp16383L, -0x1.fp16383L, 4.338197604015604524209906861060325938836e-4933L, -1.570796326794896619231321691639751442099L, UNDERFLOW_EXCEPTION);
+  TEST_c_c (catanh, -0x1.fp16383L, 0x1.fp16383L, -4.338197604015604524209906861060325938836e-4933L, 1.570796326794896619231321691639751442099L, UNDERFLOW_EXCEPTION);
+  TEST_c_c (catanh, -0x1.fp16383L, -0x1.fp16383L, -4.338197604015604524209906861060325938836e-4933L, -1.570796326794896619231321691639751442099L, UNDERFLOW_EXCEPTION);
+#endif
 
   TEST_c_c (catanh, 0.75L, 1.25L, 0.261492138795671927078652057366532140L, 0.996825126463918666098902241310446708L);
   TEST_c_c (catanh, -2, -3, -0.14694666622552975204743278515471595L, -1.3389725222944935611241935759091443L);
diff --git a/math/s_catan.c b/math/s_catan.c
index 6a1016e74c..4b7fbad30a 100644
--- a/math/s_catan.c
+++ b/math/s_catan.c
@@ -61,27 +61,44 @@ __catan (__complex__ double x)
     }
   else
     {
-      double r2, num, den, f;
+      if (fabs (__real__ x) >= 16.0 / DBL_EPSILON
+	  || fabs (__imag__ x) >= 16.0 / DBL_EPSILON)
+	{
+	  __real__ res = __copysign (M_PI_2, __real__ x);
+	  if (fabs (__real__ x) <= 1.0)
+	    __imag__ res = 1.0 / __imag__ x;
+	  else if (fabs (__imag__ x) <= 1.0)
+	    __imag__ res = __imag__ x / __real__ x / __real__ x;
+	  else
+	    {
+	      double h = __ieee754_hypot (__real__ x / 2.0, __imag__ x / 2.0);
+	      __imag__ res = __imag__ x / h / h / 4.0;
+	    }
+	}
+      else
+	{
+	  double r2, num, den, f;
 
-      r2 = __real__ x * __real__ x;
+	  r2 = __real__ x * __real__ x;
 
-      den = 1 - r2 - __imag__ x * __imag__ x;
+	  den = 1 - r2 - __imag__ x * __imag__ x;
 
-      __real__ res = 0.5 * __ieee754_atan2 (2.0 * __real__ x, den);
+	  __real__ res = 0.5 * __ieee754_atan2 (2.0 * __real__ x, den);
 
-      num = __imag__ x + 1.0;
-      num = r2 + num * num;
+	  num = __imag__ x + 1.0;
+	  num = r2 + num * num;
 
-      den = __imag__ x - 1.0;
-      den = r2 + den * den;
+	  den = __imag__ x - 1.0;
+	  den = r2 + den * den;
 
-      f = num / den;
-      if (f < 0.5)
-	__imag__ res = 0.25 * __ieee754_log (f);
-      else
-	{
-	  num = 4.0 * __imag__ x;
-	  __imag__ res = 0.25 * __log1p (num / den);
+	  f = num / den;
+	  if (f < 0.5)
+	    __imag__ res = 0.25 * __ieee754_log (f);
+	  else
+	    {
+	      num = 4.0 * __imag__ x;
+	      __imag__ res = 0.25 * __log1p (num / den);
+	    }
 	}
 
       if (fabs (__real__ res) < DBL_MIN)
diff --git a/math/s_catanf.c b/math/s_catanf.c
index d0d188ca0b..aa71e5e63a 100644
--- a/math/s_catanf.c
+++ b/math/s_catanf.c
@@ -61,27 +61,45 @@ __catanf (__complex__ float x)
     }
   else
     {
-      float r2, num, den, f;
+      if (fabsf (__real__ x) >= 16.0f / FLT_EPSILON
+	  || fabsf (__imag__ x) >= 16.0f / FLT_EPSILON)
+	{
+	  __real__ res = __copysignf ((float) M_PI_2, __real__ x);
+	  if (fabsf (__real__ x) <= 1.0f)
+	    __imag__ res = 1.0f / __imag__ x;
+	  else if (fabsf (__imag__ x) <= 1.0f)
+	    __imag__ res = __imag__ x / __real__ x / __real__ x;
+	  else
+	    {
+	      float h = __ieee754_hypotf (__real__ x / 2.0f,
+					  __imag__ x / 2.0f);
+	      __imag__ res = __imag__ x / h / h / 4.0f;
+	    }
+	}
+      else
+	{
+	  float r2, num, den, f;
 
-      r2 = __real__ x * __real__ x;
+	  r2 = __real__ x * __real__ x;
 
-      den = 1 - r2 - __imag__ x * __imag__ x;
+	  den = 1 - r2 - __imag__ x * __imag__ x;
 
-      __real__ res = 0.5f * __ieee754_atan2f (2.0f * __real__ x, den);
+	  __real__ res = 0.5f * __ieee754_atan2f (2.0f * __real__ x, den);
 
-      num = __imag__ x + 1.0f;
-      num = r2 + num * num;
+	  num = __imag__ x + 1.0f;
+	  num = r2 + num * num;
 
-      den = __imag__ x - 1.0f;
-      den = r2 + den * den;
+	  den = __imag__ x - 1.0f;
+	  den = r2 + den * den;
 
-      f = num / den;
-      if (f < 0.5f)
-	__imag__ res = 0.25f * __ieee754_logf (f);
-      else
-	{
-	  num = 4.0f * __imag__ x;
-	  __imag__ res = 0.25f * __log1pf (num / den);
+	  f = num / den;
+	  if (f < 0.5f)
+	    __imag__ res = 0.25f * __ieee754_logf (f);
+	  else
+	    {
+	      num = 4.0f * __imag__ x;
+	      __imag__ res = 0.25f * __log1pf (num / den);
+	    }
 	}
 
       if (fabsf (__real__ res) < FLT_MIN)
diff --git a/math/s_catanh.c b/math/s_catanh.c
index 5d18cd639f..55c56283b3 100644
--- a/math/s_catanh.c
+++ b/math/s_catanh.c
@@ -56,26 +56,43 @@ __catanh (__complex__ double x)
     }
   else
     {
-      double i2 = __imag__ x * __imag__ x;
+      if (fabs (__real__ x) >= 16.0 / DBL_EPSILON
+	  || fabs (__imag__ x) >= 16.0 / DBL_EPSILON)
+	{
+	  __imag__ res = __copysign (M_PI_2, __imag__ x);
+	  if (fabs (__imag__ x) <= 1.0)
+	    __real__ res = 1.0 / __real__ x;
+	  else if (fabs (__real__ x) <= 1.0)
+	    __real__ res = __real__ x / __imag__ x / __imag__ x;
+	  else
+	    {
+	      double h = __ieee754_hypot (__real__ x / 2.0, __imag__ x / 2.0);
+	      __real__ res = __real__ x / h / h / 4.0;
+	    }
+	}
+      else
+	{
+	  double i2 = __imag__ x * __imag__ x;
 
-      double num = 1.0 + __real__ x;
-      num = i2 + num * num;
+	  double num = 1.0 + __real__ x;
+	  num = i2 + num * num;
 
-      double den = 1.0 - __real__ x;
-      den = i2 + den * den;
+	  double den = 1.0 - __real__ x;
+	  den = i2 + den * den;
 
-      double f = num / den;
-      if (f < 0.5)
-	__real__ res = 0.25 * __ieee754_log (f);
-      else
-	{
-	  num = 4.0 * __real__ x;
-	  __real__ res = 0.25 * __log1p (num / den);
-	}
+	  double f = num / den;
+	  if (f < 0.5)
+	    __real__ res = 0.25 * __ieee754_log (f);
+	  else
+	    {
+	      num = 4.0 * __real__ x;
+	      __real__ res = 0.25 * __log1p (num / den);
+	    }
 
-      den = 1 - __real__ x * __real__ x - i2;
+	  den = 1 - __real__ x * __real__ x - i2;
 
-      __imag__ res = 0.5 * __ieee754_atan2 (2.0 * __imag__ x, den);
+	  __imag__ res = 0.5 * __ieee754_atan2 (2.0 * __imag__ x, den);
+	}
 
       if (fabs (__real__ res) < DBL_MIN)
 	{
diff --git a/math/s_catanhf.c b/math/s_catanhf.c
index 6376a2ccb7..3e87942369 100644
--- a/math/s_catanhf.c
+++ b/math/s_catanhf.c
@@ -56,26 +56,44 @@ __catanhf (__complex__ float x)
     }
   else
     {
-      float i2 = __imag__ x * __imag__ x;
+      if (fabsf (__real__ x) >= 16.0f / FLT_EPSILON
+	  || fabsf (__imag__ x) >= 16.0f / FLT_EPSILON)
+	{
+	  __imag__ res = __copysignf ((float) M_PI_2, __imag__ x);
+	  if (fabsf (__imag__ x) <= 1.0f)
+	    __real__ res = 1.0f / __real__ x;
+	  else if (fabsf (__real__ x) <= 1.0f)
+	    __real__ res = __real__ x / __imag__ x / __imag__ x;
+	  else
+	    {
+	      float h = __ieee754_hypotf (__real__ x / 2.0f,
+					  __imag__ x / 2.0f);
+	      __real__ res = __real__ x / h / h / 4.0f;
+	    }
+	}
+      else
+	{
+	  float i2 = __imag__ x * __imag__ x;
 
-      float num = 1.0f + __real__ x;
-      num = i2 + num * num;
+	  float num = 1.0f + __real__ x;
+	  num = i2 + num * num;
 
-      float den = 1.0f - __real__ x;
-      den = i2 + den * den;
+	  float den = 1.0f - __real__ x;
+	  den = i2 + den * den;
 
-      float f = num / den;
-      if (f < 0.5f)
-	__real__ res = 0.25f * __ieee754_logf (f);
-      else
-	{
-	  num = 4.0f * __real__ x;
-	  __real__ res = 0.25f * __log1pf (num / den);
-	}
+	  float f = num / den;
+	  if (f < 0.5f)
+	    __real__ res = 0.25f * __ieee754_logf (f);
+	  else
+	    {
+	      num = 4.0f * __real__ x;
+	      __real__ res = 0.25f * __log1pf (num / den);
+	    }
 
-      den = 1 - __real__ x * __real__ x - i2;
+	  den = 1 - __real__ x * __real__ x - i2;
 
-      __imag__ res = 0.5f * __ieee754_atan2f (2.0f * __imag__ x, den);
+	  __imag__ res = 0.5f * __ieee754_atan2f (2.0f * __imag__ x, den);
+	}
 
       if (fabsf (__real__ res) < FLT_MIN)
 	{
diff --git a/math/s_catanhl.c b/math/s_catanhl.c
index 30fd27703b..64c30b5573 100644
--- a/math/s_catanhl.c
+++ b/math/s_catanhl.c
@@ -22,6 +22,13 @@
 #include <math_private.h>
 #include <float.h>
 
+/* To avoid spurious overflows, use this definition to treat IBM long
+   double as approximating an IEEE-style format.  */
+#if LDBL_MANT_DIG == 106
+# undef LDBL_EPSILON
+# define LDBL_EPSILON 0x1p-106L
+#endif
+
 __complex__ long double
 __catanhl (__complex__ long double x)
 {
@@ -56,26 +63,44 @@ __catanhl (__complex__ long double x)
     }
   else
     {
-      long double i2 = __imag__ x * __imag__ x;
+      if (fabsl (__real__ x) >= 16.0L / LDBL_EPSILON
+	  || fabsl (__imag__ x) >= 16.0L / LDBL_EPSILON)
+	{
+	  __imag__ res = __copysignl (M_PI_2l, __imag__ x);
+	  if (fabsl (__imag__ x) <= 1.0L)
+	    __real__ res = 1.0L / __real__ x;
+	  else if (fabsl (__real__ x) <= 1.0L)
+	    __real__ res = __real__ x / __imag__ x / __imag__ x;
+	  else
+	    {
+	      long double h = __ieee754_hypotl (__real__ x / 2.0L,
+						__imag__ x / 2.0L);
+	      __real__ res = __real__ x / h / h / 4.0L;
+	    }
+	}
+      else
+	{
+	  long double i2 = __imag__ x * __imag__ x;
 
-      long double num = 1.0L + __real__ x;
-      num = i2 + num * num;
+	  long double num = 1.0L + __real__ x;
+	  num = i2 + num * num;
 
-      long double den = 1.0L - __real__ x;
-      den = i2 + den * den;
+	  long double den = 1.0L - __real__ x;
+	  den = i2 + den * den;
 
-      long double f = num / den;
-      if (f < 0.5L)
-	__real__ res = 0.25L * __ieee754_logl (f);
-      else
-	{
-	  num = 4.0L * __real__ x;
-	  __real__ res = 0.25L * __log1pl (num / den);
-	}
+	  long double f = num / den;
+	  if (f < 0.5L)
+	    __real__ res = 0.25L * __ieee754_logl (f);
+	  else
+	    {
+	      num = 4.0L * __real__ x;
+	      __real__ res = 0.25L * __log1pl (num / den);
+	    }
 
-      den = 1 - __real__ x * __real__ x - i2;
+	  den = 1 - __real__ x * __real__ x - i2;
 
-      __imag__ res = 0.5L * __ieee754_atan2l (2.0L * __imag__ x, den);
+	  __imag__ res = 0.5L * __ieee754_atan2l (2.0L * __imag__ x, den);
+	}
 
       if (fabsl (__real__ res) < LDBL_MIN)
 	{
diff --git a/math/s_catanl.c b/math/s_catanl.c
index 2440e6d151..17ab62d340 100644
--- a/math/s_catanl.c
+++ b/math/s_catanl.c
@@ -22,6 +22,13 @@
 #include <math_private.h>
 #include <float.h>
 
+/* To avoid spurious overflows, use this definition to treat IBM long
+   double as approximating an IEEE-style format.  */
+#if LDBL_MANT_DIG == 106
+# undef LDBL_EPSILON
+# define LDBL_EPSILON 0x1p-106L
+#endif
+
 __complex__ long double
 __catanl (__complex__ long double x)
 {
@@ -61,27 +68,45 @@ __catanl (__complex__ long double x)
     }
   else
     {
-      long double r2, num, den, f;
+      if (fabsl (__real__ x) >= 16.0L / LDBL_EPSILON
+	  || fabsl (__imag__ x) >= 16.0L / LDBL_EPSILON)
+	{
+	  __real__ res = __copysignl (M_PI_2l, __real__ x);
+	  if (fabsl (__real__ x) <= 1.0L)
+	    __imag__ res = 1.0L / __imag__ x;
+	  else if (fabsl (__imag__ x) <= 1.0L)
+	    __imag__ res = __imag__ x / __real__ x / __real__ x;
+	  else
+	    {
+	      long double h = __ieee754_hypotl (__real__ x / 2.0L,
+						__imag__ x / 2.0L);
+	      __imag__ res = __imag__ x / h / h / 4.0L;
+	    }
+	}
+      else
+	{
+	  long double r2, num, den, f;
 
-      r2 = __real__ x * __real__ x;
+	  r2 = __real__ x * __real__ x;
 
-      den = 1 - r2 - __imag__ x * __imag__ x;
+	  den = 1 - r2 - __imag__ x * __imag__ x;
 
-      __real__ res = 0.5L * __ieee754_atan2l (2.0L * __real__ x, den);
+	  __real__ res = 0.5L * __ieee754_atan2l (2.0L * __real__ x, den);
 
-      num = __imag__ x + 1.0L;
-      num = r2 + num * num;
+	  num = __imag__ x + 1.0L;
+	  num = r2 + num * num;
 
-      den = __imag__ x - 1.0L;
-      den = r2 + den * den;
+	  den = __imag__ x - 1.0L;
+	  den = r2 + den * den;
 
-      f = num / den;
-      if (f < 0.5L)
-	__imag__ res = 0.25L * __ieee754_logl (f);
-      else
-	{
-	  num = 4.0L * __imag__ x;
-	  __imag__ res = 0.25L * __log1pl (num / den);
+	  f = num / den;
+	  if (f < 0.5L)
+	    __imag__ res = 0.25L * __ieee754_logl (f);
+	  else
+	    {
+	      num = 4.0L * __imag__ x;
+	      __imag__ res = 0.25L * __log1pl (num / den);
+	    }
 	}
 
       if (fabsl (__real__ res) < LDBL_MIN)