# Wolstenholme numbers

The generalised harmonic number of order $$\m\$$ of $$\n\$$ is

$$H_{n,m} = \sum^n_{k=1} \frac 1 {k^m}$$

In this challenge, we'll be considering the generalised harmonic numbers of order $$\2\$$:

$$H_{n,2} = \sum^n_{k=1} \frac 1 {k^2}$$

This sequence begins $$\1, \frac 5 4, \frac {49} {36}, \frac {205} {144}, \dots\$$ and converges to $$\\frac {\pi^2} 6\$$ as $$\n \to \infty\$$.

However, we will only be considering the numerators of this sequence, which forms another sequence known as the Wolstenholme numbers (A007406):

1, 5, 49, 205, 5269, 5369, 266681, 1077749, 9778141, ...


You are to take a positive integer $$\x\$$ as input. $$\x\$$ is guaranteed to be a Wolstenholme number. You should output the next term in the sequence of Wolstenholme numbers, i.e. given the numerator of $$\H_{n,2}\$$, output the numerator of $$\H_{n+1,2}\$$

You may input and output in any convenient format. Your program should not fail due to floating point errors, meaning that attempts that calculate via the generalised harmonic numbers will be invalid if not using rational numbers. Your program should work theoretically for any valid input, but it is acceptable to fail for inputs that may be outside your language's integer maximum.

You may assume that no element is repeated in the sequence. If it turns out that there are any duplicate values in the Wolstenholme numbers sequence, any and all behaviour is acceptable when passed such values as input.

This is , so the shortest code in bytes wins

## Test cases

        x ->       out
1 ->         5
5 ->        49
49 ->       205
205 ->      5269
5269 ->      5369
5369 ->    266681
266681 ->   1077749
1077749 ->   9778141
9778141 ->   1968329
1968329 -> 239437889
239437889 -> 240505109

• Is the 9778141 -> 1968329 testcase correct? From my understanding of the problem, the input should always be smaller than the output (?) – Delfad0r Apr 25 at 21:32
• @Delfad0r My bad, I thought the sequence was entirely increasing. I've changed it to read that you should output the next term (i.e. given the numerator of $H_{n,2}$, output the numerator of $H_{n+1,2}$), which is what I originally intended – caird coinheringaahing Apr 25 at 21:35
• @LuisMendo I'm not sure if the Wolstenholme numbers repeat; answers may assume that they don't, and if it turns out that there are any repeats, any and all behaviour is acceptable. I chose "next Wolstenholme number" because I thought that the sequence was strictly increasing, and would be more interesting than the generic [sequence] I/O standards. That's (irritatingly) not the case, but changing it now would be unfair to the existing answers – caird coinheringaahing Apr 25 at 21:49
• @LuisMendo Added the rule that answerers may assume no duplicates exist, and that they don't have to handle any that do. Thanks for the feedback! – caird coinheringaahing Apr 25 at 21:54
• Related Math.SE question, asked by Jonah – caird coinheringaahing Apr 26 at 12:53

# J, 39 35 30 28 bytes

(]{~1+i.~)1*.[:+/\_2^~1+i.,]


Try it online!

-5 thanks to Bubbler's observation that the numerator is the LCM of the fraction and 1.

-2 thanks to Bubbler for a transformation.

Note: The approach below relies on the assumption that every Wolstenholme number is greater than or equal its index, which was kindly proved by Greg Martin.

• [:+/\_2^~1+i.,] Generates the sequence up to n+1, where n is the input: Scan sum +/\ of 1 over 1% the square of *: one plus 1...+ the integers from 0 to the input i.@>:.
• 1*. Numerators only.
• (]{~1+i.~) Finds index of the input in previous result, adds 1, and plucks from previous result.
• If I understand correclly, this assumes that the n-th term of the sequence is at least n. That seems to be the case, but a proof would be needed to be sure – Luis Mendo Apr 25 at 22:11
• If have a feeling that your 39-byte version and my answer do the same: compute terms until the input is found, then one more :-) – Luis Mendo Apr 25 at 22:30
• The numerator of a fraction is equal to its LCM with 1. – Bubbler Apr 26 at 3:19
• Also the part on the right of  +/\  can be golfed into _2^~1+i.,]. – Bubbler Apr 26 at 4:10
• I've added a comment below the question linking to your Math.SE question for greater visibility! – caird coinheringaahing Apr 26 at 12:53

