Skip to content

Commit

Permalink
New Bessel Algorithm, providing a higher degree of accuracy and preci…
Browse files Browse the repository at this point in the history 
…sion (#1946)

* New Bessel Algorithm, providing a higher degree of precision (12 decimal places) and still matching/exceeding MS Excel's precision across the range of values
  • Loading branch information
Mark Baker authored Mar 24, 2021
1 parent 1a7b9a4 commit 07ad800
Show file tree
Hide file tree
Showing 12 changed files with 929 additions and 492 deletions.
106 changes: 84 additions & 22 deletions src/PhpSpreadsheet/Calculation/Engineering/BesselI.php
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,6 @@
namespace PhpOffice\PhpSpreadsheet\Calculation\Engineering;

use PhpOffice\PhpSpreadsheet\Calculation\Functions;
use PhpOffice\PhpSpreadsheet\Calculation\MathTrig;

class BesselI
{
Expand All @@ -16,9 +15,12 @@ class BesselI
* Excel Function:
* BESSELI(x,ord)
*
* @param float $x The value at which to evaluate the function.
* NOTE: The MS Excel implementation of the BESSELI function is still not accurate.
* This code provides a more accurate calculation
*
* @param mixed (float) $x The value at which to evaluate the function.
* If x is nonnumeric, BESSELI returns the #VALUE! error value.
* @param int $ord The order of the Bessel function.
* @param mixed (int) $ord The order of the Bessel function.
* If ord is not an integer, it is truncated.
* If $ord is nonnumeric, BESSELI returns the #VALUE! error value.
* If $ord < 0, BESSELI returns the #NUM! error value.
Expand All @@ -28,15 +30,15 @@ class BesselI
public static function BESSELI($x, $ord)
{
$x = ($x === null) ? 0.0 : Functions::flattenSingleValue($x);
$ord = ($ord === null) ? 0.0 : Functions::flattenSingleValue($ord);
$ord = ($ord === null) ? 0 : Functions::flattenSingleValue($ord);

if ((is_numeric($x)) && (is_numeric($ord))) {
$ord = (int) floor($ord);
if ($ord < 0) {
return Functions::NAN();
}

$fResult = self::calculate($x, $ord);
$fResult = self::calculate((float) $x, $ord);

return (is_nan($fResult)) ? Functions::NAN() : $fResult;
}
Expand All @@ -46,27 +48,87 @@ public static function BESSELI($x, $ord)

private static function calculate(float $x, int $ord): float
{
if (abs($x) <= 30) {
$fResult = $fTerm = ($x / 2) ** $ord / MathTrig::FACT($ord);
$ordK = 1;
$fSqrX = ($x * $x) / 4;
do {
$fTerm *= $fSqrX;
$fTerm /= ($ordK * ($ordK + $ord));
$fResult += $fTerm;
} while ((abs($fTerm) > 1e-12) && (++$ordK < 100));

return $fResult;
// special cases
switch ($ord) {
case 0:
return self::besselI0($x);
case 1:
return self::besselI1($x);
}

return self::besselI2($x, $ord);
}

private static function besselI0(float $x): float
{
$ax = abs($x);

if ($ax < 3.75) {
$y = $x / 3.75;
$y = $y * $y;

return 1.0 + $y * (3.5156229 + $y * (3.0899424 + $y * (1.2067492
+ $y * (0.2659732 + $y * (0.360768e-1 + $y * 0.45813e-2)))));
}

$f_2_PI = 2 * M_PI;
$y = 3.75 / $ax;

return (exp($ax) / sqrt($ax)) * (0.39894228 + $y * (0.1328592e-1 + $y * (0.225319e-2 + $y * (-0.157565e-2
+ $y * (0.916281e-2 + $y * (-0.2057706e-1 + $y * (0.2635537e-1 +
$y * (-0.1647633e-1 + $y * 0.392377e-2))))))));
}

$fXAbs = abs($x);
$fResult = exp($fXAbs) / sqrt($f_2_PI * $fXAbs);
if (($ord & 1) && ($x < 0)) {
$fResult = -$fResult;
private static function besselI1(float $x): float
{
$ax = abs($x);

if ($ax < 3.75) {
$y = $x / 3.75;
$y = $y * $y;
$ans = $ax * (0.5 + $y * (0.87890594 + $y * (0.51498869 + $y * (0.15084934 + $y * (0.2658733e-1 +
$y * (0.301532e-2 + $y * 0.32411e-3))))));

return ($x < 0.0) ? -$ans : $ans;
}

return $fResult;
$y = 3.75 / $ax;
$ans = 0.2282967e-1 + $y * (-0.2895312e-1 + $y * (0.1787654e-1 - $y * 0.420059e-2));
$ans = 0.39894228 + $y * (-0.3988024e-1 + $y * (-0.362018e-2 + $y * (0.163801e-2 +
$y * (-0.1031555e-1 + $y * $ans))));
$ans *= exp($ax) / sqrt($ax);

return ($x < 0.0) ? -$ans : $ans;
}

private static function besselI2(float $x, int $ord): float
{
if ($x === 0.0) {
return 0.0;
}

$tox = 2.0 / abs($x);
$bip = 0;
$ans = 0.0;
$bi = 1.0;

for ($j = 2 * ($ord + (int) sqrt(40.0 * $ord)); $j > 0; --$j) {
$bim = $bip + $j * $tox * $bi;
$bip = $bi;
$bi = $bim;

if (abs($bi) > 1.0e+12) {
$ans *= 1.0e-12;
$bi *= 1.0e-12;
$bip *= 1.0e-12;
}

if ($j === $ord) {
$ans = $bip;
}
}

$ans *= self::besselI0($x) / $bi;

return ($x < 0.0 && (($ord % 2) === 1)) ? -$ans : $ans;
}
}
141 changes: 119 additions & 22 deletions src/PhpSpreadsheet/Calculation/Engineering/BesselJ.php
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,6 @@
namespace PhpOffice\PhpSpreadsheet\Calculation\Engineering;

use PhpOffice\PhpSpreadsheet\Calculation\Functions;
use PhpOffice\PhpSpreadsheet\Calculation\MathTrig;

class BesselJ
{
Expand All @@ -15,9 +14,12 @@ class BesselJ
* Excel Function:
* BESSELJ(x,ord)
*
* @param float $x The value at which to evaluate the function.
* NOTE: The MS Excel implementation of the BESSELJ function is still not accurate, particularly for higher order
* values with x < -8 and x > 8. This code provides a more accurate calculation
*
* @param mixed (float) $x The value at which to evaluate the function.
* If x is nonnumeric, BESSELJ returns the #VALUE! error value.
* @param int $ord The order of the Bessel function. If n is not an integer, it is truncated.
* @param mixed (int) $ord The order of the Bessel function. If n is not an integer, it is truncated.
* If $ord is nonnumeric, BESSELJ returns the #VALUE! error value.
* If $ord < 0, BESSELJ returns the #NUM! error value.
*
Expand All @@ -34,7 +36,7 @@ public static function BESSELJ($x, $ord)
return Functions::NAN();
}

$fResult = self::calculate($x, $ord);
$fResult = self::calculate((float) $x, $ord);

return (is_nan($fResult)) ? Functions::NAN() : $fResult;
}
Expand All @@ -44,28 +46,123 @@ public static function BESSELJ($x, $ord)

private static function calculate(float $x, int $ord): float
{
if (abs($x) <= 30) {
$fResult = $fTerm = ($x / 2) ** $ord / MathTrig::FACT($ord);
$ordK = 1;
$fSqrX = ($x * $x) / -4;
do {
$fTerm *= $fSqrX;
$fTerm /= ($ordK * ($ordK + $ord));
$fResult += $fTerm;
} while ((abs($fTerm) > 1e-12) && (++$ordK < 100));

return $fResult;
// special cases
switch ($ord) {
case 0:
return self::besselJ0($x);
case 1:
return self::besselJ1($x);
}

return self::besselJ2($x, $ord);
}

private static function besselJ0(float $x): float
{
$ax = abs($x);

if ($ax < 8.0) {
$y = $x * $x;
$ans1 = 57568490574.0 + $y * (-13362590354.0 + $y * (651619640.7 + $y * (-11214424.18 + $y *
(77392.33017 + $y * (-184.9052456)))));
$ans2 = 57568490411.0 + $y * (1029532985.0 + $y * (9494680.718 + $y * (59272.64853 + $y *
(267.8532712 + $y * 1.0))));

return $ans1 / $ans2;
}

$z = 8.0 / $ax;
$y = $z * $z;
$xx = $ax - 0.785398164;
$ans1 = 1.0 + $y * (-0.1098628627e-2 + $y * (0.2734510407e-4 + $y * (-0.2073370639e-5 + $y * 0.2093887211e-6)));
$ans2 = -0.1562499995e-1 + $y * (0.1430488765e-3 + $y * (-0.6911147651e-5 + $y *
(0.7621095161e-6 - $y * 0.934935152e-7)));

return sqrt(0.636619772 / $ax) * (cos($xx) * $ans1 - $z * sin($xx) * $ans2);
}

private static function besselJ1(float $x): float
{
$ax = abs($x);

if ($ax < 8.0) {
$y = $x * $x;
$ans1 = $x * (72362614232.0 + $y * (-7895059235.0 + $y * (242396853.1 + $y *
(-2972611.439 + $y * (15704.48260 + $y * (-30.16036606))))));
$ans2 = 144725228442.0 + $y * (2300535178.0 + $y * (18583304.74 + $y * (99447.43394 + $y *
(376.9991397 + $y * 1.0))));

return $ans1 / $ans2;
}

$f_PI_DIV_2 = M_PI / 2;
$f_PI_DIV_4 = M_PI / 4;
$z = 8.0 / $ax;
$y = $z * $z;
$xx = $ax - 2.356194491;

$fXAbs = abs($x);
$fResult = sqrt(Functions::M_2DIVPI / $fXAbs) * cos($fXAbs - $ord * $f_PI_DIV_2 - $f_PI_DIV_4);
if (($ord & 1) && ($x < 0)) {
$fResult = -$fResult;
$ans1 = 1.0 + $y * (0.183105e-2 + $y * (-0.3516396496e-4 + $y * (0.2457520174e-5 + $y * (-0.240337019e-6))));
$ans2 = 0.04687499995 + $y * (-0.2002690873e-3 + $y * (0.8449199096e-5 + $y *
(-0.88228987e-6 + $y * 0.105787412e-6)));
$ans = sqrt(0.636619772 / $ax) * (cos($xx) * $ans1 - $z * sin($xx) * $ans2);

return ($x < 0.0) ? -$ans : $ans;
}

private static function besselJ2(float $x, int $ord): float
{
$ax = abs($x);
if ($ax === 0.0) {
return 0.0;
}

if ($ax > $ord) {
return self::besselj2a($ax, $ord, $x);
}

return self::besselj2b($ax, $ord, $x);
}

private static function besselj2a(float $ax, int $ord, float $x)
{
$tox = 2.0 / $ax;
$bjm = self::besselJ0($ax);
$bj = self::besselJ1($ax);
for ($j = 1; $j < $ord; ++$j) {
$bjp = $j * $tox * $bj - $bjm;
$bjm = $bj;
$bj = $bjp;
}
$ans = $bj;

return ($x < 0.0 && ($ord % 2) == 1) ? -$ans : $ans;
}

private static function besselj2b(float $ax, int $ord, float $x)
{
$tox = 2.0 / $ax;
$jsum = false;
$bjp = $ans = $sum = 0.0;
$bj = 1.0;
for ($j = 2 * ($ord + (int) sqrt(40.0 * $ord)); $j > 0; --$j) {
$bjm = $j * $tox * $bj - $bjp;
$bjp = $bj;
$bj = $bjm;
if (abs($bj) > 1.0e+10) {
$bj *= 1.0e-10;
$bjp *= 1.0e-10;
$ans *= 1.0e-10;
$sum *= 1.0e-10;
}
if ($jsum === true) {
$sum += $bj;
}
$jsum = !$jsum;
if ($j === $ord) {
$ans = $bjp;
}
}
$sum = 2.0 * $sum - $bj;
$ans /= $sum;

return $fResult;
return ($x < 0.0 && ($ord % 2) === 1) ? -$ans : $ans;
}
}
Loading

0 comments on commit 07ad800

Please sign in to comment.