Advertisement
Guest User

Untitled

a guest
Feb 6th, 2020
274
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C 10.53 KB
  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
Advertisement
Advertisement
RAW Paste Data Copied
Advertisement