# JavaScript (ES6), 98 bytes

f=(n,i=0)=>(g=(p=!(j=k=++i),q=1)=>j?g(p,q*j*j--):k?g(x=p+q/k/k--,q):q?g(q,p%q):x/p)()-n?f(n,i):g()


Try it online!

### Commented

f = (                           // f is a recursive function taking:
n,                            //   n = input
i = 0                         //   i = counter, initialized to 0
) =>                            //
( g = (                         // g is a recursive function which computes the
// next Wolstenholme number in 3 steps
p = !(j = k = ++i),         // increment i and start with j = k = i, p = 0
q = 1                       // and q = 1
) =>                          //
j ?                           // [step 1] while j is not equal to 0:
g(p, q * j * j--)           //   multiply q by j² and decrement j
:                             // then:
k ?                         //   [step 2] while k is not equal to 0:
g(x = p + q / k / k--, q) //     add q / k² to p, copy p to x and decrement k
:                           //   then:
q ?                       //     [step 3] while q is not equal to 0:
g(q, p % q)             //       copy q to p and p % q to q
:                         //     then:
//       p is now equal to the GCD of x and the final
//       value of q at the end of the 1st step
x / p                   //       return x / p
)()                             // initial call to g
- n ?                           // if it's not equal to n:
f(n, i)                       //   try again until it is
:                               // else:
g()                           //   return the next Wolstenholme number


# Haskell, 78… 76 bytes

• -2 bytes thanks to kops.
f n=snd(span(/=n)$head.words.show<$>scanl1(\i j->i+1/j^2::Rational)[1..])!!1


The relevant function is f, which takes n as a string of digits, and returns the next Wolstenholme number as a string of digits.

Try it online!

• GHC.Real also provides numerator for -2 bytes (but restricted to GHC which seems fine). Separately I really thought I'd be able to get something shorter without Data.Ratio but sadly only managed to tie 78 – kops Apr 26 at 1:13
• Nevermind, words to the rescue for 76! – kops Apr 26 at 1:26
• @kops That's brilliant, thanks! Getting around import Data.Ratio is always challenging, but this show trick is very nice :) – Delfad0r Apr 26 at 7:25

# Raku, 54 bytes

