We took "east" to refer to the positive x-axis (equating the conventional directions of 2-D Cartesian axes with the principal directions in the most common cartographic convention).
The following solution was supplied by Marijke van Gans
C's floating point types don't look accurate enough (that's right, when you do it you'll see it magnify initial uncertainty by around 10^12 in the process). I reached for my copy of UBASIC, a number-theory oriented freeware BASIC, by Yuji Kida. It has integers of up to 1080 bytes long but you can transfer up to half of that accuracy to the right of the "decimal" point (i.e. fixed-point fractions).
As this was the first problem i tackled i didn't do a proper write-up at the time but put the details of the algoritmic decisions as a big initial comment in the program (reproduced below). Some more general strategy observations first:
I took two paths, one starting 100*epsilon too high (and hence bounced too far clockwise at the first hurdle), the other starting 100*epsilon too low (and hence bounced too far anticlockwise at the first hurdle). Take it from there, and see how far they diverge.
epsilon being the smallest number UBASIC can represent at the chosen level of precision (this is fixed-point reals) and 100*epsilon rather than epsilon so as not to accidentally lose it at the first bounce. I *think* this is a robust way of doing this. The result i posted at the time was calculated at "POINT 8" (some 38 digits after the decimal point) internally, and the digits i posted
> FWIW i make that 0.9952629194433541608903118 units from the origin.
were the ones that were the same (rounded, not truncated) for both paths. The program as reproduced here is a later polishing that lets a user choose the POINT size (byte pairs to the right of the point), 5 will *just* about do. Also, it now incorporates as an option the graphical version which i wrote when flummoxed that i couldn't get my photons to bounce properly, i wanted to *see* what i was doing wrong <g>.
I made the program print the divergence in x and y, as well as that in r, to avoid spuriously close r's for points far apart, and print out the t,x,y coordinates of bounce events to see the divergence marching up through the digits, i.e. avoiding spuriously close end results for wildly diverging paths intersecting near the end. Criticisms of method are invited.
--------------8<--------------------------------------------->8--------------
1 ' SIAM challenge 2002 problem 2
9 '
10 ' We'll use time steps of 0.25, each step check if the photon
11 ' should have bounced off a mirror. 0.25 is nice and exact in
12 ' binary, and smaller than the minimum mirror to mirror distance
13 ' 1/3, so no multiple collisions per step. Two things can happen.
14 '
15 ' a) The step brings us inside a mirror. Solve to find out at
16 ' what exact (non-discrete) t the bounce should have happened.
17 '
18 ' b) The step brings us outside a mirror again while the path
19 ' really intersected it. We can detect this by calculating the
20 ' distance of our line to the centre. If < 1/3, solve for when.
21 '
22 ' Either way, we must take the -sqrt(D) root always (the other
23 ' is for intersecting circle on the way out again).
24 '
25 ' But which integer grid point? Loop over all 20x20 or so
26 ' each time? No. For (a), the end point in the step can
27 ' only be inside the circle around its nearest gridpoint.
28 ' Just round x and y to integer. And for (b) we have the nice
29 ' situation that this can only happen if both endpoints of the
30 ' step have the *same* integer gridpoint as nearest (worst case
31 ' tangent, with one endpoint a whole 0.25 away from the point
32 ' of contact -- by Pythagoras, 4/12 radial and 3/12 tangential
33 ' means we're 5/12 away from the gridpoint, less than half).
34 ' Actually this means we can take (a) and (b) together. Just
35 ' look at the gridpoint nearest to the *end* of the step.
36 ' If (a), the line certainly went closer than 1/3 to *that*
37 ' gridpoint. And if (b) we only have to test our line wrt
38 ' one of the gridpoints anyway.
39 '
40 ' Calculating sin and angle from cos when cos is close to 1
41 ' is imprecise but there's no workaround, it's inherent in
42 ' the problem that we might closely graze some mirrors. So
43 ' we just have to do this to big precision, e.g. UBasic,
44 ' and then run a test with the photon starting a little bit
45 ' higher or lower. That way we can find out how sensitive
46 ' it all is to initial conditions.
49 '
50 print "Recommended screen size for text"
51 print "accuracy 2 to 6: 43-line"
52 print "accuracy 7 to 12: 50-line"
53 print "(Graphics sets own screen)"
54 print
60 print "Accuracy of fractional part to the right of the decimal point"
61 input "in 16-bit words (5 or more recommended)";W%
62 point W%
63 print
69 '
70 TWOPI=pi(2) ' 2*pi
71 MOVE=0.25
72 RADIUS=1//3 ' UBasic exact rational
73 INIX=0.5
74 INIY=0.1
75 '
80 print "In text mode we will run twice, once with the photon 100*eps"
81 print "too low down, once with it 100*eps too high up (where eps is"
82 print "the smallest representable number in this accuracy)."
83 print "In graphics mode we display the path of the bouncing photon."
84 I=input$(0) ' clear kybd buffer
85 print "Graphics display? [Y/N] ";
86 I=input$(1) ' wait for key, cursor
87 if I="y" goto 200
88 if I<>"n" goto 86
89 '
100 '''''''''''''''''''' PROGRAM FLOW for TEXT '''''''''''''''''''''''''
101 IsText=1
102 print
109 '
110 Butterfly=-100*#eps:gosub 400:gosub 900
111 GrOld=Gr
112 GxOld=Gx
113 GyOld=Gy
119 '
120 Butterfly=+100*#eps:gosub 400:gosub 900
121 print "Delta r =";Gr-GrOld
122 print "Delta x =";Gx-GxOld
123 print "Delta y =";Gy-GyOld
124 end
129 '
200 '''''''''''''''''''' PROGRAM FLOW for GRAPHICS '''''''''''''''''''''
201 IsText=0
202 Butterfly=0:gosub 300
203 J=input$(-1) ' wait for key, no cursor
204 cls:end
209 '
300 '''''''''''''''''''' MAIN LOOP entrypoint for GRAPHICS '''''''''''''
301 '
310 screen 23 ' 640x480 @ 16 colors
311 cls
312 SCALE=150
313 XORIG=320
314 YORIG=300
319 '
320 '''''''''''' draw circular mirrors
321 for I=-2 to 2:for J=-1 to 2
322 circle (XORIG+I*SCALE,YORIG-J*SCALE),RADIUS*SCALE,11 ' 11=cyan
323 next J:next I
324 circle (XORIG,YORIG),5,7 ' 7=grey, mark origin
329 '
330 '''''''''''' initial position for photon path polyline
331 Xfrom=XORIG+INIX*SCALE
332 Yfrom=YORIG-INIY*SCALE
339 '
400 '''''''''''''''''''' MAIN LOOP entrypoint for TEXT '''''''''''''''''
401 '
402 '''''''''''' initialise global time Gt and position Gx,Gy
403 Gt=0
404 Gx=INIX
405 Gy=INIY+Butterfly
406 if IsText print "Start";:gosub 800 ' Text: display coord
409 '
410 '''''''''''' velocity (cos and sin pair): Vx, Vy, angle A
411 A=0
412 Vx=1
413 Vy=0
419 '
500 '''''''''''''''''''' MAIN LOOP re-entrypoint for each step '''''''''
501 '
510 '''''''''''' global time GtS and position GxS,GyS at start of step
511 GtS=Gt
512 GxS=Gx
513 GyS=Gy
519 '
520 '''''''''''' provisional next step global position
522 Gx+=MOVE*Vx
523 Gy+=MOVE*Vy
529 '
530 '''''''''''' nearest gridpoint at end of step: Gi,Gj
531 Gi=round(Gx)
532 Gj=round(Gy)
539 '
540 '''''''''''' use local coord wrt Gi,Gj of start of step, LxS,LyS
541 LxS=GxS-Gi
542 LyS=GyS-Gj
549 '
600 '''''''''''' has our line approached Gi,Gj too closely? To
601 ' calculate the closest distance Lu, we want (p,q) dot (r,s),
602 ' where (p,q) any point on our line, we could take (LxS,LyS), and
603 ' where (r,s) a unit radius _|_ our line, we can take (Vy,-Vx)
604 '
605 Lu=Vy*LxS-Vx*LyS
606 '
607 if Lu>=+RADIUS goto 700 ' no we stayed clear
608 if Lu<=-RADIUS goto 700 ' no we stayed clear
609 '
610 ' Now we have Lu, essentially one of two orthogonal coord in a
611 ' rotated frame Lu,Lv (Lv in direction of motion, Lu constant),
612 ' and with origin Gi,Gj, just like LxS,LyS. Let's get the other.
613 ' First, the initial one @ start of step, i.e. for LxS,LyS:
614 '
615 LvS=Vx*LxS+Vy*LyS
619 '
620 '''''''''''' Next, the one at time of bouncing, via angle B wrt A:
621 B=-acos(Lu/RADIUS) ' minus because we want negative Lv solution
622 Lv=sin(B)*RADIUS
629 '
630 '''''''''''' Time of bounce
631 Lt=Lv-LvS ' local time := relative to start of step
632 Gt=GtS+Lt ' global time for display only
633 '
634 if Lt<0 goto 700 ' that's fine, we had that one
635 if Lt>=MOVE goto 700 ' fine too, we'll get it later
639 '
640 '''''''''''' x,y position of bounce
641 Gx=GxS+Vx*Lt
642 Gy=GyS+Vy*Lt
643 '
644 if IsText print "Bounce";
645 gosub 800 ' Text: display coord, Graph: draw line
649 '
650 '''''''''''' new direction
651 A+=2*B:while A<0:A+=TWOPI:wend
652 Vx=cos(A)
653 Vy=sin(A)
659 '
660 '''''''''''' post-bounce remainder of step
661 Lt-=MOVE ' this is -(remaining time)
662 Gx-=Vx*Lt
663 Gy-=Vy*Lt
669 '
700 '''''''''''''''''''' MAIN LOOP end-of-step code ''''''''''''''''''''
701 '
702 Gt=GtS+MOVE
709 '
710 '''''''''''' waypoint for debugging, commented out
711 ' if not IsText circle (XORIG+Gx*SCALE,YORIG-Gy*SCALE),2,14
712 '
720 '''''''''''' is there a next step?
721 if Gt<9.9 goto 500
729 '
790 if IsText print "End";
799 '
800 '''''''''''''''''''' DISPLAY '''''''''''''''''''''''''''''''''''''''
801 '
802 if IsText print " at t =";Gt:print "@ (x,y)=(";Gx;",";Gy;")":return
809 '
810 Xto=XORIG+Gx*SCALE
811 Yto=YORIG-Gy*SCALE
812 line (Xfrom,Yfrom)-(Xto,Yto),14 ' 14=yellow
813 Xfrom=Xto
814 Yfrom=Yto
815 return
819 '
900 '''''''''''''''''''' DISPLAY distance reached (text only) ''''''''''
901 Gr=sqrt(Gx*Gx+Gy*Gy)
902 print "Distance:";Gr;"--- Press a key... ";
903 I=input$(1):print
904 return
--------------8<--------------------------------------------->8--------------
--marijke