Mathematica Asked on March 14, 2021
The overal question I am trying to answer is:
For what $(x,y)$, which are positive integers, is the following number a perfect square number?
$$9 left(x^3 (y-2)^2+3 x^2 (y-2)-2 x (y-45) (y-2)+7 (y-1)^2right)tag1$$
Now, I am using the following code:
ParallelTable[
If[IntegerQ[
Sqrt[9 (x^3 (y - 2)^2 + 3 x^2 (y - 2) - 2 x (y - 45) (y - 2) +
7 (y - 1)^2)]], {x, y}, Nothing], {x, 1, 100}, {y, 1,
100}] /. {} -> Nothing
But that is a bit slow for bigger values of $x$ and $y$. Is there a faster way to test when the number is a perfect square.
You can do
FindInstance[n^2 == 9 (x^3 (y - 2)^2 + 3 x^2 (y - 2) - 2 x (y - 45) (y - 2) + 7 (y - 1)^2),
{x, y, n}, PositiveIntegers]
(* {{x -> 3, y -> 5, n -> 102}} *)
to find an exemplary instance.
If you want all solutions up to $x,yle s$, you can do
With[{s = 20},
Solve[{n^2 == 9 (x^3 (y - 2)^2 + 3 x^2 (y - 2) - 2 x (y - 45) (y - 2) + 7 (y - 1)^2),
1 <= x <= s && 1 <= y <= s}, {x, y, n}, PositiveIntegers]]
(* {{x -> 1, y -> 8, n -> 87},
{x -> 3, y -> 5, n -> 102},
{x -> 3, y -> 8, n -> 159},
{x -> 9, y -> 8, n -> 537}} *)
Faster: using the squareness test of this answer and a Sow
/Reap
combination, and eliminating the prefactor of 9 (see @mikado's comment):
sQ[n_] := FractionalPart@Sqrt[n + 0``1] == 0
With[{s = 1000},
Reap[Do[If[
sQ[x^3 (y - 2)^2 + 3 x^2 (y - 2) - 2 x (y - 45) (y - 2) + 7 (y - 1)^2],
Sow[{x, y}]], {x, s}, {y, s}]][[2, 1]]]
(* {{1, 8}, {1, 128}, {3, 5}, {3, 8}, {9, 8}, {11, 1}, {47, 8}} *)
Further, using the parallelization trick of this Q&A,
SetSharedFunction[ParallelSow];
ParallelSow[expr_] := Sow[expr]
With[{s = 10^4},
Reap[ParallelDo[
If[sQ[x^3 (y - 2)^2 + 3 x^2 (y - 2) - 2 x (y - 45) (y - 2) + 7 (y - 1)^2],
ParallelSow[{x, y}]], {x, s}, {y, s}]][[2, 1]]]
(* {{1, 8}, {1, 128}, {1, 1288}, {3, 5}, {3, 8}, {9, 8}, {11, 1}, {47, 8}} *)
Other than that, this kind of calculation is done much more efficiently in a low-level language like C. Here's my attempt to use pure C with 128-bit integers (because for $x=y=10^6$ we overflow 64-bit integers), going up to $s=10^6$ in about two hours:
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <stdint.h>
#include <inttypes.h>
typedef __int128 myint;
static myint
compute_isqrt(const myint x)
{
myint r = sqrt(x);
while (r*r <= x) {
if (r*r == x)
return r;
r++;
}
return -1;
}
static myint
isqrt(const myint x)
{
if (x < 0)
return -1;
switch(x & 0xf) {
case 0:
case 1:
case 4:
case 9:
return compute_isqrt(x);
default:
return -1;
}
}
#define M 1000000
int main() {
for (myint x=1; x<=M; x++)
for (myint y=1; y<=M; y++) {
myint z = x*x*x*(y-2)*(y-2)+3*x*x*(y-2)-2*x*(y-45)*(y-2)+7*(y-1)*(y-1);
myint n = isqrt(z);
if (n >= 0) {
printf("%" PRId64 " %" PRId64 " %" PRId64 "n",
(int64_t)x, (int64_t)y, (int64_t)n);
}
}
return EXIT_SUCCESS;
}
Save as perfectsquare.c
, compile with
gcc -Wall -O3 perfectsquare.c -o perfectsquare
and run with
time ./perfectsquare
Here are all the solutions ${x,y,n/3}$ up to $s=10^6$:
1 8 29
1 128 329
1 1288 3171
1 13168 32271
1 126848 310729
3 5 34
3 8 53
3 42680 225859
3 61733 326678
3 476261 2520154
3 688856 3645101
9 8 179
11 1 0
47 8 1949
15577 8 11664979
Correct answer by Roman on March 14, 2021
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP