wyj/sin
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Compare commits

base: wyj/sin:c8874275c3c6c26445889eb6e28d22c24eac29ef
Branches Tags
wyj/sin:master
...
compare: wyj/sin:master
Branches Tags
wyj/sin:master
6 Commits
c8874275c3 ... master
Author SHA1 Message Date
wyj df6d8f3689 update: update readme 2024-06-30 20:27:20 -04:00
wyj d077ea3462 update: add readme 2024-06-30 20:20:24 -04:00
wyj 76dc412053 update: use copysign() instead of d=sign(alpha) 2024-06-30 20:13:30 -04:00
wyj df5933ed5c update: fix loop to be 55, it's faster with -O2 2024-06-30 01:09:38 -04:00
wyj d03ad61916 update 2024-06-29 15:49:52 -04:00
wyj ad173aa8ae update: format 2024-06-29 15:35:48 -04:00
2 changed files with 175 additions and 93 deletions
Expand all files Collapse all files
Readme.md
+21
View File
@@ -0,0 +1,21 @@
对比 https://www.bilibili.com/video/BV1Ys3TeFEgi/ 中提到的迭代计算 $\sin$ 函数的算法和使用泰勒展开，意外地发现还是泰勒展开赢麻了……
## 当前结果
在 i5 10400ES 上，通过
```bash
gcc -march=native -O2 -o sin sin.c
```
得到结果为（用时是指循环 10000000 次总耗时，否则都太快了计时误差太大）
```bash
./sin
sin11(0.0001) = 9.999999983330736e-05, it = 55, dsin=-5.551115095370206e-17, took 1.013875 s.
sin2(0.0001) = 9.999999983333334e-05, it = 2, dsin=8.333333333333334e-23, took 0.064005 s.
sin11(0.5) = 0.4794255386042029, it = 55, dsin=-4.871561831101115e-17, took 0.977128 s.
sin2(0.5) = 0.479425538604203, it = 7, dsin=-2.333729166204778e-17, took 0.284716 s.
sin11(1) = 0.8414709848078967, it = 55, dsin=2.999280301164362e-17, took 0.981563 s.
sin2(1) = 0.8414709848078965, it = 9, dsin=-8.220635246624331e-18, took 0.421509 s.
sin11(1.3) = 0.9635581854171931, it = 55, dsin=-1.484916792996379e-17, took 0.981924 s.
sin2(1.3) = 0.963558185417193, it = 10, dsin=4.835779466409166e-18, took 0.471155 s.
sin11(1.7) = 0.9916648104524686, it = 55, dsin=-7.152306208153797e-18, took 0.991410 s.
sin2(1.7) = 0.9916648104524686, it = 10, dsin=4.239843809646965e-17, took 0.472351 s.
```
sin.c
+144 -83
View File
@@ -7,86 +7,96 @@
#define x0 0.60725293500888125616944675250493 #define x0 0.60725293500888125616944675250493
#define y0 0 #define y0 0
static double atanlist[100] = { static double atanlist[100] = {
0.78539816339744830961566084581988,0.46364760900080611621425623146121, 0.78539816339744830961566084581988, 0.46364760900080611621425623146121,
0.24497866312686415417208248121128,0.12435499454676143503135484916387, 0.24497866312686415417208248121128, 0.12435499454676143503135484916387,
0.062418809995957348473979112985505,0.031239833430268276253711744892491, 0.062418809995957348473979112985505, 0.031239833430268276253711744892491,
0.01562372862047683080280152125657,0.0078123410601011112964633918421993, 0.01562372862047683080280152125657, 0.0078123410601011112964633918421993,
0.0039062301319669718276286653114244,0.0019531225164788186851214826250767, 0.0039062301319669718276286653114244, 0.0019531225164788186851214826250767,
0.00097656218955931943040343019971729,0.00048828121119489827546923962564485, 0.00097656218955931943040343019971729, 0.00048828121119489827546923962564485,
0.00024414062014936176401672294325966,0.00012207031189367020423905864611796, 0.00024414062014936176401672294325966, 0.00012207031189367020423905864611796,
0.000061035156174208775021662569173829,0.000030517578115526096861825953438536, 0.000061035156174208775021662569173829, 0.000030517578115526096861825953438536,
0.000015258789061315762107231935812698,7.6293945311019702633884823401051e-6, 0.000015258789061315762107231935812698, 7.6293945311019702633884823401051e-6,
3.814697265606496282923075616373e-6,1.9073486328101870353653693059172e-6, 3.814697265606496282923075616373e-6, 1.9073486328101870353653693059172e-6,
9.5367431640596087942067068992311e-7,4.7683715820308885992758382144925e-7, 9.5367431640596087942067068992311e-7, 4.7683715820308885992758382144925e-7,
2.3841857910155798249094797721893e-7,1.1920928955078068531136849713792e-7, 2.3841857910155798249094797721893e-7, 1.1920928955078068531136849713792e-7,
5.9604644775390554413921062141789e-8,2.980232238769530367674013276771e-8, 5.9604644775390554413921062141789e-8, 2.980232238769530367674013276771e-8,
1.4901161193847655147092516595963e-8,7.4505805969238279871365645744954e-9, 1.4901161193847655147092516595963e-8, 7.4505805969238279871365645744954e-9,
3.7252902984619140452670705718119e-9,1.8626451492309570290958838214765e-9, 3.7252902984619140452670705718119e-9, 1.8626451492309570290958838214765e-9,
9.3132257461547851535573547768456e-10,4.6566128730773925777884193471057e-10, 9.3132257461547851535573547768456e-10, 4.6566128730773925777884193471057e-10,
2.3283064365386962890204274183882e-10,1.1641532182693481445259909272985e-10, 2.3283064365386962890204274183882e-10, 1.1641532182693481445259909272985e-10,
5.8207660913467407226496761591232e-11,2.9103830456733703613273032698904e-11, 5.8207660913467407226496761591232e-11, 2.9103830456733703613273032698904e-11,
1.4551915228366851806639597837363e-11,7.2759576141834259033201841046704e-12, 1.4551915228366851806639597837363e-11, 7.2759576141834259033201841046704e-12,
3.6379788070917129516601402005838e-12,1.818989403545856475830076118823e-12, 3.6379788070917129516601402005838e-12, 1.818989403545856475830076118823e-12,
9.0949470177292823791503881172787e-13,4.5474735088646411895751949990348e-13, 9.0949470177292823791503881172787e-13, 4.5474735088646411895751949990348e-13,
2.2737367544323205947875976170669e-13,1.1368683772161602973937988232271e-13, 2.2737367544323205947875976170669e-13, 1.1368683772161602973937988232271e-13,
5.6843418860808014869689941345026e-14,2.8421709430404007434844970695472e-14, 5.6843418860808014869689941345026e-14, 2.8421709430404007434844970695472e-14,
1.4210854715202003717422485350606e-14,7.1054273576010018587112426756617e-15, 1.4210854715202003717422485350606e-14, 7.1054273576010018587112426756617e-15,
3.5527136788005009293556213378757e-15,1.7763568394002504646778106689434e-15, 3.5527136788005009293556213378757e-15, 1.7763568394002504646778106689434e-15,
8.8817841970012523233890533447242e-16,4.440892098500626161694526672363e-16, 8.8817841970012523233890533447242e-16, 4.440892098500626161694526672363e-16,
2.2204460492503130808472633361816e-16,1.1102230246251565404236316680908e-16, 2.2204460492503130808472633361816e-16, 1.1102230246251565404236316680908e-16,
5.5511151231257827021181583404541e-17,2.7755575615628913510590791702271e-17, 5.5511151231257827021181583404541e-17, 2.7755575615628913510590791702271e-17,
1.3877787807814456755295395851135e-17,6.9388939039072283776476979255676e-18, 1.3877787807814456755295395851135e-17, 6.9388939039072283776476979255676e-18,
3.4694469519536141888238489627838e-18,1.7347234759768070944119244813919e-18, 3.4694469519536141888238489627838e-18, 1.7347234759768070944119244813919e-18,
8.6736173798840354720596224069595e-19,4.3368086899420177360298112034798e-19, 8.6736173798840354720596224069595e-19, 4.3368086899420177360298112034798e-19,
2.1684043449710088680149056017399e-19,1.0842021724855044340074528008699e-19, 2.1684043449710088680149056017399e-19, 1.0842021724855044340074528008699e-19,
5.4210108624275221700372640043497e-20,2.7105054312137610850186320021749e-20, 5.4210108624275221700372640043497e-20, 2.7105054312137610850186320021749e-20,
1.3552527156068805425093160010874e-20,6.7762635780344027125465800054371e-21, 1.3552527156068805425093160010874e-20, 6.7762635780344027125465800054371e-21,
3.3881317890172013562732900027186e-21,1.6940658945086006781366450013593e-21, 3.3881317890172013562732900027186e-21, 1.6940658945086006781366450013593e-21,
8.4703294725430033906832250067964e-22,4.2351647362715016953416125033982e-22, 8.4703294725430033906832250067964e-22, 4.2351647362715016953416125033982e-22,
2.1175823681357508476708062516991e-22,1.0587911840678754238354031258496e-22, 2.1175823681357508476708062516991e-22, 1.0587911840678754238354031258496e-22,
5.2939559203393771191770156292478e-23,2.6469779601696885595885078146239e-23, 5.2939559203393771191770156292478e-23, 2.6469779601696885595885078146239e-23,
1.3234889800848442797942539073119e-23,6.6174449004242213989712695365597e-24, 1.3234889800848442797942539073119e-23, 6.6174449004242213989712695365597e-24,
3.3087224502121106994856347682799e-24,1.6543612251060553497428173841399e-24, 3.3087224502121106994856347682799e-24, 1.6543612251060553497428173841399e-24,
8.2718061255302767487140869206996e-25,4.1359030627651383743570434603498e-25, 8.2718061255302767487140869206996e-25, 4.1359030627651383743570434603498e-25,
2.0679515313825691871785217301749e-25,1.0339757656912845935892608650875e-25, 2.0679515313825691871785217301749e-25, 1.0339757656912845935892608650875e-25,
5.1698788284564229679463043254373e-26,2.5849394142282114839731521627186e-26, 5.1698788284564229679463043254373e-26, 2.5849394142282114839731521627186e-26,
1.2924697071141057419865760813593e-26,6.4623485355705287099328804067966e-27, 1.2924697071141057419865760813593e-26, 6.4623485355705287099328804067966e-27,
3.2311742677852643549664402033983e-27,1.6155871338926321774832201016991e-27, 3.2311742677852643549664402033983e-27, 1.6155871338926321774832201016991e-27,
8.0779356694631608874161005084957e-28,4.0389678347315804437080502542479e-28, 8.0779356694631608874161005084957e-28, 4.0389678347315804437080502542479e-28,
2.0194839173657902218540251271239e-28,1.009741958682895110927012563562e-28, 2.0194839173657902218540251271239e-28, 1.009741958682895110927012563562e-28,
5.0487097934144755546350628178098e-29,2.5243548967072377773175314089049e-29, 5.0487097934144755546350628178098e-29, 2.5243548967072377773175314089049e-29,
1.2621774483536188886587657044525e-29,6.3108872417680944432938285222623e-30, 1.2621774483536188886587657044525e-29, 6.3108872417680944432938285222623e-30,
3.1554436208840472216469142611311e-30,1.5777218104420236108234571305656e-30 3.1554436208840472216469142611311e-30, 1.5777218104420236108234571305656e-30
}; };
int trig(double theta, double *xy, double *dxy){ int sign(double alpha)
{
return alpha > 0 ? 1 : (alpha < 0 ? -1 : 0);
// return (int)copysign(1.0, alpha);
}
int trig(double theta, double *xy, double *dxy)
{
int i = 0; int i = 0;
double t = 1.0; //t=2^-i double t = 1.0; // t=2^-i
int d = 0; int d = 0;
double alpha = theta; double alpha = theta;
double dnorm = 1.0;
xy[0] = x0; xy[0] = x0;
xy[1] = y0; xy[1] = y0;
while(i<100 && dnorm>1e-32){ for (i = 0; i < 56; i++)
d = alpha>0?1:(alpha<0?-1:0); {
dxy[0] = -d*xy[1]*t; d = alpha > 0 ? 1 : (alpha < 0 ? -1 : 0);
dxy[1] = d*xy[0]*t; dxy[0] = -d * xy[1] * t;
dnorm = dxy[0]*dxy[0] + dxy[1]*dxy[1]; dxy[1] = d * xy[0] * t;
xy[0] += dxy[0]; xy[0] += dxy[0];
xy[1] += dxy[1]; xy[1] += dxy[1];
alpha -= d*atanlist[i]; alpha -= d * atanlist[i];
i++; t /= 2.0;
t/=2.0;
} }
return i; return i;
} }
double sin1(double theta, int *it, double *dsin){ double sin1(double theta, int *it, double *dsin)
if(theta < 0) return -sin1(theta, it, dsin); {
if(theta > 2*PI) return sin1(theta - 2*PI, it, dsin); if (theta < 0)
if(theta > PI) return -sin1(theta - PI, it, dsin); return -sin1(theta, it, dsin);
if(theta > PI/2) return sin1(PI - theta, it, dsin); if (theta > 2 * PI)
return sin1(theta - 2 * PI, it, dsin);
if (theta > PI)
return -sin1(theta - PI, it, dsin);
if (theta > PI / 2)
return sin1(PI - theta, it, dsin);
double xy[2]; double xy[2];
double dxy[2]; double dxy[2];
@@ -96,42 +106,93 @@ double sin1(double theta, int *it, double *dsin){
return xy[1]; return xy[1];
} }
double sin2(double theta, int *it, double *dsin){ // just do sin, instead of trig()
if(theta < 0) return -sin2(theta, it, dsin); double sin11(double theta, int *it, double *dsin)
if(theta > 2*PI) return sin2(theta - 2*PI, it, dsin); {
if(theta > PI) return -sin2(theta - PI, it, dsin); if (theta < 0)
if(theta > PI/2) return sin2(PI - theta, it, dsin); return -sin11(theta, it, dsin);
if (theta > 2 * PI)
return sin11(theta - 2 * PI, it, dsin);
if (theta > PI)
return -sin11(theta - PI, it, dsin);
if (theta > PI / 2)
return sin11(PI - theta, it, dsin);
double xy[2];
double dxy[2];
int i = 0;
double t = 1.0; // t=2^-i
int d = 0;
double alpha = theta;
xy[0] = x0;
xy[1] = y0;
dxy[1] = 1;
for (i = 0; i < 55; i++)
{
// d = sign(alpha);
dxy[0] = -xy[1] * copysign(t,alpha);
dxy[1] = xy[0] * copysign(t,alpha);
xy[0] += dxy[0];
xy[1] += dxy[1];
alpha -= d * atanlist[i];
t /= 2.0;
}
*dsin = dxy[1];
*it = i;
return xy[1];
}
double sin2(double theta, int *it, double *dsin)
{
if (theta < 0)
return -sin2(theta, it, dsin);
if (theta > 2 * PI)
return sin2(theta - 2 * PI, it, dsin);
if (theta > PI)
return -sin2(theta - PI, it, dsin);
if (theta > PI / 2)
return sin2(PI - theta, it, dsin);
double result = 0; double result = 0;
double term=theta; double term = theta;
int i = 0; int i = 0;
while (i<100 && fabs(term)>1e-16){ while (i < 100 && fabs(term) > 1e-16)
{
result += term; result += term;
i += 1; i += 1;
term *= -1; term *= -1;
term /= 2*i; term /= 2 * i;
term /= (2*i+1); term /= (2 * i + 1);
term *= theta*theta; term *= theta * theta;
} }
*it = i; *it = i;
*dsin = term; *dsin = term;
return result; return result;
} }
int main(){ int main()
{
int it; int it;
double theta = 1.0; double thetalist[5] = {1e-4, 0.5, 1, 1.3, 1.7};
double sin,dsin; double sin, dsin;
clock_t start, end; clock_t start, end;
for (int i = 0; i < 5; i++)
{
start = clock(); start = clock();
for (int i=0; i<10000; i++) sin = sin1(theta, &it, &dsin); for (int j = 0; j < 10000000; j++)
sin = sin11(thetalist[i], &it, &dsin);
end = clock(); end = clock();
printf("sin1(%.16g) = %.16g, it = %d, dsin=%.16g, took %.16g s.\n", theta, sin, it, dsin, (double)(end-start)/CLOCKS_PER_SEC); printf("sin11(%.16g) = %.16g, it = %d, dsin=%.16g, took %f s.\n", thetalist[i], sin, it, dsin, (double)(end - start) / CLOCKS_PER_SEC);
start = clock(); start = clock();
for (int i=0; i<10000; i++) sin = sin2(theta, &it, &dsin); for (int j = 0; j < 10000000; j++)
sin = sin2(thetalist[i], &it, &dsin);
end = clock(); end = clock();
printf("sin2(%.16g) = %.16g, it = %d, dsin=%.16g, took %.16g s.\n", theta, sin, it, dsin, (double)(end-start)/CLOCKS_PER_SEC); printf("sin2(%.16g) = %.16g, it = %d, dsin=%.16g, took %f s.\n", thetalist[i], sin, it, dsin, (double)(end - start) / CLOCKS_PER_SEC);
}
return 0; return 0;
} }