sssms2.c (2211B)
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 const uint64_t a2b2 = a * a + b * b; 64 uint64_t rc; 65 66 for (rc = sqrt(a2b2 - 1) + 1; 67 (c = rc * rc - a2b2) <= b; rc++) { 68 if (f3 & (1 << (c & 7)) && 69 issquare(a * a + b + c * c) && 70 issquare(a + b * b + c * c)) 71 printf("%" PRIu64 " %" PRIu64 " %" PRIu64 "\n", c, b, a); 72 } 73 } 74 } 75 } 76 77 return NULL; 78 } 79 80 int main(int argc, char **argv) 81 { 82 int t = 1; 83 int i; 84 high = 5000; 85 86 if (argc >= 3) { 87 ga = strtoull(argv[1], NULL, 10); 88 high = strtoull(argv[2], NULL, 10); 89 90 if (argc >= 4) 91 t = atoi(argv[3]); 92 } 93 94 pthread_t threads[t - 1]; 95 96 for (i = 0; i < t - 1; i++) 97 pthread_create(&threads[i], NULL, worker, NULL); 98 worker(NULL); 99 100 for (i = 0; i < t - 1; i++) 101 pthread_join(threads[i], NULL); 102 103 return 0; 104 }