Wikipedia:Reference desk/Archives/Mathematics/2010 December 8

Welcome to the Wikipedia Mathematics Reference Desk Archives
The page you are currently viewing is an archive page. While you can leave answers for any questions shown below, please ask new questions on one of the current reference desk pages.

December 8

Approximating the relative value of three positive real numbers x:y:z with positive integers a:b:c

For any three positive real numbers, $x$, $y$, and $z$, I would like to find a sequence of triples $(a_{0},b_{0},c_{0}),(a_{1},b_{1},c_{1}),\ldots$ (where $a_{i}$, $b_{i}$, and $c_{i}$ are positive integers for all $i$) such that each triple approximates the relative values of the triple $(x,y,z)$ "better" than the triple before it. I am flexible on what precisely we mean by "better." Or, in other words, if I draw a ray from the origin, as $(xt,yt,zt)$ for $t>0$, how can I quickly determine what integer lattice points the ray comes closest to?

If there were only two positive real numbers, $x$ and $y$, I could use the continued fraction expansion of $x/y$. I'd set $a_{i}$ equal to the numerator of the $i$th convergent and I'd set $b_{i}$ equal to the denominator of the $i$th convergent. This would give me a sequence $(a_{0},b_{0}),(a_{1},b_{1}),\ldots$ where the pairs give better and better approximations to the relative value of $x$ and $y$. But how do I do this for three or more positive real values?

This is not a work or homework question, just idle curiosity. Thanks -- —Quantling (talk | contribs) 20:45, 8 December 2010 (UTC)

The LLL algorithm is probably overkill for that. 67.117.130.143 (talk) 22:39, 8 December 2010 (UTC)

I wrote a FORTRAN program and included subroutine to do what you want. Look through the comments to tease out the logic:

Collapsed FORTRAN source code.

program RAY ! Finds 3D grid pts closest to a 3D ray, from origin up. implicit none real*8 I,J,K ! Ray vector components. real*8 X,Y,Z ! Current coords of test point. real*8 MIN_DIST ! MINimum DIST of test points to the ray, so far. print *,"Ray V1.0" print *," " 100 print *,"Enter i,j,k values:" read (*,*) I,J,K print *," " print ("(A,F21.9,A,F21.9,A,F21.9,A)") ,"(i,j,k) = (",I,"," ,J,"," ,K," )" print *," " MIN_DIST = 999999999.9999999 ! Arbitrary large number. if ((J .gt. I) .and. (J .gt. K)) goto 300 ! If J component is greatest, ! track along Y. if ((K .gt. I) .and. (K .gt. J)) goto 400 ! If K component is greatest, ! track along Z. ! Otherwise... 200 do X = 1.0,99999999.0 ! Track along X: Y = int( (X*J)/I ) Z = int( (X*K)/I ) call LnPtD3D (I,J,K,X,Y,Z, MIN_DIST ) Y = Y 1.0 call LnPtD3D (I,J,K,X,Y,Z, MIN_DIST ) Z = Z 1.0 call LnPtD3D (I,J,K,X,Y,Z, MIN_DIST ) Y = Y - 1.0 call LnPtD3D (I,J,K,X,Y,Z, MIN_DIST ) if (MIN_DIST .lt. 0.0000001) goto 800 ! Close enough, we're done. enddo 300 do Y = 1.0,99999999.0 ! Track along Y: X = int( (Y*I)/J ) Z = int( (Y*K)/J ) call LnPtD3D (I,J,K,X,Y,Z, MIN_DIST ) X = X 1.0 call LnPtD3D (I,J,K,X,Y,Z, MIN_DIST ) Z = Z 1.0 call LnPtD3D (I,J,K,X,Y,Z, MIN_DIST ) X = X - 1.0 call LnPtD3D (I,J,K,X,Y,Z, MIN_DIST ) if (MIN_DIST .lt. 0.0000001) goto 800 ! Close enough, we're done. enddo 400 do Z = 1.0,99999999.0 ! Track along Z: X = int( (Z*I)/K ) Y = int( (Z*J)/K ) call LnPtD3D (I,J,K,X,Y,Z, MIN_DIST ) X = X 1.0 call LnPtD3D (I,J,K,X,Y,Z, MIN_DIST ) Y = Y 1.0 call LnPtD3D (I,J,K,X,Y,Z, MIN_DIST ) X = X - 1.0 call LnPtD3D (I,J,K,X,Y,Z, MIN_DIST ) if (MIN_DIST .lt. 0.0000001) goto 800 ! Close enough, we're done. enddo 800 return end ****************************************************************************** subroutine LnPtD3D (I,J,K,X,Y,Z, MIN_DIST ) ! Find LINE to ! POINT DISTANCE ! in 3D. implicit none * Inputs: real*8 I,J,K ! Ray vector components. real*8 X,Y,Z ! Current coords of test point. * Input/output: real*8 MIN_DIST ! MINimum DIST of test point to the ray, so far. * Local variables: real*8 DIST ! DISTance of test point to the ray. real*8 FAR ! how FAR test point is from the origin. real*8 PN(3) ! Projection Normal onto ray from (x,y,z) test pt. real*8 T ! T parameter along ray. integer*4 XI,YI,ZI ! X,Y,Z converted to integers. ! ! Find the parameter along the ray where the normal projection of test point ! (X,Y,Z) will fall: ! T = (I*X J*Y K*Z) / (I*I J*J K*K) ! ! Find the point coords at the T parameter: ! PN(1) = T*I PN(2) = T*J PN(3) = T*K ! ! Compute distance between the projected point and test point (X,Y,Z): ! DIST = sqrt( (X-PN(1))**2 (Y-PN(2))**2 (Z-PN(3))**2 ) if (DIST .lt. (MIN_DIST-0.00000001)) then ! We have a closer point: MIN_DIST = DIST FAR = sqrt(X*X Y*Y Z*Z) XI = int(X 0.1) ! Convert to integer for print. YI = int(Y 0.1) ! Convert to integer for print. ZI = int(Z 0.1) ! Convert to integer for print. print ("(A,I9,A,I9,A,I9,A,A,F9.8,A,F19.5)") ,"(X,Y,Z) = (" ,XI,"," ,YI,"," ,ZI,")" ,", Distance from ray = ",DIST ,", Distance from origin = ",FAR endif return end

