<< problem 501 - Eight Divisors | Divisibility of factorials - problem 549 >> |

# Problem 504: Square on the Inside

(see projecteuler.net/problem=504)

Let ABCD be a quadrilateral whose vertices are lattice points lying on the coordinate axes as follows:

A(a, 0), B(0, b), C(-c, 0), D(0, -d), where 1 <= a, b, c, d <= m and a, b, c, d, m are integers.

It can be shown that for m = 4 there are exactly 256 valid ways to construct ABCD.

Of these 256 quadrilaterals, 42 of them strictly contain a square number of lattice points.

How many quadrilaterals ABCD strictly contain a square number of lattice points for m = 100?

# My Algorithm

Pick's theorem states A = i + b/2 - 1 (see en.wikipedia.org/wiki/Pick's_theorem):

if a simple polygon has i interior integer points and b boundary integer points, then A must be a perfect square to match the conditions of the problem.

So all I have to do is finding i and b - that's what my function `countLatticePoints`

is for.

The bounding box that enclosed the quadrilateral and is parallel to the x- and y-axis has an area of

ab + bc + cd + ad = a(b + d) + c(b + d) = (a + c)(b + d)

The four triangles cover only half the area of the bounding box. My variable `inside`

equals i of the formula:

i = dfrac{1}{2} (a + c)(b + d)

The boundary integer points can be subdivided into two groups:

- 4 corners A(a, 0), B(0, b), C(-c, 0), D(0, -d)

- the edges between A and B, B and C, C and D, D and A

An edge between between two points X(x, 0) and Y(0, y) covers the same number of integer lattice points as an edge between X_2(0,0) and Y_(x,y).

And there are gcd(x,y) + 1 such integer lattice points: the origin plus all multiples of k * \left( dfrac{x}{gcd(x,y)}, dfrac{y}{gcd(x,y)} \right).

However, each corner appears twice: A in AB and DA, B in AB and BC, ...

Using just gcd(x,y) (without minus one) solves this problem.

Unfortunately, my results of `countLatticePoints`

were off by 4:

the corners A, B, C and D are not part of the original bounding-box and must be subtracted, too.

## Alternative Approaches

There are many symmetries. For example I could assume a <= b and reduce the search space accordingly.

## Note

I set up a few caches to speed up the computations considerably:

- the same GCDs have to by computed quite many times, therefore I store them in `pointsOnEdge`

- `sqrt`

is pretty slow and a simple container `isSquare`

stores a `true`

flag for each square up to the largest area of an quadrilateral (100,100,100,100) → about 20000.

# My code

… was written in C++11 and can be compiled with G++, Clang++, Visual C++. You can download it, too.

#include <iostream>
#include <vector>
// cache gcd(a,b)

std::vector<std::vector<unsigned int>> pointsOnEdge;
// greatest common divisor

template <typename T>
T gcd(T a, T b)
{
while (a != 0)
{
T c = a;
a = b % a;
b = c;
}
return b;
}
// count lattice points with Pick's theorem
// A(a, 0), B(b, 0), C(-c, 0), D(-d, 0)

unsigned int countLatticePoints(unsigned int a, unsigned int b, unsigned int c, unsigned int d)
{
// Pick's theorem: A = i + b/2 - 1
// where i is the area inside the quadrilateral
// and b is the number of lattice points on the boundary
// https://en.wikipedia.org/wiki/Pick%27s_theorem
// total area of the bounding box
auto rectangle = (a + c) * (b + d);
// only half of the area belongs to the quadrilateral
auto inside = rectangle / 2;
// points on the edges of the quadrilateral
auto boundary = pointsOnEdge[a][b] + pointsOnEdge[b][c] + pointsOnEdge[c][d] + pointsOnEdge[a][d];
// the four corners are not inside the bounding box (they are on its edges)
boundary -= 4;
// Pick's formula
return inside - (boundary / 2) - 1;
}
int main()
{
unsigned int limit = 100;
std::cin >> limit;
// precompute number of lattice points lying exactly on the sides/edges
pointsOnEdge.resize(limit + 1);
for (auto& x : pointsOnEdge)
x.resize(limit + 1);
for (unsigned int a = 1; a <= limit; a++)
for (unsigned int b = 1; b <= limit; b++)
// note: only one endpoint is included for each edge
// e.g. a=4 b=2 gcd(4,2)=2
// => has actually 3 lattice points:
// both endpoints plus point (3,1)
// but each endpoint is part of two edges therefore it would be counted twice
// that means that keeping only one endpoint is okay and much simpler
pointsOnEdge[a][b] = pointsOnEdge[b][a] = gcd(a, b);
// precompute square numbers
// note: std::vector<bool> would be sufficient but std::vector<char> is about 20% faster
std::vector<unsigned char> isSquare(countLatticePoints(limit, limit, limit, limit) + 1, false);
for (size_t i = 1; i*i < isSquare.size(); i++)
isSquare[i*i] = true;
// plain enumeration of all combinations, 100^4 = 10^8 iterations
auto count = 0;
for (unsigned int a = 1; a <= limit; a++)
for (unsigned int b = 1; b <= limit; b++)
for (unsigned int c = 1; c <= limit; c++)
for (unsigned int d = 1; d <= limit; d++)
// count only perfect squares
if (isSquare[countLatticePoints(a, b, c, d)])
count++;
// display result
std::cout << count << std::endl;
return 0;
}

This solution contains 11 empty lines, 25 comments and 2 preprocessor commands.

# Interactive test

You can submit your own input to my program and it will be instantly processed at my server:

This is equivalent to`echo 4 | ./504`

Output:

*Note:* the original problem's input `100`

__cannot__ be entered

because just copying results is a soft skill reserved for idiots.

*(this interactive test is still under development, computations will be aborted after one second)*

# Benchmark

The correct solution to the original Project Euler problem was found in 0.22 seconds on a Intel® Core™ i7-2600K CPU @ 3.40GHz.

(compiled for x86_64 / Linux, GCC flags: `-O3 -march=native -fno-exceptions -fno-rtti -std=c++11 -DORIGINAL`

)

See here for a comparison of all solutions.

**Note:** interactive tests run on a weaker (=slower) computer. Some interactive tests are compiled without `-DORIGINAL`

.

# Changelog

August 11, 2017 submitted solution

August 11, 2017 added comments

# Difficulty

Project Euler ranks this problem at **15%** (out of 100%).

# Links

projecteuler.net/thread=504 - **the** best forum on the subject (*note:* you have to submit the correct solution first)

# Heatmap

green problems solve the original Project Euler problem and have a perfect score of 100% at Hackerrank, too.

yellow problems score less than 100% at Hackerrank (but still solve the original problem).

gray problems are already solved but I haven't published my solution yet.

blue problems are solved and there wasn't a Hackerrank version of it at the time I solved it or I didn't care about it because it differed too much.

red problems are solved but exceed the time limit of one minute or the memory limit of 256 MByte.

*Please click on a problem's number to open my solution to that problem:*

1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | 19 | 20 | 21 | 22 | 23 | 24 | 25 |

26 | 27 | 28 | 29 | 30 | 31 | 32 | 33 | 34 | 35 | 36 | 37 | 38 | 39 | 40 | 41 | 42 | 43 | 44 | 45 | 46 | 47 | 48 | 49 | 50 |

51 | 52 | 53 | 54 | 55 | 56 | 57 | 58 | 59 | 60 | 61 | 62 | 63 | 64 | 65 | 66 | 67 | 68 | 69 | 70 | 71 | 72 | 73 | 74 | 75 |

76 | 77 | 78 | 79 | 80 | 81 | 82 | 83 | 84 | 85 | 86 | 87 | 88 | 89 | 90 | 91 | 92 | 93 | 94 | 95 | 96 | 97 | 98 | 99 | 100 |

101 | 102 | 103 | 104 | 105 | 106 | 107 | 108 | 109 | 110 | 111 | 112 | 113 | 114 | 115 | 116 | 117 | 118 | 119 | 120 | 121 | 122 | 123 | 124 | 125 |

126 | 127 | 128 | 129 | 130 | 131 | 132 | 133 | 134 | 135 | 136 | 137 | 138 | 139 | 140 | 141 | 142 | 143 | 144 | 145 | 146 | 147 | 148 | 149 | 150 |

151 | 152 | 153 | 154 | 155 | 156 | 157 | 158 | 159 | 160 | 161 | 162 | 163 | 164 | 165 | 166 | 167 | 168 | 169 | 170 | 171 | 172 | 173 | 174 | 175 |

176 | 177 | 178 | 179 | 180 | 181 | 182 | 183 | 184 | 185 | 186 | 187 | 188 | 189 | 190 | 191 | 192 | 193 | 194 | 195 | 196 | 197 | 198 | 199 | 200 |

201 | 202 | 203 | 204 | 205 | 206 | 207 | 208 | 209 | 210 | 211 | 212 | 213 | 214 | 215 | 216 | 217 | 218 | 219 | 220 | 221 | 222 | 223 | 224 | 225 |

226 | 227 | 228 | 229 | 230 | 231 | 232 | 233 | 234 | 235 | 236 | 237 | 238 | 239 | 240 | 241 | 242 | 243 | 244 | 245 | 246 | 247 | 248 | 249 | 250 |

251 | 252 | 253 | 254 | 255 | 256 | 257 | 258 | 259 | 260 | 261 | 262 | 263 | 264 | 265 | 266 | 267 | 268 | 269 | 270 | 271 | 272 | 273 | 274 | 275 |

276 | 277 | 278 | 279 | 280 | 281 | 282 | 283 | 284 | 285 | 286 | 287 | 288 | 289 | 290 | 291 | 292 | 293 | 294 | 295 | 296 | 297 | 298 | 299 | 300 |

301 | 302 | 303 | 304 | 305 | 306 | 307 | 308 | 309 | 310 | 311 | 312 | 313 | 314 | 315 | 316 | 317 | 318 | 319 | 320 | 321 | 322 | 323 | 324 | 325 |

326 | 327 | 328 | 329 | 330 | 331 | 332 | 333 | 334 | 335 | 336 | 337 | 338 | 339 | 340 | 341 | 342 | 343 | 344 | 345 | 346 | 347 | 348 | 349 | 350 |

351 | 352 | 353 | 354 | 355 | 356 | 357 | 358 | 359 | 360 | 361 | 362 | 363 | 364 | 365 | 366 | 367 | 368 | 369 | 370 | 371 | 372 | 373 | 374 | 375 |

376 | 377 | 378 | 379 | 380 | 381 | 382 | 383 | 384 | 385 | 386 | 387 | 388 | 389 | 390 | 391 | 392 | 393 | 394 | 395 | 396 | 397 | 398 | 399 | 400 |

401 | 402 | 403 | 404 | 405 | 406 | 407 | 408 | 409 | 410 | 411 | 412 | 413 | 414 | 415 | 416 | 417 | 418 | 419 | 420 | 421 | 422 | 423 | 424 | 425 |

426 | 427 | 428 | 429 | 430 | 431 | 432 | 433 | 434 | 435 | 436 | 437 | 438 | 439 | 440 | 441 | 442 | 443 | 444 | 445 | 446 | 447 | 448 | 449 | 450 |

451 | 452 | 453 | 454 | 455 | 456 | 457 | 458 | 459 | 460 | 461 | 462 | 463 | 464 | 465 | 466 | 467 | 468 | 469 | 470 | 471 | 472 | 473 | 474 | 475 |

476 | 477 | 478 | 479 | 480 | 481 | 482 | 483 | 484 | 485 | 486 | 487 | 488 | 489 | 490 | 491 | 492 | 493 | 494 | 495 | 496 | 497 | 498 | 499 | 500 |

501 | 502 | 503 | 504 | 505 | 506 | 507 | 508 | 509 | 510 | 511 | 512 | 513 | 514 | 515 | 516 | 517 | 518 | 519 | 520 | 521 | 522 | 523 | 524 | 525 |

526 | 527 | 528 | 529 | 530 | 531 | 532 | 533 | 534 | 535 | 536 | 537 | 538 | 539 | 540 | 541 | 542 | 543 | 544 | 545 | 546 | 547 | 548 | 549 | 550 |

I scored 13,183 points (out of 15300 possible points, top rank was 17 out of ≈60000 in August 2017) at Hackerrank's Project Euler+.

Look at my progress and performance pages to get more details.

My username at Project Euler is

**stephanbrumme**while it's stbrumme at Hackerrank.

# Copyright

I hope you enjoy my code and learn something - or give me feedback how I can improve my solutions.

All of my solutions can be used for any purpose and I am in no way liable for any damages caused.

You can even remove my name and claim it's yours. But then you shall burn in hell.

The problems and most of the problems' images were created by Project Euler.

Thanks for all their endless effort.

<< problem 501 - Eight Divisors | Divisibility of factorials - problem 549 >> |