Help Need to write a Mathematica function

[INdeWATERS]

I need help writing a Mathematica 8 function that takes a prime p as input and returns a list of the nonequivalent Pythagorean triples mod p.

For example, PythagoreanTriples[13] should return the list;
{{1,3,6},{1,4,2},{2,5,4},{2,6,1},{3,4,5},{5,6,3}}

 

[Hurkyl]

What are you having trouble with? Coming up with an algorithm? How to write a program? What the Mathematica syntax is for doing certain things?

Have you programmed before? Have you programmed in Mathematica?

 

[Dale]

How did you generate that list from the number 13? (I don't know what "nonequivalent Phthagorean triples mod 13" are).

 

[INdeWATERS]

I am very new to programming and have very little experience. But, I need help coming up with an algorithm. I know that the most general method for generating Pythagorean triples is due to Euclid. Euclid says to take positive integers s and t such that s > t, gcd(s,t) = 1, and s+t is odd. Then (s^2 - t^2,2st,s^2 + t^2) is a pythagorean triple. Furthermore, every pythagorean triple can be obtained using Euclid's formula. For example, letting s = 6 and t = 5, we obtain the PT (6^2 - 5^2,2(6)(5),6^2 + 5^2) = (11,60,61).

ANy help in writing the algorithm is much appreciated.

 

[INdeWATERS]

How did you generate that list from the number 13? (I don't know what "nonequivalent Phthagorean triples mod 13" are).

Sorry, if (a,b,c) is a PT mod p, then so is (b,a,c). We consider all of these to be "equivalent." Thus, in listing PTs mod p, we make the restriction that a should be less than or equal to b, and both b and c should be less than p/2.

 

[INdeWATERS]

So, I am interested in finding solutions to the equation

x^2 + y^2 = z^2 (mod p)

where p is a given prime. For instance, if the prime were 17, then (6,8,7) is a PT mod 17, since 6^2 mod 17 = 36 mod 17 = 2, 8^2 mod 17 = 64 mod 17 = 13, 7^2 mod 17 = 49 mod 17 = 15, and 2 + 13 = 15.

BUT! For each prime p > 5, there are at least (p - 1)/2 "different" PTs mod p. However, this bound is not always sharp. For instance, there are 16 different PTs mod 17:

(1,1,6), (1,5,3), (1,7,4), (2,2,5), (2,3,8), (2,7,6), (3,3,1), (3,4,5),
(4,4,7), (4,6,1), (5,5,4), (5,8,2), (6,6,2), (6,8,7), (7,7,8), (8,8,3)

What do I mean by "different?" Well, for any x with 0 < x < p, x^2 mod p = (p - x)^2 mod p. Thus, if (a,b,c) is a PT mod p, then so are all of the following:

(a,b,p-c), (a,p-b,c), (a,p-b,p-c), (p-a,b,c), (p-a,b,p-c), (p-a,p-b,c), (p-a,p-b,p-c).

 

[Hurkyl]

By far the easiest algorithm for finding triples to code is to simply take the list of all (a,b,c) and keep only those satisfying a^2 + b^2 = c^2 (mod p).

To deal with equivalent triples, the easiest algorithm would be to replace every triple with a canonical representative from its equivalent class. Slightly less difficult would be to find a way to iterate over a subset of (a,b,c)'s that touches each equivalence class of triples exactly once. (And once you have this list, check which ones are pythagorean)

 

[Simon_Tyler]

Hi INdeWATERS,

Before I post the code that I've written, I want to see some evidence that you've attempted the problem yourself.
What's you current best code?

My brute force Mathematica 8 code reproduces the examples for p=13 and p=17 you gave above.
And can compute the triples for the 100th prime = 541 in 5 seconds. There are 18 090 such triples.

 

[INdeWATERS]

Hi INdeWATERS,

Before I post the code that I've written, I want to see some evidence that you've attempted the problem yourself.
What's you current best code?

My brute force Mathematica 8 code reproduces the examples for p=13 and p=17 you gave above.
And can compute the triples for the 100th prime = 541 in 5 seconds. There are 18 090 such triples.

Hello,
I haven't attempted to put together a complete code just yet. I have the following so far;

PythagoreanTriples[p_] := Module[{i, j, result, s, x, y}]
s := Table[0, {j, 1, p-1)]
result = {}

I know I am missing A LOT of statements. But, like I said, I am very new to programming, and have been scratching my head for days over this.

My idea was to have s indicate the squares mod p. That is, s[[j]] = i iff i^2 = j (modp)
Then, the variable will hold my answer; initially, result is empty.

 

[GetToTheChopa]

I tried this:

PythagoreanTrip[p_] /; z>0 := With[
{result = Reduce[((Mod[x^2,p]+Mod[y^2,p]==Mod[z^2,p] && x>y) || Mod[x^2,p]+Mod[z^2,p]==Mod[y^2,p]) &&
Element[{x,y},Integers] && x>0 && y>0]},
If [result===False, {},
Sort[Map[Sort, {x,y,z} /. {ToRules[result]}]]]]

This does not work, just gives me errors. I've now been at it for hours and feel like I have accomplished nothing.

 

[Simon_Tyler]

OK, you are missing a lot of code... in any case, here's the code that I wrote over breakfast.

Implementing the simplest, brute force algorithm - but making sure not to include any of the equivalent forms.
That is, we want triples (a, b, c) such that a^2 + b^2 = c^2 mod p.
Because of the equivalence generated by the following:
(a,b,c) == (b,a,c) == (p-a,b,c) == (a,b,p-c)
1) We only need to check a,b,c <= p/2.
2) We only need to check a <= b

triples[p_?PrimeQ] := Module[{out},
  out = Flatten[Table[{k, j, i}, {i, p/2}, {j, p/2}, {k, j}], 2];
  Select[out, 0 == Mod[{1, 1, -1}.#^2, p] &]]

We can check that this reproduces the result mod 13 that you gave:

In[2]:= Sort[triples[13]] == 
        Sort[{{1, 3, 6}, {1, 4, 2}, {2, 5, 4}, {2, 6, 1}, {3, 4, 5}, {5, 6, 3}}]

Out[2]= True

But, it gets quite slow and memory inefficient for larger primes.
So here's a C-compiled version of the same thing
(but with the memory inefficient Table replaced with a Do loop):

triplesC[p_?PrimeQ] := Rest@triplesC$Helper[p]
triplesC$Helper = Compile[{{p, _Integer}},
   Module[{p2 = Floor[p/2], out = {{0, 0, 0}}, tf = False},
    Do[tf = Mod[{1, 1, -1}.{k, j, i}^2, p] == 0;
        If[tf, out = Append[out, {k, j, i}]], 
        {i, 1, p2}, {j, 1, p2}, {k, 1, j}];
    out], CompilationTarget -> "C"];

Compare timings to see that the compiled version is a lot faster:

In[5]:= p = Prime[50]
        t = triples[p]; // Timing
        tC = triplesC[p]; // Timing
        t == tC

Out[5]= 229

Out[6]= {7.43, Null}

Out[7]= {0.18, Null}

Out[8]= True

If you need to go to really large numbers, maybe split up
Mod[{1, 1, -1}.{k, j, i}^2, p]
to use PowerMod, so that you don't get any overflows.

 

[INdeWATERS]

Simon_Tyler

Wow! That totally works! But I was wondering if you could provide just a brief explanation of what each part of the code is actually doing... I am so new to programming and this code makes very little sense to me.
Thank you very much for your time.

triples[p_?PrimeQ] := Module[{out},
  out = Flatten[Table[{k, j, i}, {i, p/2}, {j, p/2}, {k, j}], 2];
  Select[out, 0 == Mod[{1, 1, -1}.#^2, p] &]]

This code is sufficient enough.

 

[Simon_Tyler]

triples[p_?PrimeQ] := ...

Defines a function called triples that takes a single argument p. This definition will only fire if p is prime.

Module[{out}, ...]

defines a local and temporary variable out. It's not really needed, but helps break the code up.

Table[{k, j, i}, {i, p/2}, {j, p/2}, {k, j}]

returns a nested list of all possible triples. Note that j >= k and that i,j,k<=p/2 so that we have a unique representation for all equivalent forms for the Pythagorean triples.
This list of list of list of triples is flattened to a plain list of triples by Flatten.

"Select" takes that list and only keeps the triples that satisfy the condition
0 == Mod[{1, 1, -1}.#^2, p] &[triple]
It is written in a pure function notation (that's the # and &), Mathematica will evaluate it something like
0 == Mod[{1, 1, -1}.#^2, p]&[{a,b,c}]
0 == Mod[{1, 1, -1}.{a,b,c}^2, p]
0 == Mod[{1, 1, -1}.{a^2,b^2,c^2}, p]
0 == Mod[a^2 + b^2 - c^2, p]

This is the list that is returned.

The compiled version works similarly. Except instead of generating all potential triples to begin with, it steps through them one by one and only keeps those that are Pythagorean.