PE 372 【类欧几里德】【Stern-Brocot Tree】

题目大意：求 $\lfloor\frac {y^2}{x^2}\rfloor$ 为奇数的个数，其中 $2*10^6<x,y <= 10^9$

设 $k = \lfloor\frac {y^2}{x^2}\rfloor$,有 $x^2k <= y^2 <x^2(k+1)$
得 $x\sqrt k <= y < x\sqrt {k+1}$
k为奇数，所以是计算在 $\sqrt kx <= y$ 区域内的整点数，奇数加，偶数减

对于直线 $y = kx$,k 为无理数，可以找到一个大于等于 k 的最小分数，答案不变，这步用 Stern-Brocot Tree 就可以实现

现在就是求 $y >= \frac {ax}b$
考虑用类欧几里德

$\sum_{i=0}^N \lfloor \frac {ax + b}c\rfloor = \lfloor\frac ac\rfloor·\frac{N(N+1)}2 + \lfloor\frac bc\rfloor·(N+1) + \sum_{i = 0}^N \lfloor\frac {(a\%c)x + b\%c}{c}\rfloor$

现在就有 $a < c,b < c$

设 $M = \lfloor \frac {aN + b}c\rfloor$
原式 = $\sum_{x = 0}^N\sum_{y=0}^M [y <\lfloor \frac {ax + b}c\rfloor]$

现在在简化这个式子
$c(y + 1) <= ax + b$
$cy + c - b <= ax$
$cy + c - b - 1 < ax$
$x >\lfloor \frac {cy + c - b - 1}a\rfloor$

现在就变成了求
$\sum_{y=0}^M\sum_{x=0}^N[x >\lfloor \frac {cy + c - b - 1}a\rfloor]$
= $\sum_{y=0}^{M-1}(N -\lfloor \frac {cy + c - b - 1}a\rfloor)$
= $N * M - \sum_{y=0}^{M-1}\lfloor \frac {cy + c - b - 1}a\rfloor)$
//这里要求的是 $x <= N$,但 y 取到 M 的时候会超出范围，而考虑到 y = M 时对答案无贡献，所以只用加到 M - 1

现在就可以在 $O(log)$ 的复杂度内计算出答案了

友情提示：这道题要用 long double

【答案】301450082318807027

#include<iostream>
#include<algorithm>
#include<cstdlib>
#include<cstdio>
#include<cstring>
#include<string>
#include<ctime>
#include<cmath>
#define D long double
#define LL long long
using namespace std;

void find(LL la,LL lb,LL ra,LL rb,LL &p,LL &q,LL lim,D aim)
{
    LL a = la + ra,b = lb + rb;
    if (a > lim || b > lim) return;
    if ((D)a / b < aim)
        find(a,b,ra,rb,p,q,lim,aim);
    else
        p = a,q = b,find(la,lb,a,b,p,q,lim,aim);
}

LL F(LL a,LL b,LL c,LL N)
{
    LL ret = N * (N + 1) / 2 * (a / c) + (N + 1) * (b / c);
    a %= c,b %= c;
    LL M = (a * N + b) / c;
    if (!M) return ret;
    LL t = a ? N * M - F(c,c - b - 1,a,M - 1) : 0;
    return ret + t;
}

LL cal(LL N,LL k)
{
    LL t = sqrt(k),a = 0,b = 0,c = 1;
    if (t * t != k) find(0,1,1,0,c,a,N,sqrt((D)k));
        else a = 1,c = t;
    return F(a,b,c,N);
}

int main()
{
    LL N = 1e9,M = 2*1e6,lim = N * N / (M + 1) / (M + 1),ans = 0;
    for (LL k = 1;k <= lim;k ++)
    {
        LL l = (LL)M * (sqrt(k)),t = cal(N,k) - cal(l,k) - (N - l) * M;
        ans += (k & 1) ? t : -t;
    }
    cout << ans << endl;

    return 0;
}