Menu

Subscribe

Project Euler problem 46

on project euler, erlang, python

Goldbach’s other conjecture

It was proposed by Christian Goldbach that every odd composite number can be written as the sum of a prime and twice a square.

        9  = 7  + 2×1^2
        15 = 7  + 2×2^2
        21 = 3  + 2×3^2
        25 = 7  + 2×3^2
        27 = 19 + 2×2^2
        33 = 31 + 2×1^2

It turns out that the conjecture was false.

What is the smallest odd composite that cannot be written as the sum of a prime and twice a square?

Link to original description
Source code examples on Github

Python version

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
#!/usr/bin/python3

import prime

MAX = 10000
prime._refresh(MAX)
squares = dict.fromkeys((x*x for x in range(1, MAX)), 1)

for x in range(35, MAX, 2):
    if not prime.isprime(x):
        is_goldbach = 0
        for p in prime.prime_list:
            if p >= x: break
            key = (x-p)/2
            if key in squares:
                is_goldbach = 1
                break
        if not is_goldbach:
            print(x)
            break

Erlang version

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
91
92
93
94
95
96
97
98
99
100
101
102
103
#!/usr/bin/env escript
%% -*- erlang -*-
%%! -smp enable -sname p46
% vim:syn=erlang

-mode(compile).
-define(MAX, 10000).

main(_) -> 
    get_primes(?MAX),
    get_squares(?MAX),
    io:format("Answer ~p ~n", [t46(35)]).

t46(N) when N < ?MAX ->
    case is_prime(N) of
        true -> t46(N+2)
        ;_   -> 
                case is_goldbach(N, get("primes")) of
                    true -> t46(N+2)
                    ;_   -> N
                end
    end;
t46(N) -> {stop, N}.

is_goldbach(N, [H|T]) when N >= H ->
   case (N-H) rem 2 of
        0 ->
           case ets:lookup(squares, (N-H) div 2) of
                [] -> is_goldbach(N, T)
                ;_ -> true
           end;
        _ -> is_goldbach(N, T) 
    end;
is_goldbach(_, _) -> false.

is_prime(N) ->
    case ets:lookup(prim, N) of
        [] -> false
        ;_ -> true
    end.

get_squares(N) ->
    ets:new(squares, [public, named_table, {write_concurrency, true} ]),
    lists:foreach(fun(E)-> ets:insert(squares, {E*E,1}) end, lists:seq(1, N)).
%----------------------------------------------- prime generator from Project Euler 10 (version 5 ---------------------------)
get_primes(N) ->
    ets:new(comp, [public, named_table, {write_concurrency, true} ]),
    ets:new(prim, [public, named_table, {write_concurrency, true}]),
    composite_mc(N),
    primes_mc(N),
    lists:sort([P || {P,_} <-ets:tab2list(prim)]).

primes_mc(N) ->
    case erlang:system_info(schedulers) of
        1 -> primes(N);
        C -> launch_primes(lists:seq(1,C), C, N, N div C)
    end.
launch_primes([1|T], C, N, R) -> P = self(), spawn(fun()-> primes(2,R), P ! {ok, prm} end), launch_primes(T, C, N, R);
launch_primes([H|[]], C, N, R)-> P = self(), spawn(fun()-> primes(R*(H-1)+1,N), P ! {ok, prm} end), wait_primes(C);
launch_primes([H|T], C, N, R) -> P = self(), spawn(fun()-> primes(R*(H-1)+1,R*H), P ! {ok, prm} end), launch_primes(T, C, N, R).

wait_primes(0) -> ok;
wait_primes(C) ->
    receive
        {ok, prm} -> wait_primes(C-1)
    after 1000    -> wait_primes(C)
    end.

primes(N) -> primes(2, N).
primes(I,N) when I =< N ->
    case ets:lookup(comp, I) of
        [] -> ets:insert(prim, {I,1})
        ;_ -> ok
    end,
    primes(I+1, N);
primes(I,N) when I > N -> ok.


composite_mc(N) -> composite_mc(N,2,round(math:sqrt(N)),erlang:system_info(schedulers)).
composite_mc(N,I,M,C) when I =< M, C > 0 ->
    C1 = case ets:lookup(comp, I) of
        [] -> comp_i_mc(I*I, I, N), C-1
        ;_ -> C
    end,
    composite_mc(N,I+1,M,C1);
composite_mc(_,I,M,_) when I > M -> ok;
composite_mc(N,I,M,0) ->
    receive
        {ok, cim} -> composite_mc(N,I,M,1)
    after 1000    -> composite_mc(N,I,M,0)
    end.

comp_i_mc(J, I, N) -> 
    Parent = self(),
    spawn(fun() ->
        comp_i(J, I, N),
        Parent ! {ok, cim}
    end).

comp_i(J, I, N) when J =< N -> ets:insert(comp, {J, 1}), comp_i(J+I, I, N);
comp_i(J, _, N) when J > N -> ok.