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

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!