Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Suppose that we have a set of 9 2-dimensional vectors $c$ by

m=3;
n=2;
c=Tuples[Range[m],n];

We want to define a subset $S$ of vectors in $c$ by the following:

$S=\{ (x,y)\in c~|~ x^2+2y^2=5\}$

Also, I want to count the cardinality of $S$.

How can I do this?

share|improve this question
    
You can use Select and Length to do it. –  Spawn1701D May 20 '13 at 18:33
    
...or Count[]. –  J. M. May 20 '13 at 18:36
    
...and of course there are no solutions in positive integers for this particular equation $x^2+2y^2=5$. For general equations of this type, please search for Cornacchia in this or other websites. –  KennyColnago May 20 '13 at 18:40
1  
One could program the Cornacchia algorithm if need be, but for the lazy: Reduce[x^2 + 2 y^2 == 5, {x, y}, Integers] –  J. M. May 20 '13 at 18:50
add comment

2 Answers

up vote 1 down vote accepted

A couple of brute force methods:

Example: $2x^2-3y^2=15$ -- there are 7 solutions up to $x,y=5\,000$.

First method.

(lim = 5*10^3;
 c = Tuples[Range[lim], 2];
 (sols = Pick[c, c^2 .{2, -3} - 15, 0])) // Timing
{2.058927,
  {{3, 1}, {9, 7}, {21, 17}, {87, 71}, {207, 169}, {861, 703}, {2049, 1673}}}
Length[sols]
7

Much faster than Select:

(c = Tuples[Range[lim], 2];
 Select[c, 2 #[[1]]^2 - 3 #[[2]]^2 == 15 &]) // Timing
{103.844565,
  {{3, 1}, {9, 7}, {21, 17}, {87, 71}, {207, 169}, {861, 703}, {2049, 1673}}}

Second method.

This one iterates through nesting. There's an iterator iterP for when the coefficients have the same sign and one iterN for when they have opposite signs. When they have opposite signs the number of solutions can be infinite, so there is a limit n on the number of iterations. This is both much faster and uses much less memory.

iterN[{x_, y_, a_, b_, eqn_}] := Switch[Sign[eqn],
   -1, {x + 1, y, a, b, eqn + a (2 x + 1)},
   1, {x, y + 1, a, b, eqn + b (2 y + 1)},
   _, Sow[{x, y}]; {x + 1, y, a, b, eqn + a (2 x + 1)}];
iterP[{x_, y_, a_, b_, eqn_}] := Switch[Sign[eqn],
   1, {x - 1, y, a, b, eqn - a (2 x - 1)},
   -1, {x, y + 1, a, b, eqn + b (2 y + 1)},
   _, Sow[{x, y}]; {x - 1, y, a, b, eqn - a (2 x - 1)}];
solveSq[a_Integer, b_Integer, c_Integer, n_Integer] := Reap[If[a b > 0,
   NestWhileList[
    iterP,
    {Ceiling[Sqrt[Abs[c/a]]], 0, Sign[a] a, Sign[a] b, Sign[a] c + a Ceiling[Sqrt[Abs[c/a]]]^2},
    First[#] > 0 &, 1, n],
   Nest[
    iterN,
    {0, 0, Sign[a] a, Sign[a] b, Sign[a] c},
    n]]]

Same example:

solveSq[2, -3, -15, 100000] // Timing
{0.057761,
  {{6606, 5394, 2, -3, -7251},
  {{{3, 1}, {9, 7}, {21, 17}, {87, 71}, {207, 169}, {861, 703}, {2049, 1673}}}}}

These are the solutions up to $x \le 6606$, $y \le 5394$.

Example: $2x^2+3y^2=105\,074$.

First method:

(m = 5*10^3;
 c = Tuples[Range[lim], 2];
 Pick[c, c^2 .{2, 3} - 105074, 0]) // Timing
{1.790960,
  {{101, 168}, {229, 8}}}

Second method:

solveSq[2, 3, -105074, 100000][[-1]] // Timing
{0.002724,
  {{{229, 8}, {101, 168}}}}
share|improve this answer
    
Thanks @Michael E2, I hadn't come across Pick before. That's useful to know that it can be so much faster than Select –  Jonathan Shock May 22 '13 at 1:12
1  
@JonathanShock Yes, if you can make the second list quickly, then Pick will be fast. In this case it's the internal, vectorized arithmetic that is much faster than testing each case, which is what Select does. –  Michael E2 May 22 '13 at 1:26
add comment

You can work this out from the comments to your question, but it's good also to have a full answer.

You can select elements of a list which fulfill a particular requirement using the following syntax as an example:

m = 3;
n = 2;
c = Tuples[Range[m], n]
Select[c, #[[1]]^2 + 2 #[[2]]^2 == 9 &]

The second argument in Select is your criterion and essentially it is a pure function which takes the first and second elements of every entry in your tuple and tests to see if they fulfill the test (in this case it is ==9, as ==5 gives you nothing). You can then ask for the Length[] of this list.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.