6502.org Forum  Projects  Code  Documents  Tools  Forum
It is currently Sun Nov 10, 2024 9:38 pm

All times are UTC




Post new topic Reply to topic  [ 87 posts ]  Go to page Previous  1, 2, 3, 4, 5, 6
Author Message
PostPosted: Sun Mar 24, 2024 1:43 am 
Offline

Joined: Wed Jan 08, 2014 3:31 pm
Posts: 578
Happy belated pi day!

My latest CS-101 Forth program computes pi using 16 bit fixed point arithmetic.
Link in github: https://github.com/Martin-H1/Forth-CS-101/blob/master/pi.fs

When run it produces this output:
Code:
pi  ok 2
.s <2> 25733 8192  ok 2

That's approximately 3.1412 which isn't half bad for 16 bits. I tested a 24 bit version and it was accurate to six decimal places. It should be possible to use a 32 bit double word arithmetic for better precision. The basic logic would remain the same. It might be a neat enhancement to perform the division and output the digits in decimal. But I will have to think that one over.

Here's the raw code if you don't want to follow the link:
Code:
\ Calculate pi using the Nilakantha infinite series. While more complicated
\ than the Leibniz formula, it is fairly easy to understand, and converges
\ on pi much more quickly.

\ The formula takes three and alternately adds and subtracts fractions with
\ a numerator of 4 and denominator that is the product of three consecutive
\ integers. So each subsequent fraction begins its set of integers with the
\ highest value used in the previous fraction.

\ Described in C syntax n starts at 2 and iterates to the desired precision.
\ pi = 3 + 4/((n)*(++n)*(++n)) - 4/((n)*(++n)*(++n)) + ...

\ Here's a three iteration example with an error slightly more than 0.0007.
\ pi = 3 + 4/(2*3*4) - 4/(4*5*6)
\        + 4/(6*7*8) - 4/(8*9*10)
\        + 4/(10*11*12) - 4/(12*13*14)
\    = 3.14088134088

\ normally this requires floating point arithmetic, but we're using fixed point
\ multiplying by a scaling factor and returned as a ratio.

1 13 lshift constant rescale
rescale 3 * constant three
rescale 4 * constant four

variable n

\ calculates (n)*(++n)*(++n)
: denominator ( -- product )
  n @ dup 1+ dup 1+ dup n ! * * ;

\ calculates a single scaled quotient term
: quotient ( -- q )
  four denominator / ;

\ calculates Qn - Qn+1
: calc_term ( -- q )
  quotient quotient - ;

\ Computes pi as a ratio of integers
: pi ( -- numerator denominator )
  2 n !
  three
  begin
    calc_term dup 0> while +
  repeat
  drop
  rescale ;


Top
 Profile  
Reply with quote  
PostPosted: Sun Mar 24, 2024 7:01 am 
Offline

Joined: Mon Jan 19, 2004 12:49 pm
Posts: 958
Location: Potsdam, DE
Youtube threw this at me a couple of days ago: https://www.youtube.com/watch?v=SPTzzSuBFlc

It calculates pi by the same mechanism (I think!) using relays...

Neil

