fix: address bug and improve C implementation of math/base/special/hyp2f1

PR-URL: #11353
Ref: #11325
Reviewed-by: Athan Reines <kgryte@gmail.com>

Author: officiallyanee

lib/node_modules/@stdlib/math/base/special/hyp2f1/manifest.json

@@ -47,6 +47,7 @@
         "@stdlib/math/base/special/ln",
         "@stdlib/math/base/special/abs",
         "@stdlib/math/base/special/pow",
+        "@stdlib/math/base/special/max",
         "@stdlib/constants/float64/pinf",
         "@stdlib/constants/float64/nan"
       ]
@@ -72,6 +73,7 @@
         "@stdlib/math/base/special/ln",
         "@stdlib/math/base/special/abs",
         "@stdlib/math/base/special/pow",
+        "@stdlib/math/base/special/max",
         "@stdlib/constants/float64/pinf",
         "@stdlib/constants/float64/nan"
       ]
@@ -97,6 +99,7 @@
         "@stdlib/math/base/special/ln",
         "@stdlib/math/base/special/abs",
         "@stdlib/math/base/special/pow",
+        "@stdlib/math/base/special/max",
         "@stdlib/constants/float64/pinf",
         "@stdlib/constants/float64/nan"
       ]

lib/node_modules/@stdlib/math/base/special/hyp2f1/src/main.c

@@ -39,6 +39,7 @@
 #include "stdlib/math/base/special/ln.h"
 #include "stdlib/math/base/special/abs.h"
 #include "stdlib/math/base/special/pow.h"
+#include "stdlib/math/base/special/max.h"
 #include "stdlib/math/base/assert/is_nan.h"
 #include "stdlib/constants/float64/pinf.h"
 #include "stdlib/constants/float64/nan.h"
@@ -235,6 +236,38 @@ static double hys2f1( double a, double b, const double c, const double x, double
 	return s;
 }
 
+/**
+* Evaluates 2F1(a, b; b; x) when `b = c` is a negative integer using AMS55 #15.4.2.
+*
+* @param a      first parameter
+* @param b      second parameter (equals c, a non-positive integer)
+* @param x      argument
+* @return       function value, or NaN if precision is insufficient
+*/
+static double hyp2f1NegCEqualBC( const double a, const double b, const double x ) {
+	double collectorMax;
+	double collector;
+	double sum;
+	double k;
+
+	collectorMax = 1.0;
+	collector = 1.0;
+	sum = 1.0;
+
+	if ( stdlib_base_abs( b ) >= 1.0e5 ) {
+		return STDLIB_CONSTANT_FLOAT64_NAN;
+	}
+	for ( k = 1.0; k <= -b; k ++ ) {
+		collector *= ( ( a + k - 1.0 ) * x ) / k;
+		collectorMax = stdlib_base_max( stdlib_base_abs( collector ), collectorMax );
+		sum += collector;
+	}
+	if ( 1.0e-16 * ( 1.0 + ( collectorMax / stdlib_base_abs( sum ) ) ) > 1.0e-7 ) {
+		return STDLIB_CONSTANT_FLOAT64_NAN;
+	}
+	return sum;
+}
+
 /**
 * Applies transformations for `|x|` near unity before performing a power series expansion.
 *
@@ -289,7 +322,7 @@ static double hyt2f1( const double a, const double b, const double c, const doub
 	d = c - a - b;
 	id = stdlib_base_round( d );
 
-	if ( x > 0.9 && !negIntA && !negIntB ) {
+	if ( x > 0.85 && !negIntA && !negIntB ) {
 		if ( isInteger( d ) == false ) {
 			// Try the power series first:
 			y = hys2f1( a, b, c, x, &err );
@@ -386,7 +419,11 @@ static double hyt2f1( const double a, const double b, const double c, const doub
 				y = -y;
 			}
 			q = stdlib_base_pow( s, id );
-			y = ( id > 0.0 ) ? y*q : y1*q;
+			if ( id > 0.0 ) {
+				y *= q;
+			} else {
+				y1 *= q;
+			}
 			y += y1;
 		}
 		*loss = err;
@@ -477,6 +514,11 @@ double stdlib_base_hyp2f1( const double a, const double b, const double c, const
 	if ( ax < 1.0 || x == -1.0 ) {
 		if ( b == c ) {
 			// 2F1(a,b;b;x) = (1-x)**(-a):
+			if ( negIntB ) {
+				// For negative integer b=c use the finite polynomial (AMS55 #15.4.2):
+				y = hyp2f1NegCEqualBC( a, b, x );
+				return y;
+			}
 			y = stdlib_base_pow( s, -a );
 			return y;
 		}

lib/node_modules/@stdlib/math/base/special/hyp2f1/test/test.native.js

@@ -43,6 +43,8 @@ var edgeCases1 = require( './fixtures/python/edge_cases1.json' );
 var edgeCases2 = require( './fixtures/python/edge_cases2.json' );
 var edgeCases3 = require( './fixtures/python/edge_cases3.json' );
 var edgeCases4 = require( './fixtures/python/edge_cases4.json' );
+var edgeCases5 = require( './fixtures/python/edge_cases5.json' );
+var edgeCases6 = require( './fixtures/python/edge_cases6.json' );
 var outliers = require( './fixtures/python/outliers.json' );
 
 
@@ -254,6 +256,50 @@ tape( 'the function correctly evaluates the hypergeometric function', opts, func
 	t.end();
 });
 
+tape( 'the function correctly evaluates the hypergeometric function', function test( t ) {
+	var expected;
+	var a;
+	var b;
+	var c;
+	var x;
+	var v;
+	var i;
+
+	a = edgeCases5.a;
+	b = edgeCases5.b;
+	c = edgeCases5.c;
+	x = edgeCases5.x;
+	expected = edgeCases5.expected;
+
+	for ( i = 0; i < x.length; i++ ) {
+		v = hyp2f1( a[ i ], b[ i ], c[ i ], x[ i ] );
+		t.strictEqual( isAlmostSameValue( v, expected[ i ], 12 ), true, 'returns expected value.' );
+	}
+	t.end();
+});
+
+tape( 'the function correctly evaluates the hypergeometric function (edge case 6)', function test( t ) {
+	var expected;
+	var a;
+	var b;
+	var c;
+	var x;
+	var v;
+	var i;
+
+	a = edgeCases6.a;
+	b = edgeCases6.b;
+	c = edgeCases6.c;
+	x = edgeCases6.x;
+	expected = edgeCases6.expected;
+
+	for ( i = 0; i < x.length; i++ ) {
+		v = hyp2f1( a[ i ], b[ i ], c[ i ], x[ i ] );
+		t.strictEqual( isAlmostSameValue( v, expected[ i ], 1 ), true, 'returns expected value.' );
+	}
+	t.end();
+});
+
 tape( 'the function correctly evaluates the hypergeometric function', opts, function test( t ) {
 	var expected;
 	var a;