Here's a run (trimmed a bit to fit on this screen):

 Ray V1.0

 Enter i,j,k values:
1.1
2.2
333

(i,j,k) = (          1.100000000,          2.200000000,        333.000000000 )

(X,Y,Z) = ( 0, 0,   1), Distance from ray = .00738621, Distance from origin =    1.00000
(X,Y,Z) = ( 1, 2, 302), Distance from ray = .00537179, Distance from origin =  302.00828
(X,Y,Z) = ( 1, 2, 303), Distance from ray = .00201442, Distance from origin =  303.00825
(X,Y,Z) = ( 3, 6, 908), Distance from ray = .00134295, Distance from origin =  908.02478
(X,Y,Z) = ( 4, 8,1211), Distance from ray = .00067147, Distance from origin = 1211.03303
(X,Y,Z) = (11,22,3330), Distance from ray = .00000000, Distance from origin = 3330.09084

Let me know if you have any questions or suggestions for ways to improve it. (One efficiency improvement might be to only do the square root required to find DIST once you've determined that you want to print that value, since SQRT is a rather CPU-intensive operation. Thus, you could get by with comparing distances squared, until then.) I also considered doing the calcs in spherical coords, since that would eliminate the need for X, Y, and Z loops in the main program, and replace them with a single R loop, but that would necessitate some conversions between Cartesian and spherical coords, so I decided against it. StuRat (talk) 07:14, 9 December 2010 (UTC)

Applied <syntaxhighlight>.—msh210℠ 07:26, 9 December 2010 (UTC)

Thanks, that makes it much more readable, although it does seem to change some text to white, which disappears against my white background, but highlighting the text fixes that prob. StuRat (talk) 07:34, 9 December 2010 (UTC)

Here's a stupid but simple way to do it: $(a_{i},b_{i},c_{i})=(\lfloor 10^{i}x\rfloor ,\lfloor 10^{i}y\rfloor ,\lfloor 10^{i}z\rfloor )$. —Bkell (talk) 08:15, 9 December 2010 (UTC)

It's been too long since I've done FORTRAN, so it is hurting me too much to try to figure out the underlying logic there. But, I am guessing that the approach presented there has a number of floating point operations that is more than linear in the number of output triples. I am hoping for something like the continued fraction approach that works for two numbers x:y; it has a count of floating point operations that is linear in the number of output pairs $(a_{i},b_{i})$. The decimal approach works in a sense but, for instance, it misses the perfect answers for 1:1/3:1/9. Any further ideas? Thanks! —Quantling (talk | contribs) 13:48, 9 December 2010 (UTC)

OK, I will list my program logic here:

1) Determine which is bigger, the i, j, or k component of the vector for the ray starting at the origin. Count up along the x, y , or z axis, accordingly. In my example, the k component was largest, so I incremented from z=1, up along the z-axis.

