-
Notifications
You must be signed in to change notification settings - Fork 0
/
euler_prob273.py
90 lines (83 loc) · 3.51 KB
/
euler_prob273.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
import numpy as np
import sympy as sp
import math
from sympy.solvers.diophantine import cornacchia
# 2032447591196869022
# Objective of this problem: To find the smaller value of (a,b) where a**2 + b**2 = given value of n
# At first glance, you might think the (a,b) pairs represent nothing.
# To understand why this problem is truly difficult, refer to the pdf document for more info.
# This code basically generates all the squarefree numbers required and solves 65535 values of n to get ALL the possible (a,b) pairs.
nums = []
sum = 0
def GenNum(num, ind):
""" Generates Squarefree Numbers from prime numbers of 4*k + 1 < 150 """
for a in range(exp[ind]+1):
if ind == len(subprimes)-1:
nums.append(num*subprimes[ind])
else:
if a+1 < exp[ind]+1:
nums.append(num*subprimes[ind])
GenNum(num, ind+1)
num *= subprimes[ind]
# Generate all prime numbers below 150
primes = list(sp.sieve.primerange(1, 150))
subprimes = []
for p in primes:
if p%4 == 1: # Take only prime numbers of 4*k+1
subprimes.append(p)
exp = [1 for i in range(len(subprimes))]
GenNum(1, 0)
nums = np.unique(np.sort(nums))
for n in nums:
sum += np.sum(np.unique(np.sort(list(cornacchia(1,1,n)), axis = 1), axis = 0)[:,0]) # Find possible a and b that sum of squares = n
print(n, sum)
print("Ans: ", sum)
# Appendix: Partial Sums of the last few values of n
"""
************************************************
* n * Partial Sum of S(n) *
************************************************
2650349212148801960854945 438361883167503705
2791544286713887969868305 443465771270435289
3171729385030533494137885 448906838435992373
3376082419233843034699789 454519564158817836
3608581638922830790317565 460323377560668066
3650480990318161191366245 466159827222344194
3671797667633887680074953 472014329472757085
3998698572860434119000545 478123723181125666
4451648499697722231595297 484570108445690793
4615011747393051487800629 491131846181094054
4718914450899086418107585 497769030201087439
4980557232335075368022461 504586158355693055
5101787834339174565621385 511487347454324292
5185941035730336207940913 518443464295241621
5229067364509798463308405 525430228425415220
5652093038942051822137849 532695148564069192
6671568706443535970427965 540587453735053936
6890907951586885098222857 548607854034277816
8246496401079387084758501 557381073878526762
9491250574827219097552237 566792021488391795
11380911322756620184847705 577097356453940261
12269177572337624687079721 587799473550659426
13595575147725476004601853 599065029208900244
16880412096169215173498945 624169561685514597
17346078636753193523112709 636896227747343573
18358988338169438400374765 663074609318915858
22258242498488611157976485 691899077434130497
23075058736965257439003145 721252430125054959
24902786161675376840112305 751744234682838224
25929705178651681039704565 782858833634592963
28260465194710259110689245 815334747202170619
29590369439167212480604033 831951366568660804
34454539757934425491114285 867814031151231117
38695098497372508628482197 886815199249933108
41232482005396935423792505 926045883265013654
47456252874136095487761185 968138618662917329
61345887861688123435398605 1015991453686706542
67977875738627380023009265 1066364919077249597
86730393183765967615563545 1123260613277207686
147951847195836062403020165 1197583793525631629
193475492486862543142410985 1282576666455437438
503036280465842612170268561 1419624364143463619
2515181402329213060851342805 2032447591196869022
"""