text:
Is q a quadratic residue of n?
Given two inputs q n determine if q is a quadratic residue of n.
That is, is there an x where x**2 == q (mod n) or is q a square mod n?
Input
Two integers q and n, where q and n are any integers 0 <= q < n.
Output
A truthy or a falsey.
Optionally, print any (or all) x that is x**2 == q (mod n)
Examples
>>> quadratic_residue(1, 5)
True
>>> quadratic_residue(3, 8)
False
>>> quadratic_residue(15, 22)
True

Rules
Your code must be a program or a function. The inputs can be in any order. This is code golf, so shortest code in bytes wins.
If anything is unclear or otherwise needs fixing, please let me know.
Bonuses

2-byte bonus if your function accepts q as any arbitrary integer.

Catalogue

var QUESTION_ID=65329;var ANSWER_FILTER="!t)IWYnsLAZle2tQ3KqrVveCRJfxcRLe";var COMMENT_FILTER="!)Q2B_A2kjfAiU78X(md6BoYk";var OVERRIDE_USER=47581;var answers=[],answers_hash,answer_ids,answer_page=1,more_answers=true,comment_page;function answersUrl(index){return"http://api.stackexchange.com/2.2/questions/"+QUESTION_ID+"/answers?page="+index+"&pagesize=100&order=desc&sort=creation&site=codegolf&filter="+ANSWER_FILTER}function commentUrl(index,answers){return"http://api.stackexchange.com/2.2/answers/"+answers.join(';')+"/comments?page="+index+"&pagesize=100&order=desc&sort=creation&site=codegolf&filter="+COMMENT_FILTER}function getAnswers(){jQuery.ajax({url:answersUrl(answer_page++),method:"get",dataType:"jsonp",crossDomain:true,success:function(data){answers.push.apply(answers,data.items);answers_hash=[];answer_ids=[];data.items.forEach(function(a){a.comments=[];var id=+a.share_link.match(/\d+/);answer_ids.push(id);answers_hash[id]=a});if(!data.has_more)more_answers=false;comment_page=1;getComments()}})}function getComments(){jQuery.ajax({url:commentUrl(comment_page++,answer_ids),method:"get",dataType:"jsonp",crossDomain:true,success:function(data){data.items.forEach(function(c){if(c.owner.user_id===OVERRIDE_USER)answers_hash[c.post_id].comments.push(c)});if(data.has_more)getComments();else if(more_answers)getAnswers();else process()}})}getAnswers();var SCORE_REG=/<h\d>\s*([^\n,<]*(?:<(?:[^\n>]*>[^\n<]*<\/[^\n>]*>)[^\n,<]*)*),.*?(\d+)(?=[^\n\d<>]*(?:<(?:s>[^\n<>]*<\/s>|[^\n<>]+>)[^\n\d<>]*)*<\/h\d>)/;var OVERRIDE_REG=/^Override\s*header:\s*/i;function getAuthorName(a){return a.owner.display_name}function process(){var valid=[];answers.forEach(function(a){var body=a.body;a.comments.forEach(function(c){if(OVERRIDE_REG.test(c.body))body='
'+c.body.replace(OVERRIDE_REG,'')+'
'});var match=body.match(SCORE_REG);if(match)valid.push({user:getAuthorName(a),size:+match[2],language:match[1],link:a.share_link,});else console.log(body)});valid.sort(function(a,b){var aB=a.size,bB=b.size;return aB-bB});var languages={};var place=1;var lastSize=null;var lastPlace=1;valid.forEach(function(a){if(a.size!=lastSize)lastPlace=place;lastSize=a.size;++place;var answer=jQuery("#answer-template").html();answer=answer.replace("{{PLACE}}",lastPlace+".").replace("{{NAME}}",a.user).replace("{{LANGUAGE}}",a.language).replace("{{SIZE}}",a.size).replace("{{LINK}}",a.link);answer=jQuery(answer);jQuery("#answers").append(answer);var lang=a.language;lang=jQuery(''+lang+'').text();languages[lang]=languages[lang]||{lang:a.language,lang_raw:lang.toLowerCase(),user:a.user,size:a.size,link:a.link}});var langs=[];for(var lang in languages)if(languages.hasOwnProperty(lang))langs.push(languages[lang]);langs.sort(function(a,b){if(a.lang_raw>b.lang_raw)return 1;if(a.lang_raw<b.lang_raw)return-1;return 0});for(var i=0;i<langs.length;++i){var language=jQuery("#language-template").html();var lang=langs[i];language=language.replace("{{LANGUAGE}}",lang.lang).replace("{{NAME}}",lang.user).replace("{{SIZE}}",lang.size).replace("{{LINK}}",lang.link);language=jQuery(language);jQuery("#languages").append(language)}}
body{text-align:left!important}#answer-list{padding:10px;width:290px;float:left}#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">

Shortest Solution by Language

Language	User	Score

Leaderboard

	Author	Language	Size

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

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

Some existing answers are assuming that 0 <= q < n. You should probably clarify whether or not this is an acceptable assumption.

I would have liked q and n to be any two integers, but so I don't break existing answers, 0 <= q < n

In this case I would have considered it reasonable to "break" the existing answers on the grounds that they weren't following the existing spec and you were just clarifying that it meant what it said rather than changing it, but it's too late now.

You  could give a small bonus for solutions accepting arbitrary q

Pyth, 9 bytes
}Em%*ddQQ

Try it online: Demonstration or Test Suite
Explanation:
}Em%*ddQQ   implicit: Q = first input number
  m     Q   map all numbers d of [0, 1, ..., Q-1] to:
    *dd       d*d
   %   Q      mod Q
            this gives the list of all quadratic residues
 E          read another input number
}           check, if it appears in the list of quadratic residues

I tried putting 7 9 as inputs, and it said "False", despite the fact that 7 is equivalent to 5^2 mod 9.

@kundor I read the integers in the reversed order. First n and than q. So try 9\n7 as input.

Mathematica, 25 bytes
AtomQ@PowerMod[#,1/2,#2]&

Mathematica, being Mathematica, naturally has a builtin for calculating modulo nth roots, via PowerMod. If a solution exists the smallest feasible solution is returned, otherwise the original expression (plus a message).
To get an actual truthy/falsy output we pass the result to AtomQ, which checks whether an expression can be broken down. Integers are atomic, returning True, whilst the non-atomic PowerMod[q,1/2,n] returns False
Thanks to @MartinBüttner for golf tips and function hunting with me.

Stupid argument ordering

What?! I never knew PowerMod could take a fractional argument!

Par, 11 9 bytes
✶X[²x%)↔,

Each character uses just one byte; see here.
Explanation
✶              ## Read two numbers
X              ## Assign second to x
[              ## Map
 ²             ## Square
 x%            ## Mod x
)              ##
↔              ## Swap
,              ## Count

Removed two bytes thanks to Jakube.

LabVIEW, 16 15 Equivalent bytes
Counted according to my meta post.

Matlab, 29
This function squares all numbers from 0 to n and checks whether a square minus q is zero mod n.
@(q,n)any(~mod((0:n).^2-q,n))

Haskell, 31 bytes
Saved 3 bytes thanks to Martin Büttner.
q#n=elem q[mod(x^2)n|x<-[1..n]]

Also 31 bytes: q#n=any(\x->mod(x*x)n==q)[0..n] and for 30 byte: q#n=[x|x<-[0..n],mod(x*x)n==q] which returns the list of x / empty list instead of True / False.

Prolog (SWI), 34 bytes
Code:
Q*N:-between(0,N,X),X*X mod N=:=Q.

Explanation:
Checks if any square between 0 and N leaves Q when divided by N.
Example:
3*8.
false

15*22.
true

Try it online here

CJam, 11 bytes
{_,2f#\f%&}

This unnamed block expects q n on the stack and leaves [q] on the stack as a truthy value or "" as a falsy value.
Test it here.
Credits to Sp3000 who also came up with this solution but "couldn't be bothered posting".
Explanation
_,  e# Duplicate n and turn into range [0 1 ... n-1]
2f# e# Square each element in the range.
\f% e# Take each element in the range modulo n.
&   e# Set intersection with q to check if any square yields q (mod n).

J, 9 bytes
e.]|i.^2:

Usage:
   1 (e.]|i.^2:) 5
1
   3 (e.]|i.^2:) 8
0
   15 (e.]|i.^2:) 22
1

Explanation:
e.]|i.^2:
    i.    [0..N-1]
      ^   to the power of
       2: 2 (constant 2 function)
  ]|      mod N
e.        contains Q? (0/1 result)

Some J mechanics trivia:
Functions are grouped by 3 iteratively from the right and if there is one left, as in our case (e. (] | (i. ^ 2:))), the grouped part is called with the right argument (N) and the left out function (e., "contains") called with the original left argument (Q) and the result of the grouped part.
(e.]|i.*i. and e.]|2^~i. also solves the problem with the same length.)
Try it online here.

Mathematica, 27 bytes
PowerModList[#,1/2,#2]!={}&

Usage:
In[1]:= PowerModList[#,1/2,#2]!={}&[1,5]

Out[1]= True

Javascript ES6, 42 bytes
(q,n)=>[...Array(n)].some((x,y)=>y*y%n==q)

Credits to @apsilers for serious bytes saved!

What is the ...Array syntax? I still don't get it.

Hope this edit is better for you.

[...Array(5)] produces array of undefined, same as Array(5) alone. I am totally confused, because removing the weird syntax and using just Array(5) breaks the code. Is there any documentation for this? my console screenshot

This might help. Although I did explain things incorrectly before :P. Basically, map won't work without the spread operator in dealing with arrays with empty slots.

@TomášZato Array(5) is an array with no own-properties except length. The array-spread [...x] fills in missing numerical properties up to length. The map function can only operate on extant properties, which Array(5) alone doesn't have. For example, try Array(5).hasOwnProperty(0) (false) versus [...Array(5)].hasOwnProperty(0) (true).

You can save 2 bytes if you do ~... instead of ...>-1. The bitwise NOT operator returns the falsey value 0 only if its input is -1; otherwise, its output is truthy.

Also, using some is shorter and (I think) equivalent: (q,n)=>[...Array(n)].some((x,y)=>y*y%n==q)

Ruby, 33 31 bytes
Saved two bytes thanks to Vasu Adari.
->q,n{(1..n).any?{|e|e*e%n==q}}

As usual Ruby's not going to beat any of the golfing languages, but it makes a good showing here.

You can lose the braces and make it ->q,n{}.

@VasuAdari Cool, I did not know that. Thank you.

Seriously, 20 bytes
,;R@,;╗@%╝`ª╜@%╛=`MΣ

Takes input on two lines: q, then n. Outputs a 0 if q is not a quadratic residue of n, else a positive number representing how many x in [1, q] (inclusive) satisfy x^2 = q (mod n).
Try it online (permalinks are having more issues, but you can copy and paste the code into a blank page in the meantime)
Explanation:
,;R      get q input, duplicate, push range(1, q+1)
@,;╗     move the list to the back of the stack, get n input, dupe, save in reg 0
@%╝      calculate q mod n and save to reg 1
`ª╜@%╛=` push this function:
  ª╜@%     square top of stack, push reg 0 value (n), swap, and mod
  ╛=       push reg 1 value (q mod n), compare equality (1 if equal else 0)
MΣ       map the function across the range, add results

Python 3, 41 40 bytes
Takes q and n and determines if q is in a list of squares from 0 squared to n-1 squared.
lambda q,n:q in[i*i%n for i in range(n)]

JavaScript (ES6), 43 bytes
(q,n)=>eval('for(x=n,r=0;x--;)r+=x*x%n==q')

Explanation
(q,n)=>
  eval(`              // eval allows us to use a for loop without {} or return
    for(x=n,r=0;x--;) // iterate over all possible values of x
      r+=x*x%n==q     // r = the number of matching x values
  `)                  // implicit: return r

Test

q = <input type="number" id="Q" />

n = <input type="number" id="N" />

<button onclick="result.innerHTML=(

(q,n)=>eval('for(x=n,r=0;x--;)r+=x*x%n==q')

)(+Q.value,+N.value)">Go</button>

This is a very interesting take on the truthy/falsey condition @user81655. Excellent work!

Japt, 10 bytes
Vo d_²%V¥U

My first official Japt golf ever! Thanks to @ETHProductions for saving a byte!
Ungolfed / Explanation
Vo d_  ²  %V¥ U
Vo dZ{Zp2 %V==U}  // implicit: U,V = inputs
Vo                // Create a range from 0 to n-1
   dZ{         }  // Check if any element Z in the range satisfies the condition:
       Zp2        // Is Z squared...
           %V     // modulo n...
             ==U  // equal to q?
                  // implicit output

Try it online!

Nice! Hint: 0oV is equivalent to Vo.

Didn't know that. Thanks!

Julia, 30 bytes
f(q,n)=q∈[i^2%n for i=0:n-1]

This is a function f that accepts two integers and returns a boolean.
Ungolfed:
function f(q::Integer, n::Integer)
    # Generate an array of quadratic residues
    x = [i^2 % n for i = 0:n-1]

    # Determine whether q is one of these values
    return q ∈ x
end

, 13 chars / 25 bytes
⟦Ѧí]ĉ⇀_²%í≔î)

