sssms2

Unnamed repository; edit this file 'description' to name the repository.
git clone https://git.inz.fi/sssms2/
Log | Files | Refs

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 }