Wolstenholme's theorem states that:

Wolstenholme's theorem

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.


Verification of Wolstenholme's theorem

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




where and the parentheses thingy is Binomial coefficient.

You can assume that 2b <= a


a b p  output
6 2 5  240360
3 1 13 3697053
7 3 13 37403621741662802118325
  • 2
    \$\begingroup\$ I feel like outputs should have a .0 on the end, to really show that there's no leftover.from the division. \$\endgroup\$ Jul 30, 2016 at 15:00
  • 3
    \$\begingroup\$ @El'endiaStarman Come on. \$\endgroup\$
    – Leaky Nun
    Jul 30, 2016 at 15:03
  • 1
    \$\begingroup\$ Would [240360] (singleton array) be an acceptable output format? \$\endgroup\$
    – Dennis
    Jul 30, 2016 at 15:16
  • 1
    \$\begingroup\$ I don't think there is one, which is why I'm asking. \$\endgroup\$
    – Dennis
    Jul 30, 2016 at 15:24
  • 2
    \$\begingroup\$ @Dennis Then make one. \$\endgroup\$
    – Leaky Nun
    Jul 30, 2016 at 15:25

14 Answers 14


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
  • 5
    \$\begingroup\$ Shouldn't the Haskell 98 version still be a valid submission? \$\endgroup\$ Jul 30, 2016 at 17:50

Jelly, 12 11 10 bytes


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

Try it online! or verify all test cases.

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, 114 109 85 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
  • \$\begingroup\$ Here is a golfed factorial lambda \$\endgroup\$ Jul 30, 2016 at 16:32

05AB1E, 11 bytes

Takes input as:

[a, b]



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

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

R, 50 48 bytes


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

  • 1
    \$\begingroup\$ How many of those spaces are necessary? \$\endgroup\$
    – Neil
    Jul 30, 2016 at 22:51
  • \$\begingroup\$ 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). \$\endgroup\$
    – rturnbull
    Jul 28, 2018 at 1:46

MATL, 13 bytes


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.


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



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

For example:

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

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

Brachylog, 52 bytes

tT;T P&t^₃D&h↰₁S&h;Pz×₎ᵐ↰₁;S-;D/

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, 85 82 75 59 bytes

import math,future

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

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


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


APL (Dyalog), 18 bytes


Try it online!


Pari/GP, 43 bytes


Try it online!


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.