Try it here (Firefox only).

I tested your code with 15 22 and it said false.

@Sherlock9 Fixed.

No custom codepage? This is not a golfing language!

There is now, but the code page was made long after the challenge.

C# 6 (.Net Framework 4.6) in LinqPad, 60 Bytes
bool b(int q,int n)=>Enumerable.Range(1,n).Any(y=>y*y%n==q);

Milky Way 1.0.2, 41 bytes
:>&{~1-:2h<:>n>;:>;<issquare(Mod(q,n))

Try it online!

JavaScript (Node.js), 33 bytes
q=>n=>f=(t=n)=>t&&f(t-1)|t*t%n==q

Try it online!
Why no one do this

C (gcc), 49 48 bytes
-1 byte thanks to ceilingcat
f(q,n,i,a){for(i=n;i--;)a=i*i%n==q?:a;return a;}

Try it online!
Does a simple for loop, brute forcing the squares.


text:
Score an asynchronous prisoner's dilemma game
In a round of the prisoner's dilemma exercise, two players each decide whether to cooperate or defect for that round. The scoring for a round is:

Player A and Player B both cooperate: 1 point for both
Player A and Player B both defect: 2 points for both
Player A cooperates and Player B defects: 3 points for cooperating Player A and 0 points for defecting Player B

You don't need to worry about strategy, though: your program will merely be tabulating the score for a game. (In case you are already familiar with prisoner's dilemma, my "points" here correspond to "years in prison.")
Your challenge is to take input that represents the players' choices over several rounds and compute their respective total scores. One player submits choices in lowercase, c and d (for cooperate and defect), and the other submits choices in uppercase, C and D. These choices are supplied to your program as a string.
Normally, players in prisoner's dilemma submit their moves simultaneously and iteratively. In this challenge, however, the players may have submitted their choices for several rounds at once. If a player's move is out of sequence, the scoring program remembers it and matches it with the next available move from the opposing player.
Here's a sample input string:
cDCddDDCcCc

To show the matches that exist in this input, I'll call out the lowercase and uppercase separately, and pair them up:
cDCddDDCcCc
c  dd   c c => cddcc
 DC  DDC C  => DCDDCC

These will be paired into the rounds:
c vs D (3 pts for lowercase-player, 0 pts for uppercase-player)
d vs C (0 pts for lowercase-player, 3 pts for uppercase-player)
d vs D (2 pts for both)
c vs D (3 pts for lowercase-player, 0 pts for uppercase-player)
c vs C (1 pt for both)

Which produces the score 9 (lowercase) to 6 (uppercase), so the output should be 9,6 (or any unambiguous delimiter).
To express it in yet another way, here's each pairing pulled out on its own row:
cDCddDDCcCc
cD
  Cd
    dD
      D c
       C  c

There is one unmatched C, because the uppercase player submitted more moves than the lowercase player. That is acceptable, and it is totally ignored for scoring purposes.
Here are the requirements:

You must write a program or function that accepts a string of the regular-expression form /[cdCD]+/, through some input mechanism (STDIN, function argument, read from file, etc.). (Your program may optionally accept input with a trailing newline.)
Your program or function must output or return the players' scores as a string. The output format must start with the lowercase player's score, followed by the uppercase player's score, separated by any non-empty, non-numeric delimiter of your choice. (A trailing newline is optional.)
If one player has more moves than the other, the excess moves are ignored.
If all moves in the input are from exclusively one player (that is, no rounds have been played at all), then each player's score is 0.
The smallest submission in bytes wins.

Test cases
Input:  cDCddDDCcCc
Output: 9,6         -- or any delimiter; I chose commas here

Input:  cccDDD
Output: 9,0

Input:  DDDDDDccc
Output: 9,0

Input:  cDcDcD
Output: 9,0

Input:  dcDDC
Output: 5,2

Input:  CcdCDDcd
Output: 6,6

Input:  Ddd
Output: 2,2

Input:  ccccccccccc
Output: 0,0

dont they normally get 2 points for cooperation and lose 1 point if both defect?

@Eumel I just copied the specs from the Wikipedia introduction, which appears to use a formulation suggested by the original authors. Also note that points here are "bad" as they correspond to years in prison. The winner is the player with the fewest points.

Is (0,0) or [0,0] OK for output?

Pyth, 23 bytes
jsMc2/L`C,@Gz-zG"cDDCdd

Test suite

Explanation:
@Gz: Lowercase letters
-zG: Uppercase letters
C,: Pair, truncate the remainder.
`: Take the string representation of the list of pairs
/L ... "cDDCdd: For each letter in "cDDCdd", count how many time it appears in the above string repr.
c2: Chop the resulting list in half.
sM: Add up each half.
j: Join on newlines and print.

