VCVRack
/
Rack
mirror of https://github.com/VCVRack/Rack.git
1
0
Fork 0
Code Issues 0 Releases 43 Wiki Activity
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
1347 Commits
4 Branches
540MB
Tree: 5b7b40e992
Branches Tags
${ item.name }
Create branch ${ searchTerm }
from '5b7b40e992'
${ noResults }
Rack/plugins/community/repos/Edge/dep/Lambert/FukushimaLambertW.cc

326 lines
7.2KB
Raw Blame History 

  1. /*

  2.   Implementation of the Fukushima method for the Lambert W function

  3. 

  4.   Copyright (C) 2015 Darko Veberic, darko.veberic@ijs.si

  5. 

  6.   This program is free software: you can redistribute it and/or modify

  7.   it under the terms of the GNU General Public License as published by

  8.   the Free Software Foundation, either version 3 of the License, or

  9.   (at your option) any later version.

  10. 

  11.   This program is distributed in the hope that it will be useful,

  12.   but WITHOUT ANY WARRANTY; without even the implied warranty of

  13.   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

  14.   GNU General Public License for more details.

  15. 

  16.   You should have received a copy of the GNU General Public License

  17.   along with this program.  If not, see <http://www.gnu.org/licenses/>.

  18. */

  19. 

  20. /*

  21.   This code is based on the following publication and its author's fortran code:

  22.   Toshio Fukushima, "Precise and fast computation of Lambert W-functions without

  23.   transcendental function evaluations", J. Comp. Appl. Math. 244 (2013) 77-89.

  24. */

  25. 

  26. #include <FukushimaLambertW.h>

  27. #include <cmath>

  28. #include <limits>

  29. #include <iostream>

  30. 

  31. using namespace std;

  32. 

  33. 

  34. namespace Fukushima {

  35. 

  36.   double

  37.   LambertWSeries(const double p)

  38.   {

  39.     static const double q[] = {

  40.       -1,

  41.       +1,

  42.       -0.333333333333333333,

  43.       +0.152777777777777778,

  44.       -0.0796296296296296296,

  45.       +0.0445023148148148148,

  46.       -0.0259847148736037625,

  47.       +0.0156356325323339212,

  48.       -0.00961689202429943171,

  49.       +0.00601454325295611786,

  50.       -0.00381129803489199923,

  51.       +0.00244087799114398267,

  52.       -0.00157693034468678425,

  53.       +0.00102626332050760715,

  54.       -0.000672061631156136204,

  55.       +0.000442473061814620910,

  56.       -0.000292677224729627445,

  57.       +0.000194387276054539318,

  58.       -0.000129574266852748819,

  59.       +0.0000866503580520812717,

  60.       -0.0000581136075044138168

  61.     };

  62.     const double ap = abs(p);

  63.     if (ap < 0.01159)

  64.       return

  65.         -1 +

  66.         p*(1 +

  67.         p*(q[2] +

  68.         p*(q[3] +

  69.         p*(q[4] +

  70.         p*(q[5] +

  71.         p*q[6]

  72.         )))));

  73.     else if (ap < 0.0766)

  74.       return

  75.         -1 +

  76.         p*(1 +

  77.         p*(q[2] +

  78.         p*(q[3] +

  79.         p*(q[4] +

  80.         p*(q[5] +

  81.         p*(q[6] +

  82.         p*(q[7] +

  83.         p*(q[8] +

  84.         p*(q[9] +

  85.         p*q[10]

  86.         )))))))));

  87.     else

  88.       return

  89.         -1 +

  90.         p*(1 +

  91.         p*(q[2] +

  92.         p*(q[3] +

  93.         p*(q[4] +

  94.         p*(q[5] +

  95.         p*(q[6] +

  96.         p*(q[7] +

  97.         p*(q[8] +

  98.         p*(q[9] +

  99.         p*(q[10] +

  100.         p*(q[11] +

  101.         p*(q[12] +

  102.         p*(q[13] +

  103.         p*(q[14] +

  104.         p*(q[15] +

  105.         p*(q[16] +

  106.         p*(q[17] +

  107.         p*(q[18] +

  108.         p*(q[19] +

  109.         p*q[20]

  110.         )))))))))))))))))));

  111.   }

  112. 

  113. 

  114.   inline

  115.   double

  116.   LambertW0ZeroSeries(const double z)

  117.   {

  118.     return

  119.       z*(1 -

  120.       z*(1 -

  121.       z*(1.5 -

  122.       z*(2.6666666666666666667 -

  123.       z*(5.2083333333333333333 -

  124.       z*(10.8 -

  125.       z*(23.343055555555555556 -

  126.       z*(52.012698412698412698 -

  127.       z*(118.62522321428571429 -

  128.       z*(275.57319223985890653 -

  129.       z*(649.78717234347442681 -

  130.       z*(1551.1605194805194805 -

  131.       z*(3741.4497029592385495 -

  132.       z*(9104.5002411580189358 -

  133.       z*(22324.308512706601434 -

  134.       z*(55103.621972903835338 -

  135.       z*136808.86090394293563

  136.       ))))))))))))))));

  137.   }

  138. 

  139. 

  140.   inline

  141.   double

  142.   FinalResult(const double w, const double y)

  143.   {

  144.     const double f0 = w - y;

  145.     const double f1 = 1 + y;

  146.     const double f00 = f0 * f0;

  147.     const double f11 = f1 * f1;

  148.     const double f0y = f0 * y;

  149.     return w - 4 * f0 * (6 * f1 * (f11 + f0y) + f00 * y) /

  150.                         (f11 * (24 * f11 + 36 * f0y) +

  151.                          f00 * (6 * y * y + 8 * f1 * y + f0y));

  152.   }

  153. 

  154. 

  155.   double

  156.   LambertW0(const double z)

  157.   {

  158.     static double e[66];

  159.     static double g[65];

  160.     static double a[12];

  161.     static double b[12];

  162. 

  163.     if (!e[0]) {

  164.       const double e1 = 1 / M_E;

  165.       double ej = 1;

  166.       e[0] = M_E;

  167.       e[1] = 1;

  168.       g[0] = 0;

  169.       for (int j = 1, jj = 2; jj < 66; ++jj) {

  170.         ej *= M_E;

  171.         e[jj] = e[j] * e1;

  172.         g[j] = j * ej;

  173.         j = jj;

  174.       }

  175.       a[0] = sqrt(e1);

  176.       b[0] = 0.5;

  177.       for (int j = 0, jj = 1; jj < 12; ++jj) {

  178.         a[jj] = sqrt(a[j]);

  179.         b[jj] = b[j] * 0.5;

  180.         j = jj;

  181.       }

  182.     }

  183.     if (abs(z) < 0.05)

  184.       return LambertW0ZeroSeries(z);

  185.     if (z < -0.35) {

  186.       const double p2 = 2 * (M_E * z + 1);

  187.       if (p2 > 0)

  188.         return LambertWSeries(sqrt(p2));

  189.       if (p2 == 0)

  190.         return -1;

  191.       cerr << "(lambertw0) Argument out of range. z=" << z << endl;

  192.       return numeric_limits<double>::quiet_NaN();

  193.     }

  194.     int n;

  195.     for (n = 0; n <= 2; ++n)

  196.       if (g[n] > z)

  197.         goto line1;

  198.     n = 2;

  199.     for (int j = 1; j <= 5; ++j) {

  200.       n *= 2;

  201.       if (g[n] > z)

  202.         goto line2;

  203.     }

  204.     cerr << "(lambertw0) Argument too large. z=" << z << endl;

  205.     return numeric_limits<double>::quiet_NaN();

  206.   line2:

  207.     {

  208.       int nh = n / 2;

  209.       for (int j = 1; j <= 5; ++j) {

  210.         nh /= 2;

  211.         if (nh <= 0)

  212.           break;

  213.         if (g[n-nh] > z)

  214.           n -= nh;

  215.       }

  216.     }

  217.   line1:

  218.     --n;

  219.     int jmax = 8;

  220.     if (z <= -0.36)

  221.       jmax = 12;

  222.     else if (z <= -0.3)

  223.       jmax = 11;

  224.     else if (n <= 0)

  225.       jmax = 10;

  226.     else if (n <= 1)

  227.       jmax = 9;

  228.     double y = z * e[n+1];

  229.     double w = n;

  230.     for (int j = 0; j < jmax; ++j) {

  231.       const double wj = w + b[j];

  232.       const double yj = y * a[j];

  233.       if (wj < yj) {

  234.         w = wj;

  235.         y = yj;

  236.       }

  237.     }

  238.     return FinalResult(w, y);

  239.   }

  240. 

  241. 

  242.   double

  243.   LambertWm1(const double z)

  244.   {

  245.     static double e[64];

  246.     static double g[64];

  247.     static double a[12];

  248.     static double b[12];

  249. 

  250.     if (!e[0]) {

  251.       const double e1 = 1 / M_E;

  252.       double ej = e1;

  253.       e[0] = M_E;

  254.       g[0] = -e1;

  255.       for (int j = 0, jj = 1; jj < 64; ++jj) {

  256.         ej *= e1;

  257.         e[jj] = e[j] * M_E;

  258.         g[jj] = -(jj+1) * ej;

  259.         j = jj;

  260.       }

  261.       a[0] = sqrt(M_E);

  262.       b[0] = 0.5;

  263.       for (int j = 0, jj = 1; jj < 12; ++jj) {

  264.         a[jj] = sqrt(a[j]);

  265.         b[jj] = b[j] * 0.5;

  266.         j = jj;

  267.       }

  268.     }

  269.     if (z >= 0) {

  270.       cerr << "(lambertwm1) Argument out of range. z=" << z << endl;

  271.       return numeric_limits<double>::quiet_NaN();

  272.     }

  273.     if (z < -0.35) {

  274.       const double p2 = 2 * (M_E * z + 1);

  275.       if (p2 > 0)

  276.         return LambertWSeries(-sqrt(p2));

  277.       if (p2 == 0)

  278.         return -1;

  279.       cerr << "(lambertwm1) Argument out of range. z=" << z << endl;

  280.       return numeric_limits<double>::quiet_NaN();

  281.     }

  282.     int n = 2;

  283.     if (g[n - 1] > z)

  284.       goto line1;

  285.     for (int j = 1; j <= 5; ++j) {

  286.       n *= 2;

  287.       if (g[n - 1] > z)

  288.         goto line2;

  289.     }

  290.     cerr << "(lambertwm1) Argument too small. z=" << z << endl;

  291.     return numeric_limits<double>::quiet_NaN();

  292.   line2:

  293.     {

  294.       int nh = n / 2;

  295.       for (int j = 1; j <= 5; ++j) {

  296.         nh /= 2;

  297.         if (nh <= 0)

  298.           break;

  299.         if (g[n-nh - 1] > z)

  300.           n -= nh;

  301.       }

  302.     }

  303.   line1:

  304.     --n;

  305.     int jmax = 11;

  306.     if (n >= 8)

  307.       jmax = 8;

  308.     else if (n >= 3)

  309.       jmax = 9;

  310.     else if (n >= 2)

  311.       jmax = 10;

  312.     double w = -n;

  313.     double y = z * e[n - 1];

  314.     for (int j = 0; j < jmax; ++j) {

  315.       const double wj = w - b[j];

  316.       const double yj = y * a[j];

  317.       if (wj < yj) {

  318.         w = wj;

  319.         y = yj;

  320.       }

  321.     }

  322.     return FinalResult(w, y);

  323.   }

  324. 

  325. }