libdspl-2.0
Digital Signal Processing Algorithm Library
math.c
1 /*
2 * Copyright (c) 2015-2020 Sergey Bakhurin
3 * Digital Signal Processing Library [http://dsplib.org]
4 *
5 * This file is part of DSPL.
6 *
7 * is free software: you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation, either version 3 of the License, or
10 * (at your option) any later version.
11 *
12 * DSPL is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with Foobar. If not, see <http://www.gnu.org/licenses/>.
19 */
20 
21 
22 #include <stdio.h>
23 #include <stdlib.h>
24 #include "dspl.h"
25 
26 
27 
28 #if DOXYGEN_ENGLISH
29 
80 #endif
81 #if DOXYGEN_RUSSIAN
82 
140 #endif
141 int DSPL_API acos_cmplx(complex_t* x, int n, complex_t *y)
142 {
143  int k, res;
144  double pi2 = 0.5 * M_PI;
145 
146  res = asin_cmplx(x, n, y);
147  if(res != RES_OK)
148  return res;
149 
150  for(k = 0; k < n; k++)
151  {
152  RE(y[k]) = pi2 - RE(y[k]);
153  IM(y[k]) = - IM(y[k]);
154  }
155  return RES_OK;
156 }
157 
158 
159 
160 
161 #ifdef DOXYGEN_ENGLISH
162 
213 #endif
214 #ifdef DOXYGEN_RUSSIAN
215 
270 #endif
271 int DSPL_API asin_cmplx(complex_t* x, int n, complex_t *y)
272 {
273  int k;
274  complex_t tmp;
275  if(!x || !y)
276  return ERROR_PTR;
277  if(n < 1)
278  return ERROR_SIZE;
279 
280  for(k = 0; k < n; k++)
281  {
282  RE(tmp) = 1.0 - CMRE(x[k], x[k]); /* 1-x[k]^2 */
283  IM(tmp) = - CMIM(x[k], x[k]); /* 1-x[k]^2 */
284  sqrt_cmplx(&tmp, 1, y+k); /* sqrt(1 - x[k]^2) */
285  RE(y[k]) -= IM(x[k]); /* j * x[k] + sqrt(1 - x[k]^2) */
286  IM(y[k]) += RE(x[k]); /* j * x[k] + sqrt(1 - x[k]^2) */
287  log_cmplx(y+k, 1, &tmp); /* log( j * x[k] + sqrt(1 - x[k]^2) ) */
288  RE(y[k]) = IM(tmp); /* -j * log( j * x[k] + sqrt(1 - x[k]^2) ) */
289  IM(y[k]) = -RE(tmp); /* -j * log( j * x[k] + sqrt(1 - x[k]^2) ) */
290  }
291  return RES_OK;
292 }
293 
294 
295 
296 
297 
298 #ifdef DOXYGEN_ENGLISH
299 
338 #endif
339 #ifdef DOXYGEN_RUSSIAN
340 
385 #endif
386 int DSPL_API bessel_i0(double* x, int n, double* y)
387 {
388  double P16[17] = { 1.0000000000000000000000801e+00,
389  2.4999999999999999999629693e-01,
390  2.7777777777777777805664954e-02,
391  1.7361111111111110294015271e-03,
392  6.9444444444444568581891535e-05,
393  1.9290123456788994104574754e-06,
394  3.9367598891475388547279760e-08,
395  6.1511873265092916275099070e-10,
396  7.5940584360755226536109511e-12,
397  7.5940582595094190098755663e-14,
398  6.2760839879536225394314453e-16,
399  4.3583591008893599099577755e-18,
400  2.5791926805873898803749321e-20,
401  1.3141332422663039834197910e-22,
402  5.9203280572170548134753422e-25,
403  2.0732014503197852176921968e-27,
404  1.1497640034400735733456400e-29};
405 
406  double P22[23] = { 3.9894228040143265335649948e-01,
407  4.9867785050353992900698488e-02,
408  2.8050628884163787533196746e-02,
409  2.9219501690198775910219311e-02,
410  4.4718622769244715693031735e-02,
411  9.4085204199017869159183831e-02,
412  -1.0699095472110916094973951e-01,
413  2.2725199603010833194037016e+01,
414  -1.0026890180180668595066918e+03,
415  3.1275740782277570164423916e+04,
416  -5.9355022509673600842060002e+05,
417  2.6092888649549172879282592e+06,
418  2.3518420447411254516178388e+08,
419  -8.9270060370015930749184222e+09,
420  1.8592340458074104721496236e+11,
421  -2.6632742974569782078420204e+12,
422  2.7752144774934763122129261e+13,
423  -2.1323049786724612220362154e+14,
424  1.1989242681178569338129044e+15,
425  -4.8049082153027457378879746e+15,
426  1.3012646806421079076251950e+16,
427  -2.1363029690365351606041265e+16,
428  1.6069467093441596329340754e+16};
429 
430  double x2;
431  int k;
432 
433  if(!x || !y)
434  return ERROR_PTR;
435  if(n < 1)
436  return ERROR_SIZE;
437 
438  for(k =0; k < n; k++)
439  {
440  if(x[k] < 0.0)
441  return ERROR_NEGATIVE;
442 
443  if(x[k] < 7.75)
444  {
445  x2 = x[k] * x[k] * 0.25;
446  polyval(P16, 16, &x2, 1, y+k);
447  y[k] = x2 * y[k] + 1.0;
448  }
449  else
450  {
451  x2 = 1.0 / x[k];
452  polyval(P22, 22, &x2, 1, y+k);
453  y[k] *= exp(x[k]) / sqrt(x[k]);
454  }
455  }
456  return RES_OK;
457 }
458 
459 
460 
461 #ifdef DOXYGEN_ENGLISH
462 
463 #endif
464 #ifdef DOXYGEN_RUSSIAN
465 
466 #endif
467 double DSPL_API dmod (double x, double y)
468 {
469  if(y == 0.0)
470  return x;
471  return x - floor(x/y) * y;
472 }
473 
474 
475 
476 
477 #ifdef DOXYGEN_ENGLISH
478 
529 #endif
530 #ifdef DOXYGEN_RUSSIAN
531 
584 #endif
585 int DSPL_API cos_cmplx(complex_t* x, int n, complex_t *y)
586 {
587  int k;
588  double ep, em, sx, cx;
589  if(!x || !y)
590  return ERROR_PTR;
591  if(n < 1)
592  return ERROR_SIZE;
593 
594  for(k = 0; k < n; k++)
595  {
596  ep = exp( IM(x[k]));
597  em = exp(-IM(x[k]));
598  sx = 0.5 * sin(RE(x[k]));
599  cx = 0.5 * cos(RE(x[k]));
600  RE(y[k]) = cx * (em + ep);
601  IM(y[k]) = sx * (em - ep);
602  }
603  return RES_OK;
604 }
605 
606 
607 
608 #ifdef DOXYGEN_ENGLISH
609 
661 #endif
662 #ifdef DOXYGEN_RUSSIAN
663 
717 #endif
718 int DSPL_API log_cmplx(complex_t* x, int n, complex_t *y)
719 {
720  int k;
721  if(!x || !y)
722  return ERROR_PTR;
723  if(n < 1)
724  return ERROR_SIZE;
725 
726  for(k = 0; k < n; k++)
727  {
728  RE(y[k]) = 0.5 * log(ABSSQR(x[k]));
729  IM(y[k]) = atan2(IM(x[k]), RE(x[k]));
730  }
731  return RES_OK;
732 }
733 
734 
735 
736 
737 
738 
739 
740 #ifdef DOXYGEN_ENGLISH
741 
793 #endif
794 #ifdef DOXYGEN_RUSSIAN
795 
848 #endif
849 int DSPL_API sin_cmplx(complex_t* x, int n, complex_t *y)
850 {
851  int k;
852  double ep, em, sx, cx;
853  if(!x || !y)
854  return ERROR_PTR;
855  if(n < 1)
856  return ERROR_SIZE;
857 
858  for(k = 0; k < n; k++)
859  {
860  ep = exp( IM(x[k]));
861  em = exp(-IM(x[k]));
862  sx = 0.5 * sin(RE(x[k]));
863  cx = 0.5 * cos(RE(x[k]));
864  RE(y[k]) = sx * (em + ep);
865  IM(y[k]) = cx * (ep - em);
866  }
867  return RES_OK;
868 }
869 
870 
871 
872 
873 
874 #ifdef DOXYGEN_ENGLISH
875 
902 #endif
903 #ifdef DOXYGEN_RUSSIAN
904 
935 #endif
936 int DSPL_API sinc(double* x, int n, double a, double* y)
937 {
938  int k;
939 
940  if(!x || !y)
941  return ERROR_PTR;
942  if(n<1)
943  return ERROR_SIZE;
944 
945  for(k = 0; k < n; k++)
946  y[k] = (x[k]==0.0) ? 1.0 : sin(a*x[k])/(a*x[k]);
947 
948  return RES_OK;
949 
950 
951 }
952 
953 
954 
955 #ifdef DOXYGEN_ENGLISH
956 
999 #endif
1000 #ifdef DOXYGEN_RUSSIAN
1001 
1042 #endif
1043 int DSPL_API sine_int(double* x, int n, double* si)
1044 {
1045  int k, sgn, p;
1046  double num, den, y, x2, x22, z, f, g;
1047 
1048  double A[8] = {+1.00000000000000000E0,
1049  -4.54393409816329991E-2,
1050  +1.15457225751016682E-3,
1051  -1.41018536821330254E-5,
1052  +9.43280809438713025E-8,
1053  -3.53201978997168357E-10,
1054  +7.08240282274875911E-13,
1055  -6.05338212010422477E-16};
1056 
1057 
1058 
1059  double B[7] = {+1.0,
1060  +1.01162145739225565E-2,
1061  +4.99175116169755106E-5,
1062  +1.55654986308745614E-7,
1063  +3.28067571055789734E-10,
1064  +4.50490975753865810E-13,
1065  +3.21107051193712168E-16};
1066 
1067 
1068 
1069  double FA[11] = {+1.000000000000000000000E0,
1070  +7.444370681619367006180E2,
1071  +1.963963728951468698010E5,
1072  +2.377503101254318340340E7,
1073  +1.430734038212746368880E9,
1074  +4.33736238870432522765E10,
1075  +6.40533830574022022911E11,
1076  +4.20968180571076940208E12,
1077  +1.00795182980368574617E13,
1078  +4.94816688199951963482E12,
1079  -4.94701168645415959931E11};
1080 
1081  double FB[10] = {+1.000000000000000000000E0,
1082  +7.464370681619276780310E2,
1083  +1.978652470315839514500E5,
1084  +2.415356701651268451440E7,
1085  +1.474789521929854649580E9,
1086  +4.58595115847765779830E10,
1087  +7.08501308149515401563E11,
1088  +5.06084464593475076774E12,
1089  +1.43468549171581016479E13,
1090  +1.11535493509914254097E13};
1091 
1092 
1093 
1094  double GA[11] = {+1.000000000000000000E0,
1095  +8.135952011516861500E2,
1096  +2.352391816264782000E5,
1097  +3.125575707957787310E7,
1098  +2.062975951467633540E9,
1099  +6.83052205423625007E10,
1100  +1.09049528450362786E12,
1101  +7.57664583257834349E12,
1102  +1.81004487464664575E13,
1103  +6.43291613143049485E12,
1104  -1.36517137670871689E12};
1105 
1106 
1107  double GB[10] = {+1.000000000000000000E0,
1108  +8.195952011514515640E2,
1109  +2.400367528355787770E5,
1110  +3.260266616470908220E7,
1111  +2.233555432780993600E9,
1112  +7.87465017341829930E10,
1113  +1.39866710696414565E12,
1114  +1.17164723371736605E13,
1115  +4.01839087307656620E13,
1116  +3.99653257887490811E13};
1117 
1118  if(!x || !si)
1119  return ERROR_PTR;
1120  if(n<1)
1121  return ERROR_SIZE;
1122 
1123 
1124  for(p = 0; p < n; p++)
1125  {
1126  sgn = x[p] > 0.0 ? 0 : 1;
1127  y = x[p] < 0.0 ? -x[p] : x[p];
1128 
1129  if(y < 4)
1130  {
1131  x2 = y * y;
1132  z = 1.0;
1133  num = 0.0;
1134  for(k = 0; k < 8; k++)
1135  {
1136  num += A[k] * z;
1137  z*=x2;
1138  }
1139  z = 1.0;
1140  den = 0.0;
1141  for(k = 0; k < 7; k++)
1142  {
1143  den += B[k]*z;
1144  z*=x2;
1145  }
1146  si[p] = x[p] * num/den;
1147  }
1148  else
1149  {
1150 
1151  x2 = 1.0/y;
1152  x22 = x2*x2;
1153  z = 1.0;
1154  num = 0.0;
1155  for(k = 0; k < 11; k++)
1156  {
1157  num += FA[k] * z;
1158  z*=x22;
1159  }
1160  z = 1.0;
1161  den = 0.0;
1162  for(k = 0; k < 10; k++)
1163  {
1164  den += FB[k]*z;
1165  z*=x22;
1166  }
1167 
1168  f = x2 * num / den;
1169 
1170  z = 1.0;
1171  num = 0.0;
1172  for(k = 0; k < 11; k++)
1173  {
1174  num += GA[k] * z;
1175  z*=x22;
1176  }
1177  z = 1.0;
1178  den = 0.0;
1179  for(k = 0; k < 10; k++)
1180  {
1181  den += GB[k]*z;
1182  z*=x22;
1183  }
1184 
1185  g = x22 * num / den;
1186 
1187  si[p] = sgn ? f * cos(y) + g * sin(y) - M_PI * 0.5 :
1188  M_PI * 0.5 - f * cos(y) - g * sin(y);
1189  }
1190  }
1191  return RES_OK;
1192 }
1193 
1194 
1195 
1196 
1197 
1198 
1199 #ifdef DOXYGEN_ENGLISH
1200 
1251 #endif
1252 #ifdef DOXYGEN_RUSSIAN
1253 
1306 #endif
1307 int DSPL_API sqrt_cmplx(complex_t* x, int n, complex_t *y)
1308 {
1309  int k;
1310  double r, zr, at;
1311  complex_t t;
1312  if(!x || !y)
1313  return ERROR_PTR;
1314  if(n < 1)
1315  return ERROR_SIZE;
1316 
1317  for(k = 0; k < n; k++)
1318  {
1319  r = ABS(x[k]);
1320  if(r == 0.0)
1321  {
1322  RE(y[k]) = 0.0;
1323  IM(y[k]) = 0.0;
1324  }
1325  else
1326  {
1327  RE(t) = RE(x[k]) + r;
1328  IM(t) = IM(x[k]);
1329  at = ABS(t);
1330  if(at == 0.0)
1331  {
1332  RE(y[k]) = 0.0;
1333  IM(y[k]) = sqrt(r);
1334  }
1335  else
1336  {
1337  zr = 1.0 / ABS(t);
1338  r = sqrt(r);
1339  RE(y[k]) = RE(t) * zr * r;
1340  IM(y[k]) = IM(t) * zr * r;
1341  }
1342  }
1343  }
1344  return RES_OK;
1345 }
int DSPL_API sine_int(double *x, int n, double *si)
Sine integral function for the real vector x.
Definition: math.c:1043
#define RE(x)
Macro sets real part of the complex number.
Definition: dspl.h:359
int DSPL_API log_cmplx(complex_t *x, int n, complex_t *y)
The logarithm function the complex vector argument x.
Definition: math.c:718
int DSPL_API acos_cmplx(complex_t *x, int n, complex_t *y)
The inverse of the cosine function the complex vector argument x.
Definition: math.c:141
#define ERROR_PTR
Pointer error. This error means that one of the required pointers (memory to be allocated for) is tra...
Definition: dspl.h:545
#define ERROR_SIZE
Error array size. This error occurs when in addition to the pointer the wrong input is passed to the ...
Definition: dspl.h:553
double complex_t[2]
Complex data type.
Definition: dspl.h:86
#define RES_OK
The function completed correctly. No errors.
Definition: dspl.h:497
#define ABSSQR(x)
The macro returns the square of the modulus of a complex number x.
Definition: dspl.h:473
int DSPL_API sqrt_cmplx(complex_t *x, int n, complex_t *y)
Square root of the complex vector argguument x.
Definition: math.c:1307
int DSPL_API sin_cmplx(complex_t *x, int n, complex_t *y)
The sine function the complex vector argument x.
Definition: math.c:849
int DSPL_API bessel_i0(double *x, int n, double *y)
Modified Bessel Function of the First Kind – [1].
Definition: math.c:386
int DSPL_API asin_cmplx(complex_t *x, int n, complex_t *y)
The inverse of the sine function the complex vector argument x.
Definition: math.c:271
int DSPL_API cos_cmplx(complex_t *x, int n, complex_t *y)
The cosine function the complex vector argument x.
Definition: math.c:585
int DSPL_API sinc(double *x, int n, double a, double *y)
Function for the real vector x.
Definition: math.c:936
#define IM(x)
Macro sets imaginary part of the complex number.
Definition: dspl.h:417
#define ERROR_NEGATIVE
Negative parameter. The function returns the given error code when it takes a negative parameter....
Definition: dspl.h:540