Author Topic: Estimating Pi  (Read 154 times)

Online FellippeHeitor

  • QB64 Developer
  • Forum Resident
  • Posts: 2209
  • LET IT = BE
    • QB64.org
Estimating Pi
« on: August 01, 2020, 09:11:11 PM »
Of course we already have _PI, this is just a cool way to visualize it.

Code: QB64: [Select]
  1. SCREEN _NEWIMAGE(600, 600, 32)
  2.  
  3. DIM inCircle AS _UNSIGNED LONG, total AS _UNSIGNED LONG
  4.  
  5.     px = RND * _WIDTH
  6.     py = RND * _HEIGHT
  7.  
  8.     total = total + 1
  9.     IF dist(_WIDTH / 2, _HEIGHT / 2, px, py) <= 300 THEN
  10.         inCircle = inCircle + 1
  11.         PSET (px, py), _RGB32(100)
  12.     ELSE
  13.         PSET (px, py), _RGB32(70)
  14.     END IF
  15.     pi = (inCircle * 4) / total
  16.  
  17.     LOCATE 1, 1
  18.     PRINT pi
  19.  
  20.     CIRCLE (_WIDTH / 2, _HEIGHT / 2), 299, _RGB32(216, 144, 22)
  21.     _DISPLAY
  22.  
  23. FUNCTION dist! (x1!, y1!, x2!, y2!)
  24.     dist! = _HYPOT((x2! - x1!), (y2! - y1!))

As seen on:
« Last Edit: August 01, 2020, 09:12:53 PM by FellippeHeitor »

Offline _vince

  • Seasoned Forum Regular
  • Posts: 251
Re: Estimating Pi
« Reply #1 on: August 02, 2020, 12:45:23 AM »
It is a great simulation! Another one is called "buffon's needles", I wanted to write a program that demonstrates the experiment graphically and then numerically computes the results to iteratively approximated pi as more needles are dropped, open challenge I guess

