Mfcttrf

For a research project I am working on currently, I need to do a very simple and straightforward calculation. Unfortunately, I do not know how to include Mathematica code here, but it is very short anyway:

b[x_] := x^2 - x + 1/2

bp[x_] := b[Mod[x, 1]]

d[n_, q_] := Sum[bp[k/n]* bp[q*k/n], k, 0, n - 1]

Now I need to compare two values of d[n,q], in particular, I need to calculate
d[1346269,514229] and d[1346269,1137064] to see which one is larger. It works perfectly fine for smaller numbers, e.g. I tried d[75025, 28657] and got the correct result in a reasonable amount of time. However, when I tried evaluating d[1346269,514229] after some time I got the result

(1/9760128332100732436)(4744910246749618660646829 - 
 4880064166050366218 Hold[$ConditionHold[$ConditionHold[
 System`Dump`AutoLoad[Hold[Sum`InfiniteSum], 
 Hold[Sum`InfiniteSum, Sum`SumInfiniteRationalSeries, 
 Sum`SumInfiniteRationalExponentialSeries, 
 Sum`SumInfiniteLogarithmicSeries, 
 Sum`SumInfiniteBernoulliSeries, 
 Sum`SumInfiniteFibonacciSeries, Sum`SumInfiniteLucasLSeries, 
 Sum`SumInfiniteArcTangentSeries, 
 Sum`SumInfiniteArcCotangentSeries, 
 Sum`SumInfiniteqArcTangentSeries, 
 Sum`SumInfiniteqArcCotangentSeries, 
 Sum`SumInfiniteStirlingNumberSeries, 
 Sum`SumInfiniteqRationalSeries], "Sum`InfiniteSum`"]]]][(1 - 
 Ceiling[(1 - Sum`FiniteSumDump`l)/1346269] + 
 Floor[(1346268 - Sum`FiniteSumDump`l)/1346269]) Mod[(
 514229 Sum`FiniteSumDump`l)/1346269, 1], Sum`FiniteSumDump`l, 
 0, 1346268, True])

Now, I am not too familiar with Mathematica, so I am not sure where the problem is exactly. However, I would need the two results of d[1346269,1137064] and d[1346269,514229] exactly (i.e. not numerically) as they are super close together, so any rounding could already change the results sufficiently much to alter the order of the two. Is there any way of computing those sums symbolically?

I find that this evaluates much faster. Mathematica is much faster at evaluating functions with a list of 1 million data points than it is at evaluate a function 1 million times with a single point each.

b[x] := x^2 - x + 1/2
bp[x_] := b[Mod[x, 1]]
d[n_, q_] := Total[bp[Range[0, n - 1]/n] bp[q Range[0, n - 1]/n]]
d[1346269, 1137064]

This takes about 17 seconds to evaluate on my machine and gives me:

$frac14599731344021534525766739760128332100732436$

Evaluating d[1346269, 514229] gives me:

$frac14599731343994714468596179760128332100732436$

The difference is:

$frac6705014292642440032083025183109$

Evaluated to 100 digits with N[difference, 100], I get

2.74792054550653169316510562953899839764312057350982296551104800194 089580643554604005134744886010933307817655482572397565420836194
$cdot 10^-7$