%I M2221 #75 Oct 29 2023 01:49:37

%S 1,3,1,5,1,5,7,5,3,5,9,1,3,7,11,7,11,13,9,7,1,15,13,15,1,13,9,5,17,13,

%T 11,9,5,17,7,17,19,1,3,15,17,7,21,19,5,11,21,19,13,1,23,5,17,19,25,13,

%U 25,23,1,5,15,27,9,19,25,17,11,5,25,27,23,29,29,25,23,19,29,13,31,31

%N a(n) is the odd member of {x,y}, where x^2 + y^2 is the n-th prime of the form 4i+1.

%C a(n)^2 + 4*A002973(n)^2 = A002144(n); A002331(n+1) = Min(a(n),2*A002973(n)) and A002330(n+1) = Max(a(n),2*A002973(n)). - _Reinhard Zumkeller_, Feb 16 2010

%C It appears that the terms in this sequence are the absolute values of the terms in A046730. - _Gerry Myerson_, Dec 02 2010

%C (a(n) - 1)/2 = A208295(n), n >= 1. - _Wolfdieter Lang_, Mar 03 2012

%C a(A267858(k)) == 1 (mod 4), k >= 1. - _Wolfdieter Lang_, Feb 18 2016

%C "the n-th prime of the form 4i+1" is A005098(n). - _Rainer Rosenthal_, Aug 24 2022

%D E. Kogbetliantz and A. Krikorian, Handbook of First Complex Prime Numbers, Gordon and Breach, NY, 1971, p. 243.

%D N. J. A. Sloane and Simon Plouffe, The Encyclopedia of Integer Sequences, Academic Press, 1995 (includes this sequence).

%H Rainer Rosenthal, <a href="/A002972/b002972.txt">Table of n, a(n) for n = 1..10000</a> (first 1000 terms from T. D. Noe)

%H S. R. Finch, <a href="https://arxiv.org/abs/math/0701251">Powers of Euler's q-Series</a>, arXiv:math/0701251 [math.NT], 2007.

%H E. Kogbetliantz and A. Krikorian, <a href="/A002970/a002970.pdf">Handbook of First Complex Prime Numbers</a>, Gordon and Breach, NY, 1971. [Annotated scans of a few pages]

%H Stan Wagon, <a href="https://doi.org/10.2307/2323912">Editor’s Corner: The Euclidean Algorithm Strikes Again</a>, The American Mathematical Monthly, vol. 97, no. 2, 1990, pp. 125-29. [Description of efficient decomposition algorithm implemented in PARI program]

%F a(n) = Min(A173330(n), A002144(n) - A173330(n)). - _Reinhard Zumkeller_, Feb 16 2010

%e The 2nd prime of the form 4i+1 is 13 = 2^2 + 3^2, so a(2)=3.

%t pmax = 1000; odd[p_] := Module[{k, m}, 2m+1 /. ToRules[Reduce[k>0 && m >= 0 && (2k)^2 + (2m+1)^2 == p, {k, m}, Integers]]]; For[n=1; p=5, p<pmax, p = NextPrime[p], If[Mod[p, 4] == 1, a[n] = odd[p]; Print["a(", n, ") = ", a[n]]; n++]]; Array[a, n-1] (* _Jean-François Alcover_, Feb 26 2016 *)

%o (PARI) decomp2sq(p) = {my (m=(p-1)/4, r, x, limit=ceil(sqrt(p))); if (p>4 && denominator(m)==1, forprime (c=2,oo, if (!issquare(Mod(c,p)), r=c; break)); x=lift (Mod(r,p)^m); until (p<limit, r=p%x; p=x; x=r); if(p^2+x^2==4*m+1, [p,x], [0,0]), [0,0])};

%o forprime (p=5, 1000, if (p%4==1, print1(select(x->x%2,decomp2sq(p))[1],", "))) \\ _Hugo Pfoertner_, Aug 27 2022

%Y Cf. A002144, A002973, A005098, A261858.

%K nonn

%O 1,2

%A _N. J. A. Sloane_

%E Better description from _Jud McCranie_, Mar 05 2003