Page 6 of 7
Re: CS-101 problems and their FORTH solutions.
Posted: Sun Mar 24, 2024 1:43 am
by Martin_H
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-1 ... ster/pi.fs
When run it produces this output:
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: Select all
\ 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 ;
Re: CS-101 problems and their FORTH solutions.
Posted: Sun Mar 24, 2024 7:01 am
by barnacle
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.)
Re: CS-101 problems and their FORTH solutions.
Posted: Sun Mar 24, 2024 7:09 am
by GARTHWILSON
(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.)
Re: CS-101 problems and their FORTH solutions.
Posted: Sun Mar 24, 2024 8:02 am
by BigEd
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)
Re: CS-101 problems and their FORTH solutions.
Posted: Sun Mar 24, 2024 12:13 pm
by BigDumbDinosaur
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. 
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. 
Re: CS-101 problems and their FORTH solutions.
Posted: Sun Mar 24, 2024 2:14 pm
by Martin_H
@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.
Re: CS-101 problems and their FORTH solutions.
Posted: Sun Mar 24, 2024 3:09 pm
by BigDumbDinosaur
@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.
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.
Re: CS-101 problems and their FORTH solutions.
Posted: Mon Mar 25, 2024 1:59 pm
by BillG
* 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.
Re: CS-101 problems and their FORTH solutions.
Posted: Tue Mar 26, 2024 2:26 am
by leepivonka
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.
Re: CS-101 problems and their FORTH solutions.
Posted: Tue Mar 26, 2024 12:52 pm
by Martin_H
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.
Re: CS-101 problems and their FORTH solutions.
Posted: Sun Jun 16, 2024 8:07 pm
by JimBoyd
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: Select all
: 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: Select all
: .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: Select all
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: Select all
: 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: Select all
: UNDER+ ( A B C -- A+C B )
ROT + SWAP ;
Re: CS-101 problems and their FORTH solutions.
Posted: Sun Jul 21, 2024 11:48 pm
by JimBoyd
The following infinite continued fraction is equal to the golden ratio.
Code: Select all
1
1 + ----------------------------------
1
1 + ------------------------------
1
1 + --------------------------
1
1 + ----------------------
1
1 + ------------------
1
1 + --------------
1
1 + ----------
1
1 + -----
1
1 + -
...
Re: CS-101 problems and their FORTH solutions.
Posted: Wed Dec 04, 2024 3:06 pm
by Martin_H
Good addition to the list of Forth examples.
Re: CS-101 problems and their FORTH solutions.
Posted: Sun Aug 31, 2025 8:54 pm
by JimBoyd
Not sure if this belongs in the thread "CS-101 problems and their FORTH solutions" or if I should have started a new Forth and mathematics thread. This example involves repetends. It was easier to program than to explain. After the basic introduction I attempt a top down explanation of the code.
My Forth does not have floating point support; however, Forth has support for scaled integer math. I believe most of us are familiar with the ratio 355 113, the approximation for Pi used in sixteen bit Forths. Using this ratio to calculate the perimeter of a circle from its diameter is straight forward.
This displays the number 78. What about displaying more precision?
A floating point result would have included a decimal point and a few digits after. The following is a proof of concept word.
Code: Select all
: TEST
25 PI */MOD 0 .R
ASCII . EMIT
3 0
DO
10 113 */MOD 0 .R
LOOP ;
TEST
This displays 78.539 .
Each digit after the decimal point is extracted by multiplying the remainder by ten and dividing by the original divisor.
Here is a generalized word to perform the same operation. It takes a number followed by the two to be used as the scale factor and how many digits after the decimal point to display.
Code: Select all
: SCALE. ( N1 N2 N3 CNT -- )
>R -ROT 2PICK */MOD
CR 0 .R ASCII . EMIT
R> 0
?DO
10 2PICK */MOD 0 .R
LOOP
2DROP ;
My Forth has the ability to log what is displayed on the screen to the printer. VICE simply writes it to a print.dump file.
Here is an excerpt from one 'session log'. The first example is the original one to show the perimeter of a circle with a diameter of twenty five.
Code: Select all
25 PI 4 SCALE.
78.5398 OK
3 30 7 6 SCALE.
12.857142 OK
The second example is "how many weeks are in three thirty-day months?"
A little more precision.
Shows the decimal expression for six sevenths is a recurring decimal.
Moving on to ratios:
That last example showed when some ratios are expressed as a decimal expression, there are repeating decimal digits. Some ratios do not have a repeating decimal expression.
For example, five fourths is just 1.25
650 divided by 64 is 10.15625 with no repeating digits.
5 divided by 6 is 0.83 with the 3 repeating.
I wrote the word RECURRING to take two unsigned numbers as its parameters and display the result of dividing the first by the second showing the decimal expression. It stops showing digits when the end of the repetend, the portion which repeats, is reached or if the decimal expression terminates. The next line displays two numbers, the first is the number of non-repeating digits and the second is the length of the repetend.
Some excerpts from a session log.
Code: Select all
5 6 RECURRING \ I typed this.
0.83 \ The computer's response
1 1 OK \
One non-repeating digit and one repeating digit.
Code: Select all
1 19 RECURRING
0.052631578947368421
0 18 OK
Zero non-repeating digits and eighteen repeating digits.
Here's a big one:
Code: Select all
1 113 RECURRING
0.0088495575221238938053097345132743362
831858407079646017699115044247787610619
469026548672566371681415929203539823
0 112 OK
Okay, it's not really that big. The biggest inverse with a maximal length repetend RECURRING can handle with a sixteen bit Forth on a Commodore 64 has a repetend length of sixty five thousand four hundred and forty six digits (65446).
Code: Select all
1 65447 RECURRING
0.0000152795391690986599844148700475193668158968325515302458477852307
.
.
.
054609072990358610784298745549834217
0 65446
On a Commodore 64 simulation (VICE) running Fleet Forth, it took just over sixteen and a half minutes from typing RECURRING to the last line displayed.
For a divisor, N , there can only be N possible remainders, including zero. If the remainder is zero, the decimal expression terminates. If the remainder is never zero, the remainder will be repeated within N iterations, therefore the decimal expression will repeat with a combined length of non-repeating and repeating digits less than N .
As I mentioned, each digit is extracted by multiplying the remainder by ten and dividing by the original divisor; however, in the word RECURRING , the remainder is not multiplied by ten, it is multiplied by the value of BASE .
At this point it might just be easier to call the decimal point the radix point.
The radix expression will be different in different number bases.
The inverse of sixty four in different bases:
Code: Select all
1 64 RECURRING
0.015625
6 0 OK
1 64 HEX RECURRING
0.04
2 0 OK
DECIMAL OK
1 64 8 BASE ! RECURRING
0.01
2 0 OK
DECIMAL OK
1 64 2 BASE ! RECURRING
0.000001
110 0 OK
110 DECIMAL . 6 OK
1 64 3 BASE ! RECURRING
0.0001021011122022
0 121 OK
121 DECIMAL . 16 OK
In base three, the inverse of sixty four is a repeating expression with no non-repeating digits and a repetend of length sixteen.
How RECURRING works.
Code: Select all
: RECURRING ( U1 U2 -- )
(RECURRING) RESULT ;
(RECURRING) takes two unsigned integers. It does not display anything. (RECURRING) sets the variable OFFSET to the start of the repetend. It returns those same two integers and a count of how many digits to display (the sum of the non-repeating digits and the repetend length) and leaves all three parameters for RESULT.
Code: Select all
: RESULT ( U1 U2 CNT -- )
-ROT TUCK 0 SWAP UM/MOD
?CR 0 U.R ASCII . EMIT
2PICK 0
?DO
BASE @ UM* 2PICK UM/MOD
1 ?LINE 0 .R
LOOP
2DROP CR OFFSET @ DUP U. - U. ;
RESULT takes two unsigned integers and the count, CNT . It returns the result of dividing the first integer by the second, displaying the unsigned quotient and a decimal point. RESULT then cycles through a ?DO LOOP for CNT times. In the loop, it multiplies the remainder by the value of BASE and divides the result by the original divisor. The quotient is displayed as the next digit and the remainder is used in the next loop cycle.
When the loop finishes, the original divisor and the final remainder are discarded. The offset is displayed (the number of non-repeating digits) followed by the difference of the count, CNT , and the offset. This is the repetend length.
Code: Select all
: (RECURRING) ( U1 U2 -- U1 U2 CNT )
OFFSET OFF
BEGIN
2DUP PROLOGUE REPETEND
IF
SWAP 0=
IF DUP OFFSET ! THEN
EXIT
THEN
1 OFFSET +!
AGAIN -;
(RECURRING) clears OFFSET before entering a loop to execute the words PROLOGUE and REPETEND . REPETEND either returns:
a FALSE flag or
the last remainder, the count of digits to display, and a TRUE flag.
If a FALSE flag is returned, OFFSET is incremented and the loop restarts.
If the last remainder, count, and TRUE flag are returned, (RECURRING) checks if the last remainder is zero. If the last remainder is zero, OFFSET is set to the count because there is no repetend in this case. Whether the last remainder is zero or not, (RECURRING) exits.
Code: Select all
VARIABLE OFFSET
: PROLOGUE ( U1 U2 -- U3 U2 )
TUCK 0 SWAP UM/MOD DROP
OFFSET @ 0
?DO
BASE @ UM* 2PICK UM/MOD DROP
LOOP
SWAP ;
PROLOGUE calculates the first remainder then runs a loop OFFSET times, each time calculating the next remainder.
Code: Select all
: REPETEND ( U1 U2 -- U3 CNT TRUE )
( U1 U2 -- FALSE )
2DUP OFFSET @
?DO
BASE @ UM* 2PICK UM/MOD DROP
2PICK OVER = OVER 0= OR
IF
I 1+ 2NIP UNLOOP
TRUE EXIT
THEN
LOOP
2DROP DROP FALSE ;
REPETEND uses the divisor as the loop limit because there will never be more digits to the right of the radix point than the size of the divisor. It uses the value of OFFSET as the initial loop index.
REPETEND saves the first remainder passed to it and compares each successive remainder to find when the sequence repeats. It also checks if the remainder is zero.
If the remainder is zero, there is no repetend and the sequence terminates.
If the remainder is not zero but it is equal to the first remainder passed to REPETEND , the repetend has been found.
Until a remainder of zero is found or a remainder which is the same as the remainder REPETEND started with is found, (RECURRING) keeps running through its BEGIN AGAIN loop until the remainder is zero or OFFSET points to the start of the repetend.
I only had to make a few changes to port this code to gforth, an ANSI Standard Forth.
I changed ?CR to CR and ASCII to [CHAR] and removed the phrase 1 ?LINE .
Re: CS-101 problems and their FORTH solutions.
Posted: Mon Sep 01, 2025 4:45 pm
by Martin_H
Jim, that's neat and shows what an integer only system is capable of.