# Definition

Wolstenholme's theorem states that:

where a and b are positive integers and p is prime, and the big parentheses thingy is Binomial coefficient.

To verify that, you will be given three inputs: a, b, p, where a and b are positive integers and p is prime.

Compute:

where a and b are positive integers and p is prime, and the parentheses thingy is Binomial coefficient.

# Specs

Since:

where and the parentheses thingy is Binomial coefficient.

You can assume that 2b <= a

# Testcases

a b p  output
6 2 5  240360
3 1 13 3697053
7 3 13 37403621741662802118325

• I feel like outputs should have a .0 on the end, to really show that there's no leftover.from the division. Jul 30 '16 at 15:00
• @El'endiaStarman Come on. Jul 30 '16 at 15:03
• Would [240360] (singleton array) be an acceptable output format? Jul 30 '16 at 15:16
• I don't think there is one, which is why I'm asking. Jul 30 '16 at 15:24
• @Dennis Then make one. Jul 30 '16 at 15:25

# Haskell, 73 71 bytes

Due to the recursion, this implementation is very slow. Unfortunately my definition of the binomial coefficient has the same length as import Math.Combinatorics.Exact.Binomial.

n#k|k<1||k>=n=1|1>0=(n-1)#(k-1)+(n-1)#k --binomial coefficient
f a b p=div((a*p)#(b*p)-a#b)p^3       --given formula


An interesting oddity is that Haskell 98 did allow for arithmetic patterns which would have shortened the same code to 64 bytes:

g a b p=div((a*p)#(b*p)-a#b)p^3
n+1#k|k<1||k>n=1|1>0=n#(k-1)+n#k

• Shouldn't the Haskell 98 version still be a valid submission? Jul 30 '16 at 17:50

# Jelly, 1211 10 bytes

ż×c/I÷S÷²}


Expects a, b and p as command-line arguments.

### How it works

ż×c/I÷S÷²}  Main link. Left argument: a, b. Right argument: p

×          Multiply; yield [pa, pb].
ż           Zipwith; yield [[a, pa], [b, pb]].
c/        Reduce columns by combinations, yielding [aCb, (pa)C(pb)].
I       Increments; yield [(pa)C(pb) - aCb].
÷      Divide; yield [((pa)C(pb) - aCb) ÷ p].
S     Sum; yield ((pa)C(pb) - aCb) ÷ p.
²}  Square right; yield p².
÷    Divide; yield  ((pa)C(pb) - aCb) ÷ p³.


# Python 2, 11410985 71 bytes

A simple implementation. Golfing suggestions welcome.

Edit: -29 bytes thanks to Leaky Nun and -14 bytes thanks to Dennis.

lambda a,b,p,f=lambda n,m:m<1or f(n-1,m-1)*n/m:(f(a*p,b*p)-f(a,b))/p**3


A simpler, same-length alternative, with thanks to Dennis, is

f=lambda n,m:m<1or f(n-1,m-1)*n/m
lambda a,b,p:(f(a*p,b*p)-f(a,b))/p**3

• Here is a golfed factorial lambda Jul 30 '16 at 16:32

# 05AB1E, 11 bytes

Takes input as:

[a, b]
p


Code:

*c¹c-²3m÷


Uses the CP-1252 encoding. Try it online!.

• Did you out-golf Dennis? Jul 30 '16 at 15:52
• If I were in Dennis' shoes I think I'd get a little tired of all these "outgolf Dennis" comments... Jul 30 '16 at 15:53
• @LuisMendo I may or may not be nuking them on a regular basis. Jul 30 '16 at 16:04
• and hes at 10. it was fun while it lasted boys Jul 31 '16 at 19:17

# R, 50 48 bytes

function(a,b,p)(choose(a*p,b*p)-choose(a,b))/p^3


As straightforward as can be... Thanks to @Neil for saving 2 bytes.

• How many of those spaces are necessary?
– Neil
Jul 30 '16 at 22:51
• 42 bytes by renaming choose and by using pryr::f to define the function: B=choose;pryr::f((B(a*p,b*p)-B(a,b))/p^3). Jul 28 '18 at 1:46

# MATL, 13 bytes

y*hZ}Xnd2G3^/


Try it online!

The last test case doesn't produce an exact integer due to numerical precision. MATL's default data type (double) can only handle exact integers up to 2^53.

