Skip to content

Useful mathematical words in Forth.

License

Notifications You must be signed in to change notification settings

frno7/forthmath

Folders and files

NameName
Last commit message
Last commit date

Latest commit

 

History

36 Commits
 
 
 
 
 
 
 
 
 
 
 
 

Repository files navigation

Mathematics in Forth

Build Status

In the examples below the top stack item is shown right-most, following Forth convention.

Primes

The file math/prime.fth defines words related to prime numbers.

The word primes gives a number of primes on the data stack. For example, 9 primes gives 23 19 17 13 11 7 5 3 2. Using traverse-primes one can for example compute the nth prime number pn by defining prime as:

: prime' ( n n -- n ) nip true ; ( True to obtain next prime number. )
: prime  ( n -- n ) 0 swap ['] prime' traverse-primes ;

Thus 9 prime gives 23 corresponding to p9. Computing the prime can be slow especially for large indices. The words prime-lower and prime-upper give quick estimates of lower and upper bounds. The inequalities prime-lowerprimeprime-upper hold for any number.

The prime-counting function denoted by π(n) is defined by the word prime-count. For example, 23 prime-count gives 9 because there are 9 primes below or equal to 23. Since this computation can be slow there are words for quick estimates of lower and upper bounds, where the inequalities prime-count-lowerprime-countprime-count-upper hold for any number.

The word primes-preceding gives on the data stack all primes below or equal to a certain number. Thus 23 primes-preceding gives 23 19 17 13 11 7 5 3 2 9 where the top integer 9 is the number of primes below or equal to 23. The word traverse-primes-preceding can be used to traverse all primes below or equal to a given number.

Prime decomposition

The file math/factor.fth defines words related to integer factorisation.

The word factors factorises a given integer into primes. For example, 12857=13·23·43 and so 12857 factors gives 43 23 13 3 on the data stack, where the top integer 3 indicates the number of prime factors and then the factors 13, 23 and 43 follow from smallest to largest.

Using traverse-factors one can for example define words for finding the smallest or largest prime factor of a given integer:

: min-factor' ( n n -- n) nip false ; ( False to ignore further prime factors. )
: min-factor  ( n -- n ) 0 swap ['] min-factor' traverse-factors ;

: max-factor' ( n n -- n) nip true ; ( True to obtain next prime factor. )
: max-factor  ( n -- n ) 0 swap ['] max-factor' traverse-factors ;

Thus 12857 min-factor gives 13 and 12857 max-factor gives 43.

The words prime? and composite? test whether a number is prime or composite. Thus 17 prime? gives true since 17 is a prime number but 18 prime? gives false since 18 is not a prime number. Numbers below 2 are neither prime nor composite.

Some integers such as 22477=7·13·13·19 have repeating primes so that 22477 factors gives 19 13 13 7 4 where 13 is repeated twice. The word factor-exponents factorises into distinct primes with corresponding exponents. Thus 22477 factor-exponents gives 19 1 13 2 7 1 3 where the top integer 3 indicates the number of distinct primes, followed by the pairs of primes and exponents corresponding to 71·132·191.

The Pollard Monte Carlo factorisation method is potentially very fast, but since it is a probabalistic method it can be slow and even fail to find factors of composite numbers. For example, pollard-rho-factors 4267640728818788929 gives the two factors 3456789019 and 1234567891 within a fraction of a second on a 64 bit Forth system.

Divisors

The files math/divisor.fth and math/gcd.fth define words related to divisor functions.

The word divisors gives the divisors of a number. For example, 52 divisors gives 52 26 13 4 2 1 6 where the top integer 6 indicates the divisor count and then the divisors follow from smallest to largest. traverse-divisors can be used to iterate over all divisors.

The sum of positive divisors function denoted by σx(n) is defined by divisor-sum. With x = 0 the function corresponds to divisor-count.

The word gcd gives the greatest common divisor. For example, 12 18 gcd gives 6. The word lcm gives the corresponding least common multiple. For example, 6 21 lcm gives 42. The word extended-gcd implements the extended Euclidean algorithm that gives two addional integers related to Bézout’s identity. For example, 240 46 extended-gcd gives -9 47 2 where 2 is the greatest common divisor and the identity holds as -9·240 + 47·46 = 2 = gcd(240, 46).

Rational numbers

The file math/rational.fth defines words related to rational numbers.

Rationals are represented as two numbers on the stack with the denominator above the cell containing the numerator. The words q+, q-, q* and q/ define rational arithmetic in a natural way. For example, 2 3 5 7 q+ gives 29 21 which corresponds to 2/3 + 5/7 = 29/21. Rational arithmetics overflow easily, especially intermediate calculations.

Exponents

The file math/exponent.fth defines words related to exponentiation.

The word ** computes integer exponentiation using efficient exponentiation by squaring. For example, 3 4 ** gives 81 which corresponds to 34 = 81.

Logarithms

The file math/log.fth defines words related to logarithms.

The words log2-floor and log2-ceiling compute log2 using integer arithmetics only, where the inequalities log2-floor ≤ log2log2-ceiling hold for any number. Similarly, ln-lower and ln-upper give quick estimates with the inequalities ln-lower ≤ ln ≤ ln-upper.

Modular arithmetic

The file math/modulo.fth defines words related to modular arithmetic.

The words +mod, -mod, *mod and **mod define modular arithmetic in a natural way. For example, the word **mod computes modular exponentiation and so 5 3 13 **mod gives 8 which corresponds to 53 = 125 = 8 (mod 13).

Matrices

The file math/matrix.fth defines words related to matrices.

Matrices are laid out in the following way on the stack: a 2×3 matrix having dimensions 2 rows and 3 columns with elements a, b and c in the first row and d, e and f in the second row is represented on the stack as a b c d e f 3 2, or stated with a different identation:

a b c
d e f  3 2

In data-space this matrix has the following layout:

Cell 0 1 2 3 4 5 6 7
Value 2 3 a b c d e f

Adding, subtracting, negating and multiplying matrices on the stack are defined by the words matrix+, matrix-, matrix-negate and matrix*. The word matrix** defines matrix exponentiation, and the words matrix0 and matrix1 give the zero and identity matrices. The word matrix-drop drops a matrix from the stack and matrix-element gives the element of particular row and column indices. The words matrix-dimensions, matrix-rows and matrix-cols give matrix dimensions. For example, multiplying a 2×3 matrix by a 3×2 matrix

     2   3   5
     7  11  13  3 2

    17  19
    23  29
    31  37  2 3

    matrix*

gives a 2×2 matrix:

   258 310
   775 933  2 2

Matrices can be moved and copied between the stack and data-space in several ways. The word matrix>allot moves a matrix off the stack to data-space appropriately reserved for by allot. The word allot>matrix is the inverse. The word allot+matrix copies a matrix from data-space to the stack, and the word matrix-allot removes a matrix from data-space by adjusting the here data-space pointer accordingly. The words matrix0allot and matrix1allot create the zero and identity matrices in data-space.

Special numbers

The files math/binomial.fth, math/bernoulli.fth and math/faulhaber.fth define words related to special numbers.

The word binomial corresponds to binomial numbers. For example, 7 3 binomial gives 35. The word bernoulli corresponds to Bernoulli numbers. For example, 12 bernoulli gives the rational number -691 2730. The word faulhaber corresponds to Faulhaber’s formula. For example, 7 2 faulhaber gives 91 which corresponds to 12+ 22+ 32+ 42+ 52+ 62=91. In general 1 faulhaber corresponds to the triangular numbers, 2 faulhaber to the square pyramidal numbers, 3 faulhaber to the squared triangular numbers, and so on.

Fibonacci numbers

The file math/fibonacci.fth defines words related to Fibonacci numbers.

The word fibonaccis gives a Fibonacci sequence with a given number of terms on the data stack. Thus 10 fibonaccis gives 34 21 13 8 5 3 2 1 1 0 starting at F0. The sequence is infinite in principle but since the terms grow at an exponential rate the number of useful terms is small due to integer overflow. The word traverse-fibonacci can be used to iterate over the Fibonacci sequence. The word fibonacci corresponds to Fn and thus 9 fibonacci gives 34, including negative indices such as -6 fibonacci giving -8.

Fibonacci numbers can also be computed in matrix form using exponentiation of a 2×2 system of linear difference equations:

: fibonacci-matrix ( -- 4 * n n n )
	1 1
	1 0 2 2 ;
: fibonacci-mod ( n n -- n )
	0 { n m r }
	fibonacci-matrix n m matrix**mod
	0 1 matrix-element to r
	matrix-drop r ;

For example, one can verify that the last two digits of F19 by 19 100 fibonacci-mod gives 81. Since matrix exponentiation is very efficient, one can compute much larger Fibonacci numbers. For example, 4267640728818788929 10000 fibonacci-mod quickly gives 4129 which corresponds to F4267640728818788929 (mod 10000).

Auxiliary

The files aux/nallot.fth, aux/reverse.fth and aux/sort.fth define auxiliary words.

Many algorithms need direct access memory to work effectively. The n>allot and nallot> words move n cells from the stack to data-space and vice versa. This makes it possible to use data space as a scratch pad, similar to >r and nr>.

The word reverse reversers a given number of items on the stack, and the word ndrop drops n items from the stack.

The generic sort function takes an execution token for comparisons. The words sort> and sort< order a given number of items on the stack, in ascending or descending order, and are implemented as:

: sort> ( n * x n -- n * x ) ['] > sort ; \ Sort in ascending order.
: sort< ( n * x n -- n * x ) ['] < sort ; \ Sort in descending order.