inz

sssms2.c (raw)

  1. /*
  2.  *         DO WHAT THE FUCK YOU WANT TO PUBLIC LICENSE
  3.  *                     Version 2, December 2004
  4.  * 
  5.  * Copyright (C) 2021 Santtu Lakkala <inz@inz.fi>
  6.  *
  7.  * Everyone is permitted to copy and distribute verbatim or modified
  8.  * copies of this license document, and changing it is allowed as long
  9.  * as the name is changed.
  10.  * 
  11.  *           DO WHAT THE FUCK YOU WANT TO PUBLIC LICENSE
  12.  *  TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
  13.  * 
  14.  *  0. You just DO WHAT THE FUCK YOU WANT TO.
  15.  *
  16.  *  compile with something like:
  17.  *    cc -std=c11 sssms2.c -o sssms2 -lm -pthread
  18.  */
  19. #include <stdint.h>
  20. #include <math.h>
  21. #include <stdio.h>
  22. #include <stdlib.h>
  23. #include <stdbool.h>
  24. #include <inttypes.h>
  25. #include <stdatomic.h>
  26. #include <pthread.h>
  27.  
  28. static const uint8_t f2[8][8] = {
  29. 	{ 19, 1, 0, 0, 17, 0, 0, 0 }, 
  30. 	{ 1, 0, 0, 0, 0, 0, 0, 0 }, 
  31. 	{ 0, 0, 0, 136, 0, 0, 0, 136 }, 
  32. 	{ 0, 0, 136, 68, 0, 0, 136, 68 }, 
  33. 	{ 17, 0, 0, 0, 17, 0, 0, 0 }, 
  34. 	{ 0, 0, 0, 0, 0, 0, 0, 0 }, 
  35. 	{ 0, 0, 0, 136, 0, 0, 0, 136 }, 
  36. 	{ 0, 0, 136, 68, 0, 0, 136, 196 }
  37. };
  38.  
  39. static inline bool issquare(uint64_t x)
  40. {
  41. 	uint64_t sqr = sqrtl(x);
  42. 	return sqr * sqr == x;
  43. }
  44.  
  45. static _Atomic uint64_t ga = ATOMIC_VAR_INIT(1);
  46. static uint64_t high;
  47.  
  48. void *worker(void *data)
  49. {
  50. 	uint64_t a;
  51. 	uint64_t b;
  52. 	uint64_t c;
  53.  
  54. 	(void)data;
  55.  
  56. 	while ((a = ga++) < high) {
  57. 		if ((a & 7) == 5)
  58. 			continue;
  59.  
  60. 		for (b = 1; b <= a; b++) {
  61. 			const uint8_t f3 = f2[a & 7][b & 7];
  62. 			if (f3) {
  63. 				uint64_t a2b2 = a * a + b * b;
  64. 				uint64_t mrc = sqrtl(a2b2 + b);
  65. 				uint64_t rc;
  66.  
  67. 				for (rc = sqrt(a2b2 - 1) + 1;
  68. 				     rc <= mrc; rc++) {
  69. 					c = rc * rc - a2b2;
  70.  
  71. 					if (f3 & (1 << (c & 7)) &&
  72. 					    issquare(a * a + b + c * c) &&
  73. 					    issquare(a + b * b + c * c))
  74. 						printf("%" PRIu64 " %" PRIu64 " %" PRIu64 "\n", c, b, a);
  75. 				}
  76. 			}
  77. 		}
  78. 	}
  79.  
  80. 	return NULL;
  81. }
  82.  
  83. int main(int argc, char **argv)
  84. {
  85. 	int t = 1;
  86. 	int i;
  87. 	high = 5000;
  88.  
  89. 	if (argc >= 3) {
  90. 		ga = strtoull(argv[1], NULL, 10);
  91. 		high = strtoull(argv[2], NULL, 10);
  92.  
  93. 		if (argc >= 4)
  94. 			t = atoi(argv[3]);
  95. 	}
  96.  
  97. 	pthread_t threads[t - 1];
  98.  
  99. 	for (i = 0; i < t - 1; i++)
  100. 		pthread_create(&threads[i], NULL, worker, NULL);
  101. 	worker(NULL);
  102.  
  103. 	for (i = 0; i < t - 1; i++)
  104. 		pthread_join(threads[i], NULL);
  105.  
  106. 	return 0;
  107. }
  108.