File indexing completed on 2024-04-06 12:26:25
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056 #include "vdt/vdtMath.h"
0057
0058 static const float SN[] = {
0059 -8.39167827910303881427E-11,
0060 4.62591714427012837309E-8,
0061 -9.75759303843632795789E-6,
0062 9.76945438170435310816E-4,
0063 -4.13470316229406538752E-2,
0064 1.00000000000000000302E0,
0065 };
0066 static const float SD[] = {
0067 2.03269266195951942049E-12,
0068 1.27997891179943299903E-9,
0069 4.41827842801218905784E-7,
0070 9.96412122043875552487E-5,
0071 1.42085239326149893930E-2,
0072 9.99999999999999996984E-1,
0073 };
0074
0075 static const float CN[] = {
0076 2.02524002389102268789E-11,
0077 -1.35249504915790756375E-8,
0078 3.59325051419993077021E-6,
0079 -4.74007206873407909465E-4,
0080 2.89159652607555242092E-2,
0081 -1.00000000000000000080E0,
0082 };
0083 static const float CD[] = {
0084 4.07746040061880559506E-12,
0085 3.06780997581887812692E-9,
0086 1.23210355685883423679E-6,
0087 3.17442024775032769882E-4,
0088 5.10028056236446052392E-2,
0089 4.00000000000000000080E0,
0090 };
0091
0092 static const float FN4[] = {
0093 4.23612862892216586994E0,
0094 5.45937717161812843388E0,
0095 1.62083287701538329132E0,
0096 1.67006611831323023771E-1,
0097 6.81020132472518137426E-3,
0098 1.08936580650328664411E-4,
0099 5.48900223421373614008E-7,
0100 };
0101 static const float FD4[] = {
0102
0103 8.16496634205391016773E0,
0104 7.30828822505564552187E0,
0105 1.86792257950184183883E0,
0106 1.78792052963149907262E-1,
0107 7.01710668322789753610E-3,
0108 1.10034357153915731354E-4,
0109 5.48900252756255700982E-7,
0110 };
0111
0112 static const float FN8[] = {
0113 4.55880873470465315206E-1,
0114 7.13715274100146711374E-1,
0115 1.60300158222319456320E-1,
0116 1.16064229408124407915E-2,
0117 3.49556442447859055605E-4,
0118 4.86215430826454749482E-6,
0119 3.20092790091004902806E-8,
0120 9.41779576128512936592E-11,
0121 9.70507110881952024631E-14,
0122 };
0123 static const float FD8[] = {
0124
0125 9.17463611873684053703E-1,
0126 1.78685545332074536321E-1,
0127 1.22253594771971293032E-2,
0128 3.58696481881851580297E-4,
0129 4.92435064317881464393E-6,
0130 3.21956939101046018377E-8,
0131 9.43720590350276732376E-11,
0132 9.70507110881952025725E-14,
0133 };
0134
0135 static const float GN4[] = {
0136 8.71001698973114191777E-2,
0137 6.11379109952219284151E-1,
0138 3.97180296392337498885E-1,
0139 7.48527737628469092119E-2,
0140 5.38868681462177273157E-3,
0141 1.61999794598934024525E-4,
0142 1.97963874140963632189E-6,
0143 7.82579040744090311069E-9,
0144 };
0145 static const float GD4[] = {
0146
0147 1.64402202413355338886E0,
0148 6.66296701268987968381E-1,
0149 9.88771761277688796203E-2,
0150 6.22396345441768420760E-3,
0151 1.73221081474177119497E-4,
0152 2.02659182086343991969E-6,
0153 7.82579218933534490868E-9,
0154 };
0155
0156 static const float GN8[] = {
0157 6.97359953443276214934E-1,
0158 3.30410979305632063225E-1,
0159 3.84878767649974295920E-2,
0160 1.71718239052347903558E-3,
0161 3.48941165502279436777E-5,
0162 3.47131167084116673800E-7,
0163 1.70404452782044526189E-9,
0164 3.85945925430276600453E-12,
0165 3.14040098946363334640E-15,
0166 };
0167 static const float GD8[] = {
0168
0169 1.68548898811011640017E0,
0170 4.87852258695304967486E-1,
0171 4.67913194259625806320E-2,
0172 1.90284426674399523638E-3,
0173 3.68475504442561108162E-5,
0174 3.57043223443740838771E-7,
0175 1.72693748966316146736E-9,
0176 3.87830166023954706752E-12,
0177 3.14040098946363335242E-15,
0178 };
0179
0180 inline float polevlf(float xx, const float *coef, int N) {
0181 float ans, x;
0182 const float *p;
0183 int i;
0184
0185 x = xx;
0186 p = coef;
0187 ans = *p++;
0188
0189 i = N;
0190 do
0191 ans = ans * x + *p++;
0192 while (--i);
0193
0194 return (ans);
0195 }
0196
0197
0198
0199
0200
0201
0202 inline float p1evlf(float xx, const float *coef, int N) {
0203 float ans, x;
0204 const float *p;
0205 int i;
0206
0207 x = xx;
0208 p = coef;
0209 ans = x + *p++;
0210 i = N - 1;
0211
0212 do
0213 ans = ans * x + *p++;
0214 while (--i);
0215
0216 return (ans);
0217 }
0218
0219 inline int sicif(float xx, float &si, float &ci) {
0220 const float MAXNUMF = 1.7014117331926442990585209174225846272e38;
0221 const float PIO2F = 1.5707963267948966192;
0222
0223 const float EUL = 0.57721566490153286061;
0224
0225 float x, z, c, s, f, g;
0226 int sign;
0227
0228 x = xx;
0229 if (x < 0.0f) {
0230 sign = -1;
0231 x = -x;
0232 } else
0233 sign = 0;
0234
0235 if (x == 0.0f) {
0236 si = 0.0;
0237 ci = -MAXNUMF;
0238 return (0);
0239 }
0240
0241 if (x > 1.0e9f) {
0242 float su, cu;
0243 vdt::fast_sincosf(x, su, cu);
0244 si = PIO2F - cu / x;
0245 ci = su / x;
0246 return (0);
0247 }
0248
0249 if (x > 4.0f)
0250 goto asympt;
0251
0252 z = x * x;
0253 s = x * polevlf(z, SN, 5) / polevlf(z, SD, 5);
0254 c = z * polevlf(z, CN, 5) / polevlf(z, CD, 5);
0255
0256 if (sign)
0257 s = -s;
0258 si = s;
0259 ci = EUL + vdt::fast_logf(x) + c;
0260 return (0);
0261
0262
0263
0264
0265
0266
0267
0268
0269
0270
0271
0272
0273
0274
0275
0276 asympt:
0277 vdt::fast_sincosf(x, s, c);
0278 z = 1.0f / (x * x);
0279 if (x < 8.0f) {
0280 f = polevlf(z, FN4, 6) / (x * p1evlf(z, FD4, 7));
0281 g = z * polevlf(z, GN4, 7) / p1evlf(z, GD4, 7);
0282 } else {
0283 f = polevlf(z, FN8, 8) / (x * p1evlf(z, FD8, 8));
0284 g = z * polevlf(z, GN8, 8) / p1evlf(z, GD8, 9);
0285 }
0286 si = PIO2F - f * c - g * s;
0287 if (sign)
0288 si = -(si);
0289 ci = f * s - g * c;
0290
0291 return (0);
0292 }