text:
High Precision Metallic Means
Background
The metallic means, starting with the famous golden mean, are defined for every natural number (positive integer), and each one is an irrational constant (it has an infinite non-recurring decimal expansion).
For a natural number , the  metallic mean is the root of a quadratic equation
The roots are always

but the metallic mean is usually given as the positive root. So for this question it will be defined by:

For  the result is the famous golden ratio:

Challenge
Your code should take 2 inputs: n and p (the order is not important as long as it is consistent)

n is a natural number indicating which metallic mean
p is a natural number indicating how many decimal places of precision

Your code should output the nth metallic mean to p decimal places precision.
Validity
Your code is valid if it works for values of n and p from 1 to 65,535.
You must output a decimal in the form
digit(s).digit(s) (without spaces)
For example, the golden mean to 9 decimal places is
1.618033988
Display the last digit without rounding, as it would appear in a longer decimal expansion. The next digit in the golden mean is a 7, but the final 8 in the example should not be rounded up to a 9.
The number of decimal digits must be p, which means any trailing zeroes must also be included.
Answers of the form

are not valid - you must use a decimal expansion.
You may output up to 1 leading newline and up to 1 trailing newline. You may not output any spaces, or any other characters besides digits and the single point/full stop/period.
Score
This is standard code golf: your score is the number of bytes in your code.