ie
(https://en.wikipedia.org/wiki/Buffon%27s_needle_problem#/media/File:Buffon_needle_experiment_compressed.gif)
« Last Edit: August 02, 2020, 12:46:39 AM by _vince »

Offline SierraKen

  • Forum Resident
  • Posts: 931
Re: Estimating Pi
« Reply #2 on: August 02, 2020, 12:55:12 AM »
Really awesome Felippe! I have no idea how you and others find Pi, was just ages ago since school. But I remembered a year ago (which turned out to be exactly a year ago), I found code to output as many digits as you want. I posted it in another thread if anyone is interested. Just funny that it's been one year, around the Sun, using Pi LOL. Oh, I also updated it tonight to make it show how many digits it has left to calculate before it shows a Notepad of it.
« Last Edit: August 02, 2020, 01:06:16 AM by SierraKen »

Online FellippeHeitor

  • QB64 Developer
  • Forum Resident
  • Posts: 2209
  • LET IT = BE
    • QB64.org
Re: Estimating Pi
« Reply #3 on: August 02, 2020, 01:04:52 AM »
Well, I guess that only proves that life is a circle, @SierraKen 😂

@_vince looking forward to your project!

Offline SierraKen

  • Forum Resident
  • Posts: 931
Re: Estimating Pi
« Reply #4 on: August 02, 2020, 01:05:25 AM »
LOL true.

Offline bplus

  • Forum Resident
  • Posts: 4548
  • B+ nots
Re: Estimating Pi
« Reply #5 on: August 02, 2020, 12:49:32 PM »
It is a great simulation! Another one is called "buffon's needles", I wanted to write a program that demonstrates the experiment graphically and then numerically computes the results to iteratively approximated pi as more needles are dropped, open challenge I guess

ie
(https://en.wikipedia.org/wiki/Buffon%27s_needle_problem#/media/File:Buffon_needle_experiment_compressed.gif)

A very slow road to Pi:
Code: QB64: [Select]
  1. _TITLE "Buffon Pi Calc " 'b+ 2020-08-2
  2. '2020-08-02 try this with Buffon Pi Calc  pi ~ 2 * nTrials / nCross  from
  3. ' pi ~ 2*L*N/(T*H) L=length of needle, T=spacing between lines N = Number of trials (or pins) H = number that cross line
  4.  
  5. CONST xmax = 800, ymax = 600
  6. SCREEN _NEWIMAGE(xmax, ymax, 32)
  7. _DELAY .25
  8. L = 100 ' both pin length and line spacing
  9. N = 0 '   number of trials
  10. C = 0 '   number of pins that cross the line
  11. SS = 0 '  pin head and point lay in Same Space
  12.  
  13.     CLS
  14.     'redraw lines
  15.     FOR y = 0 TO _HEIGHT - 1 STEP L
  16.         LINE (0, y)-(_WIDTH - 1, y), &HFFFFFFFF
  17.     NEXT
  18.  
  19.     'drop a random pin
  20.     px1 = RND * (_WIDTH - 200) + 100
  21.     py1 = RND * (_HEIGHT - 200) + 100
  22.     ra = RND * _PI(2)
  23.     px2 = px1 + L * COS(ra)
  24.     py2 = py1 + L * SIN(ra)
  25.     LINE (px1, py1)-(px2, py2), _RGB32(50 + RND * 150)
  26.     N = N + 1 '                             another trial pin
  27.  
  28.     'did our pin cross a line
  29.     FOR y = 0 TO _HEIGHT - 1 STEP L
  30.         IF py1 >= y AND py1 < y + L THEN ' one side of pin is between y and y + L
  31.             IF py2 >= y AND py2 < y + L THEN ' the pin did NOT cross a line both the head and point lay in same space
  32.                 SS = SS + 1
  33.             END IF
  34.         END IF
  35.     NEXT
  36.     C = N - SS
  37.     'Print results of trials
  38.     LOCATE 1, 1: PRINT "Trials:"; N
  39.     LOCATE 2, 1: PRINT "Crossings:"; C
  40.     LOCATE 3, 1: PRINT "Pi est:"; 2 * N / C
  41.     _DISPLAY
  42.     IF N = 1 THEN
  43.         _DELAY 2
  44.     ELSEIF N < 5 THEN
  45.         _DELAY 1
  46.     ELSEIF N < 10 THEN
  47.         _DELAY .5
  48.     ELSEIF N < 50 THEN
  49.         _DELAY .25
  50.     ELSEIF N < 1000 THEN
  51.         _DELAY .05
  52.     ELSEIF N <= 10000 THEN
  53.         IF N MOD 1000 = 0 THEN _DELAY .25
  54.     ELSE
  55.         IF N MOD 10000 = 0 THEN _DELAY .05
  56.     END IF
  57.  
  58.  
  59.  

Offline _vince

  • Seasoned Forum Regular
  • Posts: 251
Re: Estimating Pi
« Reply #6 on: August 02, 2020, 01:10:45 PM »
Yes, that is it. I've uncommented CLS to watch the needles pile up. Nice work, bplus

Offline bplus

  • Forum Resident
  • Posts: 4548
  • B+ nots
Re: Estimating Pi
« Reply #7 on: August 02, 2020, 01:25:36 PM »
Thanks @_vince  this might be slightly more efficient for check needle crossing the line, skipping a FOR loop and maybe an extra IF THEN check:
Code: QB64: [Select]
  1. _TITLE "Buffon Pi Calc " 'b+ 2020-08-2
  2. '2020-08-02 try this with Buffon Pi Calc  pi ~ 2 * nTrials / nCross  from
  3. ' pi ~ 2*L*N/(T*H) L=length of needle, T=spacing between lines N = Number of trials (or pins) H = number that cross line
  4.  
  5. CONST xmax = 800, ymax = 600
  6. SCREEN _NEWIMAGE(xmax, ymax, 32)
  7. _DELAY .25
  8. L = 100 ' both pin length and line spacing
  9. N = 0 '   number of trials
  10. C = 0 '   number of pins that cross the line
  11. 'SS = 0 '  pin head and point lay in Same Space
  12.  
  13.     CLS
  14.     'redraw lines
  15.     FOR y = 0 TO _HEIGHT - 1 STEP L
  16.         LINE (0, y)-(_WIDTH - 1, y), &HFFFFFFFF
  17.     NEXT
  18.  
  19.     'drop a random pin
  20.     px1 = RND * (_WIDTH - 200) + 100
  21.     py1 = RND * (_HEIGHT - 200) + 100
  22.     ra = RND * _PI(2)
  23.     px2 = px1 + L * COS(ra)
  24.     py2 = py1 + L * SIN(ra)
  25.     LINE (px1, py1)-(px2, py2), _RGB32(50 + RND * 150)
  26.     N = N + 1 '                             another trial pin
  27.  
  28.     'did our pin cross a line
  29.     IF py1 < py2 THEN
  30.         IF py2 >= 100 * INT(py1 / 100) + L THEN C = C + 1
  31.     ELSE
  32.         IF py1 >= 100 * INT(py2 / 100) + L THEN C = C + 1
  33.     END IF
  34.     'FOR y = 0 TO _HEIGHT - 1 STEP L
  35.     '    IF py1 >= y AND py1 < y + L THEN ' one side of pin is between y and y + L
  36.     '        IF py2 >= y AND py2 < y + L THEN ' the pin did NOT cross a line both the head and point lay in same space
  37.     '            SS = SS + 1
  38.     '            EXIT FOR ' speed up a little
  39.     '        END IF
  40.     '    END IF
  41.     'NEXT
  42.     'C = N - SS
  43.     'Print results of trials
  44.     LOCATE 1, 1: PRINT "Trials:"; N
  45.     LOCATE 2, 1: PRINT "Crossings:"; C
  46.     LOCATE 3, 1: PRINT "Pi est:"; 2 * N / C
  47.     _DISPLAY
  48.     IF N = 1 THEN
  49.         _DELAY 2
  50.     ELSEIF N < 5 THEN
  51.         _DELAY 2
  52.     ELSEIF N < 10 THEN
  53.         _DELAY 2
  54.     ELSEIF N < 50 THEN
  55.         _DELAY .25
  56.     ELSEIF N < 1000 THEN
  57.         _DELAY .05
  58.     ELSEIF N <= 10000 THEN
  59.         IF N MOD 1000 = 0 THEN _DELAY .25
  60.     ELSE
  61.         IF N MOD 10000 = 0 THEN _DELAY .05
  62.     END IF
  63.  
  64.  
  65.  

Offline bplus

  • Forum Resident
  • Posts: 4548
  • B+ nots
Re: Estimating Pi
« Reply #8 on: August 02, 2020, 02:21:40 PM »
There might be a problem with the code or the Random Number Generator as I have a definite pattern of "random" needle color "waves" along the bottom and right borders of the screen.

Anyone have an idea of the error, blunder or glitch?

 


Both versions produce this pattern.
« Last Edit: August 02, 2020, 02:23:22 PM by bplus »

Offline _vince

  • Seasoned Forum Regular
  • Posts: 251
Re: Estimating Pi
« Reply #9 on: August 02, 2020, 02:33:57 PM »
looks like the bounds of the needle coordinates are too far inside the screen

Offline bplus

  • Forum Resident
  • Posts: 4548
  • B+ nots
Re: Estimating Pi
« Reply #10 on: August 02, 2020, 04:01:50 PM »
looks like the bounds of the needle coordinates are too far inside the screen

No I don't think it was that, here is the fix of the "pattern" and even speedier calc's seems 3.140 is limit of precision.

Code: QB64: [Select]
  1. _TITLE "Buffon Pi Calc " 'b+ 2020-08-2
  2. '2020-08-02 try this with Buffon Pi Calc  pi ~ 2 * nTrials / nCross  from
  3. ' pi ~ 2*L*N/(T*H) L=length of needle, T=spacing between lines N = Number of trials (or pins) H = number that cross line
  4.  
  5. CONST xmax = 800, ymax = 600
  6. SCREEN _NEWIMAGE(xmax, ymax, 32)
  7. _DELAY .25
  8.  
  9. L = 100 ' both pin length and line spacing
  10. N = 0 '   number of trials
  11. C = 0 '   number of pins that cross the line
  12.     'CLS
  13.     'redraw lines
  14.     'drop a random pin
  15.     px1 = RND * (_WIDTH - 201) + 100
  16.     py1 = RND * (_HEIGHT - 201) + 100
  17.     ra = RND * _PI(2)
  18.     IF RND < .5 THEN dir = -1 ELSE dir = 1 ' << this fixed the pattern developing along right and bottom border of screen
  19.     px2 = px1 + dir * L * COS(ra)
  20.     py2 = py1 + dir * L * SIN(ra)
  21.  
  22.     N = N + 1 '                             another trial pin
  23.  
  24.     'did our pin cross a line
  25.     IF py1 < py2 THEN
  26.         IF py2 >= 100 * INT(py1 / 100) + L THEN C = C + 1
  27.     ELSE
  28.         IF py1 >= 100 * INT(py2 / 100) + L THEN C = C + 1
  29.     END IF
  30.     IF N < 100000 THEN 'display
  31.         FOR y = 0 TO _HEIGHT - 1 STEP L
  32.             LINE (0, y)-(_WIDTH - 1, y), &HFFFFFFFF
  33.         NEXT
  34.         LINE (px1, py1)-(px2, py2), _RGB32(50 + RND * 150)
  35.  
  36.         'Print results of trials
  37.         LOCATE 1, 1: PRINT "Trials:"; N
  38.         LOCATE 2, 1: PRINT "Crossings:"; C
  39.         LOCATE 3, 1: PRINT "Pi est:"; 2 * N / C
  40.         _DISPLAY
  41.         IF N = 1 THEN
  42.             _DELAY 2
  43.         ELSEIF N < 5 THEN
  44.             _DELAY 2
  45.         ELSEIF N < 10 THEN
  46.             _DELAY 2
  47.         ELSEIF N < 50 THEN
  48.             _DELAY .25
  49.         ELSEIF N < 1000 THEN
  50.             _DELAY .05
  51.         ELSEIF N <= 10000 THEN
  52.             IF N MOD 1000 = 0 THEN _DELAY .25
  53.         ELSE
  54.             IF N MOD 10000 = 0 THEN _DELAY .05
  55.         END IF
  56.     ELSE
  57.         IF N MOD 1000000 = 0 THEN
  58.             LOCATE 1, 1: PRINT "Trials:"; N
  59.             LOCATE 2, 1: PRINT "Crossings:"; C
  60.             LOCATE 3, 1: PRINT "Pi est:"; 2 * N / C
  61.             _DISPLAY
  62.         END IF
  63.     END IF
  64.