# Diophantine equation

# Problem 066: Diophantine equation

**Description:**

Consider quadratic Diophantine equations of the form:

`x^2 – Dy^2 = 1`

For example, when `D=13`

, the minimal solution in `x`

is `649^2 – 13×180^2 = 1`

.

It can be assumed that there are no solutions in positive integers when D is square.

By finding minimal solutions in `x`

for `D = {2, 3, 5, 6, 7}`

, we obtain the following:

```
3^2 – 2×2^2 = 1
2^2 – 3×1^2 = 1
9^2 – 5×4^2 = 1
5^2 – 6×2^2 = 1
8^2 – 7×3^2 = 1
```

Hence, by considering minimal solutions in `x`

for `D <= 7`

, the largest `x`

is obtained when `D=5`

.

Find the value of `D <= 1000`

in minimal solutions of `x`

for which the largest value of `x`

is obtained.

**Solution:**

XBXX ################

ZZZ ################

N

OOO ################

################

################

################ v 9::p4+9\g5+9::p1+9\g2+9: <

################ >+6g\9+5p:1-\!v0 <

################ v $_::9+1g\9+0p:^

v _v#-*"}"8p13:+1g13$< # <

>"@ ":*:**21p331p013p88 +>:8+0\8v >11g.@ >+8p:1-\!| >-22g/22p35*^

v p02:p010g13<$<_^#!:-1p< ^9\g2+9::<# ^*:g21g13p21-g21*g 22<

vp01g03_v#-g01:_v#- <^ <v_v#!:-1p0\0+8:p1\0 <^*53p11g13$_ ^ g

>:30p:20 g\/+2/: 30g^ >88 +>:8+0\5p:8+0\4p:8+^ > `^ >$$ $v2

> >$30g:41p:*31g-| >$164*1p164*4p012p122pv 0 >$1+:9-8- !| 03

^p02:p010g13p13+1g13 <v ># 3# 5# *# #<v ^_^#:-g8+9\g2+9::< <$

v*53p23/g22+g21g14 < $ > ^

>:::9+0g\9+1g32g*+13g+:v v 53< |!\-1:g42$< >::9+7g\9+9g- |

|\-1:p2+9\%g12 p31/g12 < # >$ 014p35 *>:9+0\ v|!-9g43-1p7g43%g12p41/g12 :+g<

>$35*:::9+4g\9+5g32g*+13g+v $ >*>2 4p35* >::24g+# 6-34p9+2g24g9+2gv>1v7

>:::9+4g\9+5g32g*+13g+v <|!\-1:g42$< |\-1:p7< ^ # < -g

|\-1:p6+9\%g12p31/g12:<> # :9+0\v| !-9 g43-1p9g43%g12p41/g12:+g9g43<v1*<|\<4

>$114p35* ^ >*>2 4p35*#9> ::2 4g+6-34p9+6g24g9+6g31g**14g+^> 4g +3^

^ _^#\-1:p< ^53$< ^ $<

**Output:**

**Stack:**

*(0)*

**Explanation:**

Okay I admit this one took me five days to complete (with two days pause in between, because of I kinda got frustrated).

I needed eight numbers that were too big for int64 so I encoded them in base-67108864 (`2^26`

).
The reason for this specific number (and I had to fail first to see the problem with bigger bases) is that the biggest calculation I do is `D_0 * D_0 * D * 1`

.
Which is maximally `2^26 * 2^26 * 2^10`

which fits *barely* in an signed 64-bit integer.

Also I needed to first write code to multiply two numbers, both in base-67108864 and both bigger than `2^63`

. Let me tell you that was no fun and long-addition is far easier to implement than long-multiplication.
Especially after first I did everything wrong and then had to redo it :/

The first running version of this program was full of debug statements (even more than normally) and had a size of 200x60. (You can look at it, it's the `Problem-066 (annotated)`

file).
But after that I managed to shrink it quite a bit :D I even *(barely)* managed to fit it in the 80x25 restriction.
I have the feeling I've gotten pretty good at compacting my programs...

But back to the interesting stuff, how did I solve this one:

I can' really explain everything in detail here, so I give you a few useful links:

- https://en.wikipedia.org/wiki/Diophantine_equation
- https://en.wikipedia.org/wiki/Pell's_equation
- https://en.wikipedia.org/wiki/Generalized_continued_fraction
- http://www.maths.surrey.ac.uk/hosted-sites/R.Knott/Fibonacci/cfINTRO.html#section9.4

If you want to have the `(x|y)`

solution for a number D:

- First you calculate the continuous fraction of
`sqrt(D)`

- For the continuous fraction you calculate the convergents
`hi / ki`

(read the link about Pell's equation) - Now you just test if
`hi`

and`ki`

are solutions, if not go on to the next convergent pair

In the end the algorithm is really not that complex (around 30 lines in C#) but all the numbers get so big - and you need to multiply and add the already big numbers together so you get even bigger immediate values. So on a last note I could say ... I wished the numbers in befunge were unlimited.

Interpreter steps: |
262 481 767 |

Execution time (BefunExec): |
55s (4.70 MHz) |

Program size: |
80 x 25 (fully conform befunge-93) |

Solution: |
661 |

Solved at: |
2015-07-14 |