2) For each plane (z=1, z=2, z=3, in my example), find the x,y,z coords of the intersection of that plane with the ray. Since only one of those coords is guaranteed to be an integer (z, in our example), the other two coords must be rounded both up and down, to create a total of 4 test points. In our example, the true intersection with the z=908 plane would be at (2.9994,6.0533,908.0000). We therefore test the points (2,6,908) (2,7,908) (3,6,908) and (3,7,908), to see if any of them is closer to the ray than our previous closest point. If so, we print it out, and store that distance as the one we need to beat for the next "close point".

3) The determination of the distance of a test point from the ray is done by first finding the T parameter (length along the ray from the origin) of the normal projection of the test point onto the ray. This parameter is then used to find the (x,y,z) coords of that normal projection, and the distance from there to the test point can then be determined.

Here's a run with the (i,j,k) values of (1,1/3,1/9):

 Ray V1.0

 Enter i,j,k values:
1
.3333333333
.1111111111

(i,j,k) = (          1.000000000,          0.333333333,          0.111111111 )

(X,Y,Z) = (1,0,0), Distance from ray = .33149677, Distance from origin = 1.00000
(X,Y,Z) = (9,3,1), Distance from ray = .00000000, Distance from origin = 9.53939

Isn't that the desired answer ? StuRat (talk) 18:37, 9 December 2010 (UTC)

The obvious (to me) definition of closeness is the angle between the two vectors, so that you want a sequence of points lying within successively thinner cones. The first algorithm that comes to mind is to simply start at (1,1,1) and at each step move to the neighbor (among the 7 in the local first octant) that minimizes that angle. This may not be linear in the output, but it explores only an asymptotically one-dimensional volume. It may be possible to (that is, I think you still get the right answer if you) optimize it by at each step incrementing the one coordinate that will then be the smallest multiple of the corresponding goal coordinate. It may also be possible to project the cone into two (or all three) planes and search the resulting infinite wedges in parallel, discarding the occasional false hit that lies within the intersection of the projections but not within the cone proper. --Tardis (talk) 15:41, 9 December 2010 (UTC)

I made a tweak to the FORTRAN program and subroutine to find some values I may have skipped before. Here's the modified source code:

Collapsed FORTRAN source code (REVISED).

