Go to the first, previous, next, last section, table of contents.


Numerics

Bit-Twiddling

(require 'logical)

The bit-twiddling functions are made available through the use of the logical package. logical is loaded by inserting (require 'logical) before the code that uses these functions.

Function: logand n1 n1
Returns the integer which is the bit-wise AND of the two integer arguments.

Example:

(number->string (logand #b1100 #b1010) 2)
   => "1000"

Function: logior n1 n2
Returns the integer which is the bit-wise OR of the two integer arguments.

Example:

(number->string (logior #b1100 #b1010) 2)
   => "1110"

Function: logxor n1 n2
Returns the integer which is the bit-wise XOR of the two integer arguments.

Example:

(number->string (logxor #b1100 #b1010) 2)
   => "110"

Function: lognot n
Returns the integer which is the 2s-complement of the integer argument.

Example:

(number->string (lognot #b10000000) 2)
   => "-10000001"
(number->string (lognot #b0) 2)
   => "-1"

Function: logtest j k
(logtest j k) == (not (zero? (logand j k)))

(logtest #b0100 #b1011) => #f
(logtest #b0100 #b0111) => #t

Function: logbit? index j
(logbit? index j) == (logtest (integer-expt 2 index) j)

(logbit? 0 #b1101) => #t
(logbit? 1 #b1101) => #f
(logbit? 2 #b1101) => #t
(logbit? 3 #b1101) => #t
(logbit? 4 #b1101) => #f

Function: ash int count
Returns an integer equivalent to (inexact->exact (floor (* int (expt 2 count)))).

Example:

(number->string (ash #b1 3) 2)
   => "1000"
(number->string (ash #b1010 -1) 2)
   => "101"

Function: logcount n
Returns the number of bits in integer n. If integer is positive, the 1-bits in its binary representation are counted. If negative, the 0-bits in its two's-complement binary representation are counted. If 0, 0 is returned.

Example:

(logcount #b10101010)
   => 4
(logcount 0)
   => 0
(logcount -2)
   => 1

Function: integer-length n
Returns the number of bits neccessary to represent n.

Example:

(integer-length #b10101010)
   => 8
(integer-length 0)
   => 0
(integer-length #b1111)
   => 4

Function: integer-expt n k
Returns n raised to the non-negative integer exponent k.

Example:

(integer-expt 2 5)
   => 32
(integer-expt -3 3)
   => -27

Function: bit-extract n start end
Returns the integer composed of the start (inclusive) through end (exclusive) bits of n. The startth bit becomes the 0-th bit in the result.

Example:

(number->string (bit-extract #b1101101010 0 4) 2)
   => "1010"
(number->string (bit-extract #b1101101010 4 9) 2)
   => "10110"

Modular Arithmetic

(require 'modular)

Function: extended-euclid n1 n2
Returns a list of 3 integers (d x y) such that d = gcd(n1, n2) = n1 * x + n2 * y.

Function: symmetric:modulus n
Returns (quotient (+ -1 n) -2) for positive odd integer n.

Function: modulus->integer modulus
Returns the non-negative integer characteristic of the ring formed when modulus is used with modular: procedures.

Function: modular:normalize modulus n
Returns the integer (modulo n (modulus->integer modulus)) in the representation specified by modulus.

The rest of these functions assume normalized arguments; That is, the arguments are constrained by the following table:

For all of these functions, if the first argument (modulus) is:

positive?
Work as before. The result is between 0 and modulus.
zero?
The arguments are treated as integers. An integer is returned.
negative?
The arguments and result are treated as members of the integers modulo (+ 1 (* -2 modulus)), but with symmetric representation; i.e. (<= (- modulus) n modulus).

If all the arguments are fixnums the computation will use only fixnums.

Function: modular:invertable? modulus k
Returns #t if there exists an integer n such that k * n == 1 mod modulus, and #f otherwise.

Function: modular:invert modulus k2
Returns an integer n such that 1 = (n * k2) mod modulus. If k2 has no inverse mod modulus an error is signaled.

Function: modular:negate modulus k2
Returns (-k2) mod modulus.

Function: modular:+ modulus k2 k3
Returns (k2 + k3) mod modulus.

Function: modular:- modulus k2 k3
Returns (k2 - k3) mod modulus.

Function: modular:* modulus k2 k3
Returns (k2 * k3) mod modulus.

The Scheme code for modular:* with negative modulus is not completed for fixnum-only implementations.

Function: modular:expt modulus k2 k3
Returns (k2 ^ k3) mod modulus.

Prime Testing and Generation

(require 'primes)

This package tests and generates prime numbers. The strategy used is as follows:

The Miller-Rabin test is a Monte-Carlo test--in other words, it's fast and it gets the right answer with high probability. For a candidate that is prime, the Miller-Rabin test is certain to report "prime"; it will never report "composite". However, for a candidate that is composite, there is a (small) probability that the Miller-Rabin test will erroneously report "prime". This probability can be made arbitarily small by adjusting the number of iterations of the Miller-Rabin test.

Function: probably-prime? candidate
Function: probably-prime? candidate iter
Returns #t if candidate is probably prime. The optional parameter iter controls the number of iterations of the Miller-Rabin test. The probability of a composite candidate being mistaken for a prime is at most (1/4)^iter. The default value of iter is 15, which makes the probability less than 1 in 10^9.

Function: primes< start count
Function: primes< start count iter
Function: primes> start count
Function: primes> start count iter
Returns a list of the first count odd probable primes less (more) than or equal to start. The optional parameter iter controls the number of iterations of the Miller-Rabin test for each candidate. The probability of a composite candidate being mistaken for a prime is at most (1/4)^iter. The default value of iter is 15, which makes the probability less than 1 in 10^9.

Theory

Rabin and Miller's result can be summarized as follows. Let p (the candidate prime) be any odd integer greater than 2. Let b (the "base") be an integer in the range 2 ... p-1. There is a fairly simple Boolean function--call it C, for "Composite"---with the following properties:

For details of C, and why it fails for at most 1/4 of the potential bases, please consult a book on number theory or cryptography such as "A Course in Number Theory and Cryptography" by Neal Koblitz, published by Springer-Verlag 1994.

There is nothing probablistic about this result. It's true for all p. If we had time to test (1/4)p + 1 different bases, we could definitively determine the primality of p. For large candidates, that would take much too long--much longer than the simple approach of dividing by all numbers up to sqrt(p). This is where probability enters the picture.

Suppose we have some candidate prime p. Pick a random integer b in the range 2 ... p-1. Compute C(p,b). If p is prime, the result will certainly be false. If p is composite, the probability is at most 1/4 that the result will be false (demonstrating that p is a strong pseudoprime to base b). The test can be repeated with other random bases. If p is prime, each test is certain to return false. If p is composite, the probability of C(p,b) returning false is at most 1/4 for each test. Since the b are chosen at random, the tests outcomes are independent. So if p is composite and the test is repeated, say, 15 times, the probability of it returning false all fifteen times is at most (1/4)^15, or about 10^-9. If the test is repeated 30 times, the probability of failure drops to at most 8.3e-25.

Rabin and Miller's result holds for all candidates p. However, if the candidate p is picked at random, the probability of the Miller-Rabin test failing is much less than the computed bound. This is because, for most composite numbers, the fraction of bases that cause the test to fail is much less than 1/4. For example, if you pick a random odd number less than 1000 and apply the Miller-Rabin test with only 3 random bases, the computed failure bound is (1/4)^3, or about 1.6e-2. However, the actual probability of failure is much less--about 7.2e-5. If you accidentally pick 703 to test for primality, the probability of failure is (161/703)^3, or about 1.2e-2, which is almost as high as the computed bound. This is because 703 is a strong pseudoprime to 161 bases. But if you pick at random there is only a small chance of picking 703, and no other number less than 1000 has that high a percentage of pseudoprime bases.

The Miller-Rabin test is sometimes used in a slightly different fashion, where it can, at least in principle, cause problems. The weaker version uses small prime bases instead of random bases. If you are picking candidates at random and testing for primality, this works well since very few composites are strong pseudo-primes to small prime bases. (For example, there is only one composite less than 2.5e10 that is a strong pseudo-prime to the bases 2, 3, 5, and 7.) The problem with this approach is that once a candidate has been picked, the test is deterministic. This distinction is subtle, but real. With the randomized test, for any candidate you pick--even if your candidate-picking procedure is strongly biased towards troublesome numbers, the test will work with high probability. With the deterministic version, for any particular candidate, the test will either work (with probability 1), or fail (with probability 1). It won't fail for very many candidates, but that won't be much consolation if your candidate-picking procedure is somehow biased toward troublesome numbers.

Prime Factorization

(require 'factor)

Function: factor k
Returns a list of the prime factors of k. The order of the factors is unspecified. In order to obtain a sorted list do (sort! (factor k) <).

Note: The rest of these procedures implement the Solovay-Strassen primality test. This test has been superseeded by the faster See section Prime Testing and Generation. However these are left here as they take up little space and may be of use to an implementation without bignums.

See Robert Solovay and Volker Strassen, A Fast Monte-Carlo Test for Primality, SIAM Journal on Computing, 1977, pp 84-85.

Function: jacobi-symbol p q
Returns the value (+1, -1, or 0) of the Jacobi-Symbol of exact non-negative integer p and exact positive odd integer q.

Function: prime? p
Returns #f if p is composite; #t if p is prime. There is a slight chance (expt 2 (- prime:trials)) that a composite will return #t.

Function: prime:trials
Is the maxinum number of iterations of Solovay-Strassen that will be done to test a number for primality.

Random Numbers

(require 'random)

Procedure: random n
Procedure: random n state
Accepts a positive integer or real n and returns a number of the same type between zero (inclusive) and n (exclusive). The values returned have a uniform distribution.

The optional argument state must be of the type produced by (make-random-state). It defaults to the value of the variable *random-state*. This object is used to maintain the state of the pseudo-random-number generator and is altered as a side effect of the random operation.

Variable: *random-state*
Holds a data structure that encodes the internal state of the random-number generator that random uses by default. The nature of this data structure is implementation-dependent. It may be printed out and successfully read back in, but may or may not function correctly as a random-number state object in another implementation.

Procedure: make-random-state
Procedure: make-random-state state
Returns a new object of type suitable for use as the value of the variable *random-state* and as a second argument to random. If argument state is given, a copy of it is returned. Otherwise a copy of *random-state* is returned.

If inexact numbers are support by the Scheme implementation, `randinex.scm' will be loaded as well. `randinex.scm' contains procedures for generating inexact distributions.

Procedure: random:uniform state
Returns an uniformly distributed inexact real random number in the range between 0 and 1.

Procedure: random:solid-sphere! vect
Procedure: random:solid-sphere! vect state
Fills vect with inexact real random numbers the sum of whose squares is less than 1.0. Thinking of vect as coordinates in space of dimension n = (vector-length vect), the coordinates are uniformly distributed within the unit n-shere. The sum of the squares of the numbers is returned.

Procedure: random:hollow-sphere! vect
Procedure: random:hollow-sphere! vect state
Fills vect with inexact real random numbers the sum of whose squares is equal to 1.0. Thinking of vect as coordinates in space of dimension n = (vector-length vect), the coordinates are uniformly distributed over the surface of the unit n-shere.

Procedure: random:normal
Procedure: random:normal state
Returns an inexact real in a normal distribution with mean 0 and standard deviation 1. For a normal distribution with mean m and standard deviation d use (+ m (* d (random:normal))).

Procedure: random:normal-vector! vect
Procedure: random:normal-vector! vect state
Fills vect with inexact real random numbers which are independent and standard normally distributed (i.e., with mean 0 and variance 1).

Procedure: random:exp
Procedure: random:exp state
Returns an inexact real in an exponential distribution with mean 1. For an exponential distribution with mean u use (* u (random:exp)).

Cyclic Checksum

(require 'make-crc)

Function: make-port-crc
Function: make-port-crc degree
Function: make-port-crc degree generator
Returns an expression for a procedure of one argument, a port. This procedure reads characters from the port until the end of file and returns the integer checksum of the bytes read.

The integer degree, if given, specifies the degree of the polynomial being computed -- which is also the number of bits computed in the checksums. The default value is 32.

The integer generator specifies the polynomial being computed. The power of 2 generating each 1 bit is the exponent of a term of the polynomial. The bit at position degree is implicit and should not be part of generator. This allows systems with numbers limited to 32 bits to calculate 32 bit checksums. The default value of generator when degree is 32 (its default) is:

(make-port-crc 32 #b00000100110000010001110110110111)

Creates a procedure to calculate the P1003.2/D11.2 (POSIX.2) 32-bit checksum from the polynomial:

     32    26    23    22    16    12    11
  ( x   + x   + x   + x   + x   + x   + x   +

      10    8    7    5    4    2    1
     x   + x  + x  + x  + x  + x  + x  + 1 )  mod 2
(require 'make-crc)
(define crc32 (slib:eval (make-port-crc)))
(define (file-check-sum file) (call-with-input-file file crc32))
(file-check-sum (in-vicinity (library-vicinity) "ratize.scm"))

=> 3553047446

Plotting on Character Devices

(require 'charplot)

The plotting procedure is made available through the use of the charplot package. charplot is loaded by inserting (require 'charplot) before the code that uses this procedure.

Variable: charplot:height
The number of rows to make the plot vertically.

Variable: charplot:width
The number of columns to make the plot horizontally.

Procedure: plot! coords x-label y-label
coords is a list of pairs of x and y coordinates. x-label and y-label are strings with which to label the x and y axes.

Example:

(require 'charplot)
(set! charplot:height 19)
(set! charplot:width 45)

(define (make-points n)
  (if (zero? n)
      '()
      (cons (cons (/ n 6) (sin (/ n 6))) (make-points (1- n)))))

(plot! (make-points 37) "x" "Sin(x)")
-|
  Sin(x)   ______________________________________________
      1.25|-                                             |
          |                                              |
         1|-       ****                                  |
          |      **    **                                |
  750.0e-3|-    *        *                               |
          |    *          *                              |
  500.0e-3|-  *            *                             |
          |  *                                           |
  250.0e-3|-                *                            |
          | *                *                           |
         0|-------------------*--------------------------|
          |                                     *        |
 -250.0e-3|-                   *               *         |
          |                     *             *          |
 -500.0e-3|-                     *                       |
          |                       *          *           |
 -750.0e-3|-                       *        *            |
          |                         **    **             |
        -1|-                          ****               |
          |____________:_____._____:_____._____:_________|
        x              2           4      

Root Finding

(require 'root)

Function: newtown:find-integer-root f df/dx x0
Given integer valued procedure f, its derivative (with respect to its argument) df/dx, and initial integer value x0 for which df/dx(x0) is non-zero, returns an integer x for which f(x) is closer to zero than either of the integers adjacent to x; or returns #f if such an integer can't be found.

To find the closest integer to a given integers square root:

(define (integer-sqrt y)
  (newton:find-integer-root
   (lambda (x) (- (* x x) y))
   (lambda (x) (* 2 x))
   (ash 1 (quotient (integer-length y) 2))))

(integer-sqrt 15) => 4

Function: integer-sqrt y
Given a non-negative integer y, returns the rounded square-root of y.

Function: newton:find-root f df/dx x0 prec
Given real valued procedures f, df/dx of one (real) argument, initial real value x0 for which df/dx(x0) is non-zero, and positive real number prec, returns a real x for which abs(f(x)) is less than prec; or returns #f if such a real can't be found.

If prec is instead a negative integer, newton:find-root returns the result of -prec iterations.

H. J. Orchard, The Laguerre Method for Finding the Zeros of Polynomials, IEEE Transactions on Circuits and Systems, Vol. 36, No. 11, November 1989, pp 1377-1381.

There are 2 errors in Orchard's Table II. Line k=2 for starting value of 1000+j0 should have Z_k of 1.0475 + j4.1036 and line k=2 for starting value of 0+j1000 should have Z_k of 1.0988 + j4.0833.

Function: laguerre:find-root f df/dz ddf/dz^2 z0 prec
Given complex valued procedure f of one (complex) argument, its derivative (with respect to its argument) df/dx, its second derivative ddf/dz^2, initial complex value z0, and positive real number prec, returns a complex number z for which magnitude(f(z)) is less than prec; or returns #f if such a number can't be found.

If prec is instead a negative integer, laguerre:find-root returns the result of -prec iterations.

Function: laguerre:find-polynomial-root deg f df/dz ddf/dz^2 z0 prec
Given polynomial procedure f of integer degree deg of one argument, its derivative (with respect to its argument) df/dx, its second derivative ddf/dz^2, initial complex value z0, and positive real number prec, returns a complex number z for which magnitude(f(z)) is less than prec; or returns #f if such a number can't be found.

If prec is instead a negative integer, laguerre:find-polynomial-root returns the result of -prec iterations.


Go to the first, previous, next, last section, table of contents.