[Question]
  [
Given a number \$n\$, we have its sum of divisors, \$\sigma(n)\ = \sum\_{d | n} {d}\$, that is, the sum of all numbers which divide \$n\$ (including \$1\$ and \$n\$). For example, \$\sigma(28) = 1 + 2 + 4 + 7 + 14 + 28 = 56\$. This is OEIS [A000203](https://oeis.org/A000203).
We can now define the sum of sum of divisors as \$S(n) = \sum\_{i=1}^{n}{\sigma(i)}\$, the sum of \$\sigma(i)\$ for all numbers from \$1\$ to \$n\$. This is OEIS [A024916](https://oeis.org/A024916).
Your task is to calculate \$S(n)\$, **in time sublinear in \$n\$, \$o(n)\$**.
# Test cases
```
10 -> 87
100 -> 8299
123 -> 12460
625 -> 321560
1000 -> 823081
1000000 (10^6) -> 822468118437
1000000000 (10^9) -> 822467034112360628
```
# Rules
* **Your complexity must be \$o(n)\$**. That is, if your code takes time \$T(n)\$ for input \$n\$, you must have \$\lim\_{n\to\infty}\frac{T(n)}n = 0\$. Examples of valid time complexities are \$O(\frac n{\log(n)})\$, \$O(\sqrt n)\$, \$O(n^\frac57 \log^4(n))\$, etc.
* You can use any reasonable I/O format.
* Note that due to the limited complexity you can't take the input in unary nor output in it (because then the I/O takes \$\Omega(n)\$ time), and the challenge might be impossible in some languages.
* Your algorithm should in theory be correct for all inputs, but it's fine if it fails for some of the big test cases (due to overflow or floating-point inaccuracies, for example).
* [Standard loopholes](https://codegolf.meta.stackexchange.com/questions/1061/loopholes-that-are-forbidden-by-default) are disallowed.
This is code golf, so the shortest answer in each language wins.
  
[Answer]
# [Ruby](https://www.ruby-lang.org/), ~~55~~ 52 bytes
Formula is adapted from [this math.SE question](https://math.stackexchange.com/q/3287652). Runs in \$\mathcal{O}(\sqrt{n})\$.
-3 bytes thanks to @dingledooper.
```
->x{(1..$.=x**0.5).sum{(v=x/_1)*_1-(v*~v-$.*~$.)/2}}
```
[Attempt This Online!](https://ato.pxeger.com/run?1=m72kqDSpcsHiNNuYpaUlaboWN0107SqqNQz19FT0bCu0tAz0TDX1iktzqzXKbCv04w01teINdTXKtOrKdFX0tOpU9DT1jWproXqdCxTSog0NYrkgNIxhZAxhmBmZwqS0tIwRTDME0zIWYtSCBRAaAA)
[Answer]
# [Python 2](https://docs.python.org/2/), 58 bytes
Runs in \$ \mathcal{O}(\sqrt{n}) \$ time.
```
f=lambda n,i=1:~i*i*(n/-~i-n/i)/2+(i<n and f(n,n/(n/-~i)))
```
[Try it online!](https://tio.run/##TY3NboMwEITvfoq9xaYQvDZ/@aGvUtGGKCvBYoFRlEtenZjSVr2sZvabnXUPfxvYLNS7YfQwPabTOvZT68f2ax4nGrijnrxEHUW5Wq511/SflwY4phqPT4ookpwmT0o4JZWaN0lnhoYvcJUcc7pBpdRyv1HXAh4FELsYhtk7qKFvnCT2MYzN/SOA2Uu1n1wXXu4geYedUgLcGCKhMHC1XS6oV1qVAvWmzOEg0NhVo8kKLQqTr8YazIMLsZ@c1RV@279NiFeIVWbL3/1/VGqbYWgudGGqFw "Python 2 – Try It Online")
## Explanation
Note: [@Neil's answer](https://codegolf.stackexchange.com/a/265132/88546) seems to use a similar approach to mine.
We first take the \$ O(n) \$-time formula \$ \sum\_{i=1}^n{i \cdot \lfloor\frac{n}{i}\rfloor} \$, noticing that there can only be \$ O(\sqrt{n}) \$ unique values of \$ \lfloor\frac{n}{i}\rfloor \$. For example, if \$ n = 50 \$, the possible values of \$ i \$ for each \$ \lfloor\frac{n}{i}\rfloor \$ are:
```
 n//i |   i
------+-------
   1  | 26-50
   2  | 17-25
   3  | 13-16
   4  | 11-12
   5  | 9-10
   6  | 8-8
   7  | 7-7
   8  | 6-6
  10  | 5-5
  12  | 4-4
  16  | 3-3
  25  | 2-2
  50  | 1-1
```
Let \$ a \$ be a sorted list of all unique values of \$ \lfloor\frac{n}{i}\rfloor \$. For this example, \$ a = [ 1, 2, 3, 4, 5, 6, 7, 8, 10, 12, 16, 25, 50 ] \$. We can then construct a new summation that solves the problem in \$ O(\sqrt{n}) \$:
$$ \sum\_{i=1}^{|a|}{\left\lfloor \frac{n}{a\_{i}} \right\rfloor \times (T(a\_{i}) - T(a\_{i-1}))} $$
where \$ T(n) \$ is the sum of the first \$ n \$ positive integers, and assuming that \$ a\_0 = 0 \$. Next, a bit of insight allows us to rearrange the summation like so:
$$ \sum\_{i=1}^{|a|}{T(a\_i) \times \left(\left\lfloor\frac{n}{a\_i}\right\rfloor - \left\lfloor\frac{n}{a\_i+1}\right\rfloor\right)} $$
To compute \$ a \$, we can apply the following recurrence:
$$ a\_i = \left\lfloor\frac{n}{\lfloor n / (1+a\_{i-1}) \rfloor}\right\rfloor $$
The proof is left as an exercise to the reader. Finally, we can compute \$ T(n) \$ using the formula \$ \frac{n(n+1)}{2} \$.
[Answer]
# [R](https://www.r-project.org), 52 bytes
```
\(x,y=1:x^.5,z=x%/%y)y%*%z+(z+1)%*%z/2-max(y)*sum(y)
```
[Attempt This Online!](https://ato.pxeger.com/run?1=m72waMEiF9ulpSVpuhY3TWI0KnQqbQ2tKuL0THWqbCtU9VUrNStVtVSrtDWqtA01QSx9I93cxAqNSk2t4tJcIAXRu684saAgp1IjWcPQQMfQAIwhhIGmjgtU0YIFEBoA)
Golfed down version of [Dominic's answer](https://codegolf.stackexchange.com/a/265126/67312) which in turn uses the same observation as ovs. Nothing crazy here; should run in \$\mathcal{O}\left(\sqrt n\right)\$ time, but saves bytes by using `1:sqrt(x)` to truncate, and `%*%` to calculate the sumproduct.
[Answer]
# [Nibbles](http://golfscript.com/nibbles/index.html), 17 bytes (34 nibbles)
```
+:.;.,;^$-2/_$-;~$/*+$~$~$;$!$,_*
```
[Attempt This Online!](https://ato.pxeger.com/run?1=m708LzMpKSe1eMGCpaUlaboWNxW1rfSs9XSs41R0jfTjVXSt61T0tbRV6oDQWkVRRSdeC6IQqn7BTQVDAy5DAzCGEFASRhlAFAIA)
Direct port of my R answer.
This feels very clunky, with half the tokens used for control-flow (`,;,;;~`) or explicit variable names (`$_$$$$$;$$_`); I suspect that one of the stack-based languages will be able to better this...
**Bit-by-bit** (in conjunction with the Nibbles [quick reference](http://golfscript.com/nibbles/quickref.html) and [tutorial](http://golfscript.com/nibbles/tutorial_basics.html)):  
y (=floor(sqrt(x)): `^$-2`  
z (=floor(x/(1..y))), saving y: `.,;^$-2/_$` (y saved as `$`)  
triangle(z), saving z, y & triangle function: `.;.,;^$-2/_$;~$/*+$~$~`  
triangle(z)-triangle(y): `.;.,;^$-2/_$-;~$/*+$~$~$;$`  
(leaves z saved as `$`, y saved as `@`, x is `_`)
So now we can easily use the saved variables to calculate D:
D = (1..y)\*z + triangle(z)-triangle(y): `+: .;.,;^$-2/_$-;~$/*+$~$~$;$ !$,_*`
[Answer]
# [R](https://www.r-project.org), ~~61~~ ~~51~~ 49 bytes
*Edit: saved ~~10~~ 12 bytes by re-arrangement, goaded (and inspired) by Giuseppe golfing-down my original version in [his answer](https://codegolf.stackexchange.com/a/265138/95126)*
```
\(x,y=1:x^.5,z=x%/%y)min((2*y+z+1)%*%z-y^3-y^2)/2
```
[Attempt This Online!](https://ato.pxeger.com/run?1=m72waMEiF9ulpSVpuhY3DWM0KnQqbQ2tKuL0THWqbCtU9VUrNXMz8zQ0jLQqtau0DTVVtVSrdCvjjIHYSFPfCKpRoTixoCCnUiNZw9BAx9AAjCEEEGjquGhC1C1YAKEB)
Golfed calculation based on the same [Math.SE question](https://math.stackexchange.com/q/3287652) used in [ovs' Ruby answer](https://codegolf.stackexchange.com/a/265116).
**Ungolfed**
```
triangle=function(x)x*(x+1)/2
D=
function(x){
    y=sqrt(x)%/%1
    k=1:y
    sum(k*x%/%k+triangle(x%/%k))-triangle(y)*y
}
```
[Answer]
# [Scala 3](https://www.scala-lang.org/), 75 bytes
Port of [@ovs's Ruby answer](https://codegolf.stackexchange.com/a/265116/110802) in Scala.
---
```
x=>{val d=sqrt(x);(BigInt(1) to d).map(k=>x-x%k-(x/k* ~(x/k)-d* ~d)/2).sum}
```
[Attempt This Online!](https://ato.pxeger.com/run?1=jZCxTsMwEIbVNU9xSyW7NG3TCoSCYgk2hJgQE2Jwm6SYpG5JXDCqwouwdIF3ok_DOU5SgWjE4vPd_9_n871_5jOe8sl21-kup4_RTME1FxI2DkAYxbDAhPBsnvtwnmX89e5GZULO76kPt1IoCEonwAqrKpUkJhdifikV8UaU0kNSizaeHNROxsdtzNYHR__QG0vh4PHMU4h96wiYjfjdj7WK3dOvKx2wjbGEQf6UKaLpWUOjoJYQ0sGCr0gSMO3qbuISPUx68GYCdUO8hXQ4poN8vSgsctdxnWrpluiDBdL60izbeGSkFZG11Afxl51IOAIBQ5AUGAPPKdvN2EsZod5MXAsZSCyjuM89LNjH-qAr48uDSCMgBEUX2QM-zYGZNlpNCCVHenXyC2JqxQ8UOnrGxlDeMxC_n6ZqKKGFU21tu7XxGw)
[Answer]
# [Charcoal](https://github.com/somebody1234/Charcoal), ~~41~~ 38 bytes
```
 Ｎθ≔θη≔⁰ζＷη«≔÷θιε≔÷θ⊕εη≧⁺×⁻ιη÷×ε⊕ε²ζ»Ｉζ
```
[Try it online!](https://tio.run/##bY7NCoMwEITP9Sly3IAF26un0l48WKT0BVJdzEJMNT8WLH321KClF/ewMHwzs1tLYeqnUCEUuvfu6rsHGhh4npyspVbDkDL5V1nKplm9JClkIDl7J7sVFdpdaKQGY4R4ynA2bsJC1wY71A4bQM6XA7tS9Iv5Rq10UClvU3anDi2UpL0FisaY/lUtEDcKj3HFRz9JZUg7OAvrYOI8D@GQrRP2o/oC "Charcoal – Try It Online") Link is to verbose version of code. Explanation: Uses the formula \$ A024916(n) = \sum\_{i=1}^n i \lfloor n/i \rfloor \$ with the linked observation that it can be calculated more efficiently by grouping terms with the same value of \$ \lfloor n/i \rfloor \$, and then applying @dingledooper's rearrangement that \$ \sum\_{j=1} \lfloor n/a\_j \rfloor (T(a\_j) - T(a\_{j-1})) = \sum\_{j=1} T(a\_j) ( \lfloor n/a\_j \rfloor - \lfloor n/a\_{j+1} \rfloor ) \$ where \$ T(n) \$ are the triangular numbers to save 3 bytes.
```
 Ｎθ
```
Input \$ n \$.
```
≔θη
```
Start counting \$ \lfloor n / a\_j \rfloor \$ down from \$ n \$.
```
≔⁰ζ
```
Start with no total.
```
Ｗη«
```
Repeat until \$ \lfloor n / a\_j \rfloor \$ is zero. This also conveniently makes a copy of the value, so I can update the original variable in the loop while still being able to access its former value.
```
≔÷θιε
```
Calculate \$ a\_j \$.
```
≔÷θ⊕εη
```
Calculate \$ \lfloor n / a\_{j+1} \rfloor = \lfloor n / ( 1 + a\_j ) \rfloor \$.
```
≧⁺×⁻ιη÷×ε⊕ε²ζ
```
Calculate \$ T(a\_j) ( \lfloor n/a\_j \rfloor - \lfloor n/a\_{j+1} \rfloor ) \$, and accumulate it to the total.
```
»Ｉζ
```
Output the final total.
On TIO, \$ 10^{11} \$ takes about \$ 10 \$ times as long as \$ 10^9 \$, so the time complexity appears to be \$ O(\sqrt n) \$.
[Answer]
# [Python](https://www.python.org), 74 bytes
```
lambda n,i=0,k=1:(k>=i)*i*k+(i<k and f(n,j:=i+1,l:=n//j)+i*(l*~l-k*~k)//2)
```
[Attempt This Online!](https://ato.pxeger.com/run?1=RU5LboMwEN37FLO0wQTbEH6qc4AqR6ha0QbUCWAQn0U2WfYS3USq2jvlNrVLo4w00rw3b-a9z5_hNL_35vJV66fvZa6D7PrYlt3roQTDUQveaFnQZqeReeg1PsWHBkpzgJoafiw0-pK3hTZheGQ-erT1zm3QeOeGhaFi_x8X7IZ-nGE6TcT2ZqrmsXpbxgl702KHM5XCFSN1P8Ie0EA_VIYKVhBA7r3wHjTsN9PQWi1znCO6cqBoZk4tZJYdRosorSkyR6zml-uHFBDsIEuJdfmbVJ4TqSI3SxUngiRq60Ck5NYiF2bVRSKTZM0mwIZ8Tti6sFeZlFkcpbf1TZHfFamIYml9EpGojKxxfgE)
### Old [Python](https://www.python.org), 88 bytes
```
lambda n,k=0:sum(i*k*(k>i)+i*(k-~(j:=n//-~i))*(k-(k:=j))//2for i in range(int(n**.5)+1))
```
[Attempt This Online!](https://ato.pxeger.com/run?1=RY5NboMwEEb3nMLLGQLBY_4MEjkBB-iiakXV0DoEgwhZdJNlL9FNpKq9U25TT0mUled78_zZX7_jx_w-2PN3Wz3-HOc21JeHfdO_vDbCBl0ly8OxB-N3PnQbgyvjzvAEu7KyURSeDCID6MpqhxhFqh0mYYSxYmrs2xaMncH6_jrFFSFe-2uWapaGcWtBYukJE_jPwSAqUa8P497MgMwY9M3INQG4iI6OE5dCCwYZLJ3nyydJEW6Ezj2Sy6SKwiMV80wqyaSXqZRDrCh1yWlXL5aa_iMTIPmU4bJwtzSRTuL8tr4Zxd3IZZyQeyeTmdLe8p0_)
Probably similar to @ovs's although found independently.
The basic idea is to regroup the given double sum by divisors and then count how many numbers in the range are a multiple of the divisor.
To achieve O(n^0.5) divisors with the same multiplicity are lumped together.
[Answer]
# JavaScript (ES6), 54 bytes
Based on the formula [found by ovs](https://codegolf.stackexchange.com/a/265116/58563).
```
x=>(g=k=>++k*k>x?--k*k*~k/2:k*(q=x/k|0)-q*~q/2+g(k))``
```
[Try it online!](https://tio.run/##bY1NDoIwGAX3noJlW1LbfuWnmBSvAkEgWkKtGMPCcPUKIRqjvM1bzCRzKR/lUN3O1zvt7an2jfajzlGrjc7D0BCTj0dK5yeTYXAwBDk9MvPkmDoyOQZhiwzGReEr2w@2q/edbVGDgmWCYxwwFqh0908Ff1PIsi0OcuUCooRvCAnEqyBBxFvGXPgkJFfix1jwlzFXlBAqkql/AQ "JavaScript (Node.js) – Try It Online")
---
# JavaScript (ES11), 45 bytes
*-3 bytes thanks to [@tsh](https://codegolf.stackexchange.com/users/44718/tsh)*
A port of [dingledooper's answer](https://codegolf.stackexchange.com/a/265137/58563).
Expects a BigInt.
```
f=(n,q=n)=>q&&(i=n/q)*(q-n/++i)*i/2n+f(n,n/i)
```
[Try it online!](https://tio.run/##bY7LDsIgEEX3fgWrBlorDPRBF/gvTW0NphmKNf4@oqmPGGZ7Tu6ZS3/v1@Fql1uJ7jSGMBmKe2@QmaPPMmoNcs9y6kvkRWFZbrnEYooOcsvC4HB183iY3ZlOlLwOkDHCOYFdgoLYqG4TGMQHy65LCVK912XViITRyHozlIQ6pcTIt6KE/v/zyX@VGNIAulJteAA "JavaScript (Node.js) – Try It Online")
[Answer]
# [K (ngn/k)](https://codeberg.org/ngn/k), ~~44~~ 42 bytes
```
{+/-[q*k;x@s]+x@q:_y%k:1+!s:_%y}{-2!x*1+x}
```
-2 thanks to @coltim
Same formula as [@ovs' answer](https://codegolf.stackexchange.com/a/265116).
[Try it online!](https://ngn.codeberg.page/k#eJxLs6rW1teNLtTKtq5wKI7VrnAotIqvVM22MtRWLLaKV62srdY1UqzQMtSuqOXiSlM3NFAwNABiI2MFMyNTEBssYIBEAwEXAB0MFQw=)
[Answer]
# [C#](http://www.mono-project.com/docs/about-mono/releases/5.0.0/#csc), ~~82~~ 80 bytes
-2 thanks to [ceilingcat](https://codegolf.stackexchange.com/users/52904/ceilingcat)
With 2 dummy parameters (`y` and `z`).
```
(x,y,z)=>Enumerable.Range(1,y=(int)Math.Sqrt(x)).Sum(i=>2*i*(z=x/i)-z*~z+y*~y)/2
```
[Try it online!](https://tio.run/##hY2xDoIwFEV3vqLje1hAanTBshiddJHBGbFoEykRaEI7@Ou1TkxKcnOXe3JP1UdVL91Bq2or1UDJVLdWX58iJzXhDkZqqEWe75VuRFf6IT6X6i4gpYaD5/FUDo@4eHUDjIhxoRuQPGehDMHyMZG4sGH0tpHxbTBhLguCXav61j9dOjmIo1QCakiXlHyDmP3Y5wC2@g9s2HpWMTncBw "C# (Visual C# Interactive Compiler) – Try It Online")
# [C#](http://www.mono-project.com/docs/about-mono/releases/5.0.0/#csc), ~~89~~ 87 bytes
Single parameter.
```
x=>{int i=0,s=0,y=(int)Math.Sqrt(x),z;for(;i++<y;)s+=i*(z=x/i)+(y*~y-z*~z)/2;return s;}
```
[Try it online!](https://tio.run/##fYy9DsIgFIV3n4IRSv@s0eWWLiZOOjk4a6XxJg1EoEnB6KsjPkAdzknO@ZKvt0VvMR4m1beoXE7uerqNsiMDEXEW3SudBEWd2xQvaJrsdHWP8vw0js4sDzBoQwE5bz0wywVmNIi5QsZpyIpPKHxqz6oGjHSTUcTCO8JqtdfK6lGWF4NOHlFJOtB1zRgskGXUbJbQrtn@Ef6M8Qs "C# (Visual C# Interactive Compiler) – Try It Online")
Both use the same formula:
$$
\begin{eqnarray\*}
D(x) &=& \Big(\sum\_{i=1}^{\lfloor\sqrt x\rfloor}i \cdot \Big\lfloor\frac x i\Big\rfloor + T(\Big\lfloor\frac x i\Big\rfloor) \ \Big) - \lfloor \sqrt x\rfloor \cdot T\big(\lfloor \sqrt x\rfloor \big)\\
&=& \Big(\sum\_{i=1}^{y}i \cdot z + T(z) \ \Big) - y \cdot T(y)\\
&=& \sum\_{i=1}^{y}i \cdot z + T(z) - T(y)\ \\
&=& \sum\_{i=1}^{y}i \cdot z + \frac{z \cdot (z-1) - y \cdot (y-1) }{2} \\
\end{eqnarray\*}
$$
where
$$
T(x) = \frac{x(x+1)}{2} \\
y = \lfloor\sqrt x\rfloor \\
z = \Big\lfloor\frac x i\Big\rfloor
$$
]