` must be used instead of s to make the case where one party never plays work.

Haskell, 139 134 bytes
g=filter
(n!m)(a,b)=(a+n,b+m)
f s=init$tail$show$foldr id(0,0)$zipWith(#)(g(>'a')s)$g(<'E')s
'c'# 'C'=1!1
'c'#_=3!0
_# 'D'=2!2
_#_=0!3

Usage example: f "cDCddDDCcCc" ->  "9,6"
15 bytes just to get the output format right, i.e. turning a pair of numbers (x,y) into a string "x,y".
How it works:
               g(>'a')s        -- extract all lowercase letters
                     g(<'E')s  -- extract all uppercase letters
         zipWith(#)            -- combine both lists element wise with function #
                               -- # calls ! depending on the combination of c/d/C/D
                               -- ! takes 2 numbers a and b and returns a function
                               -- that takes a pair (x,y) and returns (x+a,y+b)
                               -- now we have a list of such functions
    foldr id(0,0)              -- apply those functions starting with (0,0)
init$tail$show                 -- format output

Edit: @Zgarb helped me saving 5 bytes. Thanks!

LabVIEW, 77 Bytes

The code scans from tokens and uses those indicies to decide where the points go.
Counting goes like this

Python 3, 110
Saved 5 Bytes thanks to FryAmTheEggman.
Saved 7 Bytes thanks to apsillers.
Saved 26 Bytes thanks to DSM.
x=[[],[]]
a=b=0
for m in input():x[m<'E']+=m
for w,p in zip(*x):d=p>'C';c=w<'d';b+=d*2+c;a+=3-d-2*c
print(b,a)

I think it's finally all golfed out.
It scans through each character in the input and sorts it based on being uppercase or not. It then does some fancy math that abuses Python's implicit conversion of bools to ints.

JavaScript (ES6), 124 118 bytes
s=>(A=B=i=0,U=(r=x=>s.replace(/c|d/g,x))``,r(l=>U[i]&&(U[i++]<'D'?l<'d'?++A&++B:B+=3:l<'d'?A+=3:(A+=2,B+=2))),A+','+B)

Live demo
(Expanded slightly for readability.)

var f=function (s) {
    A=B=i=0;
    U=(r=function(x){return s.replace(/c|d/g,x)})("");
    r(l=>U[i]&&(U[i++]<'D'?l<'d'?++A&++B:B+=3:l<'d'?A+=3:(A+=2,B+=2)));
    return A+','+B;
}

var input = ["cDCddDDCcCc","cccDDD","DDDDDDccc","cDcDcD","dcDDC","CcdCDDcd","Ddd","ccccccccccc"];
var output = ["9,6","9,0","9,0","9,0","5,2","6,6","2,2","0,0"];
var passed = true;

for (var index=0;index<input.length;index++) {
    if (f(input[index]) !== output[index]) passed = false;
}

document.getElementById("result").textContent =
  passed ? "All tests passed." : "Some tests failed.";

Saved 6 bytes thanks to user81655.

I originally had array comprehensions, but ended up using a different method. Thanks.

Par, 49 bytes
(lW▼·L)w▼·U))t˅y])[h7%Z2*↓″4>5*-]z2↔-″0<4*+╞)t.Σ¡

One byte is used per character. See here.
Explanation
(              ## Construct array
 l             ## Read line
 W             ## Assign to w
 ▼·L)          ## Filter by immutable under lower-case
 w             ## Get w
 ▼·U)          ## Filter by immutable under upper-case
)              ##
t              ## Transpose and truncate
˅y])           ## If empty, empty 2-D matrix
[              ## Map
 h             ## Decimal to hex
 7%            ## Modulo 7
 Z             ## Assign to z
 2*↓″4>5*-     ## Score of lower case
 ]             ## Put in array
 z2↔-″0<4*+    ## Score of upper case
 ╞             ## Add to array
)              ##
t              ## Transpose and truncate
.Σ             ## Map - sum
¡              ## Empty array onto stack

Outputs in the form 9 6.

As someone who has never used (or heard of) Par, I found your explanation a pleasure to read. Thanks!

CJam, 92 83 81 bytes
This ended up longer than I thought it would...
0]K*X3tC30tG22tZ11t:L;0'a]q+{'D>}:B$_{B}%1#/z{,1>},{2<[:i:#K%L=]sY0e[{si}%}%:.+S*

Try it here.
Explanation (Dare I explain this? :O):
0]K*C3tX30tG22tZ11t:L;    e# Creates this array [0,30,0,11,0,0,0,0,0,0,0,0,3,0,0,0,22,0,0,0]
0'a]q+                    e# Creates an array that looks like [0, 'a', input string]
{'D>}:B$                  e# Sorts the array by if the int representation of each element is greater than the int value of the character 'D' (e.g. [0,C,D,a,c,d])
_{B}%1#/                  e# Finds the index of the first value in the array that is > 'D' and splits the array at that index.
z{,1>},{                  e# Zip the two sub arrays and filter for only sub arrays with more than one element. (e.g [[0,a],[C,c],[D,d]])
{2<[:i:#K%L=]s            e# For each sub array, take the first two elements, convert each to an it, calculate n=(x[0]^x[1]) mod 20, and get the nth element in the very first array, and convert it to a string
Y0e[                      e# Pad the string with 0 so it is length 2. (e.g. [["00"],["22"],["11"]])
{si}%}%:.+                e# get the numerical representation of each digit and dot sum all of them (e.g [[0,0],[2,2],[1,1] => [3,3])
S*                        e# Join with a space (e.g "3 3")

source: Stack Exchange

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

sort: votes


text:
Calculate the Kronecker symbol
Relevant links here and here, but here is the short version:
You have an input of two integers \$a\$ and \$b\$ between negative infinity and infinity (though if necessary, I can restrict the range, but the function must still accept negative inputs).
Definition of the Kronecker symbol
You must return the Kronecker symbol \$(a|b)\$ for inputs \$a\$ and \$b\$ where
$$(a|b) = (a|p_1)^{e_1} \cdot (a|p_2)^{e_2} \cdots (a|p_n)^{e_n}$$
where \$b = p_1^{e_1} \cdot p_2^{e_2} \cdots p_n^{e_n}\$, and \$p_i\$ and \$e_i\$ are the primes and exponents in the prime factorization of \$b\$.

For an odd prime \$p\$, \$(a|p)=a^{\frac{p-1}2} \mod p\$ as defined here.

For \$b = 2\$,

$$
(n|2) = \begin{cases}
0    & \text{if } n \equiv   0   \mod 2 \\
1    & \text{if } n \equiv \pm 1 \mod 8 \\
-1   & \text{if } n \equiv \mp 3 \mod 8
\end{cases}
$$

For \$b = -1\$,

$$
(n|-1) = \begin{cases}
-1 & \text{if } n  <  0 \\
1  & \text{if } n \ge 0
\end{cases}
$$

If \$a \ge b\$, \$(a|b) = (z|b)\$ where \$z = a\mod b\$. By this property, and as explained here and here, \$a\$ is a quadratic residue of \$b\$ if \$z\$ is, even though \$a \ge b\$.

\$(-1|b) = 1\$ if \$b \equiv 0,1,2 \mod 4\$ and \$-1\$ if \$b = 3 \mod 4\$.

\$(0|b)\$ is \$0\$ except for \$(0|1)\$ which is \$1\$, because \$(a|1)\$ is always \$1\$ and for negative \$a\$, \$(-a|b) = (-1|b) \times (a|b)\$.

The output of the Kronecker symbol is always \$-1, 0 \text{ or } 1\$, where the output is \$0\$ if \$a\$ and \$b\$ have any common factors. If \$b\$ is an odd prime, \$(a|b) = 1\$ if \$a\$ is a quadratic residue \$\mod b\$, and \$-1\$ if is it is not a quadratic residue.
Rules

Your code must be a program or a function.

The inputs must be in the order a b.

The output must be either \$-1, 0 \text{ or } 1\$.

This is code golf, so your code does not have to be efficient, just short.

No built-ins that directly calculate the Kronecker or the related Jacobi and Legendre symbols. Other built-ins (for prime factorization, for example) are fair game.

Examples
>>> kronecker(1, 5)
1
>>> kronecker(3, 8)
-1
>>> kronecker(15, 22)
1
>>> kronecker(21, 7)
0
>>> kronecker(5, 31)
1
>>> kronecker(31, 5)
1
>>> kronecker(7, 19)
1
>>> kronecker(19, 7)
-1
>>> kronecker(323, 455625)
1
>>> kronecker(0, 12)
0
>>> kronecker(0, 1)
1
>>> kronecker(12, 0)
0
>>> kronecker(1, 0)
1
>>> kronecker(-1, 5)
1
>>> kronecker(1, -5)
1
>>> kronecker(-1, -5)
-1
>>> kronecker(6, 7)
-1
>>> kronecker(-1, -7)
1
>>> kronecker(-6, -7)
-1

This is a complicated function, so please let me know if anything is unclear.

Are you sure you don't want to disallow built-ins? https://reference.wolfram.com/language/ref/KroneckerSymbol.html

@MartinBüttner I was editing in examples when I saw your comment. I will disallow built-ins that directly calculate the Kronecker, Jacobi or Legendre symbols, but anything else (including prime factorization functions) should be fair game.

im not entirely sure why (31|5) gives 1. There shouldnt be a qudratic residue so why isnt it -1?

also 7/19 should be 1 and 19/7 should be -1 according to the wiki you linked

@Eumel 31 is a quadratic residue mod 5, because it is congruent to 1 mod 5. I have edited the question to explain this.

@Eumel Ah, you are correct. I'll edit the examples.

If solutions have to handle negative and zero inputs correctly, you should definitely add some test cases for that.

I think kronecker(0, 12) is 0.

I don't think you've actually defined the symbol for a < 0.

Half-way there. If you now define (-1|b) I think the spec will be complete.

Only very rarely does PPCG make me feel very, very stupid, as if I were browsing Math.SE, but this is one of those rare cases

CJam (70 bytes)
{_g\zmf+f{:P2+"W>2*(
z1=
;1
7&4-z[0W0X0]=
P%P+P(2/#P%_1>P*-"N/<W=~}:*}

Online demo (test cases generated with Mathematica).
Dissection
{               e# Anonymous function. Stack: a b
  _g\zmf+       e# Factorise b, with special treatment for negatives
                e# CJam also gives special treatment to 0 and 1
                e# Stack: e.g. a [-1 2 2 5]; or a [-1 1]; or a [0 0]; or a [1 2 2 5]
  f{            e# For each "prime" factor P, find (a|P)
    :P2+        e# Extract code for P from an array formed by splitting a string
    "W>2*(      e#   a -> (a|-1)
z1=             e#   a -> (a|0)
;1              e#   a -> (a|1)
7&4-z[0W0X0]=   e#   a -> (a|2)
P%P+P(2/#P%_1>P*-" e# a -> (a|P) for odd prime P
    N/<W=~      e# Split string and select appropriate element
  }
  :*            e# Multiply the components together
}

I found several ways of evaluating (a|2) for the same character count, and have chosen to use the one with the clearest presentation.
integer array <W= is IMO a quite elegant way of doing fallbacks: if the integer is greater than the length of the array, we select the last element.
Other comments
It's disappointing that for odd prime p the direct Fermat-style (a|p) is so short, because there's a very golfy way of finding (a|n) for positive odd n which I wanted to use. The basis is Zolotarev's lemma:

If p is an odd prime and a is an integer coprime to p then the Legendre symbol (a|p) is the sign of the permutation x -> ax (mod p)

This was strengthened by Frobenius to

If a and b are coprime positive odd integers then the Jacobi symbol (a|b) is the sign of the permutation x -> ax (mod b)

and by Lerch to

If b is a positive odd integer and a is an integer coprime to b then the Jacobi symbol (a|b) is the sign of the permutation x -> ax (mod b)

See Brunyate and Clark, Extending the Zolotarev-Frobenius approach to quadratic reciprocity, The Ramanujan Journal 37.1 (2014): 25-50 for references.
And it can easily be strengthened one step further (although I haven't seen this in the literature) to

If b is a positive odd integer and a is an integer then the Jacobi symbol (a|b) is the Levi-Civita symbol of the map x -> ax (mod b).

Proof: if a is coprime to b then we use Zolotarev-Frobenius-Lerch; otherwise the map is not a permutation, and the Levi-Civita symbol is 0 as desired.
This gives the Jacobi symbol calculation
{_2m*{~>},@ff*\ff%::-:*g}

But the special treatment required for (a|-1) and (a|2) mean that I haven't found a way of calculating the Kronecker symbol which is shorter with this approach: it's shorter to factorise and treat the primes individually.

Python 3, 747 369 335 bytes
As an example answer, only slightly golfed, and to give you an idea of what an answer will look like.
And yes, the prime factorization and run-length-encoding bits are cribbed from Pyth with apologies to isaacg.
from itertools import*
def k(d,r):
 if d<0:a=-d;m=1
 else:a=d;m=0
 if r==1:return 1
 p=1;w=r;n=2;f=[]
 while n*n<=w:
  while w%n<1:w//=n;f+=n,
  n+=1
 if w>1:f+=w,
 z=[[k,len(list(g))]for k,g in groupby(f)]
 for i,j in z:
  if i==2:p*=pow(-1,(a*a-1)//8)
  x=pow(a,(i-1)//2,i)
  if x>1:x-=i
  p*=x**j
 if m:p*=pow(-1,(r-1)//2)
 return p

Apology accepted - I'm glad someone reads Pyth source code.

LabVIEW, 44 Bytes LabVIEW Primitives
Since its symetrical i swapped the inputs if a was bigger than b.
Represents the real formula now
counting like always according to

for true case

Unfortunately, (a|b) != (b|a) in all cases. In most cases, yes, but not in all of them. Though it would work if you reduced a mod b instead of swapping them.

since i have the explantion now i can edit it, give me a min

Is there any way I can test this? I don't really understand how LabView works.

that is a good question, i can think of 2 ways. First i can build an .exe and send it to you, second you can get a labview test version and i can send you the vi or you can rebuild it from the pic.

This is not 44 bytes. If you define a scoring system that is not based on the size of the file, you should call it something other than bytes.

For the moment, I will take your word for it. Send me the .exe anyway, but you have my upvote.

I thought there was a messanging system here but there isnt, how do i actually send it to you?

@Eumel Seeing as you put up a gif for some of your other LabVIEW solutions, can you put one up here? And change the heading from bytes to LabVIEW primitives?

will do when im at home

The gif appears to be for a different version to the image above. And if I'm correctly interpreting the symbols, this seems to be incapable of giving any output other than -1 or 1, which means that it must be giving incorrect outputs for some inputs.

And from looking at a reference, it seems to have some case structures for which the image only shows one case, so you haven't posted the full code.

The other case is just a 0 nothing more ill put it up though

Mathematica, 169 175 165 bytes
(1|-1)~k~0=_~k~1=1
_~k~0=0
a_~k~-1=If[a<0,-1,1]
a_~k~2=DirichletCharacter[8,2,a]
a_~k~p_/;PrimeQ@p=Mod[a^((p-1)/2),p,-1]
a_~k~b_:=1##&@@(a~k~#^#2&@@@FactorInteger@b)

Julia, 195 bytes
k(a,b)=b==0?a∈[1,-1]?1:0:b==1?1:b==2?iseven(a)?0:a%8∈[1,-1]?1:-1:b==-1?a<1?-1:1:isprime(b)&&b>2?a%b==0?0:a∈[i^2%b for i=0:b-1]?1:-1:k(a,sign(b))*prod(i->k(a,i)^factor(b)[i],keys(factor(b)))

This is a recursive function k that accepts two integers and returns an integer.
Ungolfed:
function k(a::Integer, b::Integer)
    if b == 0
        return a ∈ [1, -1] ? 1 : 0
    elseif b == 1
        return 1
    elseif b == 2
        return iseven(a) ? 0 : a % 8 ∈ [1, -1] ? 1 : -1
    elseif b == -1
        return a < 1 ? -1 : 1
    elseif isprime(b) && b > 2
        return a % b == 0 ? 0 : a ∈ [i^2 % b for i = 1:b-1] ? 1 : -1
    else
        p = factor(b)
        return k(a, sign(b)) * prod(i -> k(a, i)^p[i], keys(p))
    end
end

Haskell, 286 bytes
a#0|abs a==1=1|1<2=0
a#1=1
a#2|even a=0|mod a 8`elem`[1,7]=1|1<2=(-1)
a#b|b<0=a`div`abs a*a#(-b)|all((/=0).mod b)[2..b-1]=if elem n[0,1] then n else(-1)|1<2=product$map(a#)$f b where n=a^(div(b-1)2)`mod`b
f 1=[]
f n|n<0=(-1):f(-n)|1<2=let p=head$filter((==0).mod n)[2..n]in p:f(div n p)

Probably not completely optimized, but a valiant effort. The Kronecker symbol is defined as the infix function a#b, i.e.
*Main>323#455265
1

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

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


text:
Bernoulli Numbers
The Bernoulli numbers (specifically, the second Bernoulli numbers) are defined by the following recursive definition:

Where  denotes a combination.
Given a nonnegative integer m as input, output the decimal representation OR a reduced fraction for the mth second Bernoulli number. If you output a decimal representation, you must have at least 6 decimal places (digits after the decimal point) of precision, and it must be accurate when rounded to 6 decimal places. For example, for m = 2, 0.166666523 is acceptable because it rounds to 0.166667. 0.166666389 is not acceptable, because it rounds to 0.166666. Trailing zeroes may be omitted. Scientific notation may be used for decimal representations.
Here is the input and expected output for m up to and including 60, in scientific notation rounded to 6 decimal places, and as reduced fractions:
0 -> 1.000000e+00 (1/1)
1 -> 5.000000e-01 (1/2)
2 -> 1.666667e-01 (1/6)
3 -> 0.000000e+00 (0/1)
4 -> -3.333333e-02 (-1/30)
5 -> 0.000000e+00 (0/1)
6 -> 2.380952e-02 (1/42)
7 -> 0.000000e+00 (0/1)
8 -> -3.333333e-02 (-1/30)
9 -> 0.000000e+00 (0/1)
10 -> 7.575758e-02 (5/66)
11 -> 0.000000e+00 (0/1)
12 -> -2.531136e-01 (-691/2730)
13 -> 0.000000e+00 (0/1)
14 -> 1.166667e+00 (7/6)
15 -> 0.000000e+00 (0/1)
16 -> -7.092157e+00 (-3617/510)
17 -> 0.000000e+00 (0/1)
18 -> 5.497118e+01 (43867/798)
19 -> 0.000000e+00 (0/1)
20 -> -5.291242e+02 (-174611/330)
21 -> 0.000000e+00 (0/1)
22 -> 6.192123e+03 (854513/138)
23 -> 0.000000e+00 (0/1)
24 -> -8.658025e+04 (-236364091/2730)
25 -> 0.000000e+00 (0/1)
26 -> 1.425517e+06 (8553103/6)
27 -> 0.000000e+00 (0/1)
28 -> -2.729823e+07 (-23749461029/870)
29 -> 0.000000e+00 (0/1)
30 -> 6.015809e+08 (8615841276005/14322)
31 -> 0.000000e+00 (0/1)
32 -> -1.511632e+10 (-7709321041217/510)
33 -> 0.000000e+00 (0/1)
34 -> 4.296146e+11 (2577687858367/6)
35 -> 0.000000e+00 (0/1)
36 -> -1.371166e+13 (-26315271553053477373/1919190)
37 -> 0.000000e+00 (0/1)
38 -> 4.883323e+14 (2929993913841559/6)
39 -> 0.000000e+00 (0/1)
40 -> -1.929658e+16 (-261082718496449122051/13530)
41 -> 0.000000e+00 (0/1)
42 -> 8.416930e+17 (1520097643918070802691/1806)
43 -> 0.000000e+00 (0/1)
44 -> -4.033807e+19 (-27833269579301024235023/690)
45 -> 0.000000e+00 (0/1)
46 -> 2.115075e+21 (596451111593912163277961/282)
47 -> 0.000000e+00 (0/1)
48 -> -1.208663e+23 (-5609403368997817686249127547/46410)
49 -> 0.000000e+00 (0/1)
50 -> 7.500867e+24 (495057205241079648212477525/66)
51 -> 0.000000e+00 (0/1)
52 -> -5.038778e+26 (-801165718135489957347924991853/1590)
53 -> 0.000000e+00 (0/1)
54 -> 3.652878e+28 (29149963634884862421418123812691/798)
55 -> 0.000000e+00 (0/1)
56 -> -2.849877e+30 (-2479392929313226753685415739663229/870)
57 -> 0.000000e+00 (0/1)
58 -> 2.386543e+32 (84483613348880041862046775994036021/354)
59 -> 0.000000e+00 (0/1)
60 -> -2.139995e+34 (-1215233140483755572040304994079820246041491/56786730)

Reference implementation (in Python 3):
def factorial(n):
    if n < 1:
        return 1
    else:
        return n * factorial(n - 1)

def combination(m,k):
    if k <= m:
        return factorial(m)/(factorial(k) * factorial(m - k))
    else:
        return 0

def Bernoulli(m):
    if m == 0:
        return 1
    else:
        t = 0
        for k in range(0, m):
            t += combination(m, k) * Bernoulli(k) / (m - k + 1)
        return 1 - t

Rules

This is code-golf, so shortest code in bytes wins
You may not use any functions, either built-in or included in an external library, that compute either type of Bernoulli numbers or Bernoulli polynomials.
Your answer must give correct output for all inputs up to and including 60.

Leaderboard
The Stack Snippet at the bottom of this post generates the leaderboard from the answers a) as a list of shortest solution per language and b) as an overall leaderboard.
To make sure that your answer shows up, please start your answer with a headline, using the following Markdown template:
## Language Name, N bytes

where N is the size of your submission. If you improve your score, you can keep old scores in the headline, by striking them through. For instance:
## Ruby, <s>104</s> <s>101</s> 96 bytes

If there you want to include multiple numbers in your header (e.g. because your score is the sum of two files or you want to list interpreter flag penalties separately), make sure that the actual score is the last number in the header:
## Perl, 43 + 2 (-p flag) = 45 bytes

You can also make the language name a link which will then show up in the snippet:
## [><>](http://esolangs.org/wiki/Fish), 121 bytes

 <link rel="stylesheet" type="text/css" href="//cdn.sstatic.net/codegolf/all.css?v=83c949450c8b">

Shortest Solution by Language

Language	User	Score

Leaderboard

	Author	Language	Size

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

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

@MorganThrapp The reference implementation is only to make clearer the definition of a Bernoulli number, not to actually solve the problem.

Ah, gotcha. I thought it was a fully functional implementation.

@MorganThrapp I wrote the code, so it's far from fully functional. :P

@FryAmTheEggman `) No, because "decimal digits" refers to digits behind the decimal point. I'll clarify in the question. 2) Sure.