{first *>$_,map *.nude[0],[\+] map *⁻²,1.FatRat..*}  Try it online! • [\+] map *⁻², 1.FatRat .. * is the harmonic sequence. The members are of type FatRat, a "fat rational" type that does not degrade to floating point. (I could skip the .FatRat for seven bytes and use ordinary rationals, which give the correct answer for all of the provided test cases, but the problem states that the solution should theoretically valid for any input, so...) • map *.nude[0] maps that sequence into its numerators. *.numerator would be the more straightforward way to do this, but it's shorter to use .nude which returns the numerator and denominator as a pair. • first * >$_ finds the first element of the numerator sequence that exceeds the function argument $_. # MATL, 30 26 bytes 0@:Utptb/swyZd/ht1_)G-]0)  Try it online! ### Explanation 0 % Push 0. This initiallizes the array of numerators  % Do...while @:U % Push [1, 2, 4, ..., n^2] where n is iteration index tp % Duplicate, product. Gives denominator (not reduced) tb % Duplicate, bubble up in stack. Moves copy of [1, 2, 4, ..., n^2] to top /s % Divide element-wise and sum. Gives numerator (not reduced) wyZd % Swap, copy from below, GCD / % Divide. This gives the reduced numerator h % Concatenate to vecor of previos numerators t1_) % Duplicate, get second-to last entry G- % Subtract input ] % End. A new iteration is run if top of the stack is non-zero 0) % Get last entry. Implicit display  # Husk, 14 bytes ←t↓≠⁰mois∫m□İ\  ## Explanation ←t↓≠⁰mois∫m□İ\ İ\ infinite list [1,1/2,1/3,... m□ map to squares ∫ cumulative sum mo map to is first integer in string representation (gets numerator) ↓ drop elements till ≠⁰ input is reached t remove first item ← take the first item  # R, 112 bytes x=scan();while({d=prod((1:T)^2);n=sum(d/(1:T)^2);i=1;while((i=i+1)<n)while(!n%%i&!d%%i)n=n/i;x!=F}){T=T+1;F=n};n  Try it online! How? x=scan(); # get previous wolstenholme number x while({...do this...;x!=F}) # first 'do this', then check if x is equal to F (initially zero) { # 'do this': d=prod((1:T)^2) # denominator d = product of squares of 1..T n=sum(d/(1:T)^2) # numerator n = sum of d/1, d/4, d/9, d/16 ... d/T^2 # n/d is T-th harmonic number of order 2 # but we need to simplify the fraction: i=1;while((i=i+1)<n) # try each i from 2..(n-1): while(!n%%i&!d%%i)n=n/i # if i divides both n and d, divide n by i # at this point, n equals the simplified numerator # which is the T-th wolstenholme number } ){ # if x isn't equal to F yet T=T+1 # increase T F=n # and set F to the last wolstenholme number n } # if we got here, x was equal to F # and we've already updated n to be the next # wolstenholme number n # so we output n.  A vectorized version is shorter - 98 bytes - but quickly runs-out of memory as it tries to allocate a vector the size of the non-reduced denominator for order-2 harmonic numbers beyond the first 7... # Stax, 15 bytes ┬Z°#¡Lk╢╙°m╕%ºb  Run and debug it Reduced from 22 with help from recursive. ## Explanation {iiGx=!w^G}Z{Ju+Fr { w while result is falsy: ii push iteration index twice G GoTo unclosed curly brace } label: Z push a zero under the index { F for 1..i Ju+ square, invert, add return to GoTo r get numerator x=! is it equal to the input? ^ if so, increment G and GoTo label  # Ruby, 58 bytes ->x{r=z=1r;1until[z,z+=(r+=1)**-2]*?_=~/^#{x}.+_(\d+)/;$1}


Try it online!

# Python 2, 120 bytes

from fractions import*
x=input()
k=1
g=0
n=[]
while x not in n[:-1]:g+=Fraction(1,k*k);k+=1;n+=[g.numerator]
print n[-1]


Try it online!

fixed thanks to Jonathan Allan

• 99 bytes (initialize g as an integer, keep track of the numerator in a variable) – ovs Apr 25 at 21:18
• n<=x wont work for an input of $a(10)=1968329$ since $a(9)>a(10)$ TIO (should output 239437889). – Jonathan Allan Apr 26 at 11:36
• This would work but is probably golfable. – Jonathan Allan Apr 26 at 11:38
• @JonathanAllan Oh, thanks for catching that. – hyper-neutrino Apr 26 at 20:47

# Ruby, 61 bytes

Uses Ruby's Rational type.

->x{w=l=n=k=0r
(l=n;k+=1;w+=1/k/k;n=w.numerator)while l!=x
n}


Try it online!

73 bytes with integers:

->x{d=1
n=l=k=w=0
(l=n;k+=1;w*=k*k;w+=d;d*=k*k;n=w/w.gcd(d))while l!=x
n}


Try it online!

# Charcoal, 57 bytes

Ｎθ≔⁰η≔⁰ζＷ⁻θζ«≔ηζ⊞υＸ⊕Ｌυ²≔Πυε≔Σ÷ευδＷδ«≔﹪εκδ≔κε»≔÷Σ÷Πυυεη»Ｉη


Try it online! Link is to verbose version of code. Explanation:

Ｎθ


Input the target numerator.

≔⁰η≔⁰ζ


Initialise the next and current numerator to zero.

Ｗ⁻θζ«


Repeat until the current numerator is the target.

≔ηζ


Save the next numerator as the current numerator.

⊞υＸ⊕Ｌυ²


Add the next square to the predefined empty list.

≔Πυε≔Σ÷ευδ


Get the potential denominator and numerator of the sum by multiplying all of the fractions by the product of all the squares so far.

Ｗδ«≔﹪εκδ≔κε»


Get the GCD.

≔÷Σ÷Πυυεη


Divide the potential numerator by the GCD and save it as the next numerator.

»Ｉη


Output the next numerator.

# Jelly, (16?) 21 bytes

