Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

New Bessel Algorithm, providing a higher degree of accuracy and precision #1946

Merged
merged 3 commits into from
Mar 24, 2021
Merged

New Bessel Algorithm, providing a higher degree of accuracy and precision #1946

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
24ecbf2
New Bessel Algorithm, providing a higher degree of precision (12 deci…
Mar 23, 2021
9f4e50f
New Bessel Algorithm, providing a higher degree of precision (12+ dec…
Mar 24, 2021
a96712f
Minor refactoring
Mar 24, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
106 changes: 84 additions & 22 deletions src/PhpSpreadsheet/Calculation/Engineering/BesselI.php
View file
Open in desktop
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
Open in desktop
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