@Mego No standard float (not even quad precision) can store B_60 to quad precision. Should we then use an extended precision format if we're outputting as decimals?

@ThomasKwa A double (64-bit float) should be able to represent -2.139995e34 without much trouble. The way I read the rules, it would only really be a problem if you're not using scientific notation for the output. Getting 6 digits after the decimal point would indeed require a lot of precision if you don't write the number in scientific notation.

@ThomasKwa You should use whatever is necessary to get the output in the required precision.

I don't like the precision requirement. Some languages don't have the tools to work with floats to sufficient accuracy for B_60, and I'd rather not deal with such issues when golfing a mathematical problem. It's frustrating to write up a solution and then find that it's invalid due to what seems like a technicality.

There is a slight trap in that even if you use the exact denominators and integer rounding to reduce errors in the intermediate Bernoulli numbers to 1ulp, the naïve binomial-weighted sum has so much cancellation that things go wildly wrong.

@xnor That will always be an issue with most mathematics-based challenges. That does not mean they shouldn't be posted.

@Mego I'm not saying math challenges shouldn't be posted, just that I'd prefer more lax rules about precision.

@xnor 6 digits of accuracy seems incredibly lax already.

Julia, 23 20 bytes
Saved 3 bytes thanks to Alex A.
It uses the same formula as my Mathematica solution and PARI/GP solution.
n->n>0?-zeta(1-n)n:1

20 bytes: n->n>0?-zeta(1-n)n:1

@AlexA. I don't know why, but zeta(n) throws an error when n is a negative integer. I'm using Julia 0.2.1 on linux.

Oh my, your version of Julia is quite outdated. It works just fine for me on 0.4.1.

Mathematica, 40 28 23 22 bytes
Using the famous formula n*ζ(1−n)=−Bn, where ζ is the Riemann Zeta function.
If[#>0,-Zeta[1-#]#,1]&

The same length:
B@0=1;B@n_=-Zeta[1-n]n

Original solution, 40 bytes, using the generating function of Bernoulli numbers.
#!SeriesCoefficient[t/(1-E^-t),{t,0,#}]&

+2 if I could...

Julia, 58 bytes
B(m)=m<1?1:1-sum(k->big(binomial(m,k))*B(k)/(m-k+1),0:m-1)

This creates a recursive function B that accepts an integer and returns a BigFloat (i.e. high precision floating point).
Ungolfed:
function B(m::Integer)
    m == 0 && return 1
    return 1 - sum(k -> big(binomial(m, k)) * B(k) / (m-k+1), 0:m-1)
end

Minkolang 0.14, 97 bytes
I actually tried doing it recursively first, but my interpreter, as currently designed, actually can't do it. If you try to recurse from within a for loop, it starts a new recursion. So I went for the tabulating approach...which had precision problems. So I did the whole thing with fractions. With no built-in support for fractions. [sigh]
n1+[xxi$z0z2%1+F0c0=$1&$d4Mdm:1R:r$dz1Az0A]$:N.
11z[i0azi6M*i1azi-1+*d0c*1c2c*-1c3c*4$X]f
z1=d1+f

Try it here. Bonus: the array has all fractions for every previous Bernoulli number!
Explanation (in a bit)
n1+                 Take number from input (N) and add 1
   [                Open for loop that runs N+1 times (starts at zero)
    xx              Dump the top two values of the stack
      i$z           Store the loop counter in the register (m)
         0          Push 0
          z2%1+     Push 1 if N is even, 2 if odd
               F    Gosub; pops y,x then goes to codebox(x,y), to be returned to

    0c                                 Copy the first item on the stack
      ,                                1 if equal to 0, 0 otherwise
       $1&                             Jump 11 spaces if top of stack is not 0

                                       (If top of stack is not 0, then...)
          $d                           Duplicate whole stack
            4M                         Pop b,a and push GCD(a,b)
              dm                       Duplicate and merge (a,b,c,c -> a,c,b,c)
                :                      Divide
                 1R                    Rotate 1 item to the right (0G works too)
                   :                   Divide
                    r                  Reverse stack

                                       (In both cases...)
                     $d                Duplicate whole stack
                       z1A             Store denominator of B_m in array
                           z0A         Store numerator of B_m in array
                              ]        Close for loop
                               $:      Divide (float division)
                                 N.    Output as number and stop.

11                                           Push two 1s (a, b)
  z[                                         Open a for loop that repeats m times
    i0a                                      Retrieve numerator of B_k (p)
       zi                                    Push m, k
         6M                                  Pop k,m and push mCk (binomial) (x)
           *                                 p*x (c)
            i1a                              Retrieve denominator of B_k (q)
               zi-1+                         m-k+1 (y)
                    *                        q*y (d)
                     d                       Duplicate top of stack
                      0c*                    a*d
                         1c2c*               b*c
                              -              a*d-b*c
                               1c3c*         b*d
                                    4$X      Dump the bottom four items of stack
                                       ]f    Jump back to F

z          m
 1=        0 (if m is 1) or 1 (otherwise)
   d1+     Duplicate and add 1 (1 or 2)
      f    Jump back to F

The third line is responsible for 1/2 if m is 1 and 0/1 if m is an odd number greater than 1. The second line calculates B_m with the summation formula given in the question, and does so entirely with numerators and denominators. Otherwise it'd be a lot shorter. The first half of the first line does some bookkeeping and chooses whether to execute the second or third line, and the second half divides the numerator and denominator by their GCD (if applicable) and stores those values. And outputs the answer at the end.

Python 2, 118 bytes
Saved 6 bytes due to xsot.
Saved 6 10 more due to Peter Taylor.
n=input()
a=[n%4-1,n<2]*n;exec"a=[(a[j-1]+a[j+1])*j/2for j in range(len(a)-2)];"*~-n
print+(n<1)or-n/(2.**n-4**n)*a[1]

Uses the following identity:

where An is the nth Alternating Number, which can be formally defined as the number of alternating permutations on a set of size n, halved (see also: A000111).
The algorithm used was originally given by Knuth and Buckholtz (1967):

Let T1,k = 1 for all k = 1..n
Subsequent values of T are given by the recurrence relation:
Tn+1,k = 1/2 [ (k - 1) Tn,k-1 + (k + 1) Tn,k+1 ]
An is then given by Tn,1
(see also: A185414)

Python 2, 152 bytes
from fractions import*
n=input()
a=[n%4-1,n<2]*n
for k in range(n-1):a=[(a[j-1]+a[j+1])*j/2for j in range(n-k)]
print+(n<1)or Fraction(n*a[1],4**n-2**n)

Prints the exact fractional representation, necessary for values greater than 200 or so.

If you change range(2,n) to range(n-2) you can shorten n-k+1 to n+~k. Also, is there a reason you use >>1 instead of /2? Lastly, a trivial improvement, but you can save a few bytes by aliasing range.

Thanks for the suggestions. I originally had two expressions, when I joined them I overlooked changing >>1 with /2.

There's a one-char saving in the output line: print+(n<1)or-(-1.)**(n+n/2)*n/(4**n-2**n)*a[n%2^1%n]. And the calculation can be done for the same char count as a=[1]*(n+1);exec"a=[(a[j-1]+a[j+1])*j/2for j in range(len(a)-1)];"*(n-1)

@PeterTaylor the n+n/2 is clever; 1 doesn't need to be singled out, because all other odd values are zero anyway. The alternative calculation is actually 4 bytes shorter with bit inversions, but also considerably slower, for some reason.

@PeterTaylor It's interesting that the modified calculation even works. I got rid of Knuth's k-1 and k+1 by pre-multiplying each T_n,k by k, hence the initial vector [1,2,3,...] rather than [1,1,1,...]. I didn't notice that running an extra iteration would achieve the same.

I was working from the OEIS table, and thought that you'd found range and skipping one iteration to be a clever way to shorten the initialisation. The way you've now split out even and odd indices is very nice, and allows a further saving by pulling the sign into the definition of a: a=[(-1)**(n/2),n<2]*n. Then the return value is +(n<1)or-n/(2.**n-4**n)*a[1]. You've also got a stray semicolon at the end of line 2.

@PeterTaylor I've deviated from the table in A185414, but the first column is still the same, which is what matters. Moving (-1)**(n/2) out of the output is brilliant, and allows for the much shorter n%4-1.

CJam, 69 49 34 33 bytes
{_),:):R_:*_@f/@{_(;.-R.*}*0=\d/}

