Diophantine equation

2015-07-14

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:

vYYY     ################
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$<                              ^ $<

Start

Pause

Reset

Download

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:

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