fork(2) download 
  1. #include <stdint.h>
  2. #include <stdio.h>
  3. #include <stdlib.h>
  4. #include <math.h>
  5. #include <assert.h>
  6.  
  7. typedef uint32_t I_TYPE;
  8.  
  9. int FIG = 5;
  10. I_TYPE RADIX = 100000;
  11.  
  12. const double pi = 3.1415926535897932385;
  13.  
  14. double cos1;
  15. double sin1;
  16. double err_max;
  17.  
  18. #define max(a,b)		((a)<(b)?(b):(a))
  19.  
  20. void fr(double* s1, double* s2, double* end, size_t add, size_t t90) {
  21. 	for (; s1 < end; s1 += add, s2 += add) {
  22. 		double r0, r1, r2, r3, m0, m1;
  23.  
  24. 		r0 = s1[0];
  25. 		r1 = s2[0];
  26. 		r2 = s1[t90];
  27. 		r3 = s2[t90];
  28.  
  29. 		m0 = r2 * cos1 - r3 * sin1;
  30. 		m1 = r2 * sin1 + r3 * cos1;
  31.  
  32. 		s1[0] = r0 + m0;
  33. 		s2[0] = r0 - m0;
  34. 		s1[t90] = -r1 + m1;
  35. 		s2[t90] = r1 + m1;
  36. 	}
  37. }
  38.  
  39. void fri(double* s1, double* s2, double* end, size_t add, size_t t90) {
  40. 	for (; s1 < end; s1 += add, s2 += add) {
  41. 		double r0, r1, r2, r3, m0, m1, m2, m3;
  42.  
  43. 		r0 = s1[0];
  44. 		r1 = s2[0];
  45. 		r2 = s1[t90];
  46. 		r3 = s2[t90];
  47.  
  48. 		m0 = r0 + r1;
  49. 		m1 = -r2 + r3;
  50. 		m2 = r0 - r1;
  51. 		m3 = r2 + r3;
  52.  
  53. 		s1[0] = m0;
  54. 		s2[0] = m1;
  55. 		s1[t90] = m2 * cos1 + m3 * sin1;
  56. 		s2[t90] = -m2 * sin1 + m3 * cos1;
  57. 	}
  58. }
  59.  
  60. void fr0(double* s1, double* end, size_t add, size_t t90) {
  61. 	for (; s1 < end; s1 += add) {
  62. 		double r0, r1;
  63.  
  64. 		r0 = s1[0];
  65. 		r1 = s1[t90];
  66.  
  67. 		s1[0] = r0 + r1;
  68. 		s1[t90] = r0 - r1;
  69. 	}
  70. }
  71.  
  72. void swap_d(double* a, double* b) {
  73. 	double t = *a;
  74. 	*a = *b;
  75. 	*b = t;
  76. }
  77.  
  78. void bxchg(double* s, int p) {
  79. 	size_t t4, t5, t6, t045, t090, t180, t270;
  80.  
  81. 	t045 = (size_t)1 << p - 3;
  82. 	t090 = t045 + t045;
  83. 	t180 = t090 + t090;
  84. 	t270 = t180 + t090;
  85.  
  86. 	t4 = t090 - 4;
  87. 	t5 = t4;
  88. 	do {
  89. 		swap_d(&s[t4 + 1], &s[t180 + t5]);
  90. 		swap_d(&s[t4 + 2], &s[t090 + t5]);
  91. 		swap_d(&s[t4 + 3], &s[t270 + t5]);
  92. 		swap_d(&s[t090 + t4 + 1], &s[t180 + t5 + 2]);
  93. 		swap_d(&s[t090 + t4 + 3], &s[t270 + t5 + 2]);
  94. 		swap_d(&s[t180 + t4 + 3], &s[t270 + t5 + 1]);
  95.  
  96. 		if (t4 < t5) {
  97. 			swap_d(&s[t4], &s[t5]);
  98. 			swap_d(&s[t090 + t4 + 2], &s[t090 + t5 + 2]);
  99. 			swap_d(&s[t180 + t4 + 1], &s[t180 + t5 + 1]);
  100. 			swap_d(&s[t270 + t4 + 3], &s[t270 + t5 + 3]);
  101. 		}
  102. 		if (t4 == 0) {
  103. 			break;
  104. 		}
  105.  
  106. 		t4 -= 4;
  107. 		t6 = t045;
  108. 		for (;;) {
  109. 			t5 ^= t6;
  110. 			if ((t5 & t6) == 0) {
  111. 				break;
  112. 			}
  113. 			t6 >>= 1;
  114. 		}
  115. 	} while (1);
  116. }
  117.  
  118. void bitexchg(double* s, int p) {
  119. 	if (p <= 1) {
  120. 	}
  121. 	else if (p == 2) {
  122. 		swap_d(&s[1], &s[2]);
  123. 	}
  124. 	else if (p == 3) {
  125. 		swap_d(&s[1], &s[4]);
  126. 		swap_d(&s[3], &s[6]);
  127. 	}
  128. 	else {
  129. 		bxchg(s, p);
  130. 	}
  131. }
  132.  
  133. double* SIN_TABLE[64] = { 0 };
  134.  
  135. const double* sin_table(int p) {
  136. 	if (p >= 2 && SIN_TABLE[p] == 0) {
  137. 		size_t t090, i;
  138. 		double k;
  139.  
  140. 		t090 = (size_t)1 << p - 2;
  141. 		SIN_TABLE[p] = (double*)malloc(t090 * sizeof(double));
  142. 		if (SIN_TABLE[p] == NULL) {
  143. 			printf("malloc failed. sin_table %d\r\n", p);
  144. 			exit(0);
  145. 		}
  146. 		k = 0.5 * pi / t090;
  147. 		for (i = 0; i < t090; i++) {
  148. 			SIN_TABLE[p][i] = sin(k * i);
  149. 		}
  150. 	}
  151. 	return SIN_TABLE[p];
  152. }
  153.  
  154. void fouri(double* s, int p) {
  155. 	const double* table;
  156. 	size_t t090, t360;
  157. 	int i;
  158.  
  159. 	table = sin_table(p);
  160. 	t360 = (size_t)1 << p;
  161. 	t090 = t360 >> 2;
  162.  
  163. 	bitexchg(s, p);
  164.  
  165. 	for (i = 1; i <= p; i++) {
  166. 		size_t r090, r180, r360, k1, k2, t0, t1, t2;
  167.  
  168. 		r360 = (size_t)1 << i;
  169. 		r180 = r360 >> 1;
  170. 		r090 = r180 >> 1;
  171.  
  172. 		fr0(&s[0], &s[t360], r360, r180);
  173.  
  174. 		k1 = 1;
  175. 		k2 = r180 - 1;
  176. 		t2 = (size_t)1 << p - i;
  177. 		t0 = t2;
  178. 		t1 = t090 - t2;
  179. 		for (; k1 < r090; k1++, k2--) {
  180. 			sin1 = table[t0];
  181. 			cos1 = table[t1];
  182. 			fr(&s[k1], &s[k2], &s[t360], r360, r180);
  183. 			t0 += t2;
  184. 			t1 -= t2;
  185. 		}
  186. 	}
  187. }
  188.  
  189. void fouriinv(double* s, int p) {
  190. 	const double* table;
  191. 	size_t t090, t360;
  192. 	int i;
  193.  
  194. 	table = sin_table(p);
  195. 	t360 = (size_t)1 << p;
  196. 	t090 = t360 >> 2;
  197.  
  198. 	for (i = p; i >= 1; i--) {
  199. 		size_t r090, r180, r360, k1, k2, t0, t1, t2;
  200.  
  201. 		r360 = (size_t)1 << i;
  202. 		r180 = r360 >> 1;
  203. 		r090 = r180 >> 1;
  204.  
  205. 		fr0(&s[0], &s[t360], r360, r180);
  206.  
  207. 		k1 = 1;
  208. 		k2 = r180 - 1;
  209. 		t2 = (size_t)1 << p - i;
  210. 		t0 = t2;
  211. 		t1 = t090 - t2;
  212. 		for (; k1 < r090; k1++, k2--) {
  213. 			sin1 = table[t0];
  214. 			cos1 = table[t1];
  215. 			fri(&s[k1], &s[k2], &s[t360], r360, r180);
  216. 			t0 += t2;
  217. 			t1 -= t2;
  218. 		}
  219. 	}
  220. 	bitexchg(s, p);
  221. }
  222.  
  223. typedef struct I_ {
  224. 	I_TYPE* data;
  225. 	size_t capacity;
  226. 	size_t size;
  227. } I;
  228.  
  229. typedef struct F_ {
  230. 	double* data;
  231. 	size_t capacity;
  232. 	int p;
  233. } F;
  234.  
  235. void alloc_i(I* a, size_t cap) {
  236. 	a->data = (I_TYPE*)malloc(cap * sizeof(I_TYPE));
  237. 	if (a->data == NULL) {
  238. 		printf("malloc failed. alloc_i %zu\r\n", cap);
  239. 		exit(0);
  240. 	}
  241. 	a->capacity = cap;
  242. }
  243.  
  244. void free_i(I* a) {
  245. 	free(a->data);
  246. }
  247.  
  248. void alloc_f(F* a, int p) {
  249. 	size_t cap = (size_t)1 << p;
  250. 	a->data = (double*)malloc(cap * sizeof(double));
  251. 	if (a->data == NULL) {
  252. 		printf("malloc failed. alloc_f %d\r\n", p);
  253. 		exit(0);
  254. 	}
  255. 	a->capacity = cap;
  256. }
  257.  
  258. void free_f(F* a) {
  259. 	free(a->data);
  260. }
  261.  
  262. void swap_i(I* a, I* b) {
  263. 	I t = *a;
  264. 	*a = *b;
  265. 	*b = t;
  266. }
  267.  
  268. void i2f(F* a, const I* b, int p) {
  269. 	size_t i;
  270. 	size_t t360;
  271. 	int carry;
  272.  
  273. 	t360 = (size_t)1 << p;
  274.  
  275. 	assert(b->size <= t360);
  276. 	assert(t360 <= a->capacity);
  277.  
  278. 	a->p = p;
  279. 	carry = 0;
  280. 	for (i = 0; i < b->size - 1; i++) {
  281. 		int t = b->data[i] + carry;
  282. 		if (t < (int)RADIX / 2) {
  283. 			a->data[i] = t;
  284. 			carry = 0;
  285. 		}
  286. 		else {
  287. 			a->data[i] = (int)t - (int)RADIX;
  288. 			carry = 1;
  289. 		}
  290. 	}
  291. 	a->data[i] = b->data[i] + carry;
  292. 	i++;
  293. 	for (; i < t360; i++) {
  294. 		a->data[i] = 0.;
  295. 	}
  296. 	fouri(a->data, p);
  297. }
  298.  
  299. void mul_f2i(I* a, const F* b, const F* c, F* d) {
  300. 	int p;
  301. 	size_t i, s, t, t180, t360, k1, k2;
  302. 	int64_t carry;
  303. 	double k;
  304.  
  305. 	p = b->p;
  306. 	t180 = (size_t)1 << p - 1;
  307. 	t360 = t180 + t180;
  308.  
  309. 	assert(p == c->p);
  310. 	assert(t360 <= b->capacity);
  311. 	assert(t360 <= c->capacity);
  312. 	assert(t360 <= d->capacity);
  313.  
  314. 	d->data[0] = b->data[0] * c->data[0] * .5;
  315. 	d->data[t180] = b->data[t180] * c->data[t180] * .5;
  316.  
  317. 	k1 = 1;
  318. 	k2 = t360 - 1;
  319. 	for (; k1 < t180; k1++, k2--) {
  320. 		double re, im;
  321.  
  322. 		re = b->data[k1] * c->data[k1] - b->data[k2] * c->data[k2];
  323. 		im = b->data[k1] * c->data[k2] + b->data[k2] * c->data[k1];
  324. 		d->data[k1] = re;
  325. 		d->data[k2] = im;
  326. 	}
  327.  
  328. 	fouriinv(d->data, p);
  329.  
  330. 	s = a->capacity < t360 ? a->capacity : t360;
  331.  
  332. 	carry = 0;
  333. 	t = 1;
  334. 	k = ldexp(1., -p + 1);
  335. 	for (i = 0; i < s; i++) {
  336. 		double x;
  337. 		int64_t e, f;
  338.  
  339. 		x = d->data[i] * k;
  340. 		e = (int64_t)floor(x + .5);
  341. 		err_max = max(err_max, fabs(x - e));
  342. 		e += carry;
  343. 		carry = e / RADIX;
  344. 		f = e - carry * RADIX;
  345. 		if (f < 0) {
  346. 			f += RADIX;
  347. 			carry--;
  348. 		}
  349. 		a->data[i] = (I_TYPE)f;
  350. 		if (f != 0) {
  351. 			t = i + 1;
  352. 		}
  353. 	}
  354. 	a->size = t;
  355. 	if (carry != 0) {
  356. 		if (a->size < a->capacity) {
  357. 			assert(0 <= carry && carry < RADIX);
  358. 			a->data[a->size++] = (I_TYPE)carry;
  359. 		}
  360. 		else {
  361. 			assert(0);
  362. 		}
  363. 	}
  364. 	for (; i < t360; i++) {
  365. 		double x;
  366. 		int64_t e;
  367.  
  368. 		x = d->data[i] * k;
  369. 		e = (int64_t)floor(x + .5);
  370. 		err_max = max(err_max, fabs(x - e));
  371. 		assert(e == 0);
  372. 	}
  373. }
  374.  
  375. void add_i(I* a, const I* b, const I* c) {
  376. 	I_TYPE carry, d;
  377. 	size_t i;
  378.  
  379. 	if (b->size > c->size) {
  380. 		const I* t = b;
  381. 		b = c;
  382. 		c = t;
  383. 	}
  384. 	assert(c->size <= a->capacity);
  385.  
  386. 	carry = 0;
  387. 	for (i = 0; i < b->size; i++) {
  388. 		d = b->data[i] + c->data[i] + carry;
  389. 		carry = 0;
  390. 		if (d >= RADIX) {
  391. 			d -= RADIX;
  392. 			carry++;
  393. 		}
  394. 		a->data[i] = d;
  395. 	}
  396. 	for (; i < c->size; i++) {
  397. 		d = c->data[i] + carry;
  398. 		carry = 0;
  399. 		if (d >= RADIX) {
  400. 			d -= RADIX;
  401. 			carry++;
  402. 		}
  403. 		a->data[i] = d;
  404. 	}
  405. 	a->size = c->size;
  406. 	if (carry != 0) {
  407. 		a->data[a->size++] = carry;
  408. 		assert(a->size <= a->capacity);
  409. 	}
  410. }
  411.  
  412. void print(const char* filename, const I* a) {
  413. 	char format[10];
  414. 	size_t s;
  415. 	FILE* fp;
  416.  
  417. 	assert(a->data != NULL);
  418. 	assert(a->size > 0);
  419. 	assert(a->size <= a->capacity);
  420.  
  421. 	fp = fopen(filename, "w");
  422.  
  423. 	sprintf(format, "%%0%du", FIG);
  424.  
  425. 	s = a->size;
  426. 	fprintf(fp, "%u", a->data[--s]);
  427. 	while (s != 0) {
  428. 		fprintf(fp, format, a->data[--s]);
  429. 	}
  430. 	fclose(fp);
  431. }
  432.  
  433. int s2p(size_t s1, size_t s2) {
  434. 	int i;
  435. 	for (i = 1; i < 64; i++) {
  436. 		if (((size_t)1 << i) + 1 >= s1 + s2) {
  437. 			break;
  438. 		}
  439. 	}
  440. 	return i;
  441. }
  442.  
  443. int main(int argc, char** argv) {
  444. 	int n = 34;
  445. 	I a;
  446. 	F x;
  447.  
  448. 	int p = s2p((1llu << n) * log(2) / log(10) / FIG + 1, 1);
  449. 	size_t cap = (size_t)1 << p;
  450. 	alloc_i(&a, cap);
  451. 	alloc_f(&x, p);
  452.  
  453. 	a.data[0] = 1;
  454. 	a.size = 1;
  455.  
  456. 	for (int i = 0; i < n; i++) {
  457.  
  458. 		int p = s2p(a.size, a.size);
  459. 		i2f(&x, &a, p);
  460. 		mul_f2i(&a, &x, &x, &x);
  461. 		add_i(&a, &a, &a);
  462.  
  463. 		{
  464. 			char filename[200];
  465. 			sprintf(filename, "%d.txt", i+1);
  466. 			print(filename, &a);
  467. 		}
  468. 		printf("%d %f\n", i + 1, err_max);
  469. 	}
  470.  
  471. 	free_i(&a);
  472. 	free_f(&x);
  473. 	return 0;
  474. }
  475.  
Runtime error #stdin #stdout 0s 5464KB
comments (?)
stdin 
Standard input is empty
stdout 
Standard output is empty
https://ideone.com/xSaa6x
language:
C (gcc 8.3)
created:
3 years ago
visibility:
secret

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!