%I #58 Sep 27 2018 12:42:14
%S 2,4,8,10,14,18,28,30,38,44,58,60,64,78,80,84,94,98,120,140,144,148,
%T 164,170,198,214,218,220,228,240,248,254,270,304,318,338,340,344,350,
%U 368,408,410,430,470,480,484,494,500,504,520,528,534,578,604,630,634,644,658
%N Intersections with the x-axis of a bouncing ball on a Sophie Germain billiard table.
%C In the first quadrant of a coordinate system define a rectangular Sophie Germain billiard table with width p and length 2p+1, with vertices (0,0), (p,0), (p,2p+1) and (0,2p+1). A billiard ball (considered to be a point) starts from (0,0) at an angle of 45 degrees and hits the sides exactly p times until it hits the x-axis. The sequence gives the intersections with the x-axis of consecutive Sophie Germain prime numbers (p > 3) after p bounces.
%C The sum of all crossed lattice points (including the rectangle sides) is the sum of crossed points left under, right middle and left up respectively ((p+7)/6)^2 + (p+1)(p+4)/18 + (p+1)(p+7)/36 = ((p+4)/3)^2 (see bouncing examples).
%C The enclosed areas in the Sophie Germain billiard table also correspond to ((p+4)/3)^2.
%C The number of trajectories is a subsequence of A176045.
%C The number of trajectories with slope +1 or with slope -1 is a subsequence of A124485.
%C The sum of a term of this sequence and the corresponding Sophie Germain prime is A317510 and it appears that this is a subsequence of A179882. Checked up to and including 33295 of A317510 (Sophie Germain prime 24971).
%H Samuel King, <a href="https://nl.mathworks.com/matlabcentral/fileexchange/58354-billiard-simulator">Billiard Simulator</a>
%H Hilko Koning, <a href="http://www.hilko.net/bounce.pdf">Some bouncing examples</a>.
%F a(n) = (A005384(n)+1)/3 for n>=3. - _Michel Marcus_, Aug 25 2018
%t lst = {}; Do[If[PrimeQ[p] && PrimeQ[2 p + 1], AppendTo[lst, (p + 1)/3]], {p, 5, 2*10^3}]; lst
%t (Select[Prime@ Range[3, 300], PrimeQ[2# + 1] &] + 1)/3 (* _Robert G. Wilson v_, Aug 02 2018 *)
%o (PARI) lista(nn) = forprime(p=2, nn, if (isprime(2*p+1), print1((p+1)/3, ", "));); \\ _Michel Marcus_, Aug 25 2018
%o (GAP) a:=[];; for p in [3..2000] do if IsPrime(p) and IsPrime(2*p+1) then Add(a,(p+1)/3); fi; od; a; # _Muniru A Asiru_, Aug 28 2018
%Y Cf. A005384, A151922, A176045, A124485, A179882, A317510.
%K nonn
%O 1,1
%A _Hilko Koning_, Aug 02 2018