[math-fun] bicycle gear ratios & continued fractions
Perhaps someone has already done this. The sensor on my bike outputs both cadence (crank speed in revs/sec) and wheel speed (revs/sec). I'm trying to infer which gear the bike is in by looking at the ratio of these two numbers. My plan was to use the equivalent of Scheme's "rationalize" function (which uses a continued fraction algorithm internally), which produces an exact fraction which is within some epsilon of a given number -- in this case the wheel/crank speed. This should give me an approximation to the ratio (#front teeth)/(#rear teeth). http://www.r6rs.org/final/html/r6rs/r6rs-Z-H-14.html Unfortunately, there is enough error in the measurements of the crank speed and the wheel speed so that the output isn't stable. So far, we haven't utilized any information about which chainrings (front teeth) and cogs (rear teeth) are mounted on the bike, except for the fact that each must have an integer # of teeth. Of course, I can always make up a table of all of the gears and simply find the closest one & output it. But if the number of cogs grows large and the number of chainrings grows large (we're mathematicians, after all!), we'd rather not search a large table. Suppose now that we have upper & lower limits on the teeth on the chainrings (28 to 52) and the cogs (12 to 26), but otherwise don't know how many chainrings or cogs are mounted on the bike. How do we modify the continued fraction algorithm used within rationalize to produce fractions with numerators and denominators in the appropriate range ?
HGBaker> Perhaps someone has already done this.
The sensor on my bike outputs both cadence (crank speed in revs/sec) and wheel speed (revs/sec). I'm trying to infer which gear the bike is in by looking at the ratio of these two numbers.
My plan was to use the equivalent of Scheme's "rationalize" function (which uses a continued fraction algorithm internally), which produces an exact fraction which is within some epsilon of a given number -- in this case the wheel/crank speed. This should give me an approximation to the ratio (#front teeth)/(#rear teeth).
http://www.r6rs.org/final/html/r6rs/r6rs-Z-H-14.html
Unfortunately, there is enough error in the measurements of the crank speed and the wheel speed so that the output isn't stable.
So far, we haven't utilized any information about which chainrings (front teeth) and cogs (rear teeth) are mounted on the bike, except for the fact that each must have an integer # of teeth.
Of course, I can always make up a table of all of the gears and simply find the closest one & output it. But if the number of cogs grows large and the number of chainrings grows large (we're mathematicians, after all!), we'd rather not search a large table.
Suppose now that we have upper & lower limits on the teeth on the chainrings (28 to 52) and the cogs (12 to 26), but otherwise don't know how many chainrings or cogs are mounted on the bike.
How do we modify the continued fraction algorithm used within rationalize to produce fractions with numerators and denominators in the appropriate range ?
I sent a long msg on this, 3 MAR 05. Upshot: CF may not be best. ------------(quote from msg)-------- A trickier problem, for which the mediant algorithm seems better suited: What is the nearest fraction to x with numerator <= a and denominator <= b? E.g., (farint pi 999) 355 113 (farint pi 999 112) 333 106 But notice I said "for hand computation" For some reason (personal stupidity?), the iterated mediant is ugly to code. The following was maybe my sixth try: (defun farint (x &optional (nhi 400) (dhi nhi) (ln 0) (ld 1) (hn 1) (hd 0)) (if (> (+ ln hn) (* (+ ld hd) x)) (let* ((m (min (if (= 0 ln) nhi (floor(- nhi hn) ln)) (floor (- dhi hd) ld))) (d (- (* x ld) ln)) (k (if (= 0 d) m (ceiling (- hn (* x hd)) d)))) (if (< k m) (farint x nhi dhi (setf hn (+ (* k ln) hn)) (setf hd (+ (* k ld) hd)) (- hn ln) (- hd ld)) (let* ((n (+ (* m ln) hn)) (d (+ (* m ld) hd))) (if (< (* 2 d ld x) (+ (* ld n) (* ln d))) (values ln ld) (values n d))))) (let* ((m (min (floor (- nhi ln) hn) (if (= 0 hd) dhi (floor (- dhi ld) hd)))) (d (- hn (* x hd))) (k (if (= 0 d) m (ceiling (- (* x ld) ln) d)))) (if (< k m) (farint x nhi dhi (- (setf ln (+ (* k hn) ln)) hn) (- (setf ld (+ (* k hd) ld)) hd) ln ld) (let* ((n (+ (* m hn) ln)) (d (+ (* m hd) ld))) (if (< (* 2 d hd x) (+ (* hd n) (* hn d))) (values n d) (values hn hd))))))) Presumably the two outer branches can be obnubilated into something half as long and twice as turbid. --rwg -----------(End quote from msg)------- The text file from which I manuallly recovered this has been damaged (I suspect by an antivirus program) by the deletion of the first character of every line! For the hand algorithm, click on Rectangle Arithmetic at http://www.tweedledum.com/rwg/ --rwg
HGBaker> Perhaps someone has already done this.
The sensor on my bike outputs both cadence (crank speed in revs/sec) and wheel speed (revs/sec). I'm trying to infer which gear the bike is in by looking at the ratio of these two numbers.
My plan was to use the equivalent of Scheme's "rationalize" function (which uses a continued fraction algorithm internally), which produces an exact fraction which is within some epsilon of a given number -- in this case the wheel/crank speed. This should give me an approximation to the ratio (#front teeth)/(#rear teeth).
http://www.r6rs.org/final/html/r6rs/r6rs-Z-H-14.html
Unfortunately, there is enough error in the measurements of the crank speed and the wheel speed so that the output isn't stable.
So far, we haven't utilized any information about which chainrings (front teeth) and cogs (rear teeth) are mounted on the bike, except for the fact that each must have an integer # of teeth.
Of course, I can always make up a table of all of the gears and simply find the closest one & output it. But if the number of cogs grows large and the number of chainrings grows large (we're mathematicians, after all!), we'd rather not search a large table.
Suppose now that we have upper & lower limits on the teeth on the chainrings (28 to 52) and the cogs (12 to 26), but otherwise don't know how many chainrings or cogs are mounted on the bike.
How do we modify the continued fraction algorithm used within rationalize to produce fractions with numerators and denominators in the appropriate range ?
I sent a long msg on this, 3 MAR 05. Upshot: CF may not be best.
OTOH, (defun bstrat (x nhi &optional (c 0) (d 1)) (multiple-value-bind (f y) (floor x) (cond ((or (> f (* 2 nhi)) (and (= f (* 2 nhi)) (>= (* c y) d))) #.(SI::make-canonical-rational 1 0)) ((>= f nhi) nhi) ((= 0 (setf r (round (* nhi y) x))) f) ((+ f (/ (bstrat (/ y) (floor (- nhi r) f) (+ d (* f c)) c))))))) is a CF method and a lot shorter. But no easier to write! One hard part was to think of getting rid of the denominator bound prior to the recursion. If you want both call instead the superroutine (defun bestrat (x n &optional (d n)) (cond ((= 0 x) 0) ((< n (* d x)) (bstrat x n)) ((/ (bstrat (/ x) d))))) Another hard part is how to back up when the numerator oversteps. E.g., if x = A B C D E F G ..., and P<E is the largest acceptable term in A B C D P, then A B C D P is actually worse than A B C D unless P D C B < E-P F G ... . The *only* purpose of the c and d variables in the recursion is to compute D C B (i.e. the value of the continued fraction up to termination with its terms reversed and first term deleted!) for the borderline case P = E-P, i.e. breaking off half of an even term. E.g., pi = 3 7 15 1 292 1 1 1 2 1 3 1 14 2 1 1 2 2 2 2 ..., and 3 7 15 1 = (the famous) 355/113. The next better approximation is 3 7 15 1 146 = 52163/16604, because 146 1 15 7 > 146 1 1 1 2 1 3 1 14 ... Notice how marginal is the improvement: (c184) APPLY("/",ARGS(CFEXPAND([3,7,15,1,146]))); 52163 355 (d184) [-----, ---] 16604 113 (c185) DFLOAT(%-%PI); (d185) [- 2.66213257216208d-7, 2.66764189404967d-7] Testing the program, (bstrat pi 52163) 1 Enter BSTRAT 3.14159265358979d0 52163 | 2 Enter BSTRAT 7.06251330593105d0 16604 1 0 | 3 Enter BSTRAT 15.9965944066841d0 2351 7 1 | | 4 Enter BSTRAT 1.003417231015d0 147 106 7 | | 5 Enter BSTRAT 292.634590875012d0 146 113 106 | | 5 Exit BSTRAT 146 | | 4 Exit BSTRAT 147/146 | 3 Exit BSTRAT 2351/147 | 2 Exit BSTRAT 16604/2351 1 Exit BSTRAT 52163/16604 (bstrat pi 52162) 1 Enter BSTRAT 3.14159265358979d0 52162 | 2 Enter BSTRAT 7.06251330593105d0 16603 1 0 | 3 Enter BSTRAT 15.9965944066841d0 2350 7 1 | | 4 Enter BSTRAT 1.003417231015d0 146 106 7 | | 4 Exit BSTRAT 1 | 3 Exit BSTRAT 16 | 2 Exit BSTRAT 113/16 1 Exit BSTRAT 355/113 But the hardest part is to figure out how the reduce the numerator bound during the recursion, and the solution above just doesn't feel right. A sure fire, short, disgusting alternative is to replace (floor (- nhi (max 1 (round (* nhi y) x))) f) with (round nhi x), and then check the answer for the occasional numerator>nhi, and recurse with nhi-1! It's actually quite efficient. It can't treewalk exponentially because the errors are caught right away near the bottom, and can't propagate up the recursion. --rwg CORNERSTONE NONSECRETOR
------------(quote from msg)-------- A trickier problem, for which the mediant algorithm seems better suited: What is the nearest fraction to x with numerator <= a and denominator <= b? E.g., (farint pi 999) 355 113 (farint pi 999 112) 333 106
But notice I said "for hand computation" For some reason (personal stupidity?), the iterated mediant is ugly to code. The following was maybe my sixth try:
(defun farint (x &optional (nhi 400) (dhi nhi) (ln 0) (ld 1) (hn 1) (hd 0)) (if (> (+ ln hn) (* (+ ld hd) x)) (let* ((m (min (if (= 0 ln) nhi (floor(- nhi hn) ln)) (floor (- dhi hd) ld))) (d (- (* x ld) ln)) (k (if (= 0 d) m (ceiling (- hn (* x hd)) d)))) (if (< k m) (farint x nhi dhi (setf hn (+ (* k ln) hn)) (setf hd (+ (* k ld) hd)) (- hn ln) (- hd ld)) (let* ((n (+ (* m ln) hn)) (d (+ (* m ld) hd))) (if (< (* 2 d ld x) (+ (* ld n) (* ln d))) (values ln ld) (values n d))))) (let* ((m (min (floor (- nhi ln) hn) (if (= 0 hd) dhi (floor (- dhi ld) hd)))) (d (- hn (* x hd))) (k (if (= 0 d) m (ceiling (- (* x ld) ln) d)))) (if (< k m) (farint x nhi dhi (- (setf ln (+ (* k hn) ln)) hn) (- (setf ld (+ (* k hd) ld)) hd) ln ld) (let* ((n (+ (* m hn) ln)) (d (+ (* m hd) ld))) (if (< (* 2 d hd x) (+ (* hd n) (* hn d))) (values n d) (values hn hd)))))))
Presumably the two outer branches can be obnubilated into something half as long and twice as turbid. --rwg -----------(End quote from msg)------- The text file from which I manuallly recovered this has been damaged (I suspect by an antivirus program) by the deletion of the first character of every line!
For the hand algorithm, click on Rectangle Arithmetic at http://www.tweedledum.com/rwg/ --rwg
I haven't had a chance to play with either the mediant or the CF algorithms, but did play with searching the actual bike gear ratios where "closest" match was measured by the _log_ of the ratio instead of the ratio itself. Using logs doesn't affect things very much, but appears to be slightly closer to what is actually wanted. As you can imagine, if some of the gears have relatively close ratios, and if there is any noise and/or slippage, then there is a little instability. An extra credit problem: if you don't a priori know the wheel circumference, guess it from a long ride with lots of gear changes. Even more credit: if you know neither the wheel circumference nor the actual gear ratios, guess them from a long ride with lots & lots of gear changes. Now a real problem: from a very long & very noisy GPS track (lat,lon,altitude,time) + the cadence (but not the wheel speed), infer the circumference, chain ring teeth & cog teeth. One must decide when one is coasting and when the wheel is slipping, and somehow filter out the GPS noise to determine distance & hence speed. At 06:15 AM 4/21/2009, rwg@sdf.lonestar.org wrote:
HGBaker> Perhaps someone has already done this.
The sensor on my bike outputs both cadence (crank speed in revs/sec) and wheel speed (revs/sec). I'm trying to infer which gear the bike is in by looking at the ratio of these two numbers.
My plan was to use the equivalent of Scheme's "rationalize" function (which uses a continued fraction algorithm internally), which produces an exact fraction which is within some epsilon of a given number -- in this case the wheel/crank speed. This should give me an approximation to the ratio (#front teeth)/(#rear teeth).
http://www.r6rs.org/final/html/r6rs/r6rs-Z-H-14.html
Unfortunately, there is enough error in the measurements of the crank speed and the wheel speed so that the output isn't stable.
So far, we haven't utilized any information about which chainrings (front teeth) and cogs (rear teeth) are mounted on the bike, except for the fact that each must have an integer # of teeth.
Of course, I can always make up a table of all of the gears and simply find the closest one & output it. But if the number of cogs grows large and the number of chainrings grows large (we're mathematicians, after all!), we'd rather not search a large table.
Suppose now that we have upper & lower limits on the teeth on the chainrings (28 to 52) and the cogs (12 to 26), but otherwise don't know how many chainrings or cogs are mounted on the bike.
How do we modify the continued fraction algorithm used within rationalize to produce fractions with numerators and denominators in the appropriate range ?
I sent a long msg on this, 3 MAR 05. Upshot: CF may not be best.
OTOH, (defun bstrat (x nhi &optional (c 0) (d 1)) (multiple-value-bind (f y) (floor x) (cond ((or (> f (* 2 nhi)) (and (= f (* 2 nhi)) (>= (* c y) d))) #.(SI::make-canonical-rational 1 0)) ((>= f nhi) nhi) ((= 0 (setf r (round (* nhi y) x))) f) ((+ f (/ (bstrat (/ y) (floor (- nhi r) f) (+ d (* f c)) c))))))) is a CF method and a lot shorter. But no easier to write! One hard part was to think of getting rid of the denominator bound prior to the recursion. If you want both call instead the superroutine
(defun bestrat (x n &optional (d n)) (cond ((= 0 x) 0) ((< n (* d x)) (bstrat x n)) ((/ (bstrat (/ x) d)))))
Another hard part is how to back up when the numerator oversteps. E.g., if x = A B C D E F G ..., and P<E is the largest acceptable term in A B C D P, then A B C D P is actually worse than A B C D unless P D C B < E-P F G ... . The *only* purpose of the c and d variables in the recursion is to compute D C B (i.e. the value of the continued fraction up to termination with its terms reversed and first term deleted!) for the borderline case P = E-P, i.e. breaking off half of an even term. E.g., pi = 3 7 15 1 292 1 1 1 2 1 3 1 14 2 1 1 2 2 2 2 ..., and 3 7 15 1 = (the famous) 355/113. The next better approximation is 3 7 15 1 146 = 52163/16604, because 146 1 15 7 > 146 1 1 1 2 1 3 1 14 ... Notice how marginal is the improvement: (c184) APPLY("/",ARGS(CFEXPAND([3,7,15,1,146])));
52163 355 (d184) [-----, ---] 16604 113
(c185) DFLOAT(%-%PI);
(d185) [- 2.66213257216208d-7, 2.66764189404967d-7]
Testing the program, (bstrat pi 52163) 1 Enter BSTRAT 3.14159265358979d0 52163 | 2 Enter BSTRAT 7.06251330593105d0 16604 1 0 | 3 Enter BSTRAT 15.9965944066841d0 2351 7 1 | | 4 Enter BSTRAT 1.003417231015d0 147 106 7 | | 5 Enter BSTRAT 292.634590875012d0 146 113 106 | | 5 Exit BSTRAT 146 | | 4 Exit BSTRAT 147/146 | 3 Exit BSTRAT 2351/147 | 2 Exit BSTRAT 16604/2351 1 Exit BSTRAT 52163/16604
(bstrat pi 52162) 1 Enter BSTRAT 3.14159265358979d0 52162 | 2 Enter BSTRAT 7.06251330593105d0 16603 1 0 | 3 Enter BSTRAT 15.9965944066841d0 2350 7 1 | | 4 Enter BSTRAT 1.003417231015d0 146 106 7 | | 4 Exit BSTRAT 1 | 3 Exit BSTRAT 16 | 2 Exit BSTRAT 113/16 1 Exit BSTRAT 355/113
But the hardest part is to figure out how the reduce the numerator bound during the recursion, and the solution above just doesn't feel right. A sure fire, short, disgusting alternative is to replace (floor (- nhi (max 1 (round (* nhi y) x))) f) with (round nhi x), and then check the answer for the occasional numerator>nhi, and recurse with nhi-1! It's actually quite efficient. It can't treewalk exponentially because the errors are caught right away near the bottom, and can't propagate up the recursion. --rwg CORNERSTONE NONSECRETOR
------------(quote from msg)-------- A trickier problem, for which the mediant algorithm seems better suited: What is the nearest fraction to x with numerator <= a and denominator <= b? E.g., (farint pi 999) 355 113 (farint pi 999 112) 333 106
But notice I said "for hand computation" For some reason (personal stupidity?), the iterated mediant is ugly to code. The following was maybe my sixth try:
(defun farint (x &optional (nhi 400) (dhi nhi) (ln 0) (ld 1) (hn 1) (hd 0)) (if (> (+ ln hn) (* (+ ld hd) x)) (let* ((m (min (if (= 0 ln) nhi (floor(- nhi hn) ln)) (floor (- dhi hd) ld))) (d (- (* x ld) ln)) (k (if (= 0 d) m (ceiling (- hn (* x hd)) d)))) (if (< k m) (farint x nhi dhi (setf hn (+ (* k ln) hn)) (setf hd (+ (* k ld) hd)) (- hn ln) (- hd ld)) (let* ((n (+ (* m ln) hn)) (d (+ (* m ld) hd))) (if (< (* 2 d ld x) (+ (* ld n) (* ln d))) (values ln ld) (values n d))))) (let* ((m (min (floor (- nhi ln) hn) (if (= 0 hd) dhi (floor (- dhi ld) hd)))) (d (- hn (* x hd))) (k (if (= 0 d) m (ceiling (- (* x ld) ln) d)))) (if (< k m) (farint x nhi dhi (- (setf ln (+ (* k hn) ln)) hn) (- (setf ld (+ (* k hd) ld)) hd) ln ld) (let* ((n (+ (* m hn) ln)) (d (+ (* m hd) ld))) (if (< (* 2 d hd x) (+ (* hd n) (* hn d))) (values n d) (values hn hd)))))))
Presumably the two outer branches can be obnubilated into something half as long and twice as turbid. --rwg -----------(End quote from msg)------- The text file from which I manuallly recovered this has been damaged (I suspect by an antivirus program) by the deletion of the first character of every line!
For the hand algorithm, click on Rectangle Arithmetic at http://www.tweedledum.com/rwg/ --rwg
participants (2)
-
Henry Baker -
rwg＠sdf.lonestar.org