Leaderboard
(Using Martin's leaderboard snippet)

var QUESTION_ID=52493;function answersUrl(e){return"http://api.stackexchange.com/2.2/questions/"+QUESTION_ID+"/answers?page="+e+"&pagesize=100&order=desc&sort=creation&site=codegolf&filter="+ANSWER_FILTER}function getAnswers(){jQuery.ajax({url:answersUrl(page++),method:"get",dataType:"jsonp",crossDomain:!0,success:function(e){answers.push.apply(answers,e.items),e.has_more?getAnswers():process()}})}function shouldHaveHeading(e){var a=!1,r=e.body_markdown.split("\n");try{a|=/^#/.test(e.body_markdown),a|=["-","="].indexOf(r[1][0])>-1,a&=LANGUAGE_REG.test(e.body_markdown)}catch(n){}return a}function shouldHaveScore(e){var a=!1;try{a|=SIZE_REG.test(e.body_markdown.split("\n")[0])}catch(r){}return a}function getAuthorName(e){return e.owner.display_name}function process(){answers=answers.filter(shouldHaveScore).filter(shouldHaveHeading),answers.sort(function(e,a){var r=+(e.body_markdown.split("\n")[0].match(SIZE_REG)||[1/0])[0],n=+(a.body_markdown.split("\n")[0].match(SIZE_REG)||[1/0])[0];return r-n});var e={},a=1,r=null,n=1;answers.forEach(function(s){var t=s.body_markdown.split("\n")[0],o=jQuery("#answer-template").html(),l=(t.match(NUMBER_REG)[0],(t.match(SIZE_REG)||[0])[0]),c=t.match(LANGUAGE_REG)[1],i=getAuthorName(s);l!=r&&(n=a),r=l,++a,o=o.replace("{{PLACE}}",n+".").replace("{{NAME}}",i).replace("{{LANGUAGE}}",c).replace("{{SIZE}}",l).replace("{{LINK}}",s.share_link),o=jQuery(o),jQuery("#answers").append(o),e[c]=e[c]||{lang:c,user:i,size:l,link:s.share_link}});var s=[];for(var t in e)e.hasOwnProperty(t)&&s.push(e[t]);s.sort(function(e,a){return e.lang>a.lang?1:e.lang<a.lang?-1:0});for(var o=0;o<s.length;++o){var l=jQuery("#language-template").html(),t=s[o];l=l.replace("{{LANGUAGE}}",t.lang).replace("{{NAME}}",t.user).replace("{{SIZE}}",t.size).replace("{{LINK}}",t.link),l=jQuery(l),jQuery("#languages").append(l)}}var ANSWER_FILTER="!t)IWYnsLAZle2tQ3KqrVveCRJfxcRLe",answers=[],page=1;getAnswers();var SIZE_REG=/\d+(?=[^\d&]*(?:<(?:s>[^&]*<\/s>|[^&]+>)[^\d&]*)*$)/,NUMBER_REG=/\d+/,LANGUAGE_REG=/^#*\s*([^,]+)/;
body{text-align:left!important}#answer-list,#language-list{padding:10px;width:290px;float:left}table thead{font-weight:700}table td{padding:5px}
<link rel="stylesheet" type="text/css" href="//cdn.sstatic.net/codegolf/all.css?v=83c949450c8b">

Leaderboard

	Author	Language	Size

Winners by Language

Language	User	Score

{{PLACE}}	{{NAME}}	{{LANGUAGE}}	{{SIZE}}	Link

{{LANGUAGE}}	{{NAME}}	{{SIZE}}	Link

dc, 12
?kdd*4+v+2/p

? Push n and p onto the stack
k set precision to p
dd duplicate n twice (total three copies)
* multiply n*n
4+ add 4
v take square root
+ add n (last copy on stack)
2/ divide by 2
p print

Testcase:
$ dc -f metalmean.dc <<< "1 9"
1.618033988
$

Right tool for the job.

@Dennis its got to be the first time CJam is nearly 3 times as long as something else ;-)

R, 116 bytes
library(Rmpfr);s=scan();n=mpfr(s[1],1e6);r=(n+(4+n^2)^.5)/2;t=toString(format(r,s[2]+2));cat(substr(t,1,nchar(t)-1))

This reads two integers from STDIN and prints the result to STDOUT. You can try it online.
Ungolfed + explanation:
# Import the Rmpfr library for arbitrary precision floating point arithmetic
library(Rmpfr)

# Read two integers from STDIN
s <- scan()

# Set n equal to the first input as an mpfr object with 1e6 bits of precision
n <- mpfr(s[1], 1e6)

# Compute the result using the basic formula
r <- (n + sqrt(4 + n^2)) / 2

# Get the rounded string representation of r with 1 more digit than necessary
t <- toString(format(r, s[2] + 2))

# Print the result with p unrounded digits
cat(substr(t, 1, nchar(t) - 1))

If you don't have the Rmpfr library installed, you can install.packages("Rmpfr") and all of your dreams will come true.

Python 2, 92 Bytes
As I am now looking at the answers, it looks like the CJam answer uses the same basic method as this. It calculates the answer for n*10**p and then adds in the decimal point. It is incredibly inefficient due to the way it calculates the integer part of the square root (just adding 1 until it gets there).
n,p=input()
e=10**p;r=0
while(n*n+4)*e*e>r*r:r+=1
s=str((n*e+r-1)/2);print s[:-p]+'.'+s[-p:]

Mathematica, 50 bytes
SetAccuracy[Floor[(#+Sqrt[4+#^2])/2,10^-#2],#2+1]&

Defines an anonymous function that takes n and p in order.  I use Floor to prevent rounding with SetAccuracy, which I need in order to get decimal output.

@Arcinde I can't use machine precision numbers unfortunately, since they wouldn't be able to handle p>15.

CJam, 35 bytes
1'el+~1$*_2#2$2#4*+mQ+2/1$md@+s0'.t

Reads p first, then n.
Try it online in the CJam interpreter.
How it works
We simply compute the formula from the question for n × 10p, get the integer and fractional part of the result divided by 10p, pad the fractional part with leading zeroes to obtain p digits and print the parts separated by a dot.
1'e  e# Push 1 and 'e'.
l+   e# Read a line from STDIN and prepend the 'e'.
~    e# Evaluate. This pushes 10**p (e.g., 1e3 -> 1000) and n.
1$*  e# Copy 10**p and multiply it with n.
_2#  e# Copy n * 10**p and square it.
2$   e# Copy 10**p.
2#4* e# Square and multiply by 4.
+    e# Add (n * 10**p)**2 and 4 * 10**2p.
mQ   e# Push the integer part of the square root.
+2/  e# Add to n * 10**p and divide by 2.
1$md e# Perform modular division by 10**p.
@+s  e# Add 10**p to the fractional part and convert to string.
0'.t e# Replace the first character ('1') by a dot.

PHP, 85 78 bytes
echo bcdiv(bcadd($n=$argv[bcscale($argv[2])],bcsqrt(bcadd(4,bcpow($n,2)))),2);

It uses the BC Math mathematical extension which, on some systems, could not be available. It needs to be included on the compilation time by specifying the --enable-bcmath command line option. It is always available on Windows and it seems it is included in the PHP version bundled with OSX too.
Update:
I applied all the hacks suggested by @blackhole in their comments (thank you!) then I squeezed the initialization of $n into its first use (3 more bytes saved) and now the code fits in a single line in the code box above.

@Blackhole. 85, indeed. I have probably read 86 (did a slightly larger selection) and wrote 68 by mistake. Fixed now.

No problem :). You can have 1 byte less by the way: remove the parenthesis around the echo, just leave a space after it.

And since you expect bcscale to return true, you can use $n=$argv[bcscale($argv[2])]; and save 2 more bytes.

That's a nice hack.

Code dirtiness is an art :P. Oh, the last one: bcpow($n,2) instead of bcmul($n,$n) saves you 1 byte.

J, 27 Bytes
4 :'}:":!.(2+x)-:y+%:4+*:y'

Explanation:
4 :'                      '   | Define an explicit dyad
                       *:y    | Square y
                     4+       | Add 4
                   %:         | Square root
                 y+           | Add y
               -:             | Half
      ":!.(2+x)               | Set print precision to 2+x
    }:                        | Remove last digit, to fix rounding

Call it like this:
    9 (4 :'}:":!.(2+x)-:y+%:4+*:y') 1
1.618033988

Another, slightly cooler solution:
4 :'}:":!.(2+x){.>{:p._1,1,~-y'

Which calculates the roots of the polynomial x^2 - nx - 1. Unfortunately, the way J formats the result makes retreving the desired root slightly longer.


text:
Print an ascii spiral in O(log n) memory
You may write a program or function that receives an odd, positive integer n, where n >= 3, as either a function argument, command line arguments, or on STDIN (or equivalent for your system), and prints to STDOUT (or system equivalent) an ASCII spiral that spins inward clockwise where the top edge is exactly n characters long. The first right edge should be n+1 characters long, obviously. For instance,
Input:
11
Output:
***********
          *
********* *
*       * *
* ***** * *
* *   * * *
* * * * * *
* * *** * *
* *     * *
* ******* *
*         *
***********

The catches:

Your program must use no more than O(log n) memory.
Your program may only print the characters * (ASCII 42),  (ASCII 32), <CR> (ASCII 13) and <LF> (ASCII 10).
Your program must print the string, not return it from the function.
The Big-O restriction is only on memory, there are no restrictions on runtime.
A trailing newline is optional.
If your language does not support large integer types, you do not have to support higher than what it does support, but you may not use this as a trick to say "oh, well, I don't have to support above X so I can just make a huge array the maximum size every time"

Standard loopholes are banned, as usual.

I don't believe this is possible. One cannot store the input n in O(1) memory.

@xnor "O(1) constitutes a constant memory usage. So amount of input is inconsequential" - If input n fits in integer, then I'm sure it can be coded into a constant memory usage.

Storing the input n takes log n bits. As n gets larger, so does the space needed to store it. Are you perhaps saying to do this with a limited number of variables?

Or, alternatively, is there a limit on n?

I think he's saying you can't store the entire output at once, then just print it all at once because that will get larger. You probably have to recursively print it.

@Sp3000 @durron597 I think capping n opens another can of worms. One could just make a table hardcoding every possible spiral, though that's not golfy. Or, store the whole spiral in memory, since its size it bounded by a constant.

@xnor I think O(log n) memory is the only way to do this. Feh

@durron597 Good idea, O(log n) should do what you intend.

@xnor On the other hand, if n is not capped, languages without standard arbitrary size int types will be pretty much out of luck. Well, I guess I could include some kind of big int implementation in a C solution, but it would certainly put it at a big disadvantage against languages who are typically about equally verbose.

Do we assume a trans-dichotomous machine model?

@FUZxxl that's fine.

C, 125 121 bytes
Golfed version This has no variable k. The variable k is used in the ungolfed version just to aid readability. Also for loop conditionals are rearranged and one set of unnecessary {} removed. Another set of {} can be removed by migrating puts("") inside the brackets of the j loop in the initialization position, but this would mean a newline at the beginning of the output, so I haven't done it.
f(n){int i,j;n/=2;for(i=-n-2;i++-n-1;){if(i){for(j=-n-1;j++-n;)putchar(32+10*(n+(j*j<i*i?i:j+(i!=j|i>0))&1));puts("");}}}

Prints an n wide by n+1 high spiral like the example.
Explanation
Basically I halve the value of n (rounding down) and run two loops: an outer one i from -n/2-1 to n/2+1 to print the rows (i=0 is suppressed so we get n+1 rows) and an inner one j from (-n/2 to n/2 to print the characters.) We use expression & 1 to print stripes, and the condition j*j<i*i to decide whether to print vertical or horizontal stripes (vertical at the sides where absolute magnitude of i is larger, and horizontal at the top and bottom.) An adjustment +n is required to help with the correct termination depending on whether n/2 is odd or even.
k is normally 1, and provides an adjustment for the fact that the absolute values of i range from 1 to n/2+1 while the absolute values of j range from 0 to n/2. If k was always 1 we would get concentric rectangles, but it is inverted to 0 when i==j&i<=0 so that a diagonal row of cells is inverted, producing the spiral.
ungolfed in test program
f(n){
  int i,j,k;
  n/=2;
  for(i=-n-1;i<=n+1;i++){
    if(i){
       for(j=-n;j<=n;j++){
           k=i!=j|i>0;
           putchar(32+10*(n+(j*j<i*i?i:k+j)&1));
         }
       puts("");
     }
  }
}

int m;
main(){
  scanf("%d",&m);
  f(m);
}

Output
11
***********
          *
********* *
*       * *
* ***** * *
* *   * * *
* * * * * *
* * *** * *
* *     * *
* ******* *
*         *
***********

9
*********
        *
******* *
*     * *
* *** * *
* * * * *
* *   * *
* ***** *
*       *
*********

3
***
  *
* *
***

1
*
*

Beat me by a bit... +1 this is crazy short!

114 bytes

C, 118 bytes
m,p,x,y,d;f(n){for(m=n++/2;p<n*n;x=p%n-m,y=p++/n-m,d=y==x+1&x<0,y-=y>0,d+=x*x>y*y?x:y,putchar(x>m?10:(d+m)%2?32:42));}

Code before final golfing:
#include <stdio.h>

int m, p, x, y, d;

int f(int n) {
    for (m = n++ / 2; p < n * n; ) {
        x = p % n - m;
        y = p++ / n - m;
        d = y == x + 1 && x < 0;
        y -= y > 0;
        d += x * x > y * y ? x : y;
        if (x > m) {
            putchar(10);
        } else if ((d + m) % 2) {
            putchar(32);
        } else {
            putchar(42);
        }
    }

    return 0;
}

The key observation is that the pattern is almost a series of concentric squares. With a couple of slight wrinkles:

The y-size is one larger than the x-size. This is corrected by subtracting 1 from y for the lower half, which essentially repeats the middle row.
To turn the rectangles into a spiral, the pixels along the y = x + 1 diagonal need to be inverted up to the middle of the shape.

For the rest, the code is simply looping over all positions, calculating the Chebyshev distance from the center for each position, and emitting one of the two characters depending on the distance being even or odd. And emitting a newline for the last position of each line.
Since there are only a few scalar variables, and characters are emitted one by one, memory usage is obviously constant.

Excellent answer, but as you don't initialize p I think you fall foul of http://meta.codegolf.stackexchange.com/q/4939/15599 . I'm also not sure about declaring global variables when submitting a function. Obviously my answer would be 4 bytes shorter if I did this. I've started a meta post http://meta.codegolf.stackexchange.com/q/5532/15599

Yes, it crossed my mind that I should probably initialize p.

C++, 926 bytes
#include<iostream>
#include<string>
#include<math.h>
#define S string
using namespace std;S N(S x,int y){S z="";for(int q=0;q<y;q++){z+=x;}return z;}int main(){int n=0,t=0,g=0,fi=1;cin>>n;int t1[]={0,0,n,0};int t2[]={0,n-2,n-2,1};for(int k=0;k<n+1;k++){if((k>(n-2)/2)&&(k<(n+5)/2)){if(g==0){S d,e;if(!((n+1)%4)){cout<<N("* ",t2[0])<<"  *"<<N(" *",t2[0])<<endl<<N("* ",(n+1)/2)<<endl<<N("* ",t2[0])<<"***"<<N(" *",t2[0])<<endl;t2[2]=n-8-(n-11);t1[2]=n-4-(n-11);t1[0]--;t2[3]--;t1[3]-=2;}else{cout<<N("* ",t1[0])<<"***"<<N(" *",t2[0])<<endl<<N("* ",(n+1)/2)<<endl<<N("* ",t1[0])<<"*  "<<N(" *",t2[0])<<endl;t2[0]--;t1[2]+=2;t2[2]+=6;t1[3]--;t2[1]-=2;t2[3]-=2;}fi=0;}g=5;}else{t=1-t;int*tR;tR=t?t1:t2;cout<<N("* ",tR[0])<<N(t?"*":" ",tR[2])<<N(" *",tR[3])<<endl;if(fi){if(t){t1[0]+=k==0?0:1;t1[2]-=k==0?2:4;t1[3]++;}else{t2[0]++;t2[2]-=4;t2[3]++;}}else{if(t){t1[0]--;t1[2]+=4;t1[3]--;}else{t2[0]--;t2[2]+=4;t2[3]--;}}}}return 0;}

This is not elegant, but it doesn't take up much memory for large n.  Furthermore, there are (almost certainly) about 20 characters that can be further golfed, but I can't stand to look at it anymore.
Short Explanation:
This splits the lines in the spirals into two types: the ones with ****** in the middle, and the ones with \s\s\s\s\s in the middle.  Then it is clear that each line is composed of several "* "s, the middle, and some " *".  Figuring out exactly how many of each thing is simple if you look at the pattern for long enough. The tricky thing was printing the center of the spiral which I basically hard coded using a conditional.  This ended up being useful because the *** and \s\s\s lines switch being odd/even there.
Tests:
Input: 55 (I think the big ones look coolest)
Output:

*******************************************************
                                                      *
***************************************************** *
*                                                   * *
* ************************************************* * *
* *                                               * * *
* * ********************************************* * * *
* * *                                           * * * *
* * * ***************************************** * * * *
* * * *                                       * * * * *
* * * * ************************************* * * * * *
* * * * *                                   * * * * * *
* * * * * ********************************* * * * * * *
* * * * * *                               * * * * * * *
* * * * * * ***************************** * * * * * * *
* * * * * * *                           * * * * * * * *
* * * * * * * ************************* * * * * * * * *
* * * * * * * *                       * * * * * * * * *
* * * * * * * * ********************* * * * * * * * * *
* * * * * * * * *                   * * * * * * * * * *
* * * * * * * * * ***************** * * * * * * * * * *
* * * * * * * * * *               * * * * * * * * * * *
* * * * * * * * * * ************* * * * * * * * * * * *
* * * * * * * * * * *           * * * * * * * * * * * *
* * * * * * * * * * * ********* * * * * * * * * * * * *
* * * * * * * * * * * *       * * * * * * * * * * * * *
* * * * * * * * * * * * ***** * * * * * * * * * * * * *
* * * * * * * * * * * * *   * * * * * * * * * * * * * *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * {-- my program adds a space here btw
* * * * * * * * * * * * * *** * * * * * * * * * * * * *
* * * * * * * * * * * * *     * * * * * * * * * * * * *
* * * * * * * * * * * * ******* * * * * * * * * * * * *
* * * * * * * * * * * *         * * * * * * * * * * * *
* * * * * * * * * * * *********** * * * * * * * * * * *
* * * * * * * * * * *             * * * * * * * * * * *
* * * * * * * * * * *************** * * * * * * * * * *
* * * * * * * * * *                 * * * * * * * * * *
* * * * * * * * * ******************* * * * * * * * * *
* * * * * * * * *                     * * * * * * * * *
* * * * * * * * *********************** * * * * * * * *
* * * * * * * *                         * * * * * * * *
* * * * * * * *************************** * * * * * * *
* * * * * * *                             * * * * * * *
* * * * * * ******************************* * * * * * *
* * * * * *                                 * * * * * *
* * * * * *********************************** * * * * *
* * * * *                                     * * * * *
* * * * *************************************** * * * *
* * * *                                         * * * *
* * * ******************************************* * * *
* * *                                             * * *
* * *********************************************** * *
* *                                                 * *
* *************************************************** *
*                                                     *
*******************************************************

Input: 3
Output:

***
  *
* *
***

Note: I am not a computer scientist/CS student, and I don't know how to prove that this uses O(log n) memory.  I can only work out what to do based on the links in the question.  I would be grateful if someone could confirm/deny if this answer is valid.  My logic for this answer's validity is that it never stores any variable of size based on n except the input itself.  Instead, a for loop that runs n times computes integer values based on n.  There are the same number of those values regardless of the input.
Note2: This doesn't work for n=1 because of my method of dealing with the middle.  This would be easy to fix with conditionals, so if anyone is within a few characters of my answer, I'll fix it ;)
Play with it on ideone.

I believe it's valid, even though this much C++ code on a single line is kind of had to read. ;) Your understanding is correct. You can't use any memory with a size that depends on n. A typical example that does not meet the requirement would be some kind of string/buffer/array that holds a complete line of output.

Since this in the only answer, I have adjusted the question to not require handling n=1, as I don't consider such special casing interesting.

Haskell, 151 bytes
(#)=mod
f n=[[if y<= -(abs$x+1)||y>abs x then r$y#2/=n#2 else r$x#2==n#2|x<-[-n..n]]|y<-[-n-1..n+1],y/=0]
r b|b='*'|1<2=' '
p=putStr.unlines.f.(`div`2)

Usage example:
*Main> p 9
*********
        *
******* *
*     * *
* *** * *
* * * * *
* *   * *
* ***** *
*       *
*********

*Main> p 11
***********
          *
********* *
*       * *
* ***** * *
* *   * * *
* * * * * *
* * *** * *
* *     * *
* ******* *
*         *
***********

Thanks to Haskell's laziness this runs within constant memory. It uses the obvious approach, i.e. looping over y and x and choosing between * and , depending on

if the current position is above or below a diagonal
x resp. y is even or odd
n/2 is even or odd

Common Lisp - 346
(lambda(n &aux(d 0))(tagbody $ #6=(#7=dotimes(i n)#4=(princ"*"))#2=(#7#(i d)#5=(princ" ")#4#)#3=(terpri)#1=(#7#(i d)#4##5#)(when(> n 0)(#7#(i(1- n))#5#)#4#)#2##3#(when(> n 3)#1##4##4#(incf d)(decf n 4)(go $))(go /)@(decf d)(incf n 4)(when(> n 3)#2##5##4##3#)/ #1#(when(> n 0)#4#)(when(> n 1)(#7#(i(- n 2))#5#)#4#)#2##3##1##6#(when(> d 0)(go @))))

Iterative solution with constant memory usage. The above makes heavy uses of #n= and #n# reader variables. Even though there are more direct approaches, here I started with a recursive function and modified it to simulate recursion with goto statements: this is probably unreadable.
Output for all input values from 0 to 59.
Original recursive version, with debugging informations
(note: terpri means newline)
(defun spiral (n &optional (d 0) )
  (flet ((prefix ()
           (format t "~4d~4d | " n d)
           (dotimes (i d)
             (princ "a ")))
         (postfix ()
           (dotimes (i d)
             (princ " b"))))
    (when (= d 0) (prefix))
    (dotimes (i n) (princ "c"))
    (postfix)
    (terpri)

    (prefix)
    (when (> n 0)
      (dotimes (i (1- n)) (princ " "))
      (princ "d"))
    (postfix)
    (terpri)

    (when (> n 3)
      (prefix)
      (princ "**")
      (spiral (- n 4) (1+ d))
      (postfix)
      (princ " f")
      (terpri))

    (prefix)
    (when (> n 0)
      (princ "g"))

    (when (> n 1)
      (dotimes (i (- n 2)) (princ " "))
      (princ "h"))
    (postfix)
    (terpri)

    (prefix)
    (dotimes (i n) (princ "i"))
    ))

For example:
(spiral 8)

   8   0 | cccccccc
   8   0 |        d
   8   0 | **cccc b
   4   1 | a    d b
   4   1 | a ** b b
   0   2 | a a  b b
   0   2 | a a  b b
   0   2 | a a  b f
   4   1 | a g  h b
   4   1 | a iiii f
   8   0 | g      h
   8   0 | iiiiiiii

See also this paste with all results from 0 to 59 (not the same as above, this one is more verbose).
Iterative version, with debugging informations
(defun spiral (n &aux (d 0) )
  (flet ((prefix ()
           (format t "~4d~4d | " n d)
           (dotimes (i d)
             (princ "a ")))
         (postfix ()
           (dotimes (i d)
             (princ " b"))))
    (tagbody
     step-in
       (when (= d 0) (prefix))
       (dotimes (i n) (princ "c"))
       (postfix)
       (terpri)

       (prefix)
       (when (> n 0)
         (dotimes (i (1- n)) (princ " "))
         (princ "d"))
       (postfix)
       (terpri)

       (when (> n 3)
         (prefix)
         (princ "**")

         (incf d)
         (decf n 4)
         (go step-in))

       (go skip)

     step-out
       (decf d)
       (incf n 4)
       (when (> n 3)
         (postfix)
         (princ " f")
         (terpri))

     skip
       (prefix)
       (when (> n 0)
         (princ "g"))

       (when (> n 1)
         (dotimes (i (- n 2)) (princ " "))
         (princ "h"))
       (postfix)
       (terpri)

       (prefix)
       (dotimes (i n) (princ "i"))
       (when(> d 0)(go step-out)))))

Can you explain how this meets the memory restriction? I only see one recursion point, which is good, but can you just go into a bit more detail?

@durron597 Yes, I am working on this. This is currently O(n) because we call the function recursively a number of time proportional to n and the call stack grows accordingly, but in this case, we can simulate recursion with two loops: one with n decreasing and d increasing (until n <= 3), and another one with d decreasing to zero. I don't have much time to work on this right now, but I'll try to update the answer accordingly. Btw, there are more direct ways to print the spiral, but it was fun trying to define it recursively.

CJam, 72 bytes
li_2/:M;)__*{1$mdM-\M-_2$)=2$0<*@_*@_0>-_*e>mQ_M>2*@@+M+2%+'#S+N+N+=o}/;

This is fairly direct conversion of my C solution to CJam. Not as short as you would normally expect from a CJam solution, but this one really suffers from the memory restriction. The common benefits of building up results on the stack that gets dumped automatically at the end, and using fancy list/string operations, all go out the window. This generates and outputs the solution one character at a time. The stack only contains a few integers at runtime, and is empty at the end.
Even though it's not a great display of using a golfing language, it's still considerably shorter than the C code just because the notation is more compact.
Explanation:
li    Get input n.
_2/   Calculate n/2.
:M;   Store it in variable M
)__*  Calculate (n+1)*(n+1), which is the total number of output characters.
      Also keep a copy of n+1 on the stack.
{     Start loop over output character positions.
  1$md  Calculate divmod of position with n+1. This gives y and x of position.
  M-    Subtract M from x.
  \M-   Subtract M from y.
  _     Copy y.
  2$)   Calculate x+1.
  =     Check if y == x+1
  2$0<  Check if x < 0.
  *     Multiply the two check results. This is the result of the flip
        condition for the top-left diagonal to turn the rectangles into a spiral.
  @_*   Calculate x*x.
  @_    Get y to top of stack, and copy it.
  0>-   Subtract 1 from y if it is in the bottom half.
  _*    Calculate y*y.
  e>    Take maximum of x*x and y*y...
  mQ    ... and calculate the square root. This is the absolute value of the
        larger of the two.
  _M>   Check if the value is greater M, which means that this is the
        position of a line end.
  2*    Multiply by 2 so that we can add another condition to it later.
  @     Get result of diagonal flip condition to the stack top.
  @     Get max(x,y) to the top.
  +M+   Add the two, and add M to the whole thing. This value being even/odd
        determines if the output is a # or a space.
  2%    Check if value is odd.
  +     Add to line end condition to get a single ternary condition result.
  '#S+N+N+
        Build string "# \n\n".
  =     Use the condition result to pick the output character out of the string.
  o     Output the character.
}/    End loop over output characters.
;     Pop n+1 value off stack, to leave it empty.

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/

sort: votes


text:
Uniquely separated pixels
For an N by N image, find a set of pixels such that no separation distance is present more than once. That is, if two pixels are separated by a distance d, then they are the only two pixels that are separated by exactly d (using Euclidean distance). Note that d need not be integer.
The challenge is to find a larger such set than anyone else.
Specification
No input is required - for this contest N will be fixed at 619.
(Since people keep asking - there's nothing special about the number 619. It was chosen to be large enough to make an optimal solution unlikely, and small enough to let an N by N image be displayed without Stack Exchange automatically shrinking it. Images can be displayed full size up to 630 by 630, and I decided to go with the largest prime that doesn't exceed that.)
The output is a space separated list of integers.
Each integer in the output represents one of the pixels, numbered in English reading order from 0. For example for N = 3, the locations would be numbered in this order:
0 1 2
3 4 5
6 7 8

You may output progress information during running if you wish, as long as the final scoring output is easily available. You may output to STDOUT or to a file or whatever is easiest for pasting into the Stack Snippet Judge below.
Example
N = 3
Chosen coordinates:
(0,0)
(1,0)
(2,1)

Output:
0 1 5

Winning
The score is the number of locations in the output. Of those valid answers which have the highest score, the earliest to post output with that score wins.
Your code does not need to be deterministic. You may post your best output.

Related areas for research
(Thanks to Abulafia for the Golomb links)
While neither of these is the same as this problem, they are both similar in concept and may give you ideas on how to approach this:

Golomb ruler: the 1 dimensional case.
Golomb rectangle: a 2 dimensional extension of the Golomb ruler. A variant of the NxN (square) case known as Costas array is solved for all N.

Note that the points required for this question are not subject to the same requirements as a Golomb rectangle. A Golomb rectangle extends from the 1 dimensional case by requiring that the vector from each point to each other be unique. This means that there can be two points separated by a distance of 2 horizontally, and also two points separated by a distance of 2 vertically.
For this question, it is the scalar distance that must be unique, so there cannot be both a horizontal and a vertical separation of 2. Every solution to this question will be a Golomb rectangle, but not every Golomb rectangle will be a valid solution to this question.

Upper bounds
Dennis helpfully pointed out in chat that 487 is an upper bound on the score, and gave a proof:

According to my CJam code (619,2m*{2f#:+}%_&,), there are 118800 unique numbers that can be written as the sum of the squares of two integers between 0 and 618 (both inclusive).
  n pixels require n(n-1)/2 unique distances between each other. For n = 488, that gives 118828.

So there are 118,800 possible different lengths between all the potential pixels in the image, and placing 488 black pixels would result in 118,828 lengths, which makes it impossible for them all to be unique.
I'd be very interested to hear if anyone has a proof of a lower upper bound than this.

Leaderboard
(Best answer by each user)

Stack Snippet Judge

canvas=document.getElementById('canvas');ctx=canvas.getContext('2d');information=document.getElementById('information');problemList=document.getElementById('problemList');width=canvas.width;height=canvas.height;area=width*height;function verify(){stop();checkPoints();}function checkPoints(){valid=true;numbers=[];points=[];lengths=[];lengthDetails=[];duplicatedLengths=[];clearCanvas();rawData=pointData.value;splitData=rawData.split(' ');for(i=0;i<splitData.length;i++){datum=splitData[i];if(datum!==''){n=parseInt(datum);if(n>=area||n<0){valid=false;information.innerHTML='Out of range:'+n;break;}if(numbers.indexOf(n)>-1){valid=false;information.innerHTML='Duplicated value:'+n;break;}numbers.push(n);x=n%width;y=Math.floor(n/width);points.push([x,y]);}}if(valid){ctx.fillStyle='black';for(i=0;i<points.length-1;i++){p1=points[i];x1=p1[0];y1=p1[1];drawPoint(x1,y1);for(j=i+1;j<points.length;j++){p2=points[j];x2=p2[0];y2=p2[1];d=distance(x1,y1,x2,y2);duplicate=lengths.indexOf(d);if(duplicate>-1){duplicatedLengths.push(lengthDetails[duplicate]);duplicatedLengths.push([d,[numbers[i],numbers[j]],[x1,y1,x2,y2]]);}lengths.push(d);lengthDetails.push([d,[numbers[i],numbers[j]],[x1,y1,x2,y2]]);}}p=points[points.length-1];x=p[0];y=p[1];drawPoint(x,y);if(duplicatedLengths.length>0){ctx.strokeStyle='red';information.innerHTML='Duplicate lengths shown in red';problemTextList=[];for(i=0;i<duplicatedLengths.length;i++){data=duplicatedLengths[i];length=data[0];numberPair=data[1];coords=data[2];line(coords);specificLineText=''+length+': '+numberPair[0]+' to '+numberPair[1]+' [('+coords[0]+', '+coords[1]+') to ('+coords[2]+', '+coords[3]+')]';if(problemTextList.indexOf(specificLineText)===-1){problemTextList.push(specificLineText);}}problemText=problemTextList.join('
');problemList.innerHTML=problemText;}else{information.innerHTML='Valid: '+points.length+' points';}}}function clearCanvas(){ctx.fillStyle='white';ctx.fillRect(0,0,width,height);}function line(coords){x1=coords[0];y1=coords[1];x2=coords[2];y2=coords[3];ctx.beginPath();ctx.moveTo(x1+0.5,y1+0.5);ctx.lineTo(x2+0.5,y2+0.5);ctx.stroke();}function stop(){information.innerHTML='';problemList.innerHTML='';}function distance(x1,y1,x2,y2){xd=x2-x1;yd=y2-y1;return Math.sqrt(xd*xd+yd*yd);}function drawPoint(x,y){ctx.fillRect(x,y,1,1);}
<input id='pointData' placeholder='Paste data here' onblur='verify()'></input><button>Verify data</button>

<canvas id='canvas' width='619' height='619'></canvas>

I would have loved to see a Piet answer here

@C5H8NNaO4 the competition is open ended - no one is anywhere near an optimal solution so there's plenty of room for new answers...

Since you're offering bounties for both proved upper bound and experimental list of pixels, I assume there is some kind of application to this problem?

@Fatalize not that I'm aware of, but I'd be fascinated to hear of one. The similar problem Costas array has practical applications listed but I've found nothing on this particular problem.

(I have a great promising approach that might just stand a chance of winning, but its taking forever to run and might not be ready complete in time!)

@Will there's a 24 hour grace period after the official bounty deadline, and I'll be waiting for that if there are any answers pending. Hopefully that's enough...

Also I'll be updating the accepted answer in future if the bountied answer is later beaten - the end of the bounty isn't the end of the competition.

I'm spent 6 hours and put 77 points that invalidate less than half the space and I don't know how soon I'll get to leave the computer running longer...  Here's progress http://i.stack.imgur.com/d1qEi.png (blue squares are invalidated because of distance to point, red squares by symmetry between two points and green squares by both constraints).

@Will post the code - maybe someone will be able to run it longer...

I've been looking at this, and I believe n=487 to be the minimal upper bound on pixels. Out of curiosity, will you accept a proof that there is no lesser upper bound for the bounty?

@Mego I had guessed that there would exist an upper bound smaller than 487. If someone posts a valid answer that also contains a valid proof that 487 is minimal, then this proves the bounty cannot be assigned as described, and I will award 300 bounty to the first such proof.

@trichoplax Ok neat, I've been working on a proof, was curious if the bounty would apply, since it stated specifically an upper bound less than 487.

@Mego For any upper bound less than 487, simply a proof that it is an upper bound is sufficient to compete for the bounty. To win the bounty with 487 it will also be required to prove that it is the least possible upper bound (which is what you are working on). In either case the proof must be part of a valid answer to the main question. Hopefully that covers everything for anyone else wishing to compete with you. I didn't want to edit the question just for the bonus bounty as that would bump it to the top of the question list.

I will be very surprised and impressed if either case is proved (that there is a lower upper bound or that there is not). Hence the size of the bounty...

Python 3, 135 136 137
10 6830 20470 47750 370770 148190 306910 373250 267230 354030 30390 361470 118430 58910 197790 348450 381336 21710 183530 305050 2430 1810 365832 99038 381324 39598 262270 365886 341662 15478 9822 365950 44526 58862 24142 381150 31662 237614 118830 380846 7182 113598 306750 11950 373774 111326 272358 64310 43990 200278 381014 165310 254454 12394 382534 87894 6142 750 382478 15982 298326 70142 186478 152126 367166 1162 23426 341074 7306 76210 140770 163410 211106 207962 35282 165266 300178 120106 336110 30958 158 362758 382894 308754 88434 336918 244502 43502 54990 279910 175966 234054 196910 287284 288468 119040 275084 321268 17968 2332 86064 340044 244604 262436 111188 291868 367695 362739 370781 375723 360261 377565 383109 328689 347879 2415 319421 55707 352897 313831 302079 19051 346775 361293 328481 35445 113997 108547 309243 19439 199037 216463 62273 174471 207197 167695 296927

Found using a greedy algorithm which, at each stage, chooses the valid pixel whose set of distances to the chosen pixels overlaps the least with that of other pixels.
Specifically, the scoring is
score(P) = sum(number of pixels with D in its distance set
               for each D in P's distance set)

and the pixel with the lowest score is chosen.
The search is kicked off with the point 10 (i.e. (0, 10)). This part is adjustable, so beginning with different pixels may lead to better or worse results.
It's quite a slow algorithm, so I'm trying to add optimisations/heuristics, and maybe some backtracking. PyPy is recommended for speed.
Anyone trying to come up with an algorithm should test on N = 10, for which I've got 9 (but this took a lot of tweaking and trying different initial points):

Code
from collections import Counter, defaultdict
import sys
import time

N = 619

start_time = time.time()

def norm(p1, p2):
    return (p1//N - p2//N)**2 + (p1%N - p2%N)**2

selected = [10]
selected_dists = {norm(p1, p2) for p1 in selected for p2 in selected if p1 != p2}
pix2dist = {} # {candidate pixel: {distances to chosen}}
dist2pix = defaultdict(set)

for pixel in range(N*N):
    if pixel in selected:
        continue

    dist_list = [norm(pixel, p) for p in selected]
    dist_set = set(dist_list)

    if len(dist_set) != len(dist_list) or dist_set & selected_dists:
        continue

    pix2dist[pixel] = dist_set

    for dist in dist_set:
        dist2pix[dist].add(pixel)

while pix2dist:
    best_score = None
    best_pixel = None

    for pixel in sorted(pix2dist): # Sorting for determinism
        score = sum(len(dist2pix[d]) for d in pix2dist[pixel])

        if best_score is None or score < best_score:
            best_score = score
            best_pixel = pixel

    added_dists = pix2dist[best_pixel]
    selected_dists |= added_dists
    del pix2dist[best_pixel]
    selected.append(best_pixel)

    for d in added_dists:
        dist2pix[d].remove(best_pixel)

    to_remove = set()
    for pixel in pix2dist:
        new_dist = norm(pixel, best_pixel)

        if (new_dist in selected_dists or new_dist in pix2dist[pixel]
                or added_dists & pix2dist[pixel]):
            to_remove.add(pixel)
            continue

        pix2dist[pixel].add(new_dist)
        dist2pix[new_dist].add(pixel)

    for pixel in to_remove:
        for d in pix2dist[pixel]:
            dist2pix[d].remove(pixel)

        del pix2dist[pixel]

    print("Selected: {}, Remaining: {}, Chosen: ({}, {})".format(len(selected), len(pix2dist),
                                                                 best_pixel//N, best_pixel%N))
    sys.stdout.flush()

print(*selected)
print("Time taken:", time.time() - start_time)

I quickly brute-forced N=10 and there are many distinct layouts with 9 pts but that's the best you can do.

SWI-Prolog, score 131
Barely better than the initial answer, but I guess this will get things started a bit more. The algorithm is the same as the Python answer except for the fact that it tries pixels in an alternate way, starting with the top left pixel (pixel 0), then the bottom right pixel (pixel 383160), then pixel 1, then pixel 383159, etc.
a(Z) :-
    N = 619,
    build_list(N,Z).

build_list(N,R) :-
    M is N*N,
    get_list([M,-1],[],L),
    reverse(L,O),
    build_list(N,O,[],[],R).

get_list([A,B|C],R,Z) :-
    X is A - 1,
    Y is B + 1,
    (X =< Y,
    Z = R
    ;
    get_list([X,Y,A,B|C],[X,Y|R],Z)).

build_list(_,[],R,_,R) :- !.
build_list(N,[A|T],R,W,Z) :-
    separated_pixel(N,A,R,W,S),
    is_set(S),
    flatten([W|S],V),!,
    build_list(N,T,[A|R],V,Z)
    ;build_list(N,T,R,W,Z).

separated_pixel(N,A,L,W,R) :-
    separated_pixel(N,A,L,[],W,R).

separated_pixel(N,A,[A|T],R,W,S) :-
        separated_pixel(N,A,T,R,W,S).

separated_pixel(N,A,[B|T],R,W,S) :-
    X is (A mod N - B mod N)*(A mod N - B mod N),
    Y is (A//N - B//N)*(A//N - B//N),
    Z is X + Y,
    \+member(Z,W),
    separated_pixel(N,A,T,[Z|R],W,S).

separated_pixel(_,_,[],R,_,R).

Input:
a(A).
Output:
Z = [202089, 180052, 170398, 166825, 235399, 138306, 126354, 261759, 119490, 117393, 281623, 95521, 290446, 299681, 304310, 78491, 314776, 63618, 321423, 60433, 323679, 52092, 331836, 335753, 46989, 40402, 343753, 345805, 36352, 350309, 32701, 32470, 352329, 30256, 28089, 357859, 23290, 360097, 22534, 362132, 20985, 364217, 365098, 17311, 365995, 15965, 15156, 368487, 370980, 371251, 11713, 372078, 372337, 10316, 373699, 8893, 374417, 8313, 7849, 7586, 7289, 6922, 376588, 6121, 5831, 377399, 377639, 4941, 378494, 4490, 379179, 3848, 379453, 3521, 3420, 379963, 380033, 3017, 380409, 2579, 380636, 2450, 2221, 2006, 381235, 1875, 381369, 381442, 381682, 1422, 381784, 1268, 381918, 1087, 382144, 382260, 833, 382399, 697, 382520, 622, 382584, 382647, 382772, 384, 382806, 319, 286, 382915, 382939, 190, 172, 383005, 128, 383050, 93, 383076, 68, 383099, 52, 40, 383131, 21, 383145, 10, 383153, 4, 383158, 1, 383160, 0]

Image from Stack Snippet

Since there's a theoretical maximum of 487, even an incremental increase is significant...

Did your output as shown work with the Stack Snippet? I had specified space separated (as in my example answer) but the main reason for that was so that the Stack Snippet would work.

@trichoplax Yeah it's a typo, I start with pixel 0, I'll fix it. To get the image I selected the part of the output between the two square brackets, and removed all commas. The Stack snippet seems to work with comma-separated pixels, though.

Haskell—115 130 131 135 136
My inspiration was the Sieve of Eratosthenes and in particular The Genuine Sieve of Eratosthenes, a paper by Melissa E. O'Neill of Harvey Mudd College. My original version (which considered points in index order) sieved points extremely quickly, for some reason I can't recall I decided to shuffle the points before “sieving” them in this version (I think solely to make generating different answers easier by using a new seed in the random generator). Because the points are no longer in any sort of order there's not really any sieving going on anymore, and as a result it takes a couple minutes just to produce this single 115 point answer. A knockout Vector would probably be a better choice now.
So with this version as a checkpoint I see two branches, returning to the “Genuine Sieve” algorithm and leveraging the List monad for choice, or swapping out the Set operations for equivalents on Vector.
Edit: So for working version two I veered back toward the sieve algorithm, improved the generation of “multiples” (knocking indices out by finding points at integer coordinates on circles with radius equal to the distance between any two points, akin to generating prime multiples) and making a few constant time improvements by avoiding some needless recalculation.
For some reason I can't recompile with profiling turned on, but I believe that the major bottleneck now is backtracking. I think exploring a bit of parallelism and concurrency will produce linear speedups, but memory exhaustion will probably cap me at a 2x improvement.
Edit: Version 3 meandered a bit, I first experimented with a heuristic in taking the next 𝐧 indices (after sieving from earlier choices) and choosing the one that produced the next minimum knockout set. This ended up being far too slow, so I went back to a whole-searchspace brute force method. An idea to order the points by distance from some origin came to me, and led to an improvement by one single point (in the time my patience lasted). This version picks index 0 as the origin, it may be worth trying the center point of the plane.
Edit: I picked up 4 points by re-ordering the search space to prioritize the most distant points from the center. If you're testing my code, 135 136 is actually the second third solution found. Fast edit: This version seems the most likely to keep being productive if left running. I suspect I may tie at 137, then run out of patience waiting for 138.
One thing I noticed (that may be of help to someone) is that if you set the point ordering from the center of the plane (i.e., remove (d*d -) from originDistance) the image formed looks a bit like a sparse prime spiral.
{-# LANGUAGE RecordWildCards #-}
{-# LANGUAGE BangPatterns #-}

module Main where

import Data.Function (on)
import Data.List     (tails, sortBy)
import Data.Maybe    (fromJust)
import Data.Ratio
import Data.Set      (fromList, toList, union, difference, member)

import System.IO

sideLength :: Int
sideLength = 619

data Point = Point {  x :: !Int,  y :: !Int } deriving (Ord, Eq)
data Delta = Delta { da :: !Int, db :: !Int }

euclidean :: Delta -> Int
euclidean Delta{..} = da*da + db*db

instance Eq Delta where
  (==) = (==) `on` euclidean

instance Ord Delta where
  compare = compare `on` euclidean

delta :: Point -> Point -> Delta
delta a b = Delta (min dx dy) (max dx dy)
  where
    dx = abs (x a - x b)
    dy = abs (y a - y b)

equidistant :: Dimension -> Point -> Point -> [Point]
equidistant d a b =
  let
    (dx, dy) = (x a - x b, y a - y b)
    m = if dx == 0 then Nothing else Just (dy % dx)                    -- Slope
    w = if dy == 0 then Nothing else Just $ maybe 0 (negate . recip) m -- Negative reciprocal
    justW = fromJust w -- Moral bankruptcy
    (px, py) = ((x a + x b) % 2, (y a + y b) % 2)                      -- Midpoint
    b0 = py - (justW * px)                                             -- Y-intercept
    f q = justW * q + b0                                               -- Perpendicular bisector
  in
   maybe (if denominator px == 1 then map (Point (numerator px)) [0..d - 1] else [])
         ( map (\q -> Point q (numerator . f . fromIntegral $ q))
         . filter ((== 1) . denominator . f . fromIntegral)
         )
         (w >> return [0..d - 1])

circle :: Dimension -> Point -> Delta -> [Point]
circle d p delta' =
  let
    square = (^(2 :: Int))
    hypoteneuse = euclidean delta'
    candidates = takeWhile ((<= hypoteneuse) . square) [0..d - 1]
    candidatesSet = fromList $ map square [0..d - 1]
    legs = filter ((`member` candidatesSet) . (hypoteneuse -) . square) candidates
    pythagoreans = zipWith Delta legs
                 $ map (\l -> floor . sqrt . (fromIntegral :: Int -> Double) $ hypoteneuse - square l) legs
  in
    toList . fromList $ concatMap (knight p) pythagoreans

knight :: Point -> Delta -> [Point]
knight Point{..} Delta{..} =
    [ Point (x + da) (y - db), Point (x + da) (y + db)
    , Point (x + db) (y - da), Point (x + db) (y + da)
    , Point (x - da) (y - db), Point (x - da) (y + db)
    , Point (x - db) (y - da), Point (x - db) (y + da)
    ]

type Dimension = Int
type Index = Int

index :: Dimension -> Point -> Index
index d Point{..} = y * d + x

point :: Dimension -> Index -> Point
point d i = Point (i `rem` d) (i `div` d)

valid :: Dimension -> Point -> Bool
valid d Point{..} = 0 <= x && x < d
                 && 0 <= y && y < d

isLT :: Ordering -> Bool
isLT LT = True
isLT _  = False

sieve :: Dimension -> [[Point]]
sieve d = [i0 : sieve' is0 [i0] [] | (i0:is0) <- tails . sortBy originDistance . map (point d) $ [0..d*d - 1]]
  where
    originDistance :: Point -> Point -> Ordering
    originDistance = compare `on` ((d*d -) . euclidean . delta (point d (d*d `div` 2)))

    sieve' :: [Point] -> [Point] -> [Delta] -> [Point]
    sieve' []     _  _ = []
    sieve' (i:is) ps ds = i : sieve' is' (i:ps) ds'
      where
        ds' = map (delta i) ps ++ ds
        knockouts = fromList [k | d' <- ds
                                , k  <- circle d i d'
                                , valid d k
                                , not . isLT $ k `originDistance` i
                                ]
            `union` fromList [k | q  <- i : ps
                                , d' <- map (delta i) ps
                                , k  <- circle d q d'
                                , valid d k
                                , not . isLT $ k `originDistance` i
                                ]
            `union` fromList [e | q <- ps
                                , e <- equidistant d i q
                                , valid d e
                                , not . isLT $ e `originDistance` i
                                ]
        is' = sortBy originDistance . toList $ fromList is `difference` knockouts

main :: IO ()
main = do let answers = strictlyIncreasingLength . map (map (index sideLength)) $ sieve sideLength
          hSetBuffering stdout LineBuffering
          mapM_ (putStrLn . unwords . map show) $ answers
  where
    strictlyIncreasingLength :: [[a]] -> [[a]]
    strictlyIncreasingLength = go 0
      where
        go _ []     = []
        go n (x:xs) = if n < length x then x : go (length x) xs else go n xs

Output
1237 381923 382543 382541 1238 1857 380066 5 380687 378828 611 5571 382553 377587 375113 3705 8664 376356 602 1253 381942 370161 12376 15475 7413 383131 367691 380092 376373 362114 36 4921 368291 19180 382503 26617 3052 359029 353451 29716 382596 372674 352203 8091 25395 12959 382479 381987 35894 346031 1166 371346 336118 48276 2555 332400 46433 29675 380597 13066 382019 1138 339859 368230 29142 58174 315070 326847 56345 337940 2590 382663 320627 70553 19278 7309 82942 84804 64399 5707 461 286598 363864 292161 89126 371267 377122 270502 109556 263694 43864 382957 824 303886 248218 18417 347372 282290 144227 354820 382909 380301 382808 334361 375341 2197 260623 222212 196214 231526 177637 29884 251280 366739 39442 143568 132420 334718 160894 353132 78125 306866 140600 297272 54150 240054 98840 219257 189278 94968 226987 265881 180959 142006 218763 214475

Impressive improvements. You have 2 hours left to get to 138 before the bounty is assigned. Nice work either way...

It's looking unlikely that I'll meet that goal, I still haven't managed to generate a 137 element set. I think this method is probably tapped...

Interesting that two different answers with different approaches are hitting a maximum around the same size.

I think the upper bound is probably fairly close. Consider an infinite plane and any two points. The optimal placement of those points with any distance d minimizes the number of other points excluded from consideration by tracing circles of radius d with centers of both chosen points, where the perimeter only touches three other integer coordinates (at 90, 180, and 270 degree turns about the circle), and the perpendicular bisecting line crossing no integer coordinates. So each new point n+1 will exclude 6n other points from consideration (with optimal choice).

Python, 134 132
Here's a simple one that randomly culls some of the search space to cover a greater area.  It iterates the points in distance from an origin order.  It skips points that are same distance from origin, and early-outs if it cannot improve on the best.  It runs indefinitely.
from random import *
from bisect import *

W = H = 619
pts = []
deepest = 0
lengths = set()

def place(x, y):
    global lengths
    pos = (x, y)
    for px, py in pts:
        dist = (x-px)*(x-px) + (y-py)*(y-py)
        if dist in lengths:
            return False
    dists = set((x-px)*(x-px) + (y-py)*(y-py) for px, py in pts)
    if len(dists) != len(pts):
        return False
    lengths |= dists
    pts.append(pos)
    return True

def unplace():
    x, y = pos = pts.pop()
    for px, py in pts:
        dist = (x-px)*(x-px) + (y-py)*(y-py)
        lengths.remove(dist)

def walk(i):
    global deepest, backtrack
    depth = len(pts)
    while i < W*H:
        d, x, y, rem = order[i]
        if rem+depth <= deepest: # early out if remaining unique distances mean we can't improve
            return
        i += 1
        if place(x, y):
            j = i
            while j < W*H and order[j][0] == d: # skip those the same distance from origin
                j += 1
            walk(j)
            unplace()
            if backtrack <= depth:
                break
            if not randint(0, 5): # time to give up and explore elsewhere?
                backtrack = randint(0, len(pts))
                break
            backtrack = W*H # remove restriction
    if depth >= deepest:
        deepest = depth
        print (ox, oy), depth, "=", " ".join(str(y*W+x) for x, y in pts)

try:
    primes = (0,1,2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97)
    while True:
        backtrack = W*H
        ox, oy = choice(primes), choice(primes) # random origin coordinates
        order = sorted((float((ox-x)**2+(oy-y)**2)+random(), x, y) for x in xrange(W) for y in xrange(H))
        rem = sorted(set(int(o[0]) for o in order)) # ordered list of unique distances
        rem = {r: len(rem)-bisect_right(rem, r) for r in rem} # for each unique distance, how many remain?
        order = tuple((int(d), x, y, rem[int(d)]) for i, (d, x, y) in enumerate(order))
        walk(0)
except KeyboardInterrupt:
    print

It quickly finds solutions with 134 points:
3097 3098 2477 4333 3101 5576 1247 9 8666 8058 12381 1257 6209 15488 6837 21674 19212 26000 24783 1281 29728 33436 6863 37767 26665 14297 4402 43363 50144 52624 18651 9996 58840 42792 6295 69950 48985 34153 10644 72481 83576 63850 29233 94735 74997 60173 43477 101533 102175 24935 113313 88637 122569 11956 36098 79401 61471 135610 31796 4570 150418 57797 4581 125201 151128 115936 165898 127697 162290 33091 20098 189414 187620 186440 91290 206766 35619 69033 351 186511 129058 228458 69065 226046 210035 235925 164324 18967 254416 130970 17753 248978 57376 276798 456 283541 293423 257747 204626 298427 249115 21544 95185 231226 54354 104483 280665 518 147181 318363 1793 248609 82260 52568 365227 361603 346849 331462 69310 90988 341446 229599 277828 382837 339014 323612 365040 269883 307597 374347 316282 354625 339821 372016
For the curious, here are some brute-forced small N:
3  =  0  2  3
4  =  0  2  4  7
5  =  0  2  5 17 23
6  =  0 12 21 28 29 30
7  =  4  6 11 14 27 36 42
8  =  0  2  8 11 42 55 56
9  =  0  2  9 12 26 50 63 71
10 =  0  2  7 10 35 75 86 89  93
11 =  0 23 31 65 66 75 77 95 114 117

Have you tried running this through PyPy?

@trichoplax I always run these hobby things on both pypy and cpython, and if cpython is faster I lodge tickets on pypy.  In this particular case, pypy is quite a bit faster than cpython and that's how I got these numbers :)

I'm interested, what does "quickly" entail?

@Cain 'quickly' was about 5 mins iirc

Python 3, 130
For comparison, here's a recursive backtracker implementation:
N = 619

def norm(p1, p2):
    return (p1//N - p2//N)**2 + (p1%N - p2%N)**2

def solve(selected, dists):
    global best

    if len(selected) > best:
        print(len(selected), "|", *selected)
        best = len(selected)

    for pixel in (range(selected[-1]+1, N*N) if selected else range((N+1)//2+1)):
        # By symmetry, place first pixel in first half of top row
        added_dists = [norm(pixel, p) for p in selected]
        added_set = set(added_dists)

        if len(added_set) != len(added_dists) or added_set & dists:
            continue

        selected.append(pixel)
        dists |= added_set

        solve(selected, dists)

        selected.pop()
        dists -= added_set

print("N =", N)
best = 0
selected = []
dists = set()
solve(selected, dists)

It finds the following 130 pixel solution quickly before it starts choking:
0 1 3 7 12 20 30 44 65 80 96 122 147 181 203 251 289 360 400 474 564 592 627 660 747 890 1002 1155 1289 1417 1701 1789 1895 2101 2162 2560 2609 3085 3121 3331 3607 4009 4084 4242 4495 5374 5695 6424 6762 6808 7250 8026 8356 9001 9694 10098 11625 12881 13730 14778 15321 16091 16498 18507 19744 20163 20895 23179 25336 27397 31366 32512 33415 33949 39242 41075 46730 47394 48377 59911 61256 66285 69786 73684 79197 89530 95447 102317 107717 111751 116167 123198 126807 130541 149163 149885 154285 159655 163397 173667 173872 176305 189079 195987 206740 209329 214653 220911 230561 240814 249310 269071 274262 276855 285295 305962 306385 306515 312310 314505 324368 328071 348061 350671 351971 354092 361387 371800 376153 378169

More importantly, I'm using it to check solutions for small cases. For N <= 8, the optimal are:
1: 1 (0)
2: 2 (0 1)
3: 3 (0 1 5)
4: 4 (0 1 6 12)
5: 5 (0 1 4 11 23)
6: 6 (0 1 9 23 32 35)
7: 7 (0 2 9 20 21 40 48)
8: 7 (0 1 3 12 22 56 61)
9: 8 (0 1 3 8 15 37 62 77)
10: 9 (0 1 7 12 30 53 69 80 89)

Listed in brackets are the first lexicographical optimals.
Unconfirmed:
11: 10 (0 2 3 7 21 59 66 95 107 120)
12: 10 (0 1 3 7 33 44 78 121 130 140)

Scala, 132
Scans left-to-right and top-to-bottom like the naive solution, but tries starting at different pixel locations.
import math.pow
import math.sqrt

val height, width = 619
val area = height * width

case class Point(x: Int, y: Int)

def generate(n: Int): Set[Point] = {

  def distance(p: Point, q: Point) = {
    def square(x: Int) = x * x
    sqrt(square(q.x - p.x) + square(q.y - p.y))
  }

  def hasDuplicates(s: Seq[_]) = s.toSet.size != s.size

  def rotate(s: Vector[Point]): Vector[Point] = s.drop(n) ++ s.take(n)

  val remaining: Vector[Point] =
    rotate((for (y <- 0 until height; x <- 0 until width) yield { Point(x, y) }).toVector)
  var unique = Set.empty[Point]
  var distances = Set.empty[Double]
  for (candidate <- remaining) {
    if (!unique.exists(p => distances.contains(distance(candidate, p)))) {
      val candidateDistances = unique.toSeq.map(p => distance(candidate, p))
      if (!hasDuplicates(candidateDistances)) {
        unique = unique + candidate
        distances = distances ++ candidateDistances
      }
    }
  }
  unique
}

def print(s: Set[Point]) = {
  def toRowMajor(p: Point) = p.y*height + p.x
  println(bestPixels.map(toRowMajor).toSeq.sorted.mkString(" "))
}

var bestPixels = Set.empty[Point]
for (n <- 0 until area) {
  val pixels = generate(n)
  if (pixels.size > bestPixels.size) bestPixels = pixels
}
print(bestPixels)

Output
302 303 305 309 314 322 332 346 367 382 398 424 449 483 505 553 591 619 647 680 719 813 862 945 1014 1247 1459 1700 1740 1811 1861 1979 2301 2511 2681 2913 3114 3262 3368 4253 4483 4608 4753 5202 5522 5760 6246 6474 6579 6795 7498 8062 8573 8664 9903 10023 10567 10790 11136 12000 14153 15908 17314 17507 19331 20563 20941 22339 25131 26454 28475 31656 38328 39226 40214 50838 53240 56316 60690 61745 62374 68522 71208 78598 80204 86005 89218 93388 101623 112924 115702 118324 123874 132852 136186 139775 144948 154274 159730 182200 193642 203150 203616 213145 214149 218519 219744 226729 240795 243327 261196 262036 271094 278680 282306 289651 303297 311298 315371 318124 321962 330614 336472 343091 346698 354881 359476 361983 366972 369552 380486 382491

Just getting the ball rolling...

Python 3, score 129
This is an example answer to get things started.
Just a naive approach going through the pixels in order and choosing the first pixel that doesn't cause a duplicate separation distance, until the pixels run out.
Code
width = 619
height = 619
area = width * height
currentAttempt = 0

temporaryLengths = []
lengths = []
points = []
pixels = []
for i in range(area):
    pixels.append(0)

def generate_points():
    global lengths
    while True:
        candidate = vacantPixel()
        if isUnique(candidate):
            lengths += temporaryLengths
            pixels[candidate] = 1
            points.append(candidate)
            print(candidate)
        if currentAttempt == area:
            break
    filename = 'uniquely-separated-points.txt'
    with open(filename, 'w') as file:
        file.write(' '.join(points))

def isUnique(n):
    x = n % width
    y = int(n / width)
    temporaryLengths[:] = []
    for i in range(len(points)):
        point = points[i]
        a = point % width
        b = int(point / width)
        d = distance(x, y, a, b)
        if d in lengths or d in temporaryLengths:
            return False
        temporaryLengths.append(d)
    return True

def distance(x1, y1, x2, y2):
    xd = x2 - x1
    yd = y2 - y1
    return (xd*xd + yd*yd) ** 0.5

def vacantPixel():
    global currentAttempt
    while True:
        n = currentAttempt
        currentAttempt += 1
        if pixels[n] == 0:
            break
    return n

generate_points()

Output
0 1 3 7 12 20 30 44 65 80 96 122 147 181 203 251 289 360 400 474 564 592 627 660 747 890 1002 1155 1289 1417 1701 1789 1895 2101 2162 2560 2609 3085 3121 3331 3607 4009 4084 4242 4495 5374 5695 6424 6762 6808 7250 8026 8356 9001 9694 10098 11625 12881 13730 14778 15321 16091 16498 18507 19744 20163 20895 23179 25336 27397 31366 32512 33415 33949 39242 41075 46730 47394 48377 59911 61256 66285 69786 73684 79197 89530 95447 102317 107717 111751 116167 123198 126807 130541 149163 149885 154285 159655 163397 173667 173872 176305 189079 195987 206740 209329 214653 220911 230561 240814 249310 269071 274262 276855 285295 305962 306385 306515 312310 314505 324368 328071 348061 350671 351971 354092 361387 369933 376153

Image from Stack Snippet

Fantom 96
I used an evolution algorithm, basically add k random points at a time, do that for j different random sets, then pick the best one and repeat. Pretty terrible answer right now, but that's running it with only 2 children per generation for the sake of speed, which is nearly just random. Gonna play with the parameters a bit to see how it goes, and I probably need a better scoring function than number of free spots left.
class Pixel
{
  static const Int n := 619
  static const Int stepSize := 20
  static const Int generationSize := 5
  static const |Int, Int -> Int| d := |Int x, Int y -> Int| {
      d1 := x%n - y%n
      d2 := x/n - y/n
      return d1.pow(2) + d2.pow(2)
    }

  public static Void main(){

    //Initialize

    [Int: Int[][]] disMap := [:]
    Int[] freeSpots := (0..<n*n).toList
    Int[] pixels := [,]
    Int[] distances := [,]

    genNum := 0
    children := [,]
    while(freeSpots.size > 0){
      echo("Generation: ${genNum++} \t Spots Left: ${freeSpots.size} \t Pixels added: $pixels.size \t Distances used: $distances.size uniqueDistances: $distances.unique.size" )
      echo(distances)
      echo("Pixels: " + pixels.join(" "))
      //echo("Distances: $distances")
      //Generate children
      children = [,]
      generationSize.times{
        //echo("\tStarting child $it")
        i := Int.random(0..<freeSpots.size)
        childFreeSpots := freeSpots.dup
        childPixels := pixels.dup
        childDistances := distances.dup

        for(Int step := 0; step < stepSize; step++){

          if( i < childFreeSpots.size){
            //Choose a pixel
            pixel := childFreeSpots.removeAt(i)
            //echo("\t\tAdding pixel $pixel")

            //Remove neighbors that are the new distances away
            ///Find distances
            newDis := [,]
            childPixels.each {
              newDis.add(d(pixel, it))
            }

            //Check that there are no equal distances
            if(newDis.size != newDis.unique.size) continue

            //Remove neighbors
            childPixels.each | Int childPixel|{
              newDis.each |Int dis|{
                neighbors := getNeighbors(childPixel, dis, disMap)
                neighbors.each| Int n |{
                  index := childFreeSpots.binarySearch(n)
                  if(index >= 0) childFreeSpots.removeAt(index)
                }
              }
            }
            //echo("Removed neighbors: $test")
            //Remove all the neighbors of new pixel
            childDistances.addAll(newDis)
            childDistances.each|Int dis| {
              neighbors := getNeighbors(pixel, dis, disMap)
              childFreeSpots.removeAll(neighbors)
            }

            //Add new pixel
            childPixels.add(pixel)
          }
        }
        children.add([childPixels.dup, childDistances.dup, childFreeSpots.dup])
        echo("\tChild $it: pixels: $childPixels.size \t distances: $childDistances.size \t freeSpots: $childFreeSpots.size")
      }

      //Score children and keep best one as new parent
      Obj?[][] parent := children.max |Int[][] a, Int[][] b -> Int| { return (a.last.size  + a.first.size*10000) <=> (b.last.size + b.first.size*10000)  }
      pixels = parent.first
      distances = parent[1]
      freeSpots = parent.last

    }//End while

    //Return result
    echo("Size: " + pixels.size)
    echo(pixels.join(" "))

  }

  private static Bool checkValid(Int[] pixels){
    distances := [,]
    pixels[0..-2].each|Int p, Int i|{
      for(Int j := i + 1; j < pixels.size; j++){
        distances.add(d(p, pixels[j]))
      }
    }
    if(distances.size > distances.unique.size){
      echo("Duplicate distance found!!!!")
      echo("Pixel $pixels.last is not valid")
      return false
    }
    return true
  }

  public static Int[] getNeighbors(Int spot, Int distance, [Int : Int[][]] disMap ){
    result := [,]
    //Check hash map
    pairs := disMap.get(distance, null)

    //Find possible int pairs if not already in the map
    if(pairs == null){
      for(Int i := 0; i*i <= distance; i++ ){
        for(Int j := i; j*j + i*i <= distance; j++){
          if(i.pow(2) + j.pow(2) == distance){
            pairs.add([i, j])
          }
        }
      }
      disMap.add(distance, pairs)
    }

    pairs.each|Int[] pair|{
      //Find neighbors with pair
      x := pair.first
      y := pair.last
      2.times{
        //Positive x
        result.add(spot + x + y*n)
        result.add(spot + x - y*n)

        //negative x
        result.add(spot - x + y*n)
        result.add(spot - x - y*n)

        //Swap x and y and repeat
        temp := x
        x = y
        y = temp
      }
    }

    return result.findAll |Int i -> Bool| { i >= 0 }.unique
  }

}

Output
17595 17596 17601 17627 17670 17726 17778 17861 17956 18117 18324 18733 19145 19597 20244 21139 21857 22742 24078 25343 28577 30152 32027 34406 37008 39864 42313 44820 48049 52193 55496 59707 64551 69976 74152 79758 84392 91782 98996 104625 150212 158877 169579 178660 189201 201343 213643 225998 238177 251012 263553 276797 290790 304915 319247 332702 347266 359665 373683 125899 144678 170677 195503 220092 244336 269861 289473 308633 326736 343756 358781 374280 131880 172485 212011 245015 277131 302055 321747 347911 363717 379166 249798 284200 313870 331913 360712 378024 9704 141872 249686 293656 357038 357596 370392 381963

Oh wow, you're right, I'm sorry. Hmm, must have not copied all of it early when I tested. I will fix whatever is going on and respond with an update

Ahh, I figured it out, when adding a new pixel, I wasn't checking that it's not equidistant from two other pixels

Fixed it, but it really sucks now, I think I might be accidentally finding a worst solution instead of best one

At least it works now, so you can tweak the parameters and see if you can improve the result. Great to see another new approach. +1

Python 3, 119
I no longer remember why I named this function mc_usp, though I suspect it had something to do with Markov chains. Here I publish my code that I ran with PyPy for about 7 hours. The program tries to build up 100 different sets of pixels by randomly picking pixels until it has checked every pixel in the image, and returning one of the best sets.
On another note, at some point, we should really attempt to find an upper bound for N=619 that's better than 488, because judging from the answers on here, that number is way too high. Rowan Blush's comment about how every new point n+1 can potentially remove 6*n points with optimal choice seemed like a good idea. Unfortunately, upon inspection of the formula a(1) = 1; a(n+1) = a(n) + 6*n + 1, where a(n) is the number of points removed after adding n points to our set, this idea may not be the best fit. Checking when a(n) is greater than N**2, a(200) being larger than 619**2 seems promising, but the a(n) larger than 10**2 is a(7) and we've proved that 9 is the actual upper bound for N=10. I'll keep you posted as I try to look a better upper bound, but any suggestions are welcome.
Onto my answer. First, my set of 119 pixels.
15092 27213 294010 340676 353925 187345 127347 21039 28187 4607 23476 324112 375223 174798 246025 185935 186668 138651 273347 318338 175447 316166 158342 97442 361309 251283 29986 98029 339602 292202 304041 353401 236737 324696 42096 102574 357602 66845 40159 57866 3291 24583 254208 357748 304592 86863 19270 228963 87315 355845 55101 282039 83682 55643 292167 268632 118162 48494 378303 128634 117583 841 178939 20941 161231 247142 110205 211040 90946 170124 362592 327093 336321 291050 29880 279825 212675 138043 344012 187576 168354 28193 331713 329875 321927 129452 163450 1949 186448 50734 14422 3761 322400 318075 77824 36391 31016 33491 360713 352240 45316 79905 376004 310778 382640 383077 359178 14245 275451 362125 268047 23437 239772 299047 294065 46335 112345 382617 79986

Second, my code, which randomly picks a starting point from an octant of the 619x619 square (since the starting point is otherwise equal under rotation and reflection) and then every other point from the rest of the square.
import random
import time

start_time = time.time()
print(start_time)

def mc_usp_v3(N, z, k=100, m=1.0):
    """
    At m=1.0, it keeps randomly picking points until we've checked every point. Oh dear.
    """
    ceil = -(-N//2)
    a=random.randint(0,ceil)
    b=random.randint(a,ceil)
    r=[a*N+b]

    best_overall = r[:]
    all_best = []
    best_in_shuffle = r[:]
    num_shuffles = 0
    num_missteps = 0
    len_best = 1

    while num_shuffles < k and len(best_overall) < z:
        dist = []
        missteps = []
        points_left = list(range(N*N))
        points_left.remove(r[0])

        while len_best + num_missteps < m*N*N and len(points_left):
            index = random.randint(0, len(points_left)-1)
            point = points_left[index]
            points_left.pop(index)
            dist, better = euclid(r, point, dist, N)

            if better and len(r) + 1 > len_best:
                r.append(point)
                best_in_shuffle = r[:]
                len_best += 1
            else:
                missteps.append(point)
                num_missteps += 1

        else:
            print(num_shuffles, len(best_overall), len_best, num_missteps, time.time() - start_time)

            num_shuffles += 1
            num_missteps = 0
            missteps = []

            if len(best_in_shuffle) == len(best_overall):
                all_best.append(best_in_shuffle)
                print(best_in_shuffle)

            if len(best_in_shuffle) > len(best_overall):
                best_overall = best_in_shuffle[:]
                all_best = [best_overall]
                print(best_overall)
            a=random.randint(0,ceil)
            b=random.randint(a,ceil)
            r=[a*N+b]
            best_in_shuffle = r[:]
            len_best = 1
    return len(best_overall), all_best

def euclid(point_set, new_point, dist, N):
    new_dist = []
    unique = True
    a,b=divmod(new_point, N)
    for point in point_set:
        c,d=divmod(point, N)
        current_dist = (a-c)**2+(b-d)**2
        if current_dist in dist or current_dist in new_dist:
            unique = False
            break
        new_dist.append(current_dist)
    if unique:
        dist += new_dist
    return dist, unique

def mcusp_format(mcusp_results):
    length, all_best = mcusp_results
    return " ".join(str(i) for i in all_best[0])

print(mcusp_format(mc_usp_v3(10, 20, 100, 1.0)))
print(mcusp_format(mc_usp_v3(619, 488, 100, 1.0)))
print(time.time()-start_time)

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/


text:
Write a program that finds the most occurring paired letter in a string
The program must output the letter that is paired the most. For example, if your program was given the following string:
"Sally's friend Bobby searched for seashells."

it must output L because "ll" occurs twice, which is more frequent than the other pair "bb".
Rules:

If more than one letter has 1st place for occurrences, then output all of them in alphabetical order (e.g. "Sally's friends Jimmy and Bobby rummaged for seashells." should output both L AND M [or "LM" if you please] because they both occur more frequently than other pairs.)
Letters which are tripled, quadrupled, etc. count as one pair (e.g. "lll" in "willless" is counted as only one pair of L.)
Letter pairs must be in one word (e.g. "Sally's sociable friends Sammy and Bobby searched for fabulous seashells." should output L and not S because despite "ss" having more occurrences than "ll", they are separated by spaces.)
Count only letters from the English alphabet
Case does not matter (e.g. "Ss" is the same as "SS" or "ss", and all are counted as one pair of S.)

You may read your input from wherever you please. Shortest code wins.

Can we assume that only letters will occur in pairs or could the input contain double spaces or double ' etc?

Can we assume that at least one letter appears twice?

@MartinBüttner Yes, you can assume at least one letter pair occurs. However, other characters may appear in pairs as well. Only count letters.

Even if there is only one pair, can I still print it in a list like ['l']?

@Maltysen Yes, you may do so.

Fun challenge. My main suggestion would be to add some more tests. For example: Mix of upper/lower case for repeated letters. Tripled/quadrupled letters. Non-letter character is the most repeated.

@AlexA. It says "output all of them in alphabetical order".

@RetoKoradi: Sorry, missed that. Thanks.

Pyth, 26 25 24 16 15 bytes
.M/sfthTrrz08ZG

Try it online: Demonstration
Explanation:
.M/sfthTrrz08ZG   implicit: z = input string
         rz0      convert z to lower-char
        r   8     run-length-encoding (into tuples [count, char])
    f             filter for tuples T, which satisfy:
     thT            T[0]-1 != 0 (count > 1)
   s              join to a string
.M            G   find the elements Z of "abcd...z", which produce the highest value:
  /...........Z       count of Z in ...

eC -> s saves one byte.

Do you know any good resources which I could use to learn Pyth?

@BetaDecay You can find a tutorial about Pyth at http://pyth.readthedocs.org It doesn't cover all functionalities and tricks, but it is a good start. And if you have any questions, just ask in the chat.

Bash + GNU coreutils, 133
grep -Eo '([A-Z])\1+'<<<"${1^^}"|cut -c1|sort|uniq -c|sort -rn|while read n l
do((a-n&&a-0))&&exit||echo $l&&a=$n
done|sort|tr -d \\n

Testcases:
$ for t in "Sally's friend Bobby searched for seashells." \
> "Sally's friends Jimmy and Bobby rummaged for seashells." \
> "willless" \
> "Sally's sociable friends Sammy and Bobby searched for fabulous seashells." \
> "11ss11aa"
> do
> ./mostpaired.sh "$t"
> echo
> done
L
LM
LS
L
AS
$

Does it only count letters? (testcase: 11ss11aa -> SA)

@edc65 There I fixed it ;-).  Actually, 11ss11aa -> AS :)

I think that your sort -r needs to be sort -rn if you have 10 or more paired letters.

@TobySpeight.  Yes.  Fixed.

can make it shorter with AWK instead of while:
awk '!n{n=$1};n==$1'|grep -o .$

CJam, 29 27 bytes
leue`{2a>},s_el-$e`$z~\)-,>

Thanks to @Optimizer for golfing off 2 bytes!
Try it online in the CJam interpreter.
How it works
leu    e# Read a line from STDIN and covert to uppercase.
e`     e# Perform run-length encoding.
       e# Example: "AABBBC!!" -> [[2 'A] [3 'B] [1 'C] [2 '!]]
{2a>}, e# Filter out all pairs that are less of equal to [2].
s      e# Stringify.
       e# Example: [[2 'A] [3 'B] [2 '!]] -> "2A3B2!"
_el    e# Push a copy of the string and convert to lowercase.
       e# Example: "2A3B2!" -> "2a3b2!"
-      e# Remove all character from the second string from the first.
       e# Example: "2A3B2!" "2a3b2!" - -> "AB"
$e`$   e# Sort, perform run-length encoding and sort again.
       e# Example: "ABACABDBC" -> "AAABBBCCD"
       e#                      -> [[3 'A] [3 'B] [2 'C] [1 'D]]
                               -> [[1 'D] [2 'C] [3 'A] [3 'B]]
z~     e# Zip and dump.
       e# Example: [[1 'D] [2 'C] [3 'A] [3 'B]] -> [1 2 3 3] ['D 'C 'A 'B]
\)     e# Pop out the last element from the first array.
       e# Example: [1 2 3 3] -> [1 2 3] 3
-      e# Remove all occurrences of the popped element from the array.
       e# Example: [1 2 3] 3 -> [1 2]
,      e# Compute the length of the remainder.
>      e# Skip that many elements from the character array.

z~\)-,> should work as far as I can see.

@Optimizer: Shorter and a lot more intuitive. Thanks!

Pyth - 23 22 21 20 bytes
Uses regexp substitution to replace all of two or greater of alphabet to a temp value, and uses .Maximal to get all the have the highest occurrence. Thanks to @Jakube for pointing out the redundancy of sorting and saving a byte.
.M/:rz0+Z"{2,}"KC0KG

Takes input from stdin and outputs like ['l', 'm'] to stdout.
.M        G         Values which yield maximal amount over lowercase alphabet
 /                  Count
  :                 Regexp substitution
   rz0              Lowercased input
   +Z               String concatenate current loop var
    "{2,}"          Regexp 2 or more of previous group
   KCZ              Inline assign null byte to K and use value
  K                 Count K which is null byte

Try it online here.

C, 155
Something different, no regexps.
m=1,i=1,n[91];
main(c,a)char**a;
{
  char*t=a[1];
  for(;c=*t++;)(c&=95)>64&&c<91&&(c-(*t&95)?i=1:(c=(n[c]+=i),i=0,m=m<c?c:m));
  for(c=0;++c<91;)n[c]-m||putchar(c);
}

Python 2, 132 143 bytes
def f(x):import re;x=re.findall(r'(.)\1+',x.upper());s={l:x.count(l)for l in x};print "".join(sorted([l for l in s if s[l]==max(s.values())]))

Example run:
f("Sally's friends Jimmy and Bobby rummaged for seashells.")
LM

Probably it doesn't fulfil "Letters which are tripled, quadrupled, etc. count as one pair"

You're right! i've tried to fix it. thanks for pointing that out :)

R, 105 bytes
cat(substr(names(b<-table(regmatches(s<-toupper(readline()),gregexpr("([A-Z])\\1+",s))))[b==max(b)],1,1))

This reads a line of text from STDIN and prints a space delimited list of the most common paired letters to STDOUT.
Ungolfed + explanation:
# Read a string from STDIN, convert to uppercase
s <- toupper(readline())

# Get each match of the regex /([A-Z])\1+/
matches <- regmatches(s, gregexpr("([A-Z])\\1+", s))

# Compute the frequency of each match
freq <- table(matches)

# Get the matches with the highest frequency
highest <- names(freq)[freq == max(freq)]

# Each element of highest is the literal pair, so take the first character
first <- substr(highest, 1, 1)

# Print to STDOUT
cat(first)

Examples:
> (code)
Sally's friends Jimmy and Bobby rummaged for seashells.
L M

> (code)
Sally's friend Bobby searched for seashells.
L

> (code)
Sally's sociable friends Sammy and Bobby searched for fabulous seashells.
L

> (code)
11ss11nn
N S

You can try it online!

You can probably get rid of the toupper if you ignore case and use perl in your gregexpr.  eg cat(substr(names(b<-table(regmatches(s<-readline(),gregexpr("(\\w)\\1+",s,T,T))))[b==max(b)],1,1))

@MickyT: I had thought about that, but it looks like the OP wants the output to be uppercase, so it'd have to use toupper to ensure that anyway.

That's ashame, I missed that when I read the question.

Tried the fiddle, but it seems to run forever with no outpuy in firefox. Testcase : 11ss11nn?

@edc65 It's a problem with R-Fiddle; nothing at all works. I contacted their admin to report the issue. Fixed my regex and now your test works as expected, but it cost me 2 bytes. Thanks for pointing this stuff out, I appreciate it!

CJam, 37 bytes
leue`{~_el-\(e&},1f=$e`$_W=0=f-{,1=},

Try it online
Without regular expression support, I'm afraid it will be tough to compete with Pyth. This is the best I came up with on a first pass.
Explanation:
l     Get input.
eu    Convert it to upper case, since case does not matter.
e`    Run length encoding, to split into groups of same characters.
{     Start of block for filtering.
  ~     Unpack the length/letter pair.
  _     Copy the letter.
  el    Change copy to lower case.
  -     Subtract to compare. If result is non-zero, this is a letter.
  \     Swap count to top.
  (     Decrement to get truthy value for count > 1.
  e&    Logical and: It's a letter, and count is > 1.
},    End of filter.
1f=   Don't need the counts anymore, filter out the letters only from the RLE pairs.
$     Sort them, so that multiples of the same letter are sequential.
e`    RLE again, to count how many multiples of each letter we had.
$     And sort again, to get the count/letter pairs in order of incrementing count.
_     Copy list.
W=0=  Pick out count of last element, which is the highest count.
f-    Remove count from pairs that have the highest count. This leaves them
      as one member lists with letter only, while others still have count/letter.
{     Start block for filter.
  ,1=   Check for list length one.
},    End filter.

Q (66)
Relatively readable to boot:
{where g=max g:.Q.A#count each group y where not differ y:upper x}

Python 2, 185 159 153
i=input().lower()
d=sorted({l:len(i.split(l+l))for l in map(chr,range(97,123))}.items(),None,lambda x:x[1],1)
print sorted(c for c,k in d if k==d[0][1])

Takes input as a quoted string.

Ruby, 60
f=->s{(?a..?z).group_by{|l|s.scan(/#{l*2}+/i).size}.max[1]}

p f["Sally's friends Jimmy and Bobby rummaged for seashells."]

group_by creates a hash (dictionary) structure where the keys are the output of the block and the values are lists of letters that result in each key. In this case, the keys are counts of 2+ runs of a letter, case-insensitive. max compares each [key,value] tuple lexicographically, so it just finds the maximum key. Then [1] returns the value list part of the tuple.

C# 160 Bytes
Where s is the input:
char? d(string s){s=s.ToUpper();return s.Select((x,i)=>new{y=x,z=i==0?(char?)null:s[i-1]}).Where(x=>x.y==x.z).GroupBy(x=>x.z).OrderBy(x=>x.Count()).Last().Key;}

rs, 146 bytes
[^A-Za-z]/
*(.)\1+/\1\1
*(.)(?!\1)/
$/#
+*(.)#(?!.*?\1)/#\1
+*(.)(.*)#(.*)\1/\2#\3\1\1
#/
*(.)(\1*)/\1(_)^^((^^\1\2))
([^_])(_+)(?!_)(?=.*\2_)/
_/

Try it! Please! It took me forever to make the buttons even with the output box on that page...
Well, this was fairly...crazy. The logic here is kind of weird; I'll only post an explanation if someone asks. (Of course, I also said that for an INTERCAL answer whose explanation was requested...that I never explained... ;)

I like the interpreter, but you might want to put the debug checkbox on the same line as the buttons or something. It kinda looks weird all the way up their. Still cool! +1

Tried it (errors...) http://i.stack.imgur.com/mTioT.png

@Maltysen I'll consider that. Thanks!

@edc65 Damn...what was the complete error message? Sounds like it could be a PyPy.js bug. Or just the fact I never tested this on Firefox...

JavaScript 156 153

var x=prompt(),f={},a=0,o
x.toUpperCase().replace(/([A-Z])\1+/g,function(m,s){m=f[s]=-~f[s]
if(a==m)o.push(s)
if(a<m)a=m,o=[s]})
alert(o.sort().join(""))

f[s]?f[s]+1:1 -> -~f[s]

It fails with non-letters: Count only letters from the English alphabet

Thanks @edc65. I've added the shortcut and A-Z checking.

Your exact code, streamlned and ES6: f=x=>{x.toUpperCase(f={},a=0,o).replace(/([A-Z])\1+/g,(m,s)=>a<(m=f[s]=-~f[s])?(a=m,o=[s]):a>m?0:o.push(s));alert(o.sort().join'')} (the last 2 '' are really backticks, &#96

Bash + textutils (grep, sed), 111 chars
fold -1<<<$s|uniq -iD|sort -f|uniq -ic|sort -rn|grep -i [A-Z]|sed -n '1h;G;s/\(\s*\S\+\s\)\(.\)\n\1./\2/p'|sort

Bash + awk (instead of sed), 97 chars
fold -1<<<$s|uniq -iD|sort -f|uniq -ic|sort -rn|grep -i [A-Z]|awk '!n{n=$1};n==$1{print $2}'|sort

to test it, first assign s
s="Sally's friends Jimmy ää and Bobby rummaged ää for seashells."

R, 98 bytes
Very similar to Alex's solution, but uses a substitution rather than a match to determine consecutive letters.  Scan is used to get the input and also to split the substitution result on spaces.
cat(names(a<-table(scan(,'',t=gsub('([A-z]?)(\\1?)[^A-z]*','\\U\\2 ',scan(,''),T,T))))[a==max(a)])

A couple of tests
> cat(names(a<-table(scan(,'',t=gsub('([A-z]?)(\\1?)[^A-z]*','\\U\\2 ',scan(,''),T,T))))[a==max(a)])
1: 11 was a race horse, 22 was one too. 11 won one race and 22 one won too.
19:
Read 18 items
Read 2 items
O
> cat(names(a<-table(scan(,'',t=gsub('([A-z]?)(\\1?)[^A-z]*','\\U\\2 ',scan(,''),T,T))))[a==max(a)])
1: Sally's friends Jimmy and Bobby rummaged for seashells.
9:
Read 8 items
Read 5 items
L M
> cat(names(a<-table(scan(,'',t=gsub('([A-z]?)(\\1?)[^A-z]*','\\U\\2 ',scan(,''),T,T))))[a==max(a)])
1: 11ss11nn
2:
Read 1 item
Read 2 items
N S
>