program RAY ! Finds 3D grid pts closest to a 3D ray, from origin up. implicit none real*8 I,J,K ! Ray vector components. real*8 X,Y,Z ! Current coords of test point. real*8 MIN_DIST ! MINimum DIST of test point to the ray, so far. print *,"Ray V2.0" print *," " 100 print *,"Enter i,j,k values:" read (*,*) I,J,K print *," " print ("(A,F21.9,A,F21.9,A,F21.9,A)") ,"(i,j,k) = (",I,"," ,J,"," ,K," )" print *," " MIN_DIST = 999999999.9999999 ! Arbitrary large number. if ((J .gt. I) .and. (J .gt. K)) goto 300 ! If J component is greatest, ! track along Y. if ((K .gt. I) .and. (K .gt. J)) goto 400 ! If K component is greatest, ! track along Z. ! Otherwise... 200 if (J .lt. K) then do X = 1.0,99999999.0 ! Track along X: Y = int( (X*J)/I ) Z = int( (X*K)/I ) call LnPtD3D (I,J,K,X,Y,Z, MIN_DIST ) Y = Y 1.0 call LnPtD3D (I,J,K,X,Y,Z, MIN_DIST ) Y = Y - 1.0 Z = Z 1.0 call LnPtD3D (I,J,K,X,Y,Z, MIN_DIST ) Y = Y 1.0 call LnPtD3D (I,J,K,X,Y,Z, MIN_DIST ) if (MIN_DIST .lt. 0.0000001) goto 800 ! Close enough, we're done. enddo goto 800 ! Far enough, we're done. else do X = 1.0,99999999.0 ! Track along X: Y = int( (X*J)/I ) Z = int( (X*K)/I ) call LnPtD3D (I,J,K,X,Y,Z, MIN_DIST ) Z = Z 1.0 call LnPtD3D (I,J,K,X,Y,Z, MIN_DIST ) Z = Z - 1.0 Y = Y 1.0 call LnPtD3D (I,J,K,X,Y,Z, MIN_DIST ) Z = Z 1.0 call LnPtD3D (I,J,K,X,Y,Z, MIN_DIST ) if (MIN_DIST .lt. 0.0000001) goto 800 ! Close enough, we're done. enddo goto 800 ! Far enough, we're done. endif 300 if (I .lt. K) then do Y = 1.0,99999999.0 ! Track along Y: X = int( (Y*I)/J ) Z = int( (Y*K)/J ) call LnPtD3D (I,J,K,X,Y,Z, MIN_DIST ) X = X 1.0 call LnPtD3D (I,J,K,X,Y,Z, MIN_DIST ) X = X - 1.0 Z = Z 1.0 call LnPtD3D (I,J,K,X,Y,Z, MIN_DIST ) X = X 1.0 call LnPtD3D (I,J,K,X,Y,Z, MIN_DIST ) if (MIN_DIST .lt. 0.0000001) goto 800 ! Close enough, we're done. enddo goto 800 ! Far enough, we're done. else do Y = 1.0,99999999.0 ! Track along Y: X = int( (Y*I)/J ) Z = int( (Y*K)/J ) call LnPtD3D (I,J,K,X,Y,Z, MIN_DIST ) Z = Z 1.0 call LnPtD3D (I,J,K,X,Y,Z, MIN_DIST ) Z = Z - 1.0 X = X 1.0 call LnPtD3D (I,J,K,X,Y,Z, MIN_DIST ) Z = Z 1.0 call LnPtD3D (I,J,K,X,Y,Z, MIN_DIST ) if (MIN_DIST .lt. 0.0000001) goto 800 ! Close enough, we're done. enddo goto 800 ! Far enough, we're done. endif 400 if (I .lt. J) then do Z = 1.0,99999999.0 ! Track along Z: X = int( (Z*I)/K ) Y = int( (Z*J)/K ) call LnPtD3D (I,J,K,X,Y,Z, MIN_DIST ) X = X 1.0 call LnPtD3D (I,J,K,X,Y,Z, MIN_DIST ) X = X - 1.0 Y = Y 1.0 call LnPtD3D (I,J,K,X,Y,Z, MIN_DIST ) X = X 1.0 call LnPtD3D (I,J,K,X,Y,Z, MIN_DIST ) if (MIN_DIST .lt. 0.0000001) goto 800 ! Close enough, we're done. enddo goto 800 ! Far enough, we're done. else do Z = 1.0,99999999.0 ! Track along Z: X = int( (Z*I)/K ) Y = int( (Z*J)/K ) call LnPtD3D (I,J,K,X,Y,Z, MIN_DIST ) Y = Y 1.0 call LnPtD3D (I,J,K,X,Y,Z, MIN_DIST ) Y = Y - 1.0 X = X 1.0 call LnPtD3D (I,J,K,X,Y,Z, MIN_DIST ) Y = Y 1.0 call LnPtD3D (I,J,K,X,Y,Z, MIN_DIST ) if (MIN_DIST .lt. 0.0000001) goto 800 ! Close enough, we're done. enddo goto 800 ! Far enough, we're done. endif 800 return end ****************************************************************************** subroutine LnPtD3D (I,J,K,X,Y,Z, MIN_DIST ) ! Find LINE to ! POINT DISTANCE ! in 3D. implicit none * Inputs: real*8 I,J,K ! Ray vector components. real*8 X,Y,Z ! Current coords of test point. * Input/output: real*8 MIN_DIST ! MINimum DIST of test points to the ray, so far. * Local variables: real*8 DIST ! DISTance of test point to the ray. real*8 FAR ! how FAR test point is from the origin. real*8 PN(3) ! Projection Normal onto ray from (x,y,z) test pt. real*8 T ! T parameter along ray. integer*4 XI,YI,ZI ! X,Y,Z converted to integers. ! ! Find the parameter along the ray where the normal projection of test point ! (X,Y,Z) will fall: ! T = (I*X J*Y K*Z) / (I*I J*J K*K) ! ! Find the point coords at the T parameter: ! PN(1) = T*I PN(2) = T*J PN(3) = T*K ! ! Compute distance between the projected point and test point (X,Y,Z): ! DIST = sqrt( (X-PN(1))**2 (Y-PN(2))**2 (Z-PN(3))**2 ) if (DIST .lt. (MIN_DIST-0.00000001)) then ! We have a closer point: MIN_DIST = DIST FAR = sqrt(X*X Y*Y Z*Z) XI = int(X 0.1) ! Convert to integer for print. YI = int(Y 0.1) ! Convert to integer for print. ZI = int(Z 0.1) ! Convert to integer for print. print ("(A,I9,A,I9,A,I9,A,A,F9.8,A,F19.5)") ,"(X,Y,Z) = (" ,XI,"," ,YI,"," ,ZI,")" ,", Distance from ray = ",DIST ,", Distance from origin = ",FAR endif return end