(but on the other hand, 355 / 113 is already pi to six places: 3,14159292. How often would one need to go beyond that (it's two metres out calculating the equator!)?. I think I learned Forth's */ from Brodie using that example.)


Top
 Profile  
Reply with quote  
PostPosted: Sun Mar 24, 2024 7:09 am 
Offline
User avatar

Joined: Fri Aug 30, 2002 1:09 am
Posts: 8539
Location: Southern California
barnacle wrote:
(but on the other hand, 355 / 113 is already pi to six places: 3,14159292. How often would one need to go beyond that (it's two metres out calculating the equator!)?. I think I learned Forth's */ from Brodie using that example.)

355/113 is off by only +8.5E-8.  65298/20785 is only slightly better, at -5.1E-8.  I have a bunch of rational-number approximations for use in scaled-integer operations at http://wilsonminesco.com/16bitMathTable ... pprox.html

_________________
http://WilsonMinesCo.com/ lots of 6502 resources
The "second front page" is http://wilsonminesco.com/links.html .
What's an additional VIA among friends, anyhow?


Top
 Profile  
Reply with quote  
PostPosted: Sun Mar 24, 2024 8:02 am 
Offline
User avatar

Joined: Thu Dec 11, 2008 1:28 pm
Posts: 10977
Location: England
Thanks Martin that's a nice program implementing a very nice series! It's a new series for me - perhaps best described as an accelerated version of the Leibniz-Gregory series, but discovered much earlier. As often the case, discoveries made in different parts of the world took a long time to be communicated in either direction.

I gather that the same technique can yield a series which converges as the fifth power, and ultimately a series involving factorials which converges with over 13 digits gain from each term.

(The wonderful relay-based machine in the video linked in barnacle's post uses the much more recent BBP formula. I posted a Basic program over here to understand the approach better - I also posted a simplified version, perhaps it would be amusing to redo in Forth.)

Some readings:
Nilakantha's formula for pi (very very brief elaboration)
Nilakantha, Euler and pi by Shailesh A Shirali (1997), quite accessible, with references
The Discovery of the Series Formula for π by Leibniz, Gregory and Nilakantha by Ranjan Roy (1990)
Nilakantha's accelerated series for pi by David Brink (2015)


Top
 Profile  
Reply with quote  
PostPosted: Sun Mar 24, 2024 12:13 pm 
Offline
User avatar

Joined: Thu May 28, 2009 9:46 pm
Posts: 8481
Location: Midwestern USA
barnacle wrote:
Youtube threw this at me a couple of days ago: https://www.youtube.com/watch?v=SPTzzSuBFlc

It calculates pi by the same mechanism (I think!) using relays...

My first large-scale electronics project was a relay “computer,” built over 60 years ago using long-frame cradle relays and stepping relays, both scavenged from scrapped pinball machines.  To say it was noisy when operating was an understatement.  :D  It could add or subtract, but was massively incapable of computing pi.

In the video, the narrator laments about reliability.  Yet it appears he left out a key element of good relay circuit design.  None of his relays appears to have a coil flyback diode.  All that contact arcing has to be a killer, not to mention a significant source of heat.  The irony is having a flyback diode causes a relay to release more slowly, which means that his contraption would take even longer to compute pi.  :shock:

_________________
x86?  We ain't got no x86.  We don't NEED no stinking x86!


Top
 Profile  
Reply with quote  
PostPosted: Sun Mar 24, 2024 2:14 pm 
Offline

Joined: Wed Jan 08, 2014 3:31 pm
Posts: 578
@All, thanks for the kudos. Also, the 355 / 113 approximation is one of those better than it has any right to be ways to calculate pi. Using it a PDP-8 could accurately calculate all manner of things in only 12 bits.

@barnacle, I saw that video as well, which prompted me to write this program. I intentionally used a different algorithm because I am saving the Bailey–Borwein–Plouffe formula for an assembler program. As with the relay computer, it allows much simpler arithmetic routines.

@BigEd, thanks for the links about the various series, they look like interesting reading.

@BDD, I wondered about flyback diodes. I have never used relays for digital logic, but I have used them with microcontrollers to switch higher voltages, and the circuit I copied had a flyback diode.


Top
 Profile  
Reply with quote  
PostPosted: Sun Mar 24, 2024 3:09 pm 
Offline
User avatar

Joined: Thu May 28, 2009 9:46 pm
Posts: 8481
Location: Midwestern USA
Martin_H wrote:
@BDD, I wondered about flyback diodes. I have never used relays for digital logic, but I have used them with microcontrollers to switch higher voltages, and the circuit I copied had a flyback diode.

A flyback diode is de rigueur if a relay is being driven by a solid-state device.  Even a 5-volt relay’s coil can generate a transient of several hundred volts if not suppressed.  Needless to say, the device driving the relay won’t survive that.  :shock:

I should note that the diode’s transient current rating should be fairly high.  In most applications, I use a 1N4001, which has an IFSM rating of 30 amps over 8.3 milliseconds.  That should be more than adequate...unless your circuit is controlling traction motor contactors on a subway car.  :D

_________________
x86?  We ain't got no x86.  We don't NEED no stinking x86!


Top
 Profile  
Reply with quote  
PostPosted: Mon Mar 25, 2024 1:59 pm 
Offline

Joined: Thu Mar 12, 2020 10:04 pm
Posts: 704
Location: North Tejas
Martin_H wrote:
* While iterating a single cell I use multiply row * width nine times. Even on modern hardware integer multiply is about five times more costly than addition. Reusing the first multiply would be better.

Update: I realized that if row enumeration added width instead of 1 I would effectively by multiplying through addition. But only once per line, rather than nine times per cell. So I updated the code and it should be a performance boost.

You can multiply once to get the address of the center cell. Then add +-0 or width, +-0 or 1 to access adjacent cells.


Top
 Profile  
Reply with quote  
PostPosted: Tue Mar 26, 2024 2:26 am 
Offline

Joined: Fri Apr 15, 2016 1:03 am
Posts: 139
There is a problem with the Pi calculator program on some systems.
The word "four" evaluates to a negative number on a 16 bit signed 2s complement binary system.
Changing the divide in the word "quotient" from signed to unsigned can fix this problem.


Attachments:
F_Pi5 1.zip [1.74 KiB]
Downloaded 45 times
Top
 Profile  
Reply with quote  
PostPosted: Tue Mar 26, 2024 12:52 pm 
Offline

Joined: Wed Jan 08, 2014 3:31 pm
Posts: 578
leepivonka wrote:
There is a problem with the Pi calculator program on some systems.
The word "four" evaluates to a negative number on a 16 bit signed 2s complement binary system.
Changing the divide in the word "quotient" from signed to unsigned can fix this problem.

Thanks for the correction. I made the change and pushed it to GitHub.


Top
 Profile  
Reply with quote  
PostPosted: Sun Jun 16, 2024 8:07 pm 
Online

Joined: Fri May 05, 2017 9:27 pm
Posts: 895

I'm sure most of us are familiar with Euclid's algorithm for finding the greatest common divisor of two numbers. The following is one way of expressing this algorithm in Forth.
Code:
: GCD  ( N1 N2 -- GCD )
   BEGIN
      TUCK MOD ?DUP 0=
   UNTIL ;

What if /MOD is used instead of MOD and the quotient is displayed?
After watching a video about a fraction problem I came up with this:
Code:
: .CFTERMS  ( N1 N2 -- )
   BEGIN
      TUCK /MOD . ?DUP 0=
   UNTIL
   DROP ;

.CFTERMS takes a numerator and a denominator and displays the terms of a continued fraction. Given the numbers used in the problem from the video, 17 and 10, .CFTERMS displays 1 1 2 3.
Code:
 17               1
---- =    1 + ---------
 10                 1
              1 + -----
                      1
                  2 + -
                      3

The word RATIO takes the terms of a continued fraction and a count of the terms and returns a ratio, the numerator and denominator of a fraction equal to the continued fraction. Given the numbers 1 1 2 3 and 4 (for the count), RATIO leaves 17 and 10 on the data stack.
Code:
: RATIO  ( T1 ... TM M -- N1 N2 )
   1 0 ROT 0
   ?DO
      OVER 2SWAP * UNDER+
   LOOP ;

UNDER+ is a code word on my system. It can be defined in high level as:
Code:
: UNDER+  ( A B C -- A+C B )
   ROT + SWAP ;



Top
 Profile  
Reply with quote  
PostPosted: Sun Jul 21, 2024 11:48 pm 
Online

Joined: Fri May 05, 2017 9:27 pm
Posts: 895

The following infinite continued fraction is equal to the golden ratio.
Code:
                    1
1 + ----------------------------------
                       1
    1 + ------------------------------
                         1
        1 + --------------------------
                           1
            1 + ----------------------
                             1
                1 + ------------------
                               1
                    1 + --------------
                                 1
                        1 + ----------
                                   1
                             1 + -----
                                     1
                                 1 + -
                                     ...



Top
 Profile  
Reply with quote  
Display posts from previous:  Sort by  
Post new topic Reply to topic  [ 87 posts ]  Go to page Previous  1, 2, 3, 4, 5, 6

All times are UTC


Who is online

Users browsing this forum: No registered users and 4 guests


You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot post attachments in this forum

Search for:
Jump to: