Henry Baker has raised the matter of deterministic algorithms for square-root (mod p) . The Bumby paper seems to have been published in the late 1980's, after a burst of activity in this area. There does not appear to have been much progress since: a 2010 paper by Tsz-Wo Sze at arxiv.org/pdf/0812.2591 implies that there is no known deterministic algorithm in time polynomial in log(p). WFL

...

On 4/30/13, Fred lunnon <fred.lunnon@gmail.com> wrote:

...

Solving w^2 + x^2 + y^2 + z^2 = p , an odd prime, is rather simpler than Bumby manages to make it look. His final section mentions the Tonelli-Shanks algorithm for computing square root modulo p .

So armed with this algorithm, compute a solution of a^2 + b^2 + c^2 = 0 (mod p) in time order log(p) roughly; then set Q = w + x i + y j + z k = GCD_L ( a i + b j + c k, p) , where GCD_L denotes left-hand quaternion GCD . Now ||Q|| = p , yielding a single solution.

If all solutions are required, the simplest procedure is to continue by stepping around all p+1 points of the conic a^2 + b^2 + c^2 = 0 (mod p) in the projective plane over |F_p , starting from the solution found above, and converting each point via GCD with p as before.

The resulting set {Q} comprises one representative from each unit-association class on the right: taking products with {1, i, j, k, -1, -i, -j, -k} yields the complete solution set in time order p log(p) , along with the cardinal 8(p+1) specified by Jacobi's theorem.

Quaternion GCD is a natural extension of the Euclidean algorithm, slightly complicated by non-commutativity, and potential non-termination if both norms are even. Further details and program code are available.

Fred Lunnon

On 4/29/13, Henry Baker <hbaker1@pipeline.com> wrote:

...

Bumby! Excellent! Just exactly what I wanted. Thanks, Neil!

http://www.math.rutgers.edu/~bumby/squares1.pdf

At 06:26 AM 4/29/2013, Neil Sloane wrote:

...

Entry A000118 in the OEIS has many references. See, in particular, the Bumby reference.

On Mon, Apr 29, 2013 at 9:14 AM, Henry Baker <hbaker1@pipeline.com> wrote:

...

I'm still on the hunt for a 'ridiculously simple' way to express a positive integer as 4 squares.

Most presentations for the 4-square theorem focus on proving that it is possible, but don't provide very good algorithms for actually doing it.

At 09:47 AM 4/26/2013, Dan Asimov wrote: >How does that work with quaternions H or octonions O ? > >You can't just switch one coefficient. (Or you can, but >then you're just operating in a single complex subspace >of H or O, so it's really the same thing.) > >--Dan > >On 2013-04-26, at 5:50 AM, Fred lunnon wrote: > >> Etc. for quaternions and octonions, where available ... WFL >> >> On 4/26/13, Warut Roonguthai <warut822@gmail.com> wrote: >>> I don't know who discovered this, but I know that Pythagorean >>> triples can >>> be obtained from the fact that |z^2| = |z|^2. This is similar to >>> the way >>> you can prove that the product of numbers representable as sums of >>> two >>> squares is also representable as a sum of two squares by using the >>> fact >>> that |z1 * z2| = |z1| * |z2|. >>> >>> On Fri, Apr 26, 2013 at 3:04 PM, Bill Gosper <billgosper@gmail.com> >>> wrote: >>> >>>> In[570]:= (1+2*I)/(2+1*I) >>>> Out[570]= 4/5+(3 I)/5 >>>> >>>> In[571]:= (1+3*I)/(3+1*I) >>>> Out[571]= 3/5+(4 I)/5 >>>> >>>> In[572]:= (2+3*I)/(3+2*I) >>>> Out[572]= 12/13+(5 I)/13 >>>> >>>> In[573]:= (1+4*I)/(4+1*I) >>>> Out[573]= 8/17+(15 I)/17 >>>> >>>> In[586]:= Simplify/@ComplexExpand[(a+I*b)/(b+I*a)] >>>> Out[586]= (2 a b)/(a^2+b^2)-(I (a^2-b^2))/(a^2+b^2)

-- Dear Friends, I have now retired from AT&T. New coordinates:

Neil J. A. Sloane, President, OEIS Foundation 11 South Adelaide Avenue, Highland Park, NJ 08904, USA Phone: 732 828 6098; home page: http://NeilSloane.com Email: njasloane@gmail.com