# How to partition $[a,b]$ into $m$ equal sub-intervals and count the number of sub-intervals that intersect with a subset of $[a,b]$?

Suppose we have a subset of $$[a,b]$$, where $$a=0$$ and $$b=1$$, such as

$$A_1=\left\{\frac{1}{2^x}+\frac{1}{2^y}+\frac{1}{2^z}:x,y,z\in\mathbb{Z}\right\}\cap[0,1]$$

Which is a subset of $$[0,1]$$

If we partition $$[a,b]$$ into $$m$$ sub-intervals of equal length, where in this case $$a=0$$ and $$b=1$$, how do we count the number of sub-intervals which interesect with $$A_1$$?

How do we generalize this for any $$A_1$$? Such that the function inside the brackets in the definition of the set is arbitrary?

Take for example:

$$A_1=\left\{\frac{x^2+xy+y^2}{xy}+\frac{\sqrt{2}}{2}:x,y\in\mathbb{Z}\right\}\cap[0,1]$$

(Note the number sub-intervals that intersect with $$A_1$$ should be less than $$m$$)

# Attempt

generateA[c_Integer] :=
Select[Union@
Flatten[Table[
1/2^x + 1/2^y + 1/2^z, {x, 1, c}, {y, x, c}, {z, y, c}]],
0 <= # <= 1 &]


This generates $$A_1$$.

P[m_] := Interval /@ Partition[Subdivide[m], 2, 1]


This takes the $$m$$ sub-intervals, of equal length, of $$[a,b]$$

Total[Table[
Sign[Total[
Boole[Table[
IntervalMemberQ[P[m][[s]], generateA[100][[g]]], {g, 1,
Length[generateA[5]]}]]]], {s, 1, m}]]


This shows whether each element, for every c-value of $$A_1$$, belongs to one of the $$m$$ sub-intervals. If one element is in the $$m$$ sub-intervals, it is counted using boole; and, if multiple elements are counted, we add them up and using Sign to convert them to 1. If no elements are in the $$m$$ sub-intervals, the boole of all elements are zero and the sum of these booles is zero, hence the Sign is zero. We add the Signs of all sub-intervals to get the total number of sub-intervals intersecting with $$A_1$$.

The problem is as $$m \to \infty$$, it takes too much time to calculate all the sub-intervals that intersect with $$A_1$$. Is there a better way of doing this?

BinCounts is a bit of overkill but still efficient. Beware the bins in the partition are disjoint. I assumed the number 1 was meant to be in the last bin of the partition.

Block[{m = 1000},
With[{seq = 2^Range[-Ceiling@Log2@N@m - 2., -1.]},
Total[#] + 1 - Last[#] &@Unitize@BinCounts[
Flatten@Outer[Plus, seq, seq, seq],
{0., 1., 1./m}]
]] // RepeatedTiming

(*  {0.00024, 170}  *)


Floating-point is faster than exact integer/rational arithmetic, but you have to worry about round-off error. Once m gets to 10^15 or 16^16, we need to switch to arbitrary-precision:

Block[{m = 10^16},
With[{prec =
If[# > MachinePrecision, #, MachinePrecision] &@Log10[10. m]},
With[{seq = 2^Range[-Ceiling@Log2@N@m - SetPrecision[2, prec], -1]},
Total[#] + 1 - Last[#] &@Unitize@BinCounts[
Flatten@Outer[Plus, seq, seq, seq],
{0, SetPrecision[1, prec], N[1/m, prec]}]
]]] // RepeatedTiming

(*  {0.29, 25024}  *)


Here's an exact rational-number computation:

Block[{m = 10^16},
With[{k = Ceiling@Log2@N@m + 2},
Total[#] + 1 - Last[#] &@Unitize@BinCounts[
Total[Tuples[2^Range[-k, -1], 3], {2}],
{0, 1, 1/m}]
]] // RepeatedTiming

(*  {0.646, 25024}  *)


If the end-points are supposed to belong to each bin, then this would need adjustment (and the bins would not, strictly speaking, form a partition).

P.S. I tried it with Tuples instead of Outer, but it was slower, especially for large m.

With[{k = Ceiling@Log2@N@m + 2.},
Total[#] + 1 - Last[#] &@Unitize@BinCounts[
Total[Tuples[2^Range[-k, -1.], 3], {2}],
{0., 1., 1./m}]
]


P.P.S. Another way, with Floor:

Block[{m = 10^16},
With[{prec =
If[# > MachinePrecision, #, MachinePrecision] &@Log10[1000. m]},
With[{seq = 2^Range[-Ceiling@Log2@N@m - SetPrecision[2, prec], -1]},
Total[
UnitStep[m - #] UnitStep[#] &@DeleteDuplicates@Floor[
m*Flatten@Outer[Plus, seq, seq, seq]
]
]
]]] // RepeatedTiming

(*  {0.269, 25024}  *)

• In one of my questions I stated how we can generalize this to any set? What if we had $\left\{\frac{1}{c}:c\in\mathbb{Z}\right\}\cap[0,1]$? I wasn't able to work this out. Jul 10, 2020 at 10:42
• For example, what if the function inside the bracket in the definition of $A_1$ is arbitrary? Jul 10, 2020 at 11:38
• Suppose you have a function f[x1_, x2_,..., xn_] of n arguments. Then you have to replace seq by seq1, seq2,..., seqn, where the sequences are of the form Range[a1, b1], Range[a2, b2],..., Range[an, bn] -- and hopefully you can work out the limits a and b for each range (they will depend on f and m). Then replace Flatten@Outer[Plus,...] with Flatten@Outer[f, seq1, seq2,..., seqk]. For $1/c$, seq1 = Range[1, m+1]. For a single-variable f[x], one can replace Flatten@Outer[.. by f /@ seq1. Jul 10, 2020 at 12:07
• Is it possible to write this out in your answer. I will accept it with +15 points. Jul 10, 2020 at 12:59
• ...... For example, what would happen if we took $A_1=\left\{\frac{x^2+xy+y^2}{x+y}+\frac{\sqrt{2}}{2}:x,y\in\mathbb{Z}\right\}\cap[0,1]$.....Apologies I want to fully understand your code. Jul 10, 2020 at 13:16