### Explanation

y   % Implicitly input [a; b] (col vector) and p (number). Push another copy of [a; b]
%   Stack: [a; b], p, [a; b]
*   % Multiply the top two elements from the stack
%   Stack: [a; b], [a*p; b*p]
h   % Concatenate horizontally
%   Stack: [a, a*p; b, b*p]
Z}  % Split along first dimension
%   Stack: [a, a*p], [b, b*p]
Xn  % Vectorize nchoosek
%   Stack: [nchoosek(a,b), nchoosek(a*p,b*p)]
d   % Consecutive differences of array
%   Stack: nchoosek(a,b)-nchoosek(a*p,b*p)
2G  % Push second input again
%   Stack: nchoosek(a,b)-nchoosek(a*p,b*p), p
3^  % Raise to third power
%   Stack: nchoosek(a,b)-nchoosek(a*p,b*p), p^3
/   % Divide top two elements from the stack
%   Stack: (nchoosek(a,b)-nchoosek(a*p,b*p))/p^3
% Implicitly display


# J, 17 bytes

(!/@:*-!/@[)%]^3:


## Usage

(b,a) ( (!/@:*-!/@[)%]^3: ) p


For example:

   2 6 ( (!/@:*-!/@[)%]^3: ) 5
240360


This is just a direct implementation of the formula so far.

Note: for the 3rd testcase input numbers must be defined as extended (to handle big arithmetic):

   3x 7x ( (!/@:*-!/@[)%]^3: ) 13x
37403621741662802118325


# Brachylog, 52 bytes

tT;T P&t^₃D&h↰₁S&h;Pz×₎ᵐ↰₁;S-;D/
hḟF&⟨{-ḟ}×{tḟ}⟩;F↻/


Try it online!

Accepts input [[a, b], p].

% Predicate 1 - Given [n, r], return binomial(n, r)
hḟF              % Compute n!, set as F
&⟨               % Fork:
{-ḟ}           % (n - r)!
×              % times
{tḟ}           % r!
⟩
;F↻              % Prepend n! to that
/                % Divide n! by the product and return

% Predicate 0 (Main)
tT;T P           % Set P to the array [p, p]
&t^₃D            % Set D as p^3
&h↰₁S            % Call predicate 1 on [a, b],
%  set S as the result binomial(a, b)
&h;Pz×₎ᵐ         % Form array [ap, bp]
↰₁               % Call predicate 1 on that to get binomial(ap, bp)
;S-              % Get binomial(ap, bp) - binomial(a, b)
;D/              % Divide that by the denominator term p^3
% Implicit output


# Python 3 with SciPy, 72 bytes

from scipy.special import*
lambda a,b,p:(binom(a*p,b*p)-binom(a,b))/p**3


An anonymous function that takes input via argument and returns the result.

There's not a lot going on here; this is a direct implementation of the desired computation.

Try it on Ideone (the result is returned in exponential notation for the last test case)

# Nim, 858275 59 bytes

import math,future
(a,b,p)=>(binom(a*p,b*p)-binom(a,b))/p^3


This is an anonymous procedure; to use it, it must be passed as an argument to another procedure, which prints it. A full program that can be used for testing is given below

import math,future
proc test(x: (int, int, int) -> float) =
echo x(3, 1, 13) # substitute in your input or read from STDIN
test((a,b,p)=>(binom(a*p,b*p)-binom(a,b))/p^3)


Nim's math module's binom proc computes the binomial coefficient of its two arguments.

# Python 2, 67 bytes

def f(a,b,p):B=2<<a*p;R=(B+1)**a;print(R**p/B**(p*b)-R/B**b)%B/p**3


Try it online!

Expresses the binomial coefficients arithmetically using this method.

## JavaScript (ES6), 70 bytes

(a,b,p,c=(a,b)=>a==b|!b||c(--a,b)+c(a,--b))=>(c(a*p,b*p)-c(a,b))/p/p/p


Save 1 byte by using ES7 (/p**3 instead of /p/p/p).

# APL (Dyalog), 18 bytes

{(-/!/¨⍺1×⊂⍵)÷⍺*3}


Try it online!

# Pari/GP, 43 bytes

(a,b,p)->((c=binomial)(a*p,b*p)-c(a,b))/p^3


Try it online!