From c8874275c3c6c26445889eb6e28d22c24eac29ef Mon Sep 17 00:00:00 2001 From: Yingjie Wang Date: Sat, 29 Jun 2024 15:28:20 -0400 Subject: [PATCH] init: add sin.c --- sin.c | 137 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 137 insertions(+) create mode 100644 sin.c diff --git a/sin.c b/sin.c new file mode 100644 index 0000000..38b7191 --- /dev/null +++ b/sin.c @@ -0,0 +1,137 @@ +#include +#include +#include + +#define PI 3.14159265358979323846264338327950288419716939937 + +#define x0 0.60725293500888125616944675250493 +#define y0 0 +static double atanlist[100] = { + 0.78539816339744830961566084581988,0.46364760900080611621425623146121, + 0.24497866312686415417208248121128,0.12435499454676143503135484916387, + 0.062418809995957348473979112985505,0.031239833430268276253711744892491, + 0.01562372862047683080280152125657,0.0078123410601011112964633918421993, + 0.0039062301319669718276286653114244,0.0019531225164788186851214826250767, + 0.00097656218955931943040343019971729,0.00048828121119489827546923962564485, + 0.00024414062014936176401672294325966,0.00012207031189367020423905864611796, + 0.000061035156174208775021662569173829,0.000030517578115526096861825953438536, + 0.000015258789061315762107231935812698,7.6293945311019702633884823401051e-6, + 3.814697265606496282923075616373e-6,1.9073486328101870353653693059172e-6, + 9.5367431640596087942067068992311e-7,4.7683715820308885992758382144925e-7, + 2.3841857910155798249094797721893e-7,1.1920928955078068531136849713792e-7, + 5.9604644775390554413921062141789e-8,2.980232238769530367674013276771e-8, + 1.4901161193847655147092516595963e-8,7.4505805969238279871365645744954e-9, + 3.7252902984619140452670705718119e-9,1.8626451492309570290958838214765e-9, + 9.3132257461547851535573547768456e-10,4.6566128730773925777884193471057e-10, + 2.3283064365386962890204274183882e-10,1.1641532182693481445259909272985e-10, + 5.8207660913467407226496761591232e-11,2.9103830456733703613273032698904e-11, + 1.4551915228366851806639597837363e-11,7.2759576141834259033201841046704e-12, + 3.6379788070917129516601402005838e-12,1.818989403545856475830076118823e-12, + 9.0949470177292823791503881172787e-13,4.5474735088646411895751949990348e-13, + 2.2737367544323205947875976170669e-13,1.1368683772161602973937988232271e-13, + 5.6843418860808014869689941345026e-14,2.8421709430404007434844970695472e-14, + 1.4210854715202003717422485350606e-14,7.1054273576010018587112426756617e-15, + 3.5527136788005009293556213378757e-15,1.7763568394002504646778106689434e-15, + 8.8817841970012523233890533447242e-16,4.440892098500626161694526672363e-16, + 2.2204460492503130808472633361816e-16,1.1102230246251565404236316680908e-16, + 5.5511151231257827021181583404541e-17,2.7755575615628913510590791702271e-17, + 1.3877787807814456755295395851135e-17,6.9388939039072283776476979255676e-18, + 3.4694469519536141888238489627838e-18,1.7347234759768070944119244813919e-18, + 8.6736173798840354720596224069595e-19,4.3368086899420177360298112034798e-19, + 2.1684043449710088680149056017399e-19,1.0842021724855044340074528008699e-19, + 5.4210108624275221700372640043497e-20,2.7105054312137610850186320021749e-20, + 1.3552527156068805425093160010874e-20,6.7762635780344027125465800054371e-21, + 3.3881317890172013562732900027186e-21,1.6940658945086006781366450013593e-21, + 8.4703294725430033906832250067964e-22,4.2351647362715016953416125033982e-22, + 2.1175823681357508476708062516991e-22,1.0587911840678754238354031258496e-22, + 5.2939559203393771191770156292478e-23,2.6469779601696885595885078146239e-23, + 1.3234889800848442797942539073119e-23,6.6174449004242213989712695365597e-24, + 3.3087224502121106994856347682799e-24,1.6543612251060553497428173841399e-24, + 8.2718061255302767487140869206996e-25,4.1359030627651383743570434603498e-25, + 2.0679515313825691871785217301749e-25,1.0339757656912845935892608650875e-25, + 5.1698788284564229679463043254373e-26,2.5849394142282114839731521627186e-26, + 1.2924697071141057419865760813593e-26,6.4623485355705287099328804067966e-27, + 3.2311742677852643549664402033983e-27,1.6155871338926321774832201016991e-27, + 8.0779356694631608874161005084957e-28,4.0389678347315804437080502542479e-28, + 2.0194839173657902218540251271239e-28,1.009741958682895110927012563562e-28, + 5.0487097934144755546350628178098e-29,2.5243548967072377773175314089049e-29, + 1.2621774483536188886587657044525e-29,6.3108872417680944432938285222623e-30, + 3.1554436208840472216469142611311e-30,1.5777218104420236108234571305656e-30 +}; + +int trig(double theta, double *xy, double *dxy){ + int i = 0; + double t = 1.0; //t=2^-i + int d = 0; + double alpha = theta; + double dnorm = 1.0; + + xy[0] = x0; + xy[1] = y0; + while(i<100 && dnorm>1e-32){ + d = alpha>0?1:(alpha<0?-1:0); + dxy[0] = -d*xy[1]*t; + dxy[1] = d*xy[0]*t; + dnorm = dxy[0]*dxy[0] + dxy[1]*dxy[1]; + xy[0] += dxy[0]; + xy[1] += dxy[1]; + alpha -= d*atanlist[i]; + i++; + t/=2.0; + } + return i; +} + +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 > PI) return -sin1(theta - PI, it, dsin); + if(theta > PI/2) return sin1(PI - theta, it, dsin); + + double xy[2]; + double dxy[2]; + + *it = trig(theta, xy, dxy); + *dsin = dxy[1]; + 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 term=theta; + int i = 0; + while (i<100 && fabs(term)>1e-16){ + result += term; + i += 1; + term *= -1; + term /= 2*i; + term /= (2*i+1); + term *= theta*theta; + } + *it = i; + *dsin = term; + return result; +} + +int main(){ + int it; + double theta = 1.0; + double sin,dsin; + clock_t start, end; + + start = clock(); + for (int i=0; i<10000; i++) sin = sin1(theta, &it, &dsin); + end = clock(); + printf("sin1(%.16g) = %.16g, it = %d, dsin=%.16g, took %.16g s.\n", theta, sin, it, dsin, (double)(end-start)/CLOCKS_PER_SEC); + + start = clock(); + for (int i=0; i<10000; i++) sin = sin2(theta, &it, &dsin); + end = clock(); + printf("sin2(%.16g) = %.16g, it = %d, dsin=%.16g, took %.16g s.\n", theta, sin, it, dsin, (double)(end-start)/CLOCKS_PER_SEC); + return 0; +}