²€P©:$S:g¥® 1Ç⁼¥1#Ḣ‘Ç  Try it online! I can't prove that $$\a(n) \ge n\$$, but if so this 16 byter works too: ‘²€P©:$S:g¥®)ḟRḢ


### How?

²€P©:$S:g¥® - Link 1, get a(n): integer n ²€ - square each (of [1..n])$      - last two links as a monad - f(squares):
P         -   product (of the squares)
©        -   (copy to the register for later)
:       -   integer divide by (each square)
S     - sum
¥® - last two links as a dyad - f(sum, register)
g   - (sum) gcd (register)
:    - (sum) integer divided by (that)

1Ç⁼¥1#Ḣ‘Ç - Main Link: integer, a(z)
1         - set L=1
1#    - collect the first 1 integer(s) counting up from k=L for which:
¥      -   last two links as a dyad - f(k,a(z))
Ç        -     call Link 1 (k) -> a(k)
⁼       -     equals (a(z))?
Ḣ   - head -> z
‘  - increment -> z+1
Ç - call Link 1 (z+1) -> a(z+1)


# Vyxal, 33 30 bytes

1:{›Ȯ?≤|:λ¡²nɾ²/M∑:n¡²ġ/;†^_}Ȯ


Try it Online!

Let us try to maintain a stack with 2 elements $$\n\$$ and $$\f(n)\$$ after every step, with $$\f(n)\$$ being the required numerator.

1:     # Pushes [1, 1] to the stack.
{      # While...
›      # Increment the first value of stack -> [n+1, f(n)]
Ȯ?≤    # While condition: Is the second value of stack less than or equal to the input?
|      # If yes, then...
:      # Duplicate the top value of stack. Now it contains [n+1, n+1, f(n)]
λ      # Lambda function
¡²    # Factorial of n, squared.
nɾ²    # Inclusive range [1, n], and square every number.
/M∑    # Divide every time with first argument and return the sum.
:n¡²ġ/ # Duplicate this value, then divide it with it's gcd with factorial n, squared
;      # Return value.
†      # Call the function. Now the stack contains [f(n+1), n+1, f(n)]
^      # Reverse stack -> [f(n), n+1, f(n+1)]
_      # Pop from stack -> [n+1, f(n+1)]
}      # End while
Ȯ      # Print second item of stack.


-3 bytes thanks to Lyxal

• 1+› for one byte off – lyxal Apr 26 at 6:49
• Vyxal has ƒ(fractionify) to convert decimals to fraction arrays. may be helpful here. – Razetime Apr 26 at 6:49
• Also the first n in the lambda isn't needed because it automatically starts with its argument on the stack – lyxal Apr 26 at 6:52
• Further, 1nṡnɾ for another byte off – lyxal Apr 26 at 6:59
• Down to 29 bytes – lyxal Apr 26 at 23:59

# Python 3, 126 bytes

from math import*
def f(n,i=1,k=1):
k*=i*i
q=sum(k//j**2for j in range(1,i+1))
x=q/gcd(k,q)
return x>n and x or f(n,i+1,k)


Try it online!

Let's find the numerator if the denominator is $$\1^2 + 2^2 + ... i^2\$$. The numerator of the $$\j^{th}\$$ term will be just the denominator divided by $$\j^2\$$. So we add the numerator of all these terms and then divide it with the gcd of numerator and denominator and compare it with $$\n\$$.

# Wolfram Language (Mathematica), 54 bytes

(k=1;While[(t=Numerator[++k~HarmonicNumber~2])<=#];t)&


Try it online!

# Clojure, 90 bytes

#(loop[x 3 y 5/4](let[z(+(/ 1 x x)y)n numerator](condp = % 1 5(n y)(n z)(recur(inc x)z))))


Try it online!

Like many other languages, Clojure has native support for rational fractions, but with an annoying quirk: numerator and denominator functions throw exceptions on integer arguments. If that's not enough, even rational literals of the form x/y are auto-converted to integers wherever possible, and the only way I found to force evaluation as a rational is clojure.lang.Numbers/toRatio, which is a golfing nightmare.

This causes a problem with the initial term $$\1\$$ of the harmonic sequence, which necessitates a hack - we start the loop from the second term, and cover the first one with a special condition.