Home > Algorithms
Another Man's Floor (Posted on 2005-11-23) |
|
Let a, b, and m be positive whole numbers. Required is a fast algorithm for evaluating Σk=0..m floor(ka/b), floor(x) being the greatest integer that does not exceed the real number x.
|
Submitted by Richard
|
Rating: 4.5000 (2 votes)
|
|
Solution:
|
(Hide)
|
The title is intended to be a hint that refers to the Paul Simon song "One Man's Ceiling Is Another Man's Floor." Here the ceiling being referred to is the function ceiling, which is similar to function floor, except that non-integers are rounded up instead of down. To obtain a formula upon which a fast algorithm may be based, one may exploit the fact that 1 + floor(x) - ceiling(x) = 1 or 0 according as x is an integer or not in order to express the extra contribution that the "diagonal" makes to the lattice point count of an embedding rectangle when the lattice points of the two included triangles are added together. See "Dotty Right Triangle" for a similar analysis, although there the lattice points on the axes were also being counted, and the vertices were all lattice points. A suitable formula that can be obtained in this way is
Sum_{0<=k<=KB} floor(kA/B) + Sum_{0<=l<=KA}ceiling(lB/A)=floor(KA)floor(KB+1)
where K is any positive real number. This formula, with K fixed at K=m/b, and initially with A=a and B=b, but with A and B systematically reduced in stages as the algorithm progresses, makes possible the fast evaluation of Sum_{0<=k<=m} floor(ka/b). This is because, by extracting explicitly visible integer components from arguments of floors or ceilings for whichever of A/B or B/A is an improper fraction, it allows the computations to be organized into elementary stages that are based on, and synchronized with, the stages of Euclid's algorithm. A bc function that implements an algorithm based on this idea is as follows:
define floorsum(a,b,m) {
auto aa,bb,sum,kb,q,r,qp,rp,kr,krp
scale = 0
aa = a; bb = b; sum = 0; kb = m
while(1) {
q = aa / bb; r = aa % bb
sum = sum + (q * kb * (kb + 1)) / 2
if(r == 0) return(sum)
kr=(m * r) / b
sum = sum + kr * (kb + 1)
qp = bb / r; rp = bb % r
sum = sum - (qp * kr * (kr + 1)) / 2
if(rp == 0) return(sum)
krp =(m * rp) / b
sum = sum - kr * (krp + 1)
aa = r; bb = rp; kb = krp}
}
|
Comments: (
You must be logged in to post comments.)
|
|
Please log in:
Forums (0)
Newest Problems
Random Problem
FAQ |
About This Site
Site Statistics
New Comments (3)
Unsolved Problems
Top Rated Problems
This month's top
Most Commented On
Chatterbox:
|