Here's a run of the original program, with a new test case that shows the difference:

 Ray V1.0

 Enter i,j,k values:

93
97
101

(i,j,k) = (         93.000000000,         97.000000000,        101.000000000 )

(X,Y,Z) = ( 0, 0,  1), Distance from ray = .79938580, Distance from origin =   1.00000
(X,Y,Z) = ( 1, 1,  1), Distance from ray = .05828506, Distance from origin =   1.73205
(X,Y,Z) = (23,24, 25), Distance from ray = .01457126, Distance from origin =  41.59327
(X,Y,Z) = (93,97,101), Distance from ray = .00000000, Distance from origin = 168.10413

And here's a run of the modified version, using the same new test case, including a value I missed before:

 Ray V2.0

 Enter i,j,k values:
93
97
101

(i,j,k) = (         93.000000000,         97.000000000,        101.000000000 )

(X,Y,Z) = ( 0, 0,  1), Distance from ray = .79938580, Distance from origin =   1.00000
(X,Y,Z) = ( 0, 1,  1), Distance from ray = .78274502, Distance from origin =   1.41421
(X,Y,Z) = ( 1, 1,  1), Distance from ray = .05828506, Distance from origin =   1.73205
(X,Y,Z) = (23,24, 25), Distance from ray = .01457126, Distance from origin =  41.59327
(X,Y,Z) = (93,97,101), Distance from ray = .00000000, Distance from origin = 168.10413

The problem with the original program was in step 2. While it was correctly finding the actual intersections with the ray in an increasing order of distance from the origin, the 4 nearby test points tried for each intersection weren't necessarily tried in order of increasing distance from the origin. Now they are. StuRat (talk) 20:17, 9 December 2010 (UTC)

StuRat — yes, that's a good answer in every respect, except that I was hoping for a faster run-time.

Tardis — I like the "smaller angle is better" formulation. And I like the idea of trying to simultaneously work in the projected wedges, to somehow bring the count of floating point operations down to linear in the number of output triples. How to accomplish that though? —Quantling (talk | contribs) 20:52, 9 December 2010 (UTC)

How fast does it need to be ? I do have a check that stops further searching once an exact intersection (or one within a set distance from an exact intersection) is found, and that limits run-time considerably. My first example took between 40 and 50 milliseconds on a mid-level PC, but, as is typical for this type of program, the print statements are most of this time. If I turn off the prints, that drops to 10 milliseconds of run-time. As I said before, I could further reduce this time by only doing square roots when absolutely necessary. Would you like me to make this change, to see how much faster that will make it ? StuRat (talk) 01:17, 10 December 2010 (UTC)

By the clock, your solution is excellent; I have no practical reason to need anything faster. (Actually, I have no practical reason to need any answer at all; this is just idle curiosity on my part.) The only objection I have to your fine solution is that, as is possible with the simpler $x:y$ case, I am hopeful for an algorithm that has $O(N)$ floating point operations where $N$ is the size of the output.

I think I figured it out, actually. Let me play with it a bit more. —Quantling (talk | contribs) 13:46, 10 December 2010 (UTC)

OK, please let us know what you come up with. StuRat (talk) 22:10, 10 December 2010 (UTC)