API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Feb 6th, 2020
422
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C 10.53 KB | None | 0 0
raw download clone embed print report
  1. //
  2. //    FILE: __functions.h
  3. //
  4. //    fast and approximate transcendental functions
  5. //
  6. //    copyright (c) 2004  Robert Bristow-Johnson
  7. //
  8. //    rbj@audioimagination.com
  9. //
  10.  
  11.  
  12. #ifndef __FUNCTIONS_H
  13. #define __FUNCTIONS_H
  14.  
  15. #define TINY 1.0e-8
  16. #define HUGE 1.0e8
  17.  
  18. #define PI              (3.1415926535897932384626433832795028841972)        /* pi */
  19. #define ONE_OVER_PI     (0.3183098861837906661338147750939)
  20. #define TWOPI           (6.2831853071795864769252867665590057683943)        /* 2*pi */
  21. #define ONE_OVER_TWOPI  (0.15915494309189535682609381638)
  22. #define PI_2            (1.5707963267948966192313216916397514420986)        /* pi/2 */
  23. #define TWO_OVER_PI     (0.636619772367581332267629550188)
  24. #define LN2             (0.6931471805599453094172321214581765680755)        /* ln(2) */
  25. #define ONE_OVER_LN2    (1.44269504088896333066907387547)
  26. #define LN10            (2.3025850929940456840179914546843642076011)        /* ln(10) */
  27. #define ONE_OVER_LN10   (0.43429448190325177635683940025)
  28. #define ROOT2           (1.4142135623730950488016887242096980785697)        /* sqrt(2) */
  29. #define ONE_OVER_ROOT2  (0.707106781186547438494264988549)
  30.  
  31. #define DB_LOG2_ENERGY          (3.01029995663981154631945610163)           /* dB = DB_LOG2_ENERGY*__log2(energy) */
  32. #define DB_LOG2_AMPL            (6.02059991327962309263891220326)           /* dB = DB_LOG2_AMPL*__log2(amplitude) */
  33. #define ONE_OVER_DB_LOG2_AMPL   (0.16609640474436811218256075335)           /* amplitude = __exp2(ONE_OVER_DB_LOG2_AMPL*dB) */
  34.  
  35. #define LONG_OFFSET     4096L
  36. #define FLOAT_OFFSET    4096.0
  37.  
  38.  
  39.  
  40. float   __sqrt(float x);
  41.  
  42. float   __log2(float x);
  43. float   __exp2(float x);
  44.  
  45. float   __log(float x);
  46. float   __exp(float x);
  47.  
  48. float   __pow(float x, float y);
  49.  
  50. float   __sin_pi(float x);
  51. float   __cos_pi(float x);
  52.  
  53. float   __sin(float x);
  54. float   __cos(float x);
  55. float   __tan(float x);
  56.  
  57. float   __atan(float x);
  58. float   __asin(float x);
  59. float   __acos(float x);
  60.  
  61. float   __arg(float Imag, float Real);
  62.  
  63. float   __poly(float *a, int order, float x);
  64. float   __map(float *f, float scaler, float x);
  65. float   __discreteMap(float *f, float scaler, float x);
  66.  
  67. unsigned long __random();
  68.  
  69. #endif
  70.  
  71.  
  72.  
  73.  
  74. //
  75. //    FILE: __functions.c
  76. //
  77. //    fast and approximate transcendental functions
  78. //
  79. //    copyright (c) 2004  Robert Bristow-Johnson
  80. //
  81. //    rbj@audioimagination.com
  82. //
  83.  
  84. #define STD_MATH_LIB 0
  85.  
  86. #include "__functions.h"
  87.  
  88. #if STD_MATH_LIB
  89. #include "math.h"   // angle brackets don't work with SE markup
  90. #endif
  91.  
  92.  
  93.  
  94.  
  95. float   __sqrt(register float x)
  96. {
  97. #if STD_MATH_LIB
  98.     return (float)sqrt((double)x);
  99. #else
  100.     if (x > 5.877471754e-39)
  101.     {
  102.         register float accumulator, xPower;
  103.         register long intPart;
  104.         register union { float f; long i; } xBits;
  105.  
  106.         xBits.f = x;
  107.  
  108.         intPart = ((xBits.i) >> 23);                  /* get biased exponent */
  109.         intPart -= 127;                             /* unbias it */
  110.  
  111.         x = (float)(xBits.i & 0x007FFFFF);          /* mask off exponent leaving 0x800000*(mantissa - 1) */
  112.         x *= 1.192092895507812e-07;                 /* divide by 0x800000 */
  113.  
  114.         accumulator = 1.0 + 0.49959804148061*x;
  115.         xPower = x * x;
  116.         accumulator += -0.12047308243453*xPower;
  117.         xPower *= x;
  118.         accumulator += 0.04585425015501*xPower;
  119.         xPower *= x;
  120.         accumulator += -0.01076564682800*xPower;
  121.  
  122.         if (intPart & 0x00000001)
  123.         {
  124.             accumulator *= ROOT2;                   /* an odd input exponent means an extra sqrt(2) in the output */
  125.         }
  126.  
  127.         xBits.i = intPart >> 1;                     /* divide exponent by 2, lose LSB */
  128.         xBits.i += 127;                             /* rebias exponent */
  129.         xBits.i <<= 23;                             /* move biased exponent into exponent bits */
  130.  
  131.         return accumulator * xBits.f;
  132.     }
  133.     else
  134.     {
  135.         return 0.0;
  136.     }
  137. #endif
  138. }
  139.  
  140.  
  141.  
  142.  
  143. float   __log2(register float x)
  144. {
  145. #if STD_MATH_LIB
  146.     return (float)(ONE_OVER_LN2*log((double)x));
  147. #else
  148.     if (x > 5.877471754e-39)
  149.     {
  150.         register float accumulator, xPower;
  151.         register long intPart;
  152.  
  153.         register union { float f; long i; } xBits;
  154.  
  155.         xBits.f = x;
  156.  
  157.         intPart = ((xBits.i) >> 23);                  /* get biased exponent */
  158.         intPart -= 127;                             /* unbias it */
  159.  
  160.         x = (float)(xBits.i & 0x007FFFFF);          /* mask off exponent leaving 0x800000*(mantissa - 1) */
  161.         x *= 1.192092895507812e-07;                 /* divide by 0x800000 */
  162.  
  163.         accumulator = 1.44254494359510*x;
  164.         xPower = x * x;
  165.         accumulator += -0.71814525675041*xPower;
  166.         xPower *= x;
  167.         accumulator += 0.45754919692582*xPower;
  168.         xPower *= x;
  169.         accumulator += -0.27790534462866*xPower;
  170.         xPower *= x;
  171.         accumulator += 0.12179791068782*xPower;
  172.         xPower *= x;
  173.         accumulator += -0.02584144982967*xPower;
  174.  
  175.         return accumulator + (float)intPart;
  176.     }
  177.     else
  178.     {
  179.         return -HUGE;
  180.     }
  181. #endif
  182. }
  183.  
  184.  
  185. float   __exp2(register float x)
  186. {
  187. #if STD_MATH_LIB
  188.     return (float)exp(LN2*(double)x);
  189. #else
  190.     if (x >= -127.0)
  191.     {
  192.         register float accumulator, xPower;
  193.         register union { float f; long i; } xBits;
  194.  
  195.         xBits.i = (long)(x + FLOAT_OFFSET) - LONG_OFFSET;       /* integer part */
  196.         x -= (float)(xBits.i);                                  /* fractional part */
  197.  
  198.         accumulator = 1.0 + 0.69303212081966*x;
  199.         xPower = x * x;
  200.         accumulator += 0.24137976293709*xPower;
  201.         xPower *= x;
  202.         accumulator += 0.05203236900844*xPower;
  203.         xPower *= x;
  204.         accumulator += 0.01355574723481*xPower;
  205.  
  206.         xBits.i += 127;                                         /* bias integer part */
  207.         xBits.i <<= 23;                                         /* move biased int part into exponent bits */
  208.  
  209.         return accumulator * xBits.f;
  210.     }
  211.     else
  212.     {
  213.         return 0.0;
  214.     }
  215. #endif
  216. }
  217.  
  218.  
  219. float   __log(register float x)
  220. {
  221. #if STD_MATH_LIB
  222.     return (float)log((double)x);
  223. #else
  224.     return LN2 * __log2(x);
  225. #endif
  226. }
  227.  
  228. float   __exp(register float x)
  229. {
  230. #if STD_MATH_LIB
  231.     return (float)exp((double)x);
  232. #else
  233.     return __exp2(ONE_OVER_LN2*x);
  234. #endif
  235. }
  236.  
  237. float   __pow(float x, float y)
  238. {
  239. #if STD_MATH_LIB
  240.     return (float)pow((double)x, (double)y);
  241. #else
  242.     return __exp2(y*__log2(x));
  243. #endif
  244. }
  245.  
  246.  
  247.  
  248.  
  249. float   __sin_pi(register float x)
  250. {
  251. #if STD_MATH_LIB
  252.     return (float)sin(PI*(double)x);
  253. #else
  254.     register float accumulator, xPower, xSquared;
  255.  
  256.     register long evenIntPart = ((long)(0.5*x + 1024.5) - 1024) << 1;
  257.     x -= (float)evenIntPart;
  258.  
  259.     xSquared = x * x;
  260.     accumulator = 3.14159265358979*x;
  261.     xPower = xSquared * x;
  262.     accumulator += -5.16731953364340*xPower;
  263.     xPower *= xSquared;
  264.     accumulator += 2.54620566822659*xPower;
  265.     xPower *= xSquared;
  266.     accumulator += -0.586027023087261*xPower;
  267.     xPower *= xSquared;
  268.     accumulator += 0.06554823491427*xPower;
  269.  
  270.     return accumulator;
  271. #endif
  272. }
  273.  
  274.  
  275. float   __cos_pi(register float x)
  276. {
  277. #if STD_MATH_LIB
  278.     return (float)cos(PI*(double)x);
  279. #else
  280.     register float accumulator, xPower, xSquared;
  281.  
  282.     register long evenIntPart = ((long)(0.5*x + 1024.5) - 1024) << 1;
  283.     x -= (float)evenIntPart;
  284.  
  285.     xSquared = x * x;
  286.     accumulator = 1.57079632679490*x;                       /* series for sin(PI/2*x) */
  287.     xPower = xSquared * x;
  288.     accumulator += -0.64596406188166*xPower;
  289.     xPower *= xSquared;
  290.     accumulator += 0.07969158490912*xPower;
  291.     xPower *= xSquared;
  292.     accumulator += -0.00467687997706*xPower;
  293.     xPower *= xSquared;
  294.     accumulator += 0.00015303015470*xPower;
  295.  
  296.     return 1.0 - 2.0*accumulator*accumulator;               /* cos(w) = 1 - 2*(sin(w/2))^2 */
  297. #endif
  298. }
  299.  
  300.  
  301. float   __sin(register float x)
  302. {
  303. #if STD_MATH_LIB
  304.     return (float)sin((double)x);
  305. #else
  306.     x *= ONE_OVER_PI;
  307.     return __sin_pi(x);
  308. #endif
  309. }
  310.  
  311. float   __cos(register float x)
  312. {
  313. #if STD_MATH_LIB
  314.     return (float)cos((double)x);
  315. #else
  316.     x *= ONE_OVER_PI;
  317.     return __cos_pi(x);
  318. #endif
  319. }
  320.  
  321. float   __tan(register float x)
  322. {
  323. #if STD_MATH_LIB
  324.     return (float)tan((double)x);
  325. #else
  326.     x *= ONE_OVER_PI;
  327.     return __sin_pi(x) / __cos_pi(x);
  328. #endif
  329. }
  330.  
  331.  
  332.  
  333.  
  334. float   __atan(register float x)
  335. {
  336. #if STD_MATH_LIB
  337.     return (float)atan((double)x);
  338. #else
  339.     register float accumulator, xPower, xSquared, offset;
  340.  
  341.     offset = 0.0;
  342.  
  343.     if (x < -1.0)
  344.     {
  345.         offset = -PI_2;
  346.         x = -1.0 / x;
  347.     }
  348.     else if (x > 1.0)
  349.     {
  350.         offset = PI_2;
  351.         x = -1.0 / x;
  352.     }
  353.     xSquared = x * x;
  354.     accumulator = 1.0;
  355.     xPower = xSquared;
  356.     accumulator += 0.33288950512027*xPower;
  357.     xPower *= xSquared;
  358.     accumulator += -0.08467922817644*xPower;
  359.     xPower *= xSquared;
  360.     accumulator += 0.03252232640125*xPower;
  361.     xPower *= xSquared;
  362.     accumulator += -0.00749305860992*xPower;
  363.  
  364.     return offset + x / accumulator;
  365. #endif
  366. }
  367.  
  368.  
  369. float   __asin(register float x)
  370. {
  371. #if STD_MATH_LIB
  372.     return (float)asin((double)x);
  373. #else
  374.     return __atan(x / __sqrt(1.0 - x * x));
  375. #endif
  376. }
  377.  
  378. float   __acos(register float x)
  379. {
  380. #if STD_MATH_LIB
  381.     return (float)acos((double)x);
  382. #else
  383.     return __atan(__sqrt(1.0 - x * x) / x);
  384. #endif
  385. }
  386.  
  387.  
  388. float   __arg(float Imag, float Real)
  389. {
  390. #if STD_MATH_LIB
  391.     return (float)atan2((double)Imag, (double)Real);
  392. #else
  393.     register float accumulator, xPower, xSquared, offset, x;
  394.  
  395.     if (Imag > 0.0)
  396.     {
  397.         if (Imag <= -Real)
  398.         {
  399.             offset = PI;
  400.             x = Imag / Real;
  401.         }
  402.         else if (Imag > Real)
  403.         {
  404.             offset = PI_2;
  405.             x = -Real / Imag;
  406.         }
  407.         else
  408.         {
  409.             offset = 0.0;
  410.             x = Imag / Real;
  411.         }
  412.     }
  413.     else
  414.     {
  415.         if (Imag >= Real)
  416.         {
  417.             offset = -PI;
  418.             x = Imag / Real;
  419.         }
  420.         else if (Imag < -Real)
  421.         {
  422.             offset = -PI_2;
  423.             x = -Real / Imag;
  424.         }
  425.         else
  426.         {
  427.             offset = 0.0;
  428.             x = Imag / Real;
  429.         }
  430.     }
  431.  
  432.     xSquared = x * x;
  433.     accumulator = 1.0;
  434.     xPower = xSquared;
  435.     accumulator += 0.33288950512027*xPower;
  436.     xPower *= xSquared;
  437.     accumulator += -0.08467922817644*xPower;
  438.     xPower *= xSquared;
  439.     accumulator += 0.03252232640125*xPower;
  440.     xPower *= xSquared;
  441.     accumulator += -0.00749305860992*xPower;
  442.  
  443.     return offset + x / accumulator;
  444. #endif
  445. }
  446.  
  447.  
  448.  
  449.  
  450. float   __poly(float *a, int order, float x)
  451. {
  452.     register float accumulator = 0.0, xPower;
  453.     register int n;
  454.  
  455.     accumulator = a[0];
  456.     xPower = x;
  457.     for (n = 1; n <= order; n++)
  458.     {
  459.         accumulator += a[n] * xPower;
  460.         xPower *= x;
  461.     }
  462.  
  463.     return accumulator;
  464. }
  465.  
  466.  
  467. float   __map(float *f, float scaler, float x)
  468. {
  469.     register long i;
  470.  
  471.     x *= scaler;
  472.  
  473.     i = (long)(x + FLOAT_OFFSET) - LONG_OFFSET;         /* round down without floor() */
  474.  
  475.     return f[i] + (f[i + 1] - f[i])*(x - (float)i);       /* linear interpolate between points */
  476. }
  477.  
  478.  
  479. float   __discreteMap(float *f, float scaler, float x)
  480. {
  481.     register long i;
  482.  
  483.     x *= scaler;
  484.  
  485.     i = (long)(x + (FLOAT_OFFSET + 0.5)) - LONG_OFFSET;   /* round to nearest */
  486.  
  487.     return f[i];
  488. }
  489.  
  490.  
  491. unsigned long __random()
  492. {
  493.     static unsigned long seed0 = 0x5B7A2775, seed1 = 0x80C7169F;
  494.  
  495.     seed0 += seed1;
  496.     seed1 += seed0;
  497.  
  498.     return seed1;
  499. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!