Verify Wolstenholme's theorem

Question

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.

Task

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

Answer 1

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

Answer 2

Jelly, 12 11 10 bytes

ż×c/I÷S÷²}

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³.

Answer 3

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

Answer 4

05AB1E, 11 bytes

Takes input as:

[a, b]
p

Code:

*`c¹`c-²3m÷

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

Answer 5

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.

Answer 6

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

Answer 7

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

Answer 8

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

Answer 9

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)

Answer 10

Nim, 85 82 75 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.

Answer 11

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.

Answer 12

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).

Answer 13

APL (Dyalog), 18 bytes

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

Try it online!

Answer 14

Pari/GP, 43 bytes

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

Try it online!