abc仮説の数値検証（はい、その1つ）

 // To compile: g++ abc.cpp -O3 -fopenmp -oabc #include <string.h> #include <math.h> #include <stdbool.h> #include <stdint.h> #include <stdio.h> #include <vector> #include <set> #include <map> #include <algorithm> #include <time.h> typedef unsigned long int valType; typedef std::vector<valType> valList; typedef std::set<valType> valSet; typedef valList::iterator valListIterator; std::vector<valList> valFactors; std::vector<double> valRads; valList factors(valType n) { valList results; valType z = 2; while (z * z <= n) { if (n % z == 0) { results.push_back(z); n /= z; } else { z++; } } if (n > 1) { results.push_back(n); } return results; } valList unique_factors(valType n) { valList results = factors(n); valSet vs(results.begin(), results.end()); valList unique(vs.begin(), vs.end()); std::sort(unique.begin(), unique.end()); return unique; } double rad(valType n) { valList f = valFactors[n]; double result = 1; for (valListIterator it=f.begin(); it<f.end(); it++) { result *= *it; } return result; } bool not_mutual_primes(valType a, valType b, valType c) { valList res1 = valFactors[a], res2 = valFactors[b], res3; // = valFactors[c]; valList c12, c13, c23; set_intersection(res1.begin(),res1.end(), res2.begin(),res2.end(), back_inserter(c12)); if (c12.size() != 0) return false; res3 = valFactors[c]; set_intersection(res1.begin(),res1.end(), res3.begin(),res3.end(), back_inserter(c13)); if (c13.size() != 0) return false; set_intersection(res2.begin(),res2.end(), res3.begin(),res3.end(), back_inserter(c23)); return c23.size() == 0; } int main() { time_t start_t, end_t; time(&start_t); int cnt=0; double S = 1.2; valType N_MAX = 10000000; printf("Getting prime factors...\n"); valFactors.resize(2*N_MAX+2); valRads.resize(2*N_MAX+2); for(valType val=1; val<=2*N_MAX+1; val++) { valFactors[val] = unique_factors(val); valRads[val] = rad(val); } time(&end_t); printf("Done, T = %.2fs\n", difftime(end_t, start_t)); printf("Calculating...\n"); #pragma omp parallel for reduction(+:cnt) for(int a=1; a<=N_MAX; a++) { for(int b=a; b<=N_MAX; b++) { int c = a+b; if (c > pow(valRads[a]*valRads[b]*valRads[c], S) && not_mutual_primes(a,b,c)) { printf("%d + %d = %d\n", a,b,c); cnt += 1; } } } printf("Done, cnt=%d\n", cnt); time(&end_t); float diff_t = difftime(end_t, start_t); printf("N=%lld, T = %.2fs\n", N_MAX, diff_t); } 

 from __future__ import print_function import math import time import multiprocessing prime_factors_list = [] rad_list = [] def prime_factors(n): factors = [] # Print the number of two's that divide n while n % 2 == 0: factors.append(int(2)) n = n / 2 # n must be odd at this point so a skip of 2 ( i = i + 2) can be used for i in range(3, int(math.sqrt(n)) + 1, 2): # while i divides n , print i ad divide n while n % i == 0: factors.append(int(i)) n = n / i # Condition if n is a prime number greater than 2 if n > 2: factors.append(int(n)) return factors def rad(n): result = 1 for num in prime_factors_list[n]: result *= num return result def not_mutual_primes(a,b,c): fa, fb, fc = prime_factors_list[a], prime_factors_list[b], prime_factors_list[c] return len(fa.intersection(fb)) == 0 and len(fa.intersection(fc)) == 0 and len(fb.intersection(fc)) == 0 def calculate(N): S = 1.2 cnt = 0 for a in range(1, N): for b in range(a, N): c = a+b if c > (rad_list[a]*rad_list[b]*rad_list[c])**S and not_mutual_primes(a, b, c): print("{} + {} = {}".format(a, b, c)) cnt += 1 print("N: {}, CNT: {}".format(N, cnt)) return cnt if __name__ == '__main__': t1 = time.time() NMAX = 100000 prime_factors_list = [0]*(2*NMAX+2) rad_list = [0]*(2*NMAX+2) for p in range(1, 2*NMAX+2): prime_factors_list[p] = set(prime_factors(p)) rad_list[p] = rad(p) calculate(NMAX) print("Done", time.time() - t1)