3600 sayısını alalım. Çok özeldir, çünkü
3600 = 482 + 362
3600 = 202 + 2×402
3600 = 302 + 3×302
3600 = 452 + 7×152.
Benzer şekilde 88201 = 992 + 2802 = 2872 + 2×542 = 2832 + 3×522 = 1972 + 7×842 olduğunu görebiliriz.
1747'de Euler, hangi sayıların iki kare toplamı olarak yazılabileceğini ispatladı. Biz burada aşağıdaki dört gösterime de uyan n sayıları ile ilgileniyoruz:
n = a12 + b12
n = a22 + 2 b22
n = a32 + 3 b32
n = a72 + 7 b72,
n = a22 + 2 b22
n = a32 + 3 b32
n = a72 + 7 b72,
Burada $a_k$ ve $b_k$ pozitif tam sayılardır.
$10^7$'den küçük 75373 tane böyle sayı vardır.
$2 \times 10^9$'dan küçük böyle kaç tane sayı vardır?
Cevap: 11325263
C++:
#include <stdio.h>
#include <stdlib.h>
#define lim 2000000000
int main() {
int a,a2,b,i,j,n,si=lim/8,ct;
unsigned int *A,Bit[32],x;
A=(unsigned int*)(malloc)(si*sizeof(unsigned int));
Bit[0]=1;
for(n=1;n<32;n++) Bit[n]=2*Bit[n-1];
for(i=0;i<si;i++) A[i]=0;
for(a=1;2*a*a<lim;a++) {
a2=a*a;
for(b=a;a2+b*b<lim;b++) {
n=a2+b*b;
A[n>>3]|=Bit[(n&7)<<2];
}
}
for(a=1;a*a<lim;a++) {
a2=a*a;
for(b=1;a2+2*b*b<lim;b++) {
n=a2+2*b*b;
A[n>>3]|=Bit[1+((n&7)<<2)];
}
}
for(a=1;a*a<lim;a++) {
a2=a*a;
for(b=1;a2+3*b*b<lim;b++) {
n=a2+3*b*b;
A[n>>3]|=Bit[2+((n&7)<<2)];
}
}
for(a=1;a*a<lim;a++) {
a2=a*a;
for(b=1;a2+7*b*b<lim;b++) {
n=a2+7*b*b;
A[n>>3]|=Bit[3+((n&7)<<2)];
}
}
ct=0;
for(i=0;i<si;i++) {
x=A[i];
for(j=0;j<8;j++) {
if((x&15)==15) ct++;
x>>=4;
}
}
printf("%d\n",ct);
return 0;
}Java:
public class T229 {
public static final int solve(int n) {
byte DIV1 = 0x1;
byte DIV2 = 0x2;
byte DIV3 = 0x4;
byte DIV7 = 0x8;
byte OK = 0xf;
byte[] sieve = new byte[n];
long ul = (long) Math.floor(Math.sqrt(n)) +2 ;
//Sieves 7
for (long a = 1; a < ul; a++) {
for (long b = 1; b < ul; b++) {
long c = a*a + 7*b*b;
if(c >= n) break;
sieve[(int)c] |= DIV7;
}
}
//Sieves 3
for (long a = 1; a < ul; a++) {
for (long b = 1; b < ul; b++) {
long c = a*a + 3*b*b;
if(c >= n) break;
sieve[(int)c] |= DIV3;
}
}
//Sieves 2
for (long a = 1; a < ul; a++) {
for (long b = 1; b < ul; b++) {
long c = a*a + 2*b*b;
if(c >= n) break;
sieve[(int)c] |= DIV2;
}
}
//Sieves 1
for (long a = 1; a < ul; a++) {
for (long b = 1; b <= a; b++) {
long c = a*a + b*b;
if(c >= n) break;
sieve[(int)c] |= DIV1;
}
}
int count = 0;
for (int i = 0; i < n; i++) {
if(sieve[i] == OK) ++ count;
}
return count;
}
public static final void main(String[] args){
int n = 10000000;
System.out.println("R(" + n + ") = " + solve(n) );
n = 2000000000;
System.out.println("R(" + n + ") = " + solve(n) );
}
}Python:
import math
nmax = 2 * 10 ** 9
ns = int(math.sqrt(nmax))
primes = (ns + 1) * [True]
primes[0:2] = [False,False]
s = [1]
for x in range(2, len(primes)):
if primes[x]:
for y in range(2 * x, len(primes), x):
primes[y] = False
if x % 168 in (1, 25, 121):
for j in range(len(s)):
if x * s[j] > nmax:
break
s.append(x * s[j])
s.sort()
p = [x for x in range(len(primes)) if primes[x]]
def factor(x):
factors = {}
for q in p:
if q * q > x:
break
i = 0
while x % q == 0:
x //= q
i += 1
if i:
factors[q] = i
if x > 1:
factors[x] = 1
return factors
cc = sum(int(math.sqrt(nmax // x)) for x in s[1:])
pp = {d:[True for i in range(nmax // 168 + 1)] for d in (1, 25, 121)}
for x in p:
for k in range(0, x * 168, x):
q, r = divmod(k, 168)
if r in pp:
pq = pp[r]
for y in range(0, len(pq) - q, x):
pq[y+q] = False
sq = [int(math.sqrt(x)) for x in range(ns + 1)]
for y, pq in pp.items():
z = y
for x in range(len(pq)):
if pq[x] and ns < z <= nmax:
for j in s:
k = sq[nmax // (z * j)]
if not k:
break
cc += k
z += 168
for x in range(1, ns + 1):
f = factor(x)
if any(x for x in f if x % 4 == 1) and \
any(x for x in f if x == 3 or x % 8 in (1, 3)) and \
any(x for x in f if x == 2 or x % 6 == 1) and \
any(x for x in f if x == 2 and f[2] >= 2 or x % 14 in (1, 9, 11)):
cc += 1
print(cc)Pascal:
program bidon;
uses minimaple;
var t:array[1..1000000000] of shortint;
procedure rempli;
var a,b,k:longint;
dem:extended;
begin
dem:=2000000000;
for a:=1 to trunc(sqrt(dem/2)) do
begin
for b:=a to trunc(sqrt(dem-a*a)) do
begin
k:=a*a+b*b;
if k<=1000000000 then t[k]:=t[k] mod 10 + 10;
if k>1000000000 then
begin
k:=k-1000000000;
if k<=1000000000 then t[k]:=(t[k] div 10)*10+1;
end;
end;
end;
for a:=1 to trunc(sqrt(dem)) do
begin
for b:=1 to trunc(sqrt((dem-a*a)/2)) do
begin
k:=a*a+2*b*b;
if k<=1000000000 then if t[k] div 10=1 then t[k]:=t[k] mod 10 + 20;
if k>1000000000 then
begin
k:=k-1000000000;
if k<=1000000000 then if t[k] mod 10=1 then t[k]:=(t[k] div 10)*10+2;
end;
end;
end;
for a:=1 to trunc(sqrt(dem)) do
begin
for b:=1 to trunc(sqrt((dem-a*a)/3)) do
begin
k:=a*a+3*b*b;
if k<=1000000000 then if t[k] div 10=2 then t[k]:=t[k] mod 10 + 30;
if k>1000000000 then
begin
k:=k-1000000000;
if k<=1000000000 then if t[k] mod 10=2 then t[k]:=(t[k] div 10)*10+3;
end;
end;
end;
for a:=1 to trunc(sqrt(dem)) do
begin
for b:=1 to trunc(sqrt((dem-a*a)/7)) do
begin
k:=a*a+7*b*b;
if k<=1000000000 then if t[k] div 10=3 then t[k]:=t[k] mod 10 + 70;
if k>1000000000 then
begin
k:=k-1000000000;
if k<=1000000000 then if t[k] mod 10=3 then t[k]:=(t[k] div 10)*10+7;
end;
end;
end;
end;
procedure finalise;
var som,i:longint;
begin
som:=0;
for i:=1 to 10000000 do
if t[i] div 10=7 then som:=som+1;
writeln('vérif ',som);
som:=0;
for i:=1 to 1000000000 do
begin
if t[i] mod 10=7 then som:=som+1;
if t[i] div 10=7 then som:=som+1;
end;
writeln('résultat ',som);
end;
begin
rempli;
finalise;
end.