[Question]
  [
Your task is to implement integer sequence [A130826](https://oeis.org/A130826):
> 
> **an** is the smallest positive integer such that **an - n** is an entire multiple of **3** and twice the number of divisors of **(an - n) / 3** gives the **n**th term in the first differences of the sequence produced by the Flavius Josephus sieve.
> 
> 
> 
Lost yet? Well, it's actually quite easy.
The [Flavius Josephus sieve](https://oeis.org/A000960) defines an integer sequence as follows.
1. Start with the sequence of positive integers and set **k = 2**.
2. Remove every **k**th integer of the sequence, starting with the **k**th.
3. Increment **k** and go back to step 2.
**fn** is the **n**th integer (1-indexed) that never gets removed.
If – as usual – **σ0(k)** denotes the number of positive divisors of the integer **k**, we can define **an** as the smallest positive integer such that **2σ0((an - n) / 3) = fn+1 - fn**.
### Challenge
Write a program or function that takes a positive integer **n** as input and prints or returns **an**.
Standard [code-golf](/questions/tagged/code-golf "show questions tagged 'code-golf'") rules apply. May the shortest code win!
### Worked examples
If we remove every second element of the positive integers, we are left with
```
 1  3  5  7  9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 ...
```
After removing every third element of the remainder, we get
```
 1  3  7  9 13 15 19 21 25 27 31 33 37 39 ...
```
Now, removing every fourth, then fifth, then sixth element gets us
```
 1  3  7 13 15 19 25 27 31 37 39 ...
 1  3  7 13 19 25 27 31 39 ...
 1  3  7 13 19 27 31 39 ...
 1  3  7 13 19 27 39 ...
```
The last row shows the terms **f1** to **f7**.
The differences of the consecutive elements of the these terms are
```
 2  4  6  6  8 12
```
Dividing these forward differences by **2**, we get
```
 1  2  3  3  4  6 
```
These are the target divisor counts.
* **4** is the first integer **k** such that **σ0((k - 1) / 3) = 1**. In fact, **σ0(1) = 1**.
* **8** is the first integer **k** such that **σ0((k - 2) / 3) = 2**. In fact, **σ0(2) = 2**.
* **15** is the first integer **k** such that **σ0((k - 3) / 3) = 3**. In fact, **σ0(4) = 3**.
* **16** is the first integer **k** such that **σ0((k - 4) / 3) = 3**. In fact, **σ0(4) = 3**.
* **23** is the first integer **k** such that **σ0((k - 5) / 3) = 4**. In fact, **σ0(6) = 4**.
* **42** is the first integer **k** such that **σ0((k - 6) / 3) = 6**. In fact, **σ0(12) = 6**.
### Test cases
```
   n     a(n)
   1        4
   2        8
   3       15
   4       16
   5       23
   6       42
   7       55
   8      200
   9       81
  10       46
  11      119
  12      192
  13      205
  14   196622
  15    12303
  16       88
  17      449
  18      558
  19      127
  20     1748
  21   786453
  22       58
  23     2183
  24     3096
  25     1105
  26   786458
  27 12582939
  28      568
  29     2189
  30     2730
```
  
[Answer]
# Jelly, ~~30~~ ~~29~~ ~~27~~ 25 bytes
*Saved 2 bytes thanks to @Dennis notifying me about `Æd`, and another 2 bytes for combining the two chains.*
```
RUð÷‘Ċ×µ/
‘Ç_ÇH0Æd=¥1#×3+
```
[Try it online!](https://tio.run/nexus/jelly#@x8UenjD4e2PGmYc6To8/dBWfS4g83B7/OF2D4PDbSm2h5YaKh@ebqz9//9/MwA)
This was probably the most fun I've ever had with Jelly. I started from the second line, which calculates **fn** from **n** using the formula on OEIS, and is quite beautiful.
### Explanation
```
RUð÷‘Ċ×µ/    Helper link to calculate Fn. Argument: n
R            Get numbers [1..n]
 U           Reverse
        /    Reduce by "round up to next 2 multiples":
   ÷           Divide by the next number
    ‘          Increment to skip a multiple
     Ċ         Ceil (round up)
      ×        Multiply by the next number
‘Ç_ÇH0Æd=¥1#×3+    Main link. Argument: n
‘                  Increment n
 Ç                 Calculate Fn+1
   Ç               Calculate Fn
  _                Subtract
    H              Divide by 2
     0    1#       Starting from 0, find the first candidate for (an-n)/3
                   that satisfies...
      Æd             σ0((an-n)/3)
        =            = (Fn+1-Fn)/2
            ×3     Multiply by 3 to turn (an-n)/3 into an-n
              +    Add n to turn an-n into an
```
[Answer]
# [Python 2](https://docs.python.org/2/), ~~121~~ ~~119~~ 118 bytes
```
n=input();r=range(1,4**n);d=s,=r*1,
for k in r:del s[k::k+1];d+=sum(k%j<1for j in r)*2,
print d.index(s[n]-s[n-1])*3+n
```
Run time should be roughly **O(16n)** with **O(4n)** memory usage. Replacing `4**n` with `5<<n` – which I *think* is sufficient – would improve this dramatically, but I'm not convinced that it works for arbitrarily large values of **n**.
[Try it online!](https://tio.run/nexus/python2#HcpLCsIwEADQfU8xGyG/ClFxkTgnKVkIEyWNjmXSgrePn81bvc5YeNlWpaOgXPmelXcnY1hHwuZQjHfD7SVQoTBIoPyANtUQqvUpksW2PVXdzRf/S/M/aXNwwyKFV6B9Ycpv1SZO45fRJ22Olns/fwA "Python 2 – TIO Nexus")
### Asymptotic behavior and upper bounds of an
Define **bn** as **(an - n)/3**, i.e., the smallest positive integer **k** such that **σ0(k) = ½(fn+1 - fn)**.
As noted on the OEIS page, **fn ~ ¼πn2**, so **fn+1 - fn ~ ¼π(n + 1)2 - ¼πn2 = ¼π(2n + 1) ~ ½πn**.
This way, **½(fn+1 - fn) ~ ¼πn**. If the actual number is a prime **p**, the smallest positive integer with **p** divisors is **2p-1**, so **bn** can be approximated by **2cn**, where **cn ~ ¼πn**.
Therefore **bn < 4n** will hold for sufficiently large **n**, and given that **2¼πn < 2n << (2n)2 = 4n**, I'm confident there are no counterexamples.
### How it works
```
n=input();r=range(1,4**n);d=s,=r*1,
```
This sets up a few references for our iterative process.
* **n** is the user input: a positive integer.
* **r** is the list **[1, ..., 4n - 1]**.
* **s** is a copy of **r**.
Repeating the list once with `r*1` creates a shallow copy, so modifying **s** won't modify **r**.
* **d** is initialized as the tuple **(s)**.
This first value is not important. All others will hold divisor counts of positive integers.
```
for k in r:del s[k::k+1];d+=sum(k%j<1for j in r)*2,
```
For each integer **k** from **1** to **4n - 1**, we do the following.
* `del s[k::k+1]` takes every **(k + 1)**th integer in **s** – starting with the **(k + 1)**th – and deletes that slice from **s**.
This is a straightforward way of storing an initial interval of the Flavius Josephus sieve in **s**. It will compute much more than the required **n + 1** initial terms, but using a single `for` loop to update both **s** and **d** saves some bytes.
* `d+=sum(k%j<1for j in r)*2,` counts how many elements of **r** divide **k** evenly and appends **2σ0(k)** to **d**.
Since **d** was initialized as a singleton tuple, **2σ0(k)** is stored at index **k**.
```
print d.index(s[n]-s[n-1])*3+n
```
This finds the first index of **fn+1 - fn** in **d**, which is the smallest **k** such that **2σ0(k) = fn+1 - fn**, then computes **an** as **3k + 1** and prints the result.
[Answer]
# Java 8, 336, 305, 303, 287, 283 279 bytes
57 bytes removed thanks to Kritixi Lithos
Golfed
```
class f{static int g(int s,int N){return s<1?N+1:g(s-1,N+N/s);}static int h(int k){int u=0,t=1,i;for(;u!=(g(k,k)-g(k,k-1))/2;t++)for(i=1,u=0;i<=t;)if(t%i++<1)u++;return 3*t-3+k;}public static void main(String[]a){System.out.print(h(new java.util.Scanner(System.in).nextInt()));}}
```
Ungolfed
```
class f {
    static int g(int s,int N){return s < 1 ? N + 1 : g(s - 1, N + N / s);}
    static int h(int k) {
        int u = 0, t = 1, i;
        // get the first number with v divisors
        while(u != (g(k, k) - g(k, k - 1))/2){
            u = 0;
            for (i = 1; i <= t; i++)
                if (t % i < 1) u++;
            t++;
        }
        // 3*(t-1)+k = 3*t+k-3
        return 3 * t + k - 3;
    }
    public static void main(String[] a) {
        System.out.print(h(new java.util.Scanner(System.in).nextInt()));
    }
}
```
[Answer]
# Mathematica, ~~130~~ ~~116~~ ~~106~~ 103 bytes
```
3Catch@Do[f=#2⌈#/#2+1⌉&~Fold~Reverse@Range@#&;If[Tr[2+0Divisors@k]==f[#+1]-f@#,Throw@k],{k,∞}]+#&
```
or
```
3Catch@Do[f=#2⌈#/#2+1⌉&~Fold~Reverse@Range@#&;If[2DivisorSum[k,1&]==f[#+1]-f@#,Throw@k],{k,∞}]+#&
```
Ended up being almost identical to @Pietu1998's Jelly code...
**Explanation**
```
Catch@
```
`Catch` whatever is `Throw`-ed (thrown).
```
Do[ ... ,{k,∞}]
```
Infinite loop; `k` starts from `1` and increments every iteration.
```
f= ...
```
Assign `f`:
```
Reverse@Range@#
```
Find `{1, 2, ... , n}`. Reverse it.
```
#2⌈#/#2+1⌉&
```
A function that outputs ceil(n1/n2 + 1) \* n2
```
f= ... ~Fold~ ... &
```
Assign `f` a function that recursively applies the above function to the list from two steps above, using each output as the first input and each element of the list as the second input. The initial "output" (first input) is the first element of the list.
```
Tr[2+0Divisors@k]==f[#+1]-f@#
```
Check whether twice the number of divisors of `k` is equal to f(n + 1) - f(n).
```
If[ ... ,Throw@k]
```
If the condition is `True`, `Throw` the value of `k`. If not, continue looping.
```
3 ... +#&
```
Multiply the output by 3 and add n.
**130 byte version**
```
Catch@Do[s=#+1;a=k-#;If[3∣a&&2DivisorSigma[0,a/3]==Differences[Nest[i=1;Drop[#,++i;;;;i]&,Range[s^2],s]][[#]],Throw@k],{k,∞}]&
```
[Answer]
# [Perl 6](http://perl6.org), ~~154~~ ~~149~~ ~~136~~ 107 bytes
```
->\n{n+3*first ->\o{([-] ->\m{m??&?BLOCK(m-1).rotor(m+0=>1).flat!!1..*}(n)[n,n-1])/2==grep o%%*,1..o},^Inf}
```
Ungolfed:
```
-> \n {                    # Anonymous sub taking argument n
  n + 3 * first -> \o {    # n plus thrice the first integer satisfying:
    (                      #
      [-]                  #
      -> \m {              # Compute nth sieve iteration:
        m                  # If m is nonzero,
          ?? &?BLOCK(m-1).rotor(m+0=>1).flat # then recurse and remove every (m+1)-th element;
          !! 1..*          # the base case is all of the positive integers
      }                    #
      (n)                  # Get the nth sieve
      [n,n-1]              # Get the difference between the nth and (n-1)th elements (via the [-] reduction operator above)
    ) / 2                  # and divide by 2;
    ==                     # We want the number that equals
    grep o %% *, 1..o      # the number of divisors of o.
  }
  ,^Inf
}
```
[Answer]
# [05AB1E](https://github.com/Adriandmen/05AB1E), ~~35~~ ~~34~~ 39 bytes
```
1Qi4ë[N3*¹+NÑg·¹D>‚vyy<LRvy/>îy*}}‚Æ(Q#
```
It looks awful, so is runtime performance. It takes several seconds for input that yield small values. Don't try numbers like 14; it will eventually find the result but it would take ages.
## Explanation
It works as 2 sequentially called programs. The first one computes **Fn+1 - Fn** and the second one determines **an** based on its definition, using a bruteforce approach.
**Fn+1 - Fn** is evaluated for each iteration even though it is loop invariant. It makes the code time inefficient, but it makes the code shorter.
[Try it online!](https://tio.run/nexus/05ab1e#ATMAzP//xb5ITMK5wrXCvERnTMK@PiUww4rDj33CpcK5PMOoVVtOMyrCuStOw5FnwrdYUSP//zg "05AB1E – TIO Nexus")
]