Online demo
Thanks to Cabbie407, whose answer made me aware of the Akiyama–Tanigawa algorithm.
Dissection
{           e# Function: takes n on the stack
  _),:)     e# Stack: n [1 2 3 ... n+1]
  :R        e# Store that array in R
  _:*       e# Stack: n [1 2 3 ... n+1] (n+1)!
  _@f/      e# Stack: n (n+1)! [(n+1)!/1 (n+1)!/2 (n+1)!/3 ... (n+1)!/(n+1)]
            e#   representing [1/1 1/2 ... 1/(n+1)] but avoiding floating point
  @{        e# Repeat n times:
    _(;.-   e#   Take pairwise differences
    R.*     e#   Pointwise multiply by 1-based indices
  }*        e#   Note that the tail of the array accumulates junk, but we don't care
  0=\d/     e# Take the first element and divide by (n+1)!
}

Multiplying through by n! to prevent precision loss is clever, if not slightly ridiculous. I wonder if the algorithm couldn't be refactored slightly to avoid this.

I don't think refactoring could avoid the need to scale for the simple reason that since we know that every other Bernoulli number is 0 there's obviously a lot of subtraction of similar values going on, so there's a lot of places where catastrophic loss of significance can occur.

Python 3, 112 bytes
Edit: I cleaned up this answer. If you want to see all the other ways I thought of to answer this question in Python 2 and 3, look in the revisions.
If I don't use the lookup table (and I use memoization instead), I manage to get the recursive definition to 112 bytes! WOO! Note that b(m) returns a Fraction. As usual, the byte count and a link for testing.
from fractions import*
def b(m):
 s=k=0;p=1
 while k<m:a=m-k;s+=Fraction(p*b(k))/-~a;p=p*a//-~k;k+=1
 return 1-s

And a function that does use a lookup table, and returns the whole table of fractions from b(0) to b(m), inclusive.
from fractions import*
def b(m,r=[]):
 s=k=0;p=1
 while k<m:
  if k>=len(r):r=b(k,r)
  a=m-k;s+=Fraction(p*r[k])/-~a;p=p*a//-~k;k+=1
 return r+[1-s]

I think you can omit the trailing zero on float literals, e.g. 1. instead of 1.0.

@AlexA. Done. Removed .0 from s entirely, because it will quickly become a float later.

Can you use p=v=1;exec('[...];p+=1'*k) instead of your innermost loop?

PARI/GP, 52 23 bytes
Using the famous formula n*ζ(1−n)=−Bn, where ζ is the Riemann Zeta function.
n->if(n,-n*zeta(1-n),1)

Original solution, 52 bytes, using the generating function of Bernoulli numbers.
n->n!*polcoeff(-x/sum(i=1,n+1,(-x)^i/i!)+O(x^n*x),n)

Can only upvote once. It is a shame it's not exact, though.

According to the documentation the zeta function is calculated using Bernoulli numbers, in fact.

@primo, yes, I regard all of the answers which use built-in zeta as cheating.

Even easier, bernfrac and bernreal are 8 bytes each and they're already functions, so no need for the n->. But +1 for a good solution.

Mathematica, 52 48 42 bytes
1-Sum[#~Binomial~k#0@k/(#-k+1),{k,0,#-1}]&

Unnamed function that uses the literal definition.

Is the Sign@# necessary?

I tested it on my computer. After removing the Sign@#, it still returns the correct answer for 0 .

PARI/GP, 45 bytes
n->if(n,2*n/(2^n-4^n)*real(polylog(1-n,I)),1)

Using the same formula as my Python answer, with An generated via polylog.

Test Script
Run gp, at the prompt paste the following:
n->if(n,2*n/(2^n-4^n)*real(polylog(1-n,I)),1)
for(i=0, 60, print(i, ": ", %(i)))

Thank you for providing a test script - it made testing this a lot easier!

@Mego for you and me both ;)

Python 2, 132 130 bytes
import math,fractions
f=math.factorial
B=lambda m:~-m*m%2or 1+sum(B(k)*f(m)/f(k)/f(m-k)/fractions.Fraction(k+~m)for k in range(m))

This is just a golfed version of the reference implementation.
This is a little slow in practice, but can be sped up significantly with memoization:
import math,fractions
f=math.factorial

def memoize(f):
 memo = {}
 def helper(x):
  if x not in memo:
   memo[x] = f(x)
  return memo[x]
 return helper

@memoize
def B(m):
 return~-m*m%2or 1+sum(B(k)*f(m)/f(k)/f(m-k)/fractions.Fraction(k+~m)for k in range(m))

for m in range(61):
 print(B(m))

You can try this version online on Ideone.

gawk4, 79 bytes
77 bytes code + 2 bytes for -M flag
PREC^=2{for(n=$0;m++<=n;)for($(j=m)=1/m;j>1;)$j=(-$j+$--j)*j;printf"%.6f",$1}

It's an implementation of the Akiyama–Tanigawa algorithm from the Wikipedia page.
Had some trouble with the "6-decimal-digits-rule", because this prints the whole number and then 6 digits, but there is no list here to compare the results.
A flaw is that this prints a minus sign in front of the 0.000000 al lot of times, but I don't think that's wrong.
Usage example
echo 58 | awk -M 'PREC^=2{for(n=$0;m++<=n;)for($(j=m)=1/m;j>1;)$j=(-$j+$--j)*j;printf"%.6f",$1}'

Output from 0 to 60

0 -> 1.000000
1 -> 0.500000
2 -> 0.166667
3 -> -0.000000
4 -> -0.033333
5 -> 0.000000
6 -> 0.023810
7 -> 0.000000
8 -> -0.033333
9 -> 0.000000
10 -> 0.075758
11 -> -0.000000
12 -> -0.253114
13 -> -0.000000
14 -> 1.166667
15 -> -0.000000
16 -> -7.092157
17 -> -0.000000
18 -> 54.971178
19 -> -0.000000
20 -> -529.124242
21 -> -0.000000
22 -> 6192.123188
23 -> 0.000000
24 -> -86580.253114
25 -> 0.000000
26 -> 1425517.166667
27 -> 0.000000
28 -> -27298231.067816
29 -> 0.000000
30 -> 601580873.900642
31 -> 0.000000
32 -> -15116315767.092157
33 -> 0.000000
34 -> 429614643061.166667
35 -> 0.000000
36 -> -13711655205088.332772
37 -> 0.000000
38 -> 488332318973593.166667
39 -> -0.000000
40 -> -19296579341940068.148633
41 -> -0.000000
42 -> 841693047573682615.000554
43 -> -0.000000
44 -> -40338071854059455413.076812
45 -> -0.000000
46 -> 2115074863808199160560.145390
47 -> -0.000000
48 -> -120866265222965259346027.311937
49 -> -0.000000
50 -> 7500866746076964366855720.075758
51 -> -0.000000
52 -> -503877810148106891413789303.052201
53 -> -0.000000
54 -> 36528776484818123335110430842.971178
55 -> -0.000000
56 -> -2849876930245088222626914643291.067816
57 -> -0.000000
58 -> 238654274996836276446459819192192.149718
59 -> -0.000000
60 -> -21399949257225333665810744765191097.392674

Would printf"%e" work?

No, it wouldn't, because the 0.00000s are only very small and not really zero.

R, 93 bytes
function(m){if(m==0){1}else{v=c();for(k in 0:(m-1))v=c(v,choose(m,k)*f(k)/(m-k+1));1-sum(v)}}

Not really original as a solution. If any comment, please feel free !
Ungolfed :
function(m)
    if(m==0){1}
    else
         v=c()
         for(k in 0:(m-1))
            v=c(v,choose(m,k)*f(k)/(m-k+1))

1-sum(v)

I know this is a bit late now but you can save 3 bytes by changing the if/else statement order and use m>0 as well as looping over 1:m-1 instead.

GolfScript, 63 bytes
~:i.!+.[3i&(2i>]*i(,{i\-,{1$1$(=2$2$)=+*2/}%\;}/~\2i?.(*\--1?**

Online Demo.
Using the same formula as my Python answer.

Test Script
61,{[.`
  ~:i.!+.[3i&(2i>]*i(,{i\-,{1$1$(=2$2$)=+*2/}%\;}/~\2i?.(*\--1?**
]p}/

The apphb link will time out on this. If you don't have GolfScript installed locally, I recommend using the anarchy golf interpreter (use form, select GolfScript, paste, submit).

Perl, 101 bytes
#!perl -p
@a=($_%4-1,$_<2)x$_;
@a=map$_*($a[$_-1]+$a[$_+1])/2,0..@a-3for 2..$_;
$_=!$_||$_/(4**$_-2**$_)*$a[1]

Counting the shebang as three, input is taken from stdin.
Using the same formula as my Python answer.

Sample Usage
$ echo 60 | perl bernoulli.pl
-2.13999492572253e+034

Online Demo.

Ruby, 66 61 bytes
This is a Ruby version of my Python answer.
b=->m{s,p=0r,1;m.times{|k|a=m-k;s+=p*b[k]/-~a;p=p*a/-~k};1-s}

Since this uses Rational in its answers, I'm fairly sure this works up to 60, but I had trouble running even b[24], so I implemented the lookup table again for 86 81 80 bytes.
t=->m{s,p,r=0r,1,m>0?t[m-1]:[];m.times{|k|a=m-k;s+=p*r[k]/-~a;p=p*a/-~k};r<<1-s}

Actually, 46 45 bytes
I've been meaning to do a Seriously/Actually answer for months and now I can. Golfing suggestions welcome. Try it online!
Edit: In February of 2017, there was an update to Actually that changed which function literals are which. Enjoy.
This uses the explicit definition of the Bernoulli numbers on Wikipedia.
;╖ur⌠;╝ur⌠;;0~ⁿ(╛█*╜(uⁿ*⌡MΣ╛uk⌡M┬i@;π;)♀\*@k▼

Ungolfing
;╖     Duplicate and save m in register 0.
ur     Range [0..m]
  ⌠      Start first for loop
  ;╝     Duplicate and save k in register 1.
  ur     Range [0..k]
    ⌠      Start second for loop (as string).
    ;;     Duplicate v twice.
    0~ⁿ    Push -1, and pow() to get (-1)**v.
    (╛█    Rotate a duplicate v to TOS, push k, and binom(k, v).
    *      Multiply (-1)**v by binom(k, v).
    ╜(uⁿ   Push m, rotate last duplicate v to TOS, increment, and pow() to get (v+1)**m.
    *      (-1)**v * binom(k, v) * (v+1)**m
    ⌡      End second for loop (string turned to function).
  MΣ     Map over range [0..v] and sum
  ╛u     Push k and increment (the denominator)
           (Note: second for loop does numerators only as denominator only depends on k)
  k      Push fraction in list form [numerator, denominator]
  ⌡      End first for loop
M      Map over range [0..k]
┬i@    Transpose all of the fractions, flatten and swap.
         Stack: [denominators] [numerators]
;π     Duplicate and take product of denominators.
;)     Duplicate product and move to bottom of stack.
         Stack: product [denominators] [numerators] product
♀\     For all items in denominators, integer divide product by item.
         Return a list of these scaled-up denominators.
*      Dot product of numerators and the scaled-up denominators as new numerator.
         (In effect, getting the fractions to the same denominator and summing them)
@k     Swap new numerator and product (new denominator) and turn into a list (fraction).
▼      Divide fraction by gcd(numerator, denominator) (Simplify fraction).

Uses commands Seriously didn't have in November 2015? Man, this uses an entirely new language that didn't exist in November 2015! I'm so proud...

J, 10 bytes
(%1-^@-)t:

Computes the nth Bernoulli number by finding the nth coefficient of the exponential generating function of x/(1 - e-x).
Usage
If the input is given integer or floats as an argument, it will output a float. If given an extended integer, marked with a suffix x, it will output either an extended integer or a rational, two extended integers separated by r.
   f =: (%1-^@-)t:
   f 1
0.5
   f 1x
1r2
   (,.f"0) i. 10x
0     1
1   1r2
2   1r6
3     0
4 _1r30
5     0
6  1r42
7     0
8 _1r30
9     0

Explanation
(%1-^@-)t: Input: n
(      )t: Takes a monad and creates a new monad that
           computes the coefficients of its egf
(      )   A monad that operates on x
      -      Negate x
    ^@       Computes its exponential, e^-x
  1-         Subtract it from 1
 %           Divide x by it, x/(1 - e^-x)

Axiom, 134 147 bytes
b(n:NNI):FRAC INT==(v:=[1/1];k:=1;repeat(k>n=>break;r:=1-reduce(+,[binomial(k,j)*v.(j+1)/(k-j+1)for j in 0..k-1]);v:=append(v,[r]);k:=k+1);v.(n+1))

ungolf and test
(23) -> b
   (23)
   b n ==
           1
     v := [-]
           1
     k := 1
     repeat
       if n < k
         then break
         else
                               binomial(k,j)v(j + 1)
           r := 1 - reduce(+,[[--------------------- for j in 0..(k - 1)]])
                                     k - j + 1
           v := append(v,[r])
           k := k + 1
     v(n + 1)
                                                   Type: FunctionCalled b
(50) -> [[i,b(i)]  for i in [0,1,2,3,4,5,6,7,8,9,10]]
   (50)
             1     1              1            1              1             5
   [[0,1],[1,-],[2,-],[3,0],[4,- --],[5,0],[6,--],[7,0],[8,- --],[9,0],[10,--]]
             2     6             30           42             30            66
                                         Type: List List Fraction Integer

(51) -> b 1000
   (51)
   -
   18243104738661887254572640256857788879338336867042906052197158157641126_
    2572624911158657472577321069709615489924627495522908087488299539455188_
    7918567582241551668492697244184914012242579830955617098629924652251740_
    9791915637226361428342780548971002281045465308441161372350696920220116_
    2441791760680262602019620260255790058416539271332852806000966628467639_
    0683434226380702951226108116666172815817157023611889303668166839919156_
    3797683877845690114843122753427426880591799883780255338278664578660218_
    5045895962670442011443630321460259486764674312436994856054301765557425_
    1371150213401051058408679874766352952749178734973676859834707623881634_
    6251471489942512878190574323531299070406930309477389251738705417680653_
    1183648189451892725726445949589759600705334767585389769924857630972963_
    9976364832442643512622073858780110731539833099817555775136008111170797_
    6250597322951308884900670113339167641953793994512377610306198429310933_
    1214632141683542607746641232089854815064629129596536997380608256428801_
    9784909897301658268809203555030692846151917069465607257641149187197651_
    0905515966840312411845543650593021402849221691341852819791233589301994_
    1012291773441794027493574651881059432274494354092231954894280742068472_
    7146192942133436054611475404867886313250114399681532753236429290625909_
    3411000391368336312138915621701535954814084208794241665492294270773347_
    6055878415765927582014214726584822236443691314366097570085473354584000_
    9985915190584047337934331297339403392719579093995842312746836871169674_
    9786460913411872527166990047126222109345933847358924230951718379883743_
    2563465487604170316077418754242710065269818190591271690695446633836120_
    3745255515267088218383996330164203403732365333352120338272021319718003_
    5994220458994876460018350270385634117807768745161622933834063145505621_
    9106004731529642292049578901
     /
    342999030
                                                   Type: Fraction Integer

(52) -> b 60
           1215233140483755572040304994079820246041491
   (52)  - -------------------------------------------
                             56786730
                                                   Type: Fraction Integer

APL(NARS), 83 chars, 166 bytes
r←B w;i
r←,1⋄i←0x⋄w+←1
→3×⍳w≤i+←1⋄r←r,1-+/{(1+i-⍵)÷⍨(⍵!i)×r[⍵+1]}¨0..i-1⋄→2
r←r[i]

Input as integer output as big rational
  B 0
1
  B 1
1r2
  B 2
1r6
  B 3
0
  B 4
¯1r30
  B 10
5r66
  B 100
¯94598037819122125295227433069493721872702841533066936133385696204311395415197247711r33330
  B 1000
¯1824310473866188725457264025685778887933833686704290605219715815764112625726249111586574725773210697096154899246
  27495522908087488299539455188791856758224155166849269724418491401224257983095561709862992465225174097919156
  37226361428342780548971002281045465308441161372350696920220116244179176068026260201962026025579005841653927
  13328528060009666284676390683434226380702951226108116666172815817157023611889303668166839919156379768387784
  56901148431227534274268805917998837802553382786645786602185045895962670442011443630321460259486764674312436
  99485605430176555742513711502134010510584086798747663529527491787349736768598347076238816346251471489942512
  87819057432353129907040693030947738925173870541768065311836481894518927257264459495897596007053347675853897
  69924857630972963997636483244264351262207385878011073153983309981755577513600811117079762505973229513088849
  00670113339167641953793994512377610306198429310933121463214168354260774664123208985481506462912959653699738
  06082564288019784909897301658268809203555030692846151917069465607257641149187197651090551596684031241184554
  36505930214028492216913418528197912335893019941012291773441794027493574651881059432274494354092231954894280
  74206847271461929421334360546114754048678863132501143996815327532364292906259093411000391368336312138915621
  70153595481408420879424166549229427077334760558784157659275820142147265848222364436913143660975700854733545
  84000998591519058404733793433129733940339271957909399584231274683687116967497864609134118725271669900471262
  22109345933847358924230951718379883743256346548760417031607741875424271006526981819059127169069544663383612
  03745255515267088218383996330164203403732365333352120338272021319718003599422045899487646001835027038563411
  78077687451616229338340631455056219106004731529642292049578901r342999030

Haskell, 95 bytes
import Data.Ratio
p=product
b m=sum[p[k-v+1..k]*(v+1)^m%(p[-v..0-1]*(k+1))|k<-[0..m],v<-[0..k]]

This implements the explicit definition of Bernoulli numbers outlined on the Wikipedia page.

Perl 6, 83 bytes
my &B={$^m??1-[+] (^$m).map: {combinations($m,$_)*B($_)/($m+1-$_)}!!1};say B slurp;

A faster, 114-byte solution:
my @b=1;for 1..+slurp() {@b.push: 1-[+] (^$^m).map: {([*] $m+1-$_..$m)*@b[$_]/($m+1-$_)/([*] 1..$_)}};say @b[*-1];

Your code for a code golf challenge needs to be as short as possible, even if it takes a few lifetimes of the universe to terminate for certain inputs.

Javascript, 168 bytes
function h(b,a){return a?h(a,b%a):b}for(var c=[],a=[],e=0,b,d,f,g;e<=k;)for(c[b=d=e]=1,a[e]=++e;d;)f=c[--d]*a[b]-(c[b]*=g=a[d]),r=h(f*=b,g=a[b]*=g),c[d]=f/r,a[--b]=g/r;

Set the 'k' variable to the Bernoulli number you want, and the result is c[0] over a[0]. (numerator & denominator)
Sample Usage
k = 2;
console.log(c[0] + "/" + a[0]);

Not as small as the others, but the only one I've written that comes close. See https://marquisdegeek.com/code_ada99 for my other (non-golf) attempts.

Axiom, 57 bytes
g(n)==factorial(n)*coefficient(taylor(t*%e^t/(%e^t-1)),n)

code for test and results
(18) -> [[i, g(i)]  for i in 0..29]
   (18)
              1      1                1              1                1
   [[0,1], [1,-], [2,-], [3,0], [4,- --], [5,0], [6,--], [7,0], [8,- --],
              2      6               30             42               30
                5                  691               7                 3617
    [9,0], [10,--], [11,0], [12,- ----], [13,0], [14,-], [15,0], [16,- ----],
               66                 2730               6                  510
                43867                 174611               854513
    [17,0], [18,-----], [19,0], [20,- ------], [21,0], [22,------], [23,0],
                 798                    330                  138
          236364091               8553103                 23749461029
    [24,- ---------], [25,0], [26,-------], [27,0], [28,- -----------], [29,0]]
             2730                    6                        870
                                       Type: List List Expression Integer

(19) -> g 60
           1215233140483755572040304994079820246041491
   (19)  - -------------------------------------------
                             56786730
                                                 Type: Expression Integer

one has to note that the function is not the one someone wrote above
but t*%e^t/(%e^t-1)) with %e Euler costant

Pyth, 22 bytes
L?b-1sm*.cbdcyd-btdUb1

Try it online!
Defines a function that is called as y<number>, e.g. yQ.
L                      # y=lambda b:
 ?b                  1 # ... if b else 1
   -1                  # 1 -
     s                 #     sum(
      m            Ub  #         map(lambda d: ... , range(b))
       *.cbd           #           combinations(b, d) *
            cyd        #             y(d) /      (float division)
               -btd    #                    b - (d - 1)

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


text:
Guess the hidden communities
Background
This challenge is about the stochastic block model.
Basically, you are given an undirected graph, where the nodes represent people, and the edges represent social connections between them.
The people are divided into two "communities", and two people are more likely to have a connection if they are in the same community.
Your task is to guess these hidden communities based on the connections.
Input
Your input is a 2D array of bits of size 100×100.
It represents the adjacency matrix of an undirected graph of size 100.
The graph is generated by the following random process:

Each node is placed in community 0 with probability 0.5, otherwise it's placed in community 1. The communities are not visible in the input.
For each pair of distinct nodes that are in the same community, put an edge between them with probability 0.3.
For each pair of distinct nodes that are in different communities, put an edge between them with probability 0.1.

All random choices are independent.
Output
Your output is a 1D array of bits of length 100.
It represents your program's guess about the underlying communities.
You are not required to correctly guess which nodes belong to community 0 and which to community 1 (which would be impossible), just the partition into two communities.
Scoring
In this repository, you'll find a text file containing 30 lines.
Each line represents a graph generated by the above process, and it consists of the array of the hidden communities, a semicolon, and a comma-separated list of rows of the adjacency matrix.
You should read this text file, call your submission program on the adjacency matrix of each line, and compare its output to the underlying community structure.
Your submission's score on a particular graph is
min(sum_i(V[i] == W[i]), sum_i(V[i] == 1-W[i]))

where V is the array of hidden communities, W is your program's guess based on the adjacency matrix, and i ranges over their indices.
Basically, this measures the number of mis-classified nodes, taking into account that you're not required to guess which community is which.
The lowest average score over all 30 test cases wins, with byte count used as a tie breaker.
Example
Consider the input matrix (6×6 for clarity)
011000
101100
110000
010011
000101
000110

It corresponds to the 6-node graph
a-b-d
|/ /|
c e-f

consisting of two triangles with one edge between them.
Your program might guess that each triangle is a separate community, and output 000111, where the first three nodes are in community 0 and the rest in community 1.
The output 111000 is equivalent, since the order of the communities doesn't matter.
If the true community structure was 001111, then your score would be 1, because one node is mis-classified.
Additional rules
Your program should take no more than 5 seconds on any single test case on a reasonably modern computer (I'm unable to test most answers, so I'll just trust you here).
You can write a full program or a function, and standard loopholes are disallowed.
You are allowed to optimize your answer toward the parameters used here (size 100, probabilities 0.3 and 0.1), but not the actual test cases.
I keep the right to generate a new set of graphs if I suspect foul play.
Test program
Here's a Python 3 function that computes the score of a Python 3 function submission.
import time

def test_from_file(submission, filename):
    f = open(filename)
    num_graphs = 0
    total_score = 0
    for line in f:
        num_graphs += 1
        a, _, c = line.partition(';')
        V = [int(node) for node in a]
        E = [[int(node) for node in row] for row in c[:-1].split(',')]
        start = time.perf_counter()
        W = submission(E)
        end = time.perf_counter()
        x = sum(V[i]==W[i] for i in range(len(V)))
        score = min(x, len(V) - x)
        print("Scored {0:>2d} in {1:f} seconds".format(score, end - start))
        total_score += score
    f.close()
    return total_score / num_graphs

R, 16/30 = 0.53333...

Edit:
Here is a version that has been golfed down to 259 bytes. I am sure a much better bytecount is possible. This version takes advantage of the fact that T is a global variable in R which can be treated as having a value of 1. If T has been reassigned in your R environment, results may differ.
The code produces some warnings, but they can be ignored.
f=function(A){L=function(p){q=1-p;z=.1*p%o%q;x=z+t(z)+.3*(p%o%p+q%o%q);c(sum(log(x^A*(1-x)^(1-A))),(A/x-(1-A)/(1-x))%*%(2*p-1))};x=tanh(rnorm(100))/2+.5;for(i in x){t=.5+tanh(atanh(2*x-1)+T*2*(L(x)[2:101]>0)-T)/2;if(L(t)>L(x)){x=t;T=2*T}else{T=T/2}};round(x)}

If the adjacency matrix of your graph is saved as a 100x100 matrix A in R, then the vector of communities is f(A). It will probably make less than one mistake!

Explanation
It treats the problem as fitting a statistical model (it helps that it knows the values 0.3 and 0.1) and fits it using a sort of gradient ascent algorithm. The community membership vector is treated as containing values in [0,1], rather than just zeroes and ones. This enables the model to be fitted using a sort of gradient ascent. (Ascent and descent is actually by a fixed amount delta in each coordinate, to avoid various problems. If a step of gradient ascent fails, delta is adjusted down.) It's quite fast and does pretty well.
I used today's date as the random seed.
LL <- function(p, A){
# log-lieklihood of the model
  q <- 1-p
  x <- 0.1*(p %o% q) + 0.1*(q %o% p) + 0.3*(p %o% p) + 0.3*(q %o% q)
  sum(A*log(x) + (1-A)*log(1-x))
}

expit <- function(x) exp(x)/(1+exp(x))

gradLL <- function(p, A){
# (scalar multiple of) the gradient of the log-likelihood
  q <- 1-p
  x <- 0.1*(p %o% q) + 0.1*(q %o% p) + 0.3*(p %o% p) + 0.3*(q %o% q)
  drop((A/x - (1-A)/(1-x)) %*% (0.4*p - 0.2))
}

fit_model <- function(A, N=100){
# gradient ascent
  x <- rnorm(100)
  score <- LL(expit(x), A)
  delta <- 1
  for (i in 1:N){
    grad <- gradLL(expit(x), A)
    test <- x + delta*((grad > 0) - (grad < 0))

    test_score <- LL(expit(test), A)
    if (test_score > score){
      score <- test_score
      x <- test
      #print(score)
      delta <- (1/0.7)*delta
    } else {
      delta <- 0.7*delta
    }
  }
round(expit(x))
}

generate_A <- function(){
# generate a matrix for testing
  comms <- sample(0:1, 100, replace=T)
  A <- matrix(rep(0, 100^2), nc=100)
  for (i in 1:(nrow(A)-1)){
    for (j in (i+1):ncol(A)){
      prob <- ifelse(comms[i]==comms[j], 0.3, 0.1)
      if (runif(1) < prob) A[i,j] <- 1
    }
  }
  A <- (A + t(A))/2
  list(A=A, comms=comms)
}

test <- function(){
# test the fitting function once
  dat <- generate_A()
  fit <- fit_model(dat$A)
  sum(fit==dat$comms)
}

### performance ###

strings <- readLines("community-graphs.txt", 30)

total <- 0

set.seed(011215) # note that output is random!

for (i in 1:30){
  str_data <- strsplit(strings[i], ";")[[1]]
  comms <- as.numeric(strsplit(str_data[1], "")[[1]])
  A_str <- strsplit(str_data[2], ",")[[1]]
  A <- matrix(rep(0, 100^2), nc=100)
  for (r in 1:100){
    A[r,] <- as.numeric(strsplit(A_str[r], "")[[1]])
  }
  out <- fit_model(A)
  total <- total + min(sum(out==comms), sum(out!=comms))
}

total # 16

JavaScript (ES6), Score 16.83
Just a quick and simple solution that doesn't get a very good score. I know this isn't code-golf but here's my (lightly golfed) solution in 134 bytes:
m=>(l=m.split`,`,n=l.map((_,i)=>[!i,0]),l.map((g,a)=>[...g].map((h,b)=>n[b][n[a][0]<n[a][1]==h|0]++)),n.map(g=>o+=g[0]<g[1]|0,o=""),o)

Explanation
m=>(                                    // m = comma delimited input matrix string

  // Initialisation
  l=m.split`,`,                         // l = input matrix as array
  n=l.map((_,i)=>                       // n = array of each node with a count of their
                                        //     graphs to nodes of each community
    [!i,0]                              // initialise the first node to community 0
  ),

  // Find graphs
  l.map(                                // compute each node
    (g,a)=>                             // a = index of node A, g = graphs of node A
      [...g].map(                       // compute each of the node's graphs
        (h,b)=>                         // b = index of node B
                                        // h = 1 if node A has a graph to node B
          n[b][n[a][0]<n[a][1]==h|0]++  // increment node B's number of graphs to node A's
      )                                 //     community if they have a graph, else
  ),                                    //     assume and increment the other community

  // Output the results
  n.map(                                // iterate over the graphs
    g=>                                 // g = current node's graphs
      o+=g[0]<g[1]|0,                   // add the number of the predominant community to
                                        //     the output
      o=""                              // o = output string
  ),
  o                                     // return o
)

Possible Improvements

I could implement two passes in the "Find graphs" step, with the first pass getting the probable community of each node and reinitialising each node to that community

Test

var solution = m=>(l=m.split`,`,n=l.map((_,i)=>[!i,0]),l.map((g,a)=>[...g].map((h,b)=>n[b][n[a][0]<n[a][1]==h|0]++)),n.map(g=>o+=g[0]<g[1]|0,o=""),o);

function test() {
  var lines = input.value.split("\n"),
      output = [],
      totalScore = 0;
  for(var i = 0; i < lines.length; i++) {
    var parts = lines[i].split(";"),
        correct = parts[0],
        matrix = parts[1],
        result = solution(matrix);
    var sum = 0,
        sumInverse= 0;
    for(var c = 0; c < correct.length; c++) {
      if(correct[c] == result[c]) sum++;
    }
    for(var c = 0; c < correct.length; c++) {
      if(correct[c] == 1 - result[c]) sumInverse++;
    }
    var score = Math.min(sum, sumInverse);
    totalScore += score;
    output.push("Test " + (i + 1) + ": Score = " + score + ", Result = " + result);
  }
  output.splice(0, 0, "Average Score: " + (totalScore / lines.length));
  results.innerHTML = output.join("\n");
}
Enter the correct result followed by ";" then the input matrix delimited by ",".

Put each test on a separate line.

<textarea id="input" rows="5" cols="60">111000;011000,101100,110000,010011,000101,000110</textarea>

<button onclick="test()">Go</button>

Mathematica, score 31/15 = 2.066666...
Range[100]/.Thread[Alternatives@@@FindGraphCommunities[AdjacencyGraph[#,Method->{"Modularity","StopTest"->(Length[#2]==2&)}]]->{0,1}]&

Please don't make me recalculate the score.

We meet again, Mathematica builtin. :P


text:
Is it a factor of a polynomial?
A polynomial is divisible by a factor (x-n) if f(n)=0 for a function f.  Your job: to determine if a polynomial function f(x) is divisible by (x-n).
The input
The input is in the form of (x-n), (Polynomial).  Remember, if n is negative, (x-n) will be in the input form of (x+n).  For the polynomial, all exponents will be put in as ^.  Coefficients will be written next to the variable x.  An example polynomial could be 2x^2 + x^1.  There will be no spaces between anything. The term x will be inputed as x^1.  So what would "normally" look like (x - 1) will be (x^1-1).  The coefficients and powers will always be integers.  The coefficient one will be implicit if it is just x. I.e., x can be interpreted as 1x
The output
A boolean value.  Truthy, or Falsey.
Thanks to @AlexA. For helping me clarify this!
Examples
Input:(x^1-1),(x^1-1)
Output: True

Input: (x^1+2),(2x^2+4x^1+2)
Output: False

Input: (x^1+7),(x^2-49)
Output: True

Rules

This is code-golf, so shortest code in bytes wins

Unfortunately, I don't know how to implement the snippet leaderboard.  If anyone knows how, feel free to edit the post.

Will the input be a string with that exact form, i.e. parens around the divisor candidate, a comma with zero or one space, and parens around the polynomial?

Possible duplicate

Definitely not a duplicate of that.

@intrepidcoder This isn't a duplicate because the question isn't to factor a polynomial.  It is to see if a polynomial can be divided by a linear factor.

Will polynomial coefficients always be integers?

That is correct.

In the "Input" section you say that terms with a coefficient of 1 will be written as 1x^..., but in the examples you leave out the 1. Which is it?

The input will be x^....  The coefficient of 1 is implicit if it is just x

Then you should change where it says An example polynomial could be 2x^2 + 1x^1.

Nice catch! Fixed it!

There are still some inconsistencies in the described notation and the actually used notation: the input is really in the form (x^1-n),(Polynomial); it's not true that "all exponents" are given explicitly, because the constant term is e.g. -1 rather than -1x^0 (or -x^0); the sentence about the implicit 1 should use x^a or a specific value of a rather than just x. It would also be worth stating explicitly that terms with 0 coefficient are (or only may be?) omitted, and saying whether the terms are always ordered by exponent, descending.

I suggest using lists/array of integers as inputs insted of strings one has to parse.

@flawr keep in mind that if the input was an array of integers, this challenge would be 100x easier.

@intboolstring Yes it would be, and I suggested this because I thought you wanted to do a challenge that was about checking divisibility and not mainly about parsing, as many other code-golf challenges...

@flawr.  Ok, makes sense.

Pyth - 39 bytes
This is a monstrous combination of regexp and eval. I like the approach, but will try to improve the implementation.
It uses the Polynomial Remainder Theorem.
K_sPe:z"-|\+"3!v.ssXPtw,\^\x,"**""*K"\*

Doesn't work online because of eval use.

Casio Basic, 19 bytes
judge(mod(b,a)=0

As it turns out, the fx-CP400 can do mod on algebraic expressions!
Polynomial and factor should be entered as expressions. 16 bytes for the code, 3 bytes to enter a,b into the parameter value box.

VBScript, 118 116 bytes
a=inputbox(""):for i=0 to 9:a=replace(a,i&"x",i&"*x"):next:b=split(a,","):x=-eval(b(0)):msgbox not cbool(eval(b(1)))

Since we know that the first part of the input is a linear polynomial, we only need to check if its root matches that of the second polynomial; and we need to prepare the term for eval by inserting * as needed.

MATLAB, 103 99 97 95 93 bytes
I'm trying some different things, and got this to work to save a couple of bytes:
eval([regexprep(input(''),{'.+?1(.+)\),','(\d)x'},{'x=str2num(''$1'');disp(~','$1\*x'}) 41]);

If I can reduce that down even further, I will post an explanation.

Old code an explanation
t=sscanf(input(''),'(x^1%d),%s')';x=-t(1);disp(~eval(regexprep([t(2:end) ''],'(\d)x','$1\*x')))

This also works with Octave. You can try it online. I've save the program as a script named isFactor.m, so you can just enter isFactor at the prompt. [Note: in Octave spits out a warning while running - MATLAB doesn't generate this].
The input must be in the format '(x^1+7),(x^2-49)' as per the question. The quotes marks are added so MATLAB/Octave knows it is a string.
The output is either a 0 or a 1 depending on whether or not it is true or false.

So, the code works as follows. First we request an input, and then parse it. The parsing string extracts the signed number after the first (x^1 in the string - this is our value of n. Then it continues to extract the string (%s) after the ), in the input - this is our expression.
t=sscanf(input(''),'(x^1%d),%s')';

Next, we extract the value of n and set x equal to it - we are going to evaluate whether the expression equals zero when n==x, so this is why we store the value to x. Also we negate the extracted number, because of the minus sign when parsing.
x=-t(1);

We will then display the output which is a boolean
disp(

The output is basically the logical negation of our evaluated equation. If f(x) is zero, this will return 1, otherwise it will result in zero.
     ~eval(

We are evaluating the input expression, but to do this, we need to reformat it slightly so that MATLAB can understand. When we read the string, it is actually an array of double type, so we need to convert that to a character array. Before the conversion we also get rid of the first element as that is what we used for n. We then need to replace any occurrence of x which is preceded by a number (e.g 4x) by the same thing but with a multiplication (*) sign in between so MATLAB can calculate it.
           regexprep(char([t(2:end) ''],'(\d)x','$1\*x')
     )
)

Axiom 77 180 Bytes
f(a:UP(x,INT),b:UP(x,INT)):Boolean==(ground?(a)or ground?(b)=>false;p:=b;r:=a;if degree(a::POLY INT,x)>degree(b::POLY INT,x)then(p:=a;r:=b);(p rem r)$UP(x,FRAC INT)~=0=>false;true)

the previous solution
v(a,b)==(ground?(a) or ground?(b) or (b rem a)$UP(x,FRAC INT)~=0=>false;true)

was wrong because it assume degree(b)>=degree(a) one bug i wrote...
test and results
(3) -> f(x^1-1,x^1-1)
   (3)  true
                                                            Type: Boolean
(4) -> f(x^1+1,2*x^2+4*x^1+2)
   (4)  true
                                                            Type: Boolean
(5) -> f(x^1+2,2*x^2+4*x^1+2)
   (5)  false
                                                            Type: Boolean
(6) -> f(x^1+7,x^2-49)
   (6)  true
                                                            Type: Boolean
(7) -> f(1, 1)
   (7)  false
                                                            Type: Boolean
(8) -> f(1, x^2+1)
   (8)  false
                                                            Type: Boolean
(9) -> f(x^8-1, x^2-1)
   (9)  true

Mathematica, 30 bytes
Try it online!
PolynomialRemainder[#2,#,x]<1&


text:
Guess my password! (Code Golf)
So, like the title says, today you are going to guess my password (sort of). You need to make a big long list of every possible password I might be using on my 4 digit online vault.
Specifications:

It could contain any 26 letters of the alphabet in UPPER CASE.
It could contain any numbers 0-9.
So, basically, each digit has 36 different possibilities.
There are four digits.
Repeat characters are allowed! (Hint: Sample space is 1,679,616 or
36^4)

Program Requirements/Rules:

Program must output all possible answers to a text file or to the console, line by line. (One answer per line)
You may not make use of any premade tools for this (Idk if any exist, but just in case)
It must generate all 36^4 possible answers

Remember, this is a Code-golf. So shortest answer in bytes wins. Please use the standard title format of ##Language Name, Byte Size. Even if there is a winner, feel free to post your answer. (It's always fun to see the different solutions in different languages).
Edit: Not a duplicate of sugested, mines has different context, and different paremeters. (I require numbers, duplicate characters, and only UPPERCASE)

I'm not sure if this should be closed as a duplicate of that question; it requires digits as well as letters, so answers can't just be copied over. Then again, it is pretty similar...

@ETHproductions yeah, IDK. I think its sort of different. I mean, the practice/methods to solve are similar/same but then the application is different.

Both questions are "Count in base B" for some B. There may be superficial differences, but they're fundamentally the same question.

If it didn't have the word "password" in it, there is now way it could be thought of as a copy of that question. Answer to this question in Perl 6 .say for [X~] (0..35)».base(36)xx 4, answer to the other question is say 'a'.../password/ ( or print "$_ " for 'a'.../password/ if you want it to start printing immediately. ) There is no easy way to answer either question starting out with the answer for the other one. I say remove all mention of the word password, and try to re-open.

@BradGilbertb2gills thank you. Thats what I'm trying to say here...

@BradGilbertb2gills you think I should re open?

Pyth - 8 bytes
^+jkUTG4

Try smaller password length online here.

I'm not entirely sure I understand, the output I got when I ran your program is in pairs of two? I said the password is 4 digit.

@AshwinGupta yes, because the online interpreter would reject a program with 36^4 permutations.

Ah wait nevermind, I see, you changed the 4 to 2 for running purposes. Sorry my bad.

I'm afraid to upvote when your score is 7 777 :)

You are currently the winner unless someone can beat you. Where are the Japt and Seriously guys? They both are short-syntax golfing languages with a chance.

@AshwinGupta Sadly, the best Japt answer I can come up with is 22 bytes: 36p4 o £'0³+Xs36 u)t-4 (to try this at home, change 36p4 to 36p2)

@ETHproductions ok, still a great answer! I love seeing these tiny golf language responses. It makes me really want to learn one. I think I might try and learn some Pyth. (I know a bit of python already so it might help). Im glad I only made you guys do 4 digits though. I was initially considering 8 digits with lower and upper case. Luckily, I thought against it because of the sheer amount of answers (62^8 = 218,340,105,584,896 or 218 trillion answers. YIKES!)

@Ashwin Yeah, the code actually finishes in ~45 seconds on my browser, as opposed to the 2.5 years it would take to do 8 digits ;)

@ETHproductions LOL. What? You don't want to leave your pc on for 2.5 years to run a golf program? (jk)

Haskell, 38 Bytes
mapM(\a->['A'..'Z']++['0'..'9'])"abcd"

yeah, it's not gonna beat Pyth. Oh well...

Don't worry, not much beats Pyth. from what little I've seen here on CodeGolf SE in the last week or so. Still a great answer.

:p too true... I'll take beating another non-golf language for now...

Well, you beat Ruby! @quartata although it was close. 15 Bytes...

Ruby, 64 53 bytes
[*'A'..'Z',*'0'..'9'].repeated_permutation(4){|a|p a}

Probably golfable further, I'm not that good at Ruby.


text:
Implement the Method of Finite Differences
The Method of Finite Differences is a technique used to find the next term in a sequence of numbers, given that the sequence is governed by consecutive values of a certain polynomial. Given a list of n terms, the Method will be able to determine the degree of the polynomial, which is a maximum of n+1.
The Method works by arranging the numbers in order, like so.
4 8 15 16 23 42

Then the difference between any number and its neighbor to the right (given that it has a neighbor to the right) is written in the row below. Note that entries must be written one space (ASCII 20) to the right of the longest entry in the previous column.
4 8 15 16 23 42
4 7 1  7  19

The process is repeated until there is one entry left in the bottom row.
4   8  15 16 23 42
4   7  1  7  19
3   -6 6  12
-9  12 6
21  -6
-27

Then the bottom row's entry is duplicated.
4   8   15 16 23 42
4   7   1  7  19
3   -6  6  12
-9  12  6
21  -6
-27 -27

Then the above process for finding differences is reversed until there is a new entry in the first row, which is separated from the other entries by  |
4   8   15  16  23  42 | 46
4   7   1   7   19  4
3   -6  6   12  -15
-9  12  6   -27
21  -6  -33
-27 -27

The Method moves down rows until all of the entries in a row are equal. Take the starting numbers 2 5 10 17 26. This is the finished table.
2 5 10 17 26 | 37
3 5 7  9  11
2 2 2  2

The degree of the polynomial is the number of rows minus 1. As you might tell from the previous example, the polynomial is x^2+1.

Your program will take as its input a list of integers separated by whitespace. You may assume that the list contains only integers and has at least 2 elements (they need not be distinct; the resulting chart would have only 1 row). The program will output the entire table generated by the Method of Finite Differences, as defined above.

This is code golf. Standard rules apply. Shortest code in bytes wins.

So is the output the entire table, the top row, or just the rightmost number? How flexible is the formatting?

@Zgarb The entire table. I'd like the formatting to be as I described it, with the numbers in even columns separated by one space.

Related.

Also related

Are trailing spaces allowed?

@MartinBüttner Sure.

CJam, 74 73 71 68 bytes
q~]{_2ew::-Wf*_)-}g0]W%{_W=@W=++_}*;)"|"\++]W%z{:s_:,$W=)f{Se]}}%zN*

Try it online
The latest version uses the middle section of @MartinBüttner's proposed solution, which is 3 bytes shorter than what I had. The first part for calculating the difference rows was identical. The part for the final formatting was the same length, so I kept my original version. Martin's savings came from making the sum part of the list while the new values are calculated, and using a * instead of a % loop. That also simplified the insertion of the "|".
Explanation:
l       Get one line of input.
~]      Interpret input and wrap in list.
{       While loop to construct additional lines.
  _2ew    Copy and build pairs of subsequent numbers.
  ::-     Reduce all number pairs by subtraction.
  Wf*     Multiply all of them by -1, since the reduction calculated 1st minus 2nd,
          but we need 2nd minus 1st.
  _)-     Check if all numbers in list are equal by popping off last and removing
          value from the list. This is truthy if all numbers are not equal.
}g      End of while loop to construct additional lines.
0       Add a 0 at the end, which will be the starting point of the sum for
        calculating the new values.
]       Wrap all lines in a list.
W%      Revert order of list, since adding the additional numbers in each row
        goes in bottom-up order.
{       Start of loop over rows.
  _W=     Get last value of row.
  @       Rotate previous row to top.
  W=      Extract last value as new increment. For the 0 that was added as an
          additional row, this is a comparison with -1, which gives the desired
          value of 0.
  +       Add increment to last value.
  +       Append value to list.
  _       Copy row for use in next loop iteration.
}*      End of loop over rows.
;       Remove extra copy of last row.
)       Pop off last value of last row.
"|"     Push the vertical bar.
\       Swap vertical bar and last value.
++      Concatenate the 3 parts.
W%      Revert order of list, to get the rows back in original order.
z       Transpose list for handling column alignment.
{       Loop over columns.
  :s      Convert values to strings.
  _:,     Copy and calculate the length of all values.
  $W=     Get maximum length by sorting and getting last value.
  )       Increment to add space between columns.
  f{      Apply to all values.
    Se]     Right pad to column width with spaces.
  }       End apply to all values.
}%      End loop over columns.
z       Transpose back to row order.
N*      Join rows with newlines.

I was going to post my own but it turned out almost exactly the same as yours. After using some more of your ideas I got it down to 70. Feel free to use it: q~]{_2ew::-Wf*_)-}g0]W%{_W=@W=++_}*;)"|"\++]W%z{:s_:,:e>f{Se]}}%zSf*N*

Also, both your and my code occasionally print trailing spaces. If that's allowed (I've already asked), then you can save two more bytes by ditching the Sf* and incrementing before Se].

@MartinBüttner Thanks, I'll check it out later today. I was trying to avoid trailing spaces, but I see now that it wasn't entirely successful in the case where the last number in a row is not the longest. Allowing them will simplify the formatting.

APL (Dyalog Unicode), 93 bytes
{↑⍕¨¨{(⊂1⌽'|',⍨¯1⌽⊃⍵),(¯1↓1↓⍵),⊂(⊢,⊢/)⊃⌽⍵}⊃{(⊂⍺,(⊃⌽⍺)+⊃⌽⊃⍵),⍵}/⌽⊂@1⌽{⍵,⊂-2-/⊃⌽⍵}⍣{1≥≢∪⊃⌽⍺}⊂⍵}

Try it online!
The output format is a bit annoying, and takes up most of the space here.