5. Mathematical Packages

5.1 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. These functions behave as though operating on integers in two's-complement representation.

5.1.1 Bitwise Operations

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: bitwise-if mask n0 n1

Returns an integer composed of some bits from integer n0 and some from integer n1. A bit of the result is taken from n0 if the corresponding bit of integer mask is 1 and from n1 if that bit of mask is 0.

Function: logtest j k

(logtest j k) == (not (zero? (logand j k))) (logtest #b0100 #b1011) => #f (logtest #b0100 #b0111) => #t

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

5.1.2 Bit Within Word

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: copy-bit index from bit

Returns an integer the same as from except in the indexth bit, which is 1 if bit is #t and 0 if bit is #f.

Example:

(number->string (copy-bit 0 0 #t) 2) => "1" (number->string (copy-bit 2 0 #t) 2) => "100" (number->string (copy-bit 2 #b1111 #f) 2) => "1011"

5.1.3 Fields of Bits

Function: logical:ones n

Returns the smallest non-negative integer having n binary ones.

Function: bit-field 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-field #b1101101010 0 4) 2) => "1010" (number->string (bit-field #b1101101010 4 9) 2) => "10110"

Function: copy-bit-field to start end from

Returns an integer the same as to except possibly in the start (inclusive) through end (exclusive) bits, which are the same as those of from. The 0-th bit of from becomes the startth bit of the result.

Example:

(number->string (copy-bit-field #b1101101010 0 4 0) 2) => "1101100000" (number->string (copy-bit-field #b1101101010 0 4 -1) 2) => "1101101111"

Function: ash n count

Returns an integer equivalent to (inexact->exact (floor (* n (expt 2 count)))).

Example:

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

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

5.1.4 Bit order and Lamination

Function: logical:rotate k count len

Returns the low-order len bits of k cyclically permuted count bits towards high-order.

Example:

(number->string (logical:rotate #b0100 3 4) 2) => "10" (number->string (logical:rotate #b0100 -1 4) 2) => "10"

Function: bit-reverse k n

Returns the low-order k bits of n with the bit order reversed. The low-order bit of n is the high order bit of the returned value.

(number->string (bit-reverse 8 #xa7) 16) => "e5"

Function: integer->list k len

Function: integer->list k

integer->list returns a list of len booleans corresponding to each bit of the given integer. #t is coded for each 1; #f for 0. The len argument defaults to (integer-length k).

Function: list->integer list

list->integer returns an integer formed from the booleans in the list list, which must be a list of booleans. A 1 bit is coded for each #t; a 0 bit for #f.

integer->list and list->integer are inverses so far as equal? is concerned.

Function: booleans->integer bool1 ...

Returns the integer coded by the bool1 ... arguments.

Function: bitwise:laminate k1 ...

Returns an integer composed of the bits of k1 ... interlaced in argument order. Given k1, ... kn, the n low-order bits of the returned value will be the lowest-order bit of each argument.

Function: bitwise:delaminate count k

Returns a list of count integers comprised of every counth bit of the integer k.

For any non-negative integers k and count:

(eqv? k (bitwise:laminate (bitwise:delaminate count k)))

5.1.5 Gray code

A Gray code is an ordering of non-negative integers in which exactly one bit differs between each pair of successive elements. There are multiple Gray codings. An n-bit Gray code corresponds to a Hamiltonian cycle on an n-dimensional hypercube.

Gray codes find use communicating incrementally changing values between asynchronous agents. De-laminated Gray codes comprise the coordinates of Peano-Hilbert space-filling curves.

Function: integer->gray-code k

Converts k to a Gray code of the same integer-length as k.

Function: gray-code->integer k

Converts the Gray code k to an integer of the same integer-length as k.

For any non-negative integer k,

(eqv? k (gray-code->integer (integer->gray-code k)))

Function: = k1 k2

Function: gray-code<? k1 k2

Function: gray-code>? k1 k2

Function: gray-code<=? k1 k2

Function: gray-code>=? k1 k2

These procedures return #t if their Gray code arguments are (respectively): equal, monotonically increasing, monotonically decreasing, monotonically nondecreasing, or monotonically nonincreasing.

For any non-negative integers k1 and k2, the Gray code predicate of (integer->gray-code k1) and (integer->gray-code k2) will return the same value as the corresponding predicate of k1 and k2.

5.2 Modular Arithmetic

(require 'modular)

Function: mod x1 x2

Function: rem x1 x2

These procedures implement the Common-Lisp functions of the same names. The real number x2 must be non-zero. mod returns (- x1 (* x2 (floor (/ x1 x2)))). rem returns (- x1 (* x2 (truncate (/ x1 x2)))).

If x1 and x2 are integers, then mod behaves like modulo and rem behaves like remainder.

(mod -90 360) => 270 (rem -90 180) => -90 (mod 540 360) => 180 (rem 540 360) => 180 (mod (* 5/2 pi) (* 2 pi)) => 1.5707963267948965 (rem (* -5/2 pi) (* 2 pi)) => -1.5707963267948965

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.

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 n2

Returns an integer n such that 1 = (n * n2) mod modulus. If n2 has no inverse mod modulus an error is signaled.

Function: modular:negate modulus n2

Returns (-n2) mod modulus.

Function: modular:+ modulus n2 n3

Returns (n2 + n3) mod modulus.

Function: modular:- modulus n2 n3

Returns (n2 - n3) mod modulus.

Function: modular:* modulus n2 n3

Returns (n2 * n3) mod modulus.

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

Function: modular:expt modulus n2 n3

Returns (n2 ^ n3) mod modulus.

5.3 Prime Numbers

(require 'factor)

Variable: prime:prngs

prime:prngs is the random-state (see section 5.4 Random Numbers) used by these procedures. If you call these procedures from more than one thread (or from interrupt), random may complain about reentrant calls.

• 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.

Variable: prime:trials

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

Function: prime? n

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

Function: primes< start count

Returns a list of the first count prime numbers less than start. If there are fewer than count prime numbers less than start, then the returned list will have fewer than start elements.

Function: primes> start count

Returns a list of the first count prime numbers greater than start.

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) <).

5.4 Random Numbers

A pseudo-random number generator is only as good as the tests it passes. George Marsaglia of Florida State University developed a battery of tests named DIEHARD (http://stat.fsu.edu/~geo/diehard.html). `diehard.c' has a bug which the patch http://swissnet.ai.mit.edu/ftpdir/users/jaffer/diehard.c.pat corrects.

SLIB's PRNG generates 8 bits at a time. With the degenerate seed `0', the numbers generated pass DIEHARD; but when bits are combined from sequential bytes, tests fail. With the seed `http://swissnet.ai.mit.edu/~jaffer/SLIB.html', all of those tests pass.

5.4.1 Exact Random Numbers

(require 'random)

Function: random n state

Function: random n

n must be an exact positive integer. random returns an exact integer between zero (inclusive) and n (exclusive). The values returned by random are uniformly distributed from 0 to n.

The optional argument state must be of the type returned by (seed->random-state) or (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 calls to random.

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.

Function: copy-random-state state

Returns a new copy of argument state.

Function: copy-random-state

Returns a new copy of *random-state*.

Function: seed->random-state seed

Returns a new object of type suitable for use as the value of the variable *random-state* or as a second argument to random. The number or string seed is used to initialize the state. If seed->random-state is called twice with arguments which are equal?, then the returned data structures will be equal?. Calling seed->random-state with unequal arguments will nearly always return unequal states.

Function: make-random-state

Function: make-random-state obj

Returns a new object of type suitable for use as the value of the variable *random-state* or as a second argument to random. If the optional argument obj is given, it should be a printable Scheme object; the first 50 characters of its printed representation will be used as the seed. Otherwise the value of *random-state* is used as the seed.

5.4.2 Inexact Random Numbers

(require 'random-inexact)

Function: random:uniform

Function: random:uniform state

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

Function: random:exp

Function: 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)).

Function: random:normal

Function: 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: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: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.

5.5 Fast Fourier Transform

(require 'fft)

Function: fft array

array is an array of (expt 2 n) numbers. fft returns an array of complex numbers comprising the Discrete Fourier Transform of array.

Function: fft-1 array

fft-1 returns an array of complex numbers comprising the inverse Discrete Fourier Transform of array.

(fft '#(1 0+i -1 0-i 1 0+i -1 0-i)) =>

#(0.0 0.0 0.0+628.0783185208527e-18i 0.0
  0.0 0.0 8.0-628.0783185208527e-18i 0.0)

(fft-1 '#(0 0 0 0 0 0 8 0)) =>

#(1.0 -61.23031769111886e-18+1.0i -1.0 61.23031769111886e-18-1.0i
  1.0 -61.23031769111886e-18+1.0i -1.0 61.23031769111886e-18-1.0i)

5.6 Cyclic Checksum

(require 'crc) Cyclic Redundancy Checks using Galois field GF(2) polynomial arithmetic are used for error detection in many data transmission and storage applications.

The generator polynomials for various CRC protocols are availble from many sources. But the polynomial is just one of many parameters which must match in order for a CRC implementation to interoperate with existing systems:

• the byte-order and bit-order of the data stream;

• whether the CRC or its inverse is being calculated;

• the initial CRC value; and

• whether and where the CRC value is appended (inverted or non-inverted) to the data stream.

There is even some controversy over the polynomials themselves.

Constant: crc-32-polynomial

For CRC-32, http://www2.sis.pitt.edu/~jkabara/tele-2100/lect08.html gives x^32+x^26+x^23+x^16+x^12+x^11+x^10+x^8+x^7+x^5+x^4+x^2+x^1+1.

But http://www.cs.ncl.ac.uk/people/harry.whitfield/home.formal/CRCs.html, http://duchon.umuc.edu/Web_Pages/duchon/99_f_cm435/ShiftRegister.htm, http://spinroot.com/spin/Doc/Book91_PDF/ch3.pdf, http://www.erg.abdn.ac.uk/users/gorry/course/dl-pages/crc.html, http://www.rad.com/networks/1994/err_con/crc_most.htm, and http://www.gpfn.sk.ca/~rhg/csc8550s02/crc.html, http://www.nobugconsulting.ro/crc.php give x^32+x^26+x^23+x^22+x^16+x^12+x^11+x^10+x^8+x^7+x^5+x^4+x^2+x+1.

SLIB crc-32-polynomial uses the latter definition.

Constant: crc-ccitt-polynomial

http://www.math.grin.edu/~rebelsky/Courses/CS364/2000S/Outlines/outline.12.html, http://duchon.umuc.edu/Web_Pages/duchon/99_f_cm435/ShiftRegister.htm, http://www.cs.ncl.ac.uk/people/harry.whitfield/home.formal/CRCs.html, http://www2.sis.pitt.edu/~jkabara/tele-2100/lect08.html, and http://www.gpfn.sk.ca/~rhg/csc8550s02/crc.html give CRC-CCITT: x^16+x^12+x^5+1.

Constant: crc-16-polynomial

http://www.math.grin.edu/~rebelsky/Courses/CS364/2000S/Outlines/outline.12.html, http://duchon.umuc.edu/Web_Pages/duchon/99_f_cm435/ShiftRegister.htm, http://www.cs.ncl.ac.uk/people/harry.whitfield/home.formal/CRCs.html, http://www.gpfn.sk.ca/~rhg/csc8550s02/crc.html, and http://www.usb.org/developers/data/crcdes.pdf give CRC-16: x^16+x^15+x^2+1.

Constant: crc-12-polynomial

http://www.math.grin.edu/~rebelsky/Courses/CS364/2000S/Outlines/outline.12.html, http://www.cs.ncl.ac.uk/people/harry.whitfield/home.formal/CRCs.html, http://www.it.iitb.ac.in/it605/lectures/link/node4.html, and http://spinroot.com/spin/Doc/Book91_PDF/ch3.pdf give CRC-12: x^12+x^11+x^3+x^2+1.

But http://www.ffldusoe.edu/Faculty/Denenberg/Topics/Networks/Error_Detection_Correction/crc.html, http://duchon.umuc.edu/Web_Pages/duchon/99_f_cm435/ShiftRegister.htm, http://www.eng.uwi.tt/depts/elec/staff/kimal/errorcc.html, http://www.ee.uwa.edu.au/~roberto/teach/itc314/java/CRC/, http://www.gpfn.sk.ca/~rhg/csc8550s02/crc.html, and http://www.efg2.com/Lab/Mathematics/CRC.htm give CRC-12: x^12+x^11+x^3+x^2+x+1.

These differ in bit 1 and calculations using them return different values. With citations near evenly split, it is hard to know which is correct.

Constant: crc-10-polynomial

http://www.math.grin.edu/~rebelsky/Courses/CS364/2000S/Outlines/outline.12.html gives CRC-10: x^10+x^9+x^5+x^4+1; but http://cell-relay.indiana.edu/cell-relay/publications/software/CRC/crc10.html, http://www.it.iitb.ac.in/it605/lectures/link/node4.html, http://www.gpfn.sk.ca/~rhg/csc8550s02/crc.html, http://www.techfest.com/networking/atm/atm.htm, http://www.protocols.com/pbook/atmcell2.htm, and http://www.nobugconsulting.ro/crc.php give CRC-10: x^10+x^9+x^5+x^4+x+1.

Constant: crc-08-polynomial

http://www.math.grin.edu/~rebelsky/Courses/CS364/2000S/Outlines/outline.12.html, http://www.cs.ncl.ac.uk/people/harry.whitfield/home.formal/CRCs.html, http://www.it.iitb.ac.in/it605/lectures/link/node4.html, and http://www.nobugconsulting.ro/crc.php give CRC-8: x^8+x^2+x^1+1

Constant: atm-hec-polynomial

http://cell-relay.indiana.edu/cell-relay/publications/software/CRC/32bitCRC.tutorial.html and http://www.gpfn.sk.ca/~rhg/csc8550s02/crc.html give ATM HEC: x^8+x^2+x+1.

Constant: dowcrc-polynomial

http://www.cs.ncl.ac.uk/people/harry.whitfield/home.formal/CRCs.html gives DOWCRC: x^8+x^5+x^4+1.

Constant: usb-token-polynomial

http://www.usb.org/developers/data/crcdes.pdf and http://www.nobugconsulting.ro/crc.php give USB-token: x^5+x^2+1.

Each of these polynomial constants is a string of `1's and `0's, the exponent of each power of x in descending order.

Function: crc:make-table poly

poly must be string of `1's and `0's beginning with `1' and having length greater than 8. crc:make-table returns a vector of 256 integers, such that:

(set! crc (logxor (ash (logand (+ -1 (ash 1 (- deg 8))) crc) 8) (vector-ref crc-table (logxor (ash crc (- 8 deg)) byte))))

will compute the crc with the 8 additional bits in byte; where crc is the previous accumulated CRC value, deg is the degree of poly, and crc-table is the vector returned by crc:make-table.

If the implementation does not support deg-bit integers, then crc:make-table returns #f.

Function: cksum file

Computes the P1003.2/D11.2 (POSIX.2) 32-bit checksum of file.

(require 'crc) (cksum (in-vicinity (library-vicinity) "ratize.scm")) => 157103930

Function: cksum port

Computes the checksum of the bytes read from port until the end-of-file.

cksum-string, which returns the P1003.2/D11.2 (POSIX.2) 32-bit checksum of the bytes in str, can be defined as follows:

(require 'string-port)
(define (cksum-string str) (call-with-input-string str cksum))

Function: crc16 file

Computes the USB data-packet (16-bit) CRC of file.

Function: crc16 port

Computes the USB data-packet (16-bit) CRC of the bytes read from port until the end-of-file.

crc16 calculates the same values as the crc16.pl program given in http://www.usb.org/developers/data/crcdes.pdf.

Function: crc5 file

Computes the USB token (5-bit) CRC of file.

Function: crc5 port

Computes the USB token (5-bit) CRC of the bytes read from port until the end-of-file.

crc5 calculates the same values as the crc5.pl program given in http://www.usb.org/developers/data/crcdes.pdf.

5.7 Graphing

5.7.1 Character Plotting

(require 'charplot)

Variable: charplot:dimensions

A list of the maximum height (number of lines) and maximum width (number of columns) for the graph, its scales, and labels.

The default value for charplot:dimensions is the output-port-height and output-port-width of current-output-port.

Procedure: plot coords x-label y-label

coords is a list or vector of coordinates, lists 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:dimensions '(20 55)) (define (make-points n) (if (zero? n) '() (cons (list (/ n 6) (sin (/ n 6))) (make-points (1- n))))) (plot (make-points 40) "x" "Sin(x)") -| Sin(x) _________________________________________ 1|- **** | | ** ** | 0.75|- * * | | * * | 0.5|- * * | | * *| 0.25|- * * | | * * | 0|-------------------*------------------*--| | * | -0.25|- * * | | * * | -0.5|- * | | * * | -0.75|- * * | | ** ** | -1|- **** | |:_____._____:_____._____:_____._____:____| x 2 4 6

Procedure: plot func x1 x2

Procedure: plot func x1 x2 npts

Plots the function of one argument func over the range x1 to x2. If the optional integer argument npts is supplied, it specifies the number of points to evaluate func at.

(plot sin 0 (* 2 pi)) -| _________________________________________ 1|-: **** | | : ** ** | 0.75|-: * * | | : * * | 0.5|-: ** ** | | : * * | 0.25|-:** ** | | :* * | 0|-*------------------*--------------------| | : * * | -0.25|-: ** ** | | : * * | -0.5|-: * ** | | : * * | -0.75|-: * ** | | : ** ** | -1|-: **** | |_:_____._____:_____._____:_____._____:___| 0 2 4 6

Procedure: histograph data label

Creates and displays a histogram of the numerical values contained in vector or list data

(require 'random-inexact) (histograph (do ((idx 99 (+ -1 idx)) (lst '() (cons (* .02 (random:normal)) lst))) ((negative? idx) lst)) "normal") -| _________________________________________ 8|- : I | | : I | 7|- I I : I | | I I : I | 6|- III I :I I | | III I :I I | 5|- IIIIIIIIII I | | IIIIIIIIII I | 4|- IIIIIIIIIIII | | IIIIIIIIIIII | 3|-I I I IIIIIIIIIIII II I | | I I I IIIIIIIIIIII II I | 2|-I I I IIIIIIIIIIIIIIIII I | | I I I IIIIIIIIIIIIIIIII I | 1|-II I I IIIIIIIIIIIIIIIIIIIII I I I | | II I I IIIIIIIIIIIIIIIIIIIII I I I | 0|-IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII----| |__.____:____.____:____.____:____.____:___| normal -0.025 0 0.025 0.05

5.7.2 PostScript Graphing

(require 'eps-graph)

This is a graphing package creating encapsulated-PostScript files. Its motivations and design choice are described in http://swissnet.ai.mit.edu/~jaffer/Docupage/grapheps

A dataset to be plotted is taken from a 2-dimensional array. Corresponding coordinates are in rows. Coordinates from any pair of columns can be plotted.

Function: create-postscript-graph filename.eps size elt1 ...

filename.eps should be a string naming an output file to be created. size should be an exact integer, a list of two exact integers, or #f. elt1, ... are values returned by graphing primitives described here.

create-postscript-graph creates an Encapsulated-PostScript file named filename.eps containing graphs as directed by the elt1, ... arguments.

The size of the graph is determined by the size argument. If a list of two integers, they specify the width and height. If one integer, then that integer is the width and the height is 3/4 of the width. If #f, the graph will be 800 by 600.

Function: whole-page

Pushes a rectangle for the whole encapsulated page onto the PostScript stack. This pushed rectangle is an implicit argument to partition-page or setup-plot.

5.7.2.1 Column Ranges

A range is a list of two numbers, the minimum and the maximum. Ranges can be given explicity or computed in PostScript by column-range.

Function: column-range array k

Returns the range of values in 2-dimensional array column k.

Function: pad-range range p

Expands range by p/100 on each end.

Function: snap-range range

Expands range to round number of ticks.

Function: combine-ranges range1 range2 ...

Returns the minimal range covering all range1, range2, ...

Function: setup-plot x-range y-range pagerect

Function: setup-plot x-range y-range

x-range and y-range should each be a list of two numbers or the value returned by pad-range, snap-range, or combine-range. pagerect is the rectangle bounding the graph to be drawn; if missing, the rectangle from the top of the PostScript stack is popped and used.

Based on the given ranges, setup-plot sets up scaling and margins for making a graph. The margins are sized proportional to the fontheight value at the time of the call to setup-plot. setup-plot sets two variables:

plotrect

The region where data points will be plotted.

graphrect

The pagerect argument to setup-plot. Includes plotrect, legends, etc.

5.7.2.2 Drawing the Graph

Function: plot-column array x-column y-column proc3s

Plots points with x coordinate in x-column of array and y coordinate y-column of array. The symbol proc3s specifies the type of glyph or drawing style for presenting these coordinates.

line

Draws line connecting points in order.

mountain

Fill area below line connecting points.

cloud

Fill area above line connecting points.

impulse

Draw line from x-axis to each point.

bargraph

Draw rectangle from x-axis to each point.

disc

Solid round dot.

point

Minimal point -- invisible if linewidth is 0.

square

Square box.

diamond

Square box at 45.o

plus

Plus sign.

cross

X sign.

triup

Triangle pointing upward

tridown

Triangle pointing downward

pentagon

Five sided polygon

circle

Hollow circle

5.7.2.3 Graphics Context

Function: in-graphic-context arg ...

Saves the current graphics state, executes args, then restores to saved graphics state.

Function: set-color color

color should be a string naming a Resene color, a saturate color, or a number between 0 and 100.

set-color sets the PostScript color to the color of the given string, or a grey value between black (0) and white (100).

Function: set-font name fontheight

name should be a (case-sensitive) string naming a PostScript font. fontheight should be a positive real number.

set-font Changes the current PostScript font to name with height equal to fontheight. The default font is Helvetica (12pt).

Line parameters do no affect fonts; they do effect glyphs.

Function: set-linewidth w

The default linewidth is 1. Setting it to 0 makes the lines drawn as skinny as possible. Linewidth must be much smaller than glyphsize for readable glyphs.

Function: set-linedash j k

Lines are drawn j-on k-off.

Function: set-linedash j

Lines are drawn j-on j-off.

Function: set-linedash

Turns off dashing.

Function: set-glyphsize w

Sets the (PostScript) variable glyphsize to w. The default glyphsize is 6.

5.7.2.4 Rectangles

A rectangle is a list of 4 numbers; the first two elements are the x and y coordinates of lower left corner of the rectangle. The other two elements are the width and height of the rectangle.

Function: whole-page

Pushes a rectangle for the whole encapsulated page onto the PostScript stack. This pushed rectangle is an implicit argument to partition-page or setup-plot.

Function: partition-page xparts yparts

Pops the rectangle currently on top of the stack and pushes xparts * yparts sub-rectangles onto the stack in decreasing y and increasing x order. If you are drawing just one graph, then you don't need partition-page.

Variable: plotrect

The rectangle where data points should be plotted. plotrect is set by setup-plot.

Variable: graphrect

The pagerect argument of the most recent call to setup-plot. Includes plotrect, legends, etc.

Function: fill-rect rect

fills rect with the current color.

Function: outline-rect rect

Draws the perimiter of rect in the current color.

Function: clip-to-rect rect

Modifies the current graphics-state so that nothing will be drawn outside of the rectangle rect. Use in-graphic-context to limit the extent of clip-to-rect.

5.7.2.5 Legending

Function: title-top title subtitle

Function: title-top title

Puts a title line and an optional subtitle line above the graphrect.

Function: title-bottom title subtitle

Function: title-bottom title

Puts a title line and an optional subtitle line below the graphrect.

Variable: topedge

Variable: bottomedge

These edge coordinates of graphrect are suitable for passing as the first argument to rule-horizontal.

Variable: leftedge

Variable: rightedge

These edge coordinates of graphrect are suitable for passing as the first argument to rule-vertical.

Function: rule-vertical x-coord text tick-width

Draws a vertical ruler with X coordinate x-coord and labeled with string text. If tick-width is positive, then the ticks are tick-width long on the right side of x-coord; and text and numeric legends are on the left. If tick-width is negative, then the ticks are -tick-width long on the left side of x-coord; and text and numeric legends are on the right.

Function: rule-horizontal x-coord text tick-height

Draws a horizontal ruler with X coordinate x-coord and labeled with string text. If tick-height is positive, then the ticks are tick-height long on the right side of x-coord; and text and numeric legends are on the left. If tick-height is negative, then the ticks are -tick-height long on the left side of x-coord; and text and numeric legends are on the right.

Function: y-axis

Draws the y-axis.

Function: x-axis

Draws the x-axis.

Function: grid-verticals

Draws vertical lines through graphrect at each tick on the vertical ruler.

Function: grid-horizontals

Draws horizontal lines through graphrect at each tick on the horizontal ruler.

5.7.2.6 Legacy Plotting

Variable: graph:dimensions

A list of the width and height of the graph to be plotted using plot.

Function: plot func x1 x2 npts

Function: plot func x1 x2

Creates and displays using (system "gv tmp.eps") an encapsulated PostScript graph of the function of one argument func over the range x1 to x2. If the optional integer argument npts is supplied, it specifies the number of points to evaluate func at.

Function: plot coords x-label y-label

coords is a list or vector of coordinates, lists of x and y coordinates. x-label and y-label are strings with which to label the x and y axes.

5.7.2.7 Example Graph

The file `am1.5.html', a table of solar irradiance, is fetched with `wget' if it isn't already in the working directory. The file is read and stored into an array, irradiance.

create-postscript-graph is then called to create an encapsulated-PostScript file, `solarad.eps'. The size of the page is set to 600 by 300. whole-page is called and leaves the rectangle on the PostScript stack. setup-plot is called with a literal range for x and computes the range for column 1.

Two calls to top-title are made so a different font can be used for the lower half. in-graphic-context is used to limit the scope of the font change. The graphing area is outlined and a rule drawn on the left side.

Because the X range was intentionally reduced, in-graphic-context is called and clip-to-rect limits drawing to the plotting area. A black line is drawn from data column 1. That line is then overlayed with a mountain plot of the same column colored "Bright Sun".

After returning from the in-graphic-context, the bottom ruler is drawn. Had it been drawn earlier, all its ticks would have been painted over by the mountain plot.

The color is then changed to `seagreen' and the same graphrect is setup again, this time with a different Y scale, 0 to 1000. The graphic context is again clipped to plotrect, linedash is set, and column 2 is plotted as a dashed line. Finally the rightedge is ruled. Having the line and its scale both in green helps disambiguate the scales.

(require 'eps-graph)
(require 'line-i/o)
(require 'string-port)

(define irradiance
  (let ((url "http://www.pv.unsw.edu.au/am1.5.html")
        (file "am1.5.html"))
    (define (read->list line)
      (define elts '())
      (call-with-input-string line
        (lambda (iprt) (do ((elt (read iprt) (read iprt)))
                           ((eof-object? elt) elts)
                         (set! elts (cons elt elts))))))
    (if (not (file-exists? file))
        (system (string-append "wget -c -O" file " " url)))
    (call-with-input-file file
      (lambda (iprt)
        (define lines '())
        (do ((line (read-line iprt) (read-line iprt)))
            ((eof-object? line)
             (let ((nra (create-array (Ar64)
                                      (length lines)
                                      (length (car lines)))))
               (do ((lns lines (cdr lns))
                    (idx (+ -1 (length lines)) (+ -1 idx)))
                   ((null? lns) nra)
                 (do ((kdx (+ -1 (length (car lines))) (+ -1 kdx))
                      (lst (car lns) (cdr lst)))
                     ((null? lst))
                   (array-set! nra (car lst) idx kdx)))))
          (if (and (positive? (string-length line))
                   (char-numeric? (string-ref line 0)))
              (set! lines (cons (read->list line) lines))))))))

(let ((xrange '(.25 2.5)))
  (create-postscript-graph
   "solarad.eps" '(600 300)
   (whole-page)
   (setup-plot xrange (column-range irradiance 1))
   (title-top
    "Solar Irradiance   http://www.pv.unsw.edu.au/am1.5.html")
   (in-graphic-context
    (set-font "Helvetica-Oblique" 12)
    (title-top
     ""
     "Key Centre for Photovoltaic Engineering UNSW - Air Mass 1.5 Global Spectrum"))
   (outline-rect plotrect)
   (rule-vertical leftedge "W/(m^2.um)" 10)
   (in-graphic-context (clip-to-rect plotrect)
                       (plot-column irradiance 0 1 'line)
                       (set-color "Bright Sun")
                       (plot-column irradiance 0 1 'mountain)
                       )
   (rule-horizontal bottomedge "Wavelength in .um" 5)
   (set-color 'seagreen)

   (setup-plot xrange '(0 1000) graphrect)
   (in-graphic-context (clip-to-rect plotrect)
                       (set-linedash 5 2)
                       (plot-column irradiance 0 2 'line))
   (rule-vertical rightedge "Integrated .W/(m^2)" -10)
   ))

(system "gv solarad.eps")

5.8 Solid Modeling

<A NAME="Solid"> (require 'solid) </A>

http://swissnet.ai.mit.edu/~jaffer/Solid/#Example gives an example use of this package.

Function: vrml node ...

Returns the VRML97 string (including header) of the concatenation of strings nodes, ....

Function: vrml-append node1 node2 ...

Returns the concatenation with interdigitated newlines of strings node1, node2, ....

Function: vrml-to-file file node ...

Writes to file named file the VRML97 string (including header) of the concatenation of strings nodes, ....

Function: world:info title info ...

Returns a VRML97 string setting the title of the file in which it appears to title. Additional strings info, ... are comments.

VRML97 strings passed to vrml and vrml-to-file as arguments will appear in the resulting VRML code. This string turns off the headlight at the viewpoint:

" NavigationInfo {headlight FALSE}"

Function: scene:panorama front right back left top bottom

Specifies the distant images on the inside faces of the cube enclosing the virtual world.

Function: scene:sphere colors angles

colors is a list of color objects. Each may be of type color, a 24-bit sRGB integer, or a list of 3 numbers between 0.0 and 1.0.

angles is a list of non-increasing angles the same length as colors. Each angle is between 90 and -90 degrees. If 90 or -90 are not elements of angles, then the color at the zenith and nadir are taken from the colors paired with the angles nearest them.

scene:sphere fills horizontal bands with interpolated colors on the backgroud sphere encasing the world.

Function: scene:sky-and-dirt

Returns a blue and brown backgroud sphere encasing the world.

Function: scene:sky-and-grass

Returns a blue and green backgroud sphere encasing the world.

Function: scene:sun latitude julian-day hour turbidity strength

Function: scene:sun latitude julian-day hour turbidity

latitude is the virtual place's latitude in degrees. julian-day is an integer from 0 to 366, the day of the year. hour is a real number from 0 to 24 for the time of day; 12 is noon. turbidity is the degree of fogginess described in See section turbidity.

scene:sun returns a bright yellow, distant sphere where the sun would be at hour on julian-day at latitude. If strength is positive, included is a light source of strength (default 1).

Function: scene:overcast latitude julian-day hour turbidity strength

Function: scene:overcast latitude julian-day hour turbidity

latitude is the virtual place's latitude in degrees. julian-day is an integer from 0 to 366, the day of the year. hour is a real number from 0 to 24 for the time of day; 12 is noon. turbidity is the degree of cloudiness described in See section turbidity.

scene:overcast returns an overcast sky as it might look at hour on julian-day at latitude. If strength is positive, included is an ambient light source of strength (default 1).

Function: scene:viewpoint name distance compass pitch

Function: scene:viewpoint name distance compass

Returns a viewpoint named name facing the origin and placed distance from it. compass is a number from 0 to 360 giving the compass heading. pitch is a number from -90 to 90, defaulting to 0, specifying the angle from the horizontal.

Function: scene:viewpoints proximity

Returns 6 viewpoints, one at the center of each face of a cube with sides 2 * proximity, centered on the origin.

Light Sources

In VRML97, lights shine only on objects within the same children node and descendants of that node. Although it would have been convenient to let light direction be rotated by solid:rotation, this restricts a rotated light's visibility to objects rotated with it.

To workaround this limitation, these directional light source procedures accept either Cartesian or spherical coordinates for direction. A spherical coordinate is a list (theta azimuth); where theta is the angle in degrees from the zenith, and azimuth is the angle in degrees due west of south.

It is sometimes useful for light sources to be brighter than `1'. When intensity arguments are greater than 1, these functions gang multiple sources to reach the desired strength.

Function: light:ambient color intensity

Function: light:ambient color

Ambient light shines on all surfaces with which it is grouped.

color is a an object of type color, a 24-bit sRGB integer, or a list of 3 numbers between 0.0 and 1.0. If color is #f, then the default color will be used. intensity is a real non-negative number defaulting to `1'.

light:ambient returns a light source or sources of color with total strength of intensity (or 1 if omitted).

Function: light:directional color direction intensity

Function: light:directional color direction

Function: light:directional color

Directional light shines parallel rays with uniform intensity on all objects with which it is grouped.

color is a an object of type color, a 24-bit sRGB integer, or a list of 3 numbers between 0.0 and 1.0. If color is #f, then the default color will be used.

direction must be a list or vector of 2 or 3 numbers specifying the direction to this light. If direction has 2 numbers, then these numbers are the angle from zenith and the azimuth in degrees; if direction has 3 numbers, then these are taken as a Cartesian vector specifying the direction to the light source. The default direction is upwards; thus its light will shine down.

intensity is a real non-negative number defaulting to `1'.

light:directional returns a light source or sources of color with total strength of intensity, shining from direction.

Function: light:beam attenuation radius aperture peak

Function: light:beam attenuation radius aperture

Function: light:beam attenuation radius

Function: light:beam attenuation

attenuation is a list or vector of three nonnegative real numbers specifying the reduction of intensity, the reduction of intensity with distance, and the reduction of intensity as the square of distance. radius is the distance beyond which the light does not shine. radius defaults to `100'.

aperture is a real number between 0 and 180, the angle centered on the light's axis through which it sheds some light. peak is a real number between 0 and 90, the angle of greatest illumination.

Function: light:point location color intensity beam

Function: light:point location color intensity

Function: light:point location color

Function: light:point location

Point light radiates from location, intensity decreasing with distance, towards all objects with which it is grouped.

color is a an object of type color, a 24-bit sRGB integer, or a list of 3 numbers between 0.0 and 1.0. If color is #f, then the default color will be used. intensity is a real non-negative number defaulting to `1'. beam is a structure returned by light:beam or #f.

light:point returns a light source or sources at location of color with total strength intensity and beam properties. Note that the pointlight itself is not visible. To make it so, place an object with emissive appearance at location.

Function: light:spot location direction color intensity beam

Function: light:spot location direction color intensity

Function: light:spot location direction color

Function: light:spot location direction

Function: light:spot location

Spot light radiates from location towards direction, intensity decreasing with distance, illuminating objects with which it is grouped.

direction must be a list or vector of 2 or 3 numbers specifying the direction to this light. If direction has 2 numbers, then these numbers are the angle from zenith and the azimuth in degrees; if direction has 3 numbers, then these are taken as a Cartesian vector specifying the direction to the light source. The default direction is upwards; thus its light will shine down.

color is a an object of type color, a 24-bit sRGB integer, or a list of 3 numbers between 0.0 and 1.0. If color is #f, then the default color will be used.

intensity is a real non-negative number defaulting to `1'.

light:spot returns a light source or sources at location of direction with total strength color. Note that the spotlight itself is not visible. To make it so, place an object with emissive appearance at location.

Object Primitives

Function: solid:box geometry appearance

Function: solid:box geometry

geometry must be a number or a list or vector of three numbers. If geometry is a number, the solid:box returns a cube with sides of length geometry centered on the origin. Otherwise, solid:box returns a rectangular box with dimensions geometry centered on the origin. appearance determines the surface properties of the returned object.

Function: solid:cylinder radius height appearance

Function: solid:cylinder radius height

Returns a right cylinder with dimensions radius and (abs height) centered on the origin. If height is positive, then the cylinder ends will be capped. appearance determines the surface properties of the returned object.

Function: solid:disk radius thickness appearance

Function: solid:disk radius thickness

thickness must be a positive real number. solid:disk returns a circular disk with dimensions radius and thickness centered on the origin. appearance determines the surface properties of the returned object.

Function: solid:cone radius height appearance

Function: solid:cone radius height

Returns an isosceles cone with dimensions radius and height centered on the origin. appearance determines the surface properties of the returned object.

Function: solid:pyramid side height appearance

Function: solid:pyramid side height

Returns an isosceles pyramid with dimensions side and height centered on the origin. appearance determines the surface properties of the returned object.

Function: solid:sphere radius appearance

Function: solid:sphere radius

Returns a sphere of radius radius centered on the origin. appearance determines the surface properties of the returned object.

Function: solid:ellipsoid geometry appearance

Function: solid:ellipsoid geometry

geometry must be a number or a list or vector of three numbers. If geometry is a number, the solid:ellipsoid returns a sphere of diameter geometry centered on the origin. Otherwise, solid:ellipsoid returns an ellipsoid with diameters geometry centered on the origin. appearance determines the surface properties of the returned object.

Function: solid:basrelief width height depth colorray appearance

Function: solid:basrelief width height depth appearance

Function: solid:basrelief width height depth

One of width, height, or depth must be a 2-dimensional array; the others must be real numbers giving the length of the basrelief in those dimensions. The rest of this description assumes that height is an array of heights.

solid:basrelief returns a width by depth basrelief solid with heights per array height with the buttom surface centered on the origin.

If present, appearance determines the surface properties of the returned object. If present, colorray must be an array of objects of type color, 24-bit sRGB integers or lists of 3 numbers between 0.0 and 1.0.

If colorray's dimensions match height, then each element of colorray paints its corresponding vertex of height. If colorray has all dimensions one smaller than height, then each element of colorray paints the corresponding face of height. Other dimensions for colorray are in error.

Surface Attributes

Function: solid:color diffuseColor ambientIntensity specularColor shininess emissiveColor transparency

Function: solid:color diffuseColor ambientIntensity specularColor shininess emissiveColor

Function: solid:color diffuseColor ambientIntensity specularColor shininess

Function: solid:color diffuseColor ambientIntensity specularColor

Function: solid:color diffuseColor ambientIntensity

Function: solid:color diffuseColor

Returns an appearance, the optical properties of the objects with which it is associated. ambientIntensity, shininess, and transparency must be numbers between 0 and 1. diffuseColor, specularColor, and emissiveColor are objects of type color, 24-bit sRGB integers or lists of 3 numbers between 0.0 and 1.0. If a color argument is omitted or #f, then the default color will be used.

Function: solid:texture image color scale rotation center translation

Function: solid:texture image color scale rotation center

Function: solid:texture image color scale rotation

Function: solid:texture image color scale

Function: solid:texture image color

Function: solid:texture image

Returns an appearance, the optical properties of the objects with which it is associated. image is a string naming a JPEG or PNG image resource. color is #f, a color, or the string returned by solid:color. The rest of the optional arguments specify 2-dimensional transforms applying to the image.

scale must be #f, a number, or list or vector of 2 numbers specifying the scale to apply to image. rotation must be #f or the number of degrees to rotate image. center must be #f or a list or vector of 2 numbers specifying the center of image relative to the image dimensions. translation must be #f or a list or vector of 2 numbers specifying the translation to apply to image.

Aggregating Objects

Function: solid:center-row-of number solid spacing

Returns a row of number solid objects spaced evenly spacing apart.

Function: solid:center-array-of number-a number-b solid spacing-a spacing-b

Returns number-b rows, spacing-b apart, of number-a solid objects spacing-a apart.

Function: solid:center-pile-of number-a number-b number-c solid spacing-a spacing-b spacing-c

Returns number-c planes, spacing-c apart, of number-b rows, spacing-b apart, of number-a solid objects spacing-a apart.

Function: solid:arrow center

center must be a list or vector of three numbers. Returns an upward pointing metallic arrow centered at center.

Function: solid:arrow

Returns an upward pointing metallic arrow centered at the origin.

Spatial Transformations

Function: solid:translation center solid ...

center must be a list or vector of three numbers. solid:translation Returns an aggregate of solids, ... with their origin moved to center.

Function: solid:scale scale solid ...

scale must be a number or a list or vector of three numbers. solid:scale Returns an aggregate of solids, ... scaled per scale.

Function: solid:rotation axis angle solid ...

axis must be a list or vector of three numbers. solid:rotation Returns an aggregate of solids, ... rotated angle degrees around the axis axis.

5.9 Color

<A NAME="Color"></A>

http://swissnet.ai.mit.edu/~jaffer/Color

The goals of this package are to provide methods to specify, compute, and transform colors in a core set of additive color spaces. The color spaces supported should be sufficient for working with the color data encountered in practice and the literature.

5.9.1 Color Data-Type

<A NAME="Color_Data-Type"></A> (require 'color)

Function: color? obj

Returns #t if obj is a color.

Function: color? obj typ

Returns #t if obj is a color of color-space typ. The symbol typ must be one of:

• CIEXYZ
• RGB709
• L*a*b*
• L*u*v*
• sRGB
• e-sRGB
• L*C*h

Function: make-color space arg ...

Returns a color of type space.

Function: color-space color

Returns the symbol for the color-space in which color is embedded.

Function: color-precision color

For colors in digital color-spaces, color-precision returns the number of bits used for each of the R, G, and B channels of the encoding. Otherwise, color-precision returns #f

Function: color-white-point color

Returns the white-point of color in all color-spaces except CIEXYZ.

Function: convert-color color space white-point

Function: convert-color color space

Function: convert-color color e-sRGB precision

Converts color into space at optional white-point.

5.9.1.1 External Representation

Each color encoding has an external, case-insensitive representation. To ensure portability, the white-point for all color strings is D65. (4)

The X, Y, Z, L, a, b, u, v, C, h, R, G, and B fields are (Scheme) real numbers within the appropriate ranges.

The R, G, and B, fields are non-negative exact decimal integers within the appropriate ranges.

Several additional syntaxes are supported by string->color:

Where RRGGBB is a non-negative six-digit hexadecimal number.

Function: color->string color

Returns a string representation of color.

Function: string->color string

Returns the color represented by string. If string is not a syntactically valid notation for a color, then string->color returns #f.

5.9.1.2 White

We experience color relative to the illumination around us. CIEXYZ coordinates, although subject to uniform scaling, are objective. Thus other color spaces are specified relative to a white point in CIEXYZ coordinates.

The white point for digital color spaces is set to D65. For the other spaces a white-point argument can be specified. The default if none is specified is the white-point with which the color was created or last converted; and D65 if none has been specified.

Constant: D65

Is the color of 6500.K (blackbody) illumination. D65 is close to the average color of daylight.

Constant: D50

Is the color of 5000.K (blackbody) illumination. D50 is the color of indoor lighting by incandescent bulbs, whose filaments have temperatures around 5000.K.

5.9.2 Color Spaces

<A NAME="Color_Spaces"></A>

Measurement-based Color Spaces

The tristimulus color spaces are those whose component values are proportional measurements of light intensity. The CIEXYZ(1931) system provides 3 sets of spectra to convolve with a spectrum of interest. The result of those convolutions is coordinates in CIEXYZ space. All tristimuls color spaces are related to CIEXYZ by linear transforms, namely matrix multiplication. Of the color spaces listed here, CIEXYZ and RGB709 are tristimulus spaces.

Color Space: CIEXYZ

The CIEXYZ color space covers the full gamut. It is the basis for color-space conversions.

CIEXYZ is a list of three inexact numbers between 0 and 1.1. '(0. 0. 0.) is black; '(1. 1. 1.) is white.

Function: ciexyz->color xyz

xyz must be a list of 3 numbers. If xyz is valid CIEXYZ coordinates, then ciexyz->color returns the color specified by xyz; otherwise returns #f.

Function: color:ciexyz x y z

Returns the CIEXYZ color composed of x, y, z. If the coordinates do not encode a valid CIEXYZ color, then an error is signaled.

Function: color->ciexyz color

Returns the list of 3 numbers encoding color in CIEXYZ.

Color Space: RGB709

BT.709-4 (03/00) Parameter values for the HDTV standards for production and international programme exchange specifies parameter values for chromaticity, sampling, signal format, frame rates, etc., of high definition television signals.

An RGB709 color is represented by a list of three inexact numbers between 0 and 1. '(0. 0. 0.) is black '(1. 1. 1.) is white.

Function: rgb709->color rgb

rgb must be a list of 3 numbers. If rgb is valid RGB709 coordinates, then rgb709->color returns the color specified by rgb; otherwise returns #f.

Function: color:rgb709 r g b

Returns the RGB709 color composed of r, g, b. If the coordinates do not encode a valid RGB709 color, then an error is signaled.

Function: color->rgb709 color

Returns the list of 3 numbers encoding color in RGB709.

Perceptual Uniformity

Although properly encoding the chromaticity, tristimulus spaces do not match the logarithmic response of human visual systems to intensity. Minimum detectable differences between colors correspond to a smaller range of distances (6:1) in the L*a*b* and L*u*v* spaces than in tristimulus spaces (80:1). For this reason, color distances are computed in L*a*b* (or L*C*h).

Color Space: L*a*b*

Is a CIE color space which better matches the human visual system's perception of color. It is a list of three numbers:

• 0 <= L* <= 100 (CIE Lightness)

• -500 <= a* <= 500
• -200 <= b* <= 200

Function: l*a*b*->color L*a*b* white-point

L*a*b* must be a list of 3 numbers. If L*a*b* is valid L*a*b* coordinates, then l*a*b*->color returns the color specified by L*a*b*; otherwise returns #f.

Function: color:l*a*b* L* a* b* white-point

Returns the L*a*b* color composed of L*, a*, b* with white-point.

Function: color:l*a*b* L* a* b*

Returns the L*a*b* color composed of L*, a*, b*. If the coordinates do not encode a valid L*a*b* color, then an error is signaled.

Function: color->l*a*b* color white-point

Returns the list of 3 numbers encoding color in L*a*b* with white-point.

Function: color->l*a*b* color

Returns the list of 3 numbers encoding color in L*a*b*.

Color Space: L*u*v*

Is another CIE encoding designed to better match the human visual system's perception of color.

Function: l*u*v*->color L*u*v* white-point

L*u*v* must be a list of 3 numbers. If L*u*v* is valid L*u*v* coordinates, then l*u*v*->color returns the color specified by L*u*v*; otherwise returns #f.

Function: color:l*u*v* L* u* v* white-point

Returns the L*u*v* color composed of L*, u*, v* with white-point.

Function: color:l*u*v* L* u* v*

Returns the L*u*v* color composed of L*, u*, v*. If the coordinates do not encode a valid L*u*v* color, then an error is signaled.

Function: color->l*u*v* color white-point

Returns the list of 3 numbers encoding color in L*u*v* with white-point.

Function: color->l*u*v* color

Returns the list of 3 numbers encoding color in L*u*v*.

Cylindrical Coordinates

HSL (Hue Saturation Lightness), HSV (Hue Saturation Value), HSI (Hue Saturation Intensity) and HCI (Hue Chroma Intensity) are cylindrical color spaces (with angle hue). But these spaces are all defined in terms device-dependent RGB spaces.

One might wonder if there is some fundamental reason why intuitive specification of color must be device-dependent. But take heart! A cylindrical system can be based on L*a*b* and is used for predicting how close colors seem to observers.

Color Space: L*C*h

Expresses the *a and b* of L*a*b* in polar coordinates. It is a list of three numbers:

• 0 <= L* <= 100 (CIE Lightness)

• C* (CIE Chroma) is the distance from the neutral (gray) axis.
• 0 <= h <= 360 (CIE Hue) is the angle.

The colors by quadrant of h are:

0	red, orange, yellow	90
90	yellow, yellow-green, green	180
180	green, cyan (blue-green), blue	270
270	blue, purple, magenta	360

Function: l*c*h->color L*C*h white-point

L*C*h must be a list of 3 numbers. If L*C*h is valid L*C*h coordinates, then l*c*h->color returns the color specified by L*C*h; otherwise returns #f.

Function: color:l*c*h L* C* h white-point

Returns the L*C*h color composed of L*, C*, h with white-point.

Function: color:l*c*h L* C* h

Returns the L*C*h color composed of L*, C*, h. If the coordinates do not encode a valid L*C*h color, then an error is signaled.

Function: color->l*c*h color white-point

Returns the list of 3 numbers encoding color in L*C*h with white-point.

Function: color->l*c*h color

Returns the list of 3 numbers encoding color in L*C*h.

Digital Color Spaces

The color spaces discussed so far are impractical for image data because of numerical precision and computational requirements. In 1998 the IEC adopted A Standard Default Color Space for the Internet - sRGB (http://www.w3.org/Graphics/Color/sRGB). sRGB was cleverly designed to employ the 24-bit (256x256x256) color encoding already in widespread use; and the 2.2 gamma intrinsic to CRT monitors.

Conversion from CIEXYZ to digital (sRGB) color spaces is accomplished by conversion first to a RGB709 tristimulus space with D65 white-point; then each coordinate is individually subjected to the same non-linear mapping. Inverse operations in the reverse order create the inverse transform.

Color Space: sRGB

Is "A Standard Default Color Space for the Internet". Most display monitors will work fairly well with sRGB directly. Systems using ICC profiles (5) should work very well with sRGB.

Function: srgb->color rgb

rgb must be a list of 3 numbers. If rgb is valid sRGB coordinates, then srgb->color returns the color specified by rgb; otherwise returns #f.

Function: color:srgb r g b

Returns the sRGB color composed of r, g, b. If the coordinates do not encode a valid sRGB color, then an error is signaled.

Color Space: xRGB

Represents the equivalent sRGB color with a single 24-bit integer. The most significant 8 bits encode red, the middle 8 bits blue, and the least significant 8 bits green.

Function: color->srgb color

Returns the list of 3 integers encoding color in sRGB.

Function: color->xrgb color

Returns the 24-bit integer encoding color in sRGB.

Function: xrgb->color k

Returns the sRGB color composed of the 24-bit integer k.

Color Space: e-sRGB

Is "Photography - Electronic still picture imaging - Extended sRGB color encoding" (PIMA 7667:2001). It extends the gamut of sRGB; and its higher precision numbers provide a larger dynamic range.

A triplet of integers represent e-sRGB colors. Three precisions are supported:

e-sRGB10

0 to 1023

e-sRGB12

0 to 4095

e-sRGB16

0 to 65535

Function: e-srgb->color precision rgb

precision must be the integer 10, 12, or 16. rgb must be a list of 3 numbers. If rgb is valid e-sRGB coordinates, then e-srgb->color returns the color specified by rgb; otherwise returns #f.

Function: color:e-srgb 10 r g b

Returns the e-sRGB10 color composed of integers r, g, b.

Function: color:e-srgb 12 r g b

Returns the e-sRGB12 color composed of integers r, g, b.

Function: color:e-srgb 16 r g b

Returns the e-sRGB16 color composed of integers r, g, b. If the coordinates do not encode a valid e-sRGB color, then an error is signaled.

Function: color->e-srgb precision color

precision must be the integer 10, 12, or 16. color->e-srgb returns the list of 3 integers encoding color in sRGB10, sRGB12, or sRGB16.

5.9.3 Spectra

<A NAME="Spectra"></A> The following functions compute colors from spectra, scale color luminance, and extract chromaticity. XYZ is used in the names of procedures for unnormalized colors; the coordinates of CIEXYZ colors are constrained as described in 5.9.2 Color Spaces.

(require 'color-space)

A spectrum may be represented as:

• A procedure of one argument accepting real numbers from 380e-9 to 780e-9, the wavelength in meters; or
• A vector of real numbers representing intensity samples evenly spaced over some range of wavelengths overlapping the range 380e-9 to 780e-9.

CIEXYZ values are calculated as dot-product with the X, Y (Luminance), and Z Spectral Tristimulus Values. The files `cie1931.xyz' and `cie1964.xyz' in the distribution contain these CIE-defined values.

Feature: cie1964

Loads the Spectral Tristimulus Values defining CIE 1964 Supplementary Standard Colorimetric Observer.

Feature: cie1931

Loads the Spectral Tristimulus Values defining CIE 1931 Supplementary Standard Colorimetric Observer.

Feature: ciexyz

Requires Spectral Tristimulus Values, defaulting to cie1931.

(require 'cie1964) or (require 'cie1931) will load-ciexyz specific values used by the following spectrum conversion procedures. The spectrum conversion procedures (require 'ciexyz) to assure that a set is loaded.

Function: spectrum->XYZ proc

proc must be a function of one argument. spectrum->XYZ computes the CIEXYZ(1931) values for the spectrum returned by proc when called with arguments from 380e-9 to 780e-9, the wavelength in meters.

Function: spectrum->XYZ spectrum x1 x2

x1 and x2 must be positive real numbers specifying the wavelengths (in meters) corresponding to the zeroth and last elements of vector or list spectrum. spectrum->XYZ returns the CIEXYZ(1931) values for a light source with spectral values proportional to the elements of spectrum at evenly spaced wavelengths between x1 and x2.

Compute the colors of 6500.K and 5000.K blackbody radiation:

(require 'color-space) (define xyz (spectrum->XYZ (blackbody-spectrum 6500))) (define y_n (cadr xyz)) (map (lambda (x) (/ x y_n)) xyz) => (0.9687111145512467 1.0 1.1210875945303613) (define xyz (spectrum->XYZ (blackbody-spectrum 5000))) (map (lambda (x) (/ x y_n)) xyz) => (0.2933441826889158 0.2988931825387761 0.25783646831201573)

Function: spectrum->CIEXYZ proc

Function: spectrum->CIEXYZ spectrum x1 x2

spectrum->CIEXYZ computes the CIEXYZ(1931) values for the spectrum, scaled so their sum is 1.

Function: spectrum->chromaticity proc

Function: spectrum->chromaticity spectrum x1 x2

Computes the chromaticity for the given spectrum.

Function: wavelength->XYZ w

Function: wavelength->chromaticity w

Function: wavelength->CIEXYZ w

w must be a number between 380e-9 to 780e-9. wavelength->XYZ returns (unnormalized) XYZ values for a monochromatic light source with wavelength w. wavelength->chromaticity returns the chromaticity for a monochromatic light source with wavelength w. wavelength->CIEXYZ returns XYZ values for the saturated color having chromaticity of a monochromatic light source with wavelength w.

Function: blackbody-spectrum temp

Function: blackbody-spectrum temp span

Returns a procedure of one argument (wavelength in meters), which returns the radiance of a black body at temp.

The optional argument span is the wavelength analog of bandwidth. With the default span of 1.nm (1e-9.m), the values returned by the procedure correspond to the power of the photons with wavelengths w to w+1e-9.

Function: temperature->XYZ x

The positive number x is a temperature in degrees kelvin. temperature->XYZ computes the CIEXYZ(1931) values for the spectrum of a black body at temperature x.

Compute the chromaticities of 6500.K and 5000.K blackbody radiation:

(require 'color-space) (XYZ->chromaticity (temperature->XYZ 6500)) => (0.3135191660557008 0.3236456786200268) (XYZ->chromaticity (temperature->XYZ 5000)) => (0.34508082841161052 0.3516084965163377)

Function: temperature->CIEXYZ x

The positive number x is a temperature in degrees kelvin. temperature->CIEXYZ computes the CIEXYZ(1931) values for the spectrum of a black body at temperature x, scaled to be just inside the RGB709 gamut.

Function: temperature->chromaticity x

Function: XYZ:normalize xyz

xyz is a list of three non-negative real numbers. XYZ:normalize returns a list of numbers proportional to xyz; scaled so their sum is 1.

Function: XYZ:normalize-colors colors ...

colors is a list of XYZ triples. XYZ:normalize-colors scales all the triples by a common factor such that the maximum sum of numbers in a scaled triple is 1.

Function: XYZ->chromaticity xyz

Returns a two element list: the x and y components of xyz normalized to 1 (= x + y + z).

Function: chromaticity->CIEXYZ x y

Returns the list of x, and y, 1 - y - x.

Function: chromaticity->whitepoint x y

Returns the CIEXYZ(1931) values having luminosity 1 and chromaticity x and y.

Many color datasets are expressed in xyY format; chromaticity with CIE luminance (Y). But xyY is not a CIE standard like CIEXYZ, CIELAB, and CIELUV. Although chrominance is well defined, the luminance component is sometimes scaled to 1, sometimes to 100, but usually has no obvious range. With no given whitepoint, the only reasonable course is to ascertain the luminance range of a dataset and normalize the values to lie from 0 to 1.

Function: XYZ->xyY xyz

Returns a three element list: the x and y components of XYZ normalized to 1, and CIE luminance Y.

Function: xyY->XYZ xyY

Function: xyY:normalize-colors colors

colors is a list of xyY triples. xyY:normalize-colors scales each chromaticity so it sums to 1 or less; and divides the Y values by the maximum Y in the dataset, so all lie between 0 and 1.

Function: xyY:normalize-colors colors n

If n is positive real, then xyY:normalize-colors divides the Y values by n times the maximum Y in the dataset.

If n is an exact non-positive integer, then xyY:normalize-colors divides the Y values by the maximum of the Ys in the dataset excepting the -n largest Y values.

In all cases, returned Y values are limited to lie from 0 to 1.

Why would one want to normalize to other than 1? If the sun or its reflection is the brightest object in a scene, then normalizing to its luminance will tend to make the rest of the scene very dark. As with photographs, limiting the specular highlights looks better than darkening everything else.

The results of measurements being what they are, xyY:normalize-colors is extremely tolerant. Negative numbers are replaced with zero, and chromaticities with sums greater than one are scaled to sum to one.

5.9.4 Color Difference Metrics

<A NAME="Color_Difference_Metrics"></A>

(require 'color-space)

The low-level metric functions operate on lists of 3 numbers, lab1, lab2, lch1, or lch2.

(require 'color)

The wrapped functions operate on objects of type color, color1 and color2 in the function entries.

Function: L*a*b*:DE* lab1 lab2

Returns the Euclidean distance between lab1 and lab2.

Function: CIE:DE* color1 color2 white-point

Function: CIE:DE* color1 color2

Returns the Euclidean distance in L*a*b* space between color1 and color2.

Function: L*C*h:DE*94 lch1 lch2 parametric-factors

Function: L*C*h:DE*94 lch1 lch2

Function: CIE:DE*94 color1 color2 parametric-factors

Function: CIE:DE*94 color1 color2

Measures distance in the L*C*h cylindrical color-space. The three axes are individually scaled (depending on C*) in their contributions to the total distance.

The CIE has defined reference conditions under which the metric with default parameters can be expected to perform well. These are:

• The specimens are homogeneous in colour.
• The colour difference (CIELAB) is <= 5 units.
• They are placed in direct edge contact.
• Each specimen subtends an angle of >4 degrees to the assessor, whose colour vision is normal.
• They are illuminated at 1000 lux, and viewed against a background of uniform grey, with L* of 50, under illumination simulating D65.

The parametric-factors argument is a list of 3 quantities kL, kC and kH. parametric-factors independently adjust each colour-difference term to account for any deviations from the reference viewing conditions. Under the reference conditions explained above, the default is kL = kC = kH = 1.

The Color Measurement Committee of The Society of Dyers and Colorists in Great Britain created a more sophisticated color-distance function for use in judging the consistency of dye lots. With CMC:DE* it is possible to use a single value pass/fail tolerance for all shades.

Function: CMC-DE lch1 lch2 parametric-factors

Function: CMC-DE lch1 lch2 l c

Function: CMC-DE lch1 lch2 l

Function: CMC-DE lch1 lch2

Function: CMC:DE* color1 color2 l c

Function: CMC:DE* color1 color2

CMC:DE is a L*C*h metric. The parametric-factors argument is a list of 2 numbers l and c. l and c parameterize this metric. 1 and 1 are recommended for perceptibility; the default, 2 and 1, for acceptability.

5.9.5 Color Conversions

<A NAME="Color_Conversions"></A>

This package contains the low-level color conversion and color metric routines operating on lists of 3 numbers. There is no type or range checking.

(require 'color-space)

Constant: CIEXYZ:D65

Is the color of 6500.K (blackbody) illumination. D65 is close to the average color of daylight.

Constant: CIEXYZ:D50

Is the color of 5000.K (blackbody) illumination. D50 is the color of indoor lighting by incandescent bulbs.

Constant: CIEXYZ:A

Constant: CIEXYZ:B

Constant: CIEXYZ:C

Constant: CIEXYZ:E

CIE 1931 illuminants normalized to 1 = y.

Function: color:linear-transform matrix row

Function: CIEXYZ->RGB709 xyz

Function: RGB709->CIEXYZ srgb

Function: CIEXYZ->L*u*v* xyz white-point

Function: CIEXYZ->L*u*v* xyz

Function: L*u*v*->CIEXYZ L*u*v* white-point

Function: L*u*v*->CIEXYZ L*u*v*

The white-point defaults to CIEXYZ:D65.

Function: CIEXYZ->L*a*b* xyz white-point

Function: CIEXYZ->L*a*b* xyz

Function: L*a*b*->CIEXYZ L*a*b* white-point

Function: L*a*b*->CIEXYZ L*a*b*

The XYZ white-point defaults to CIEXYZ:D65.

Function: L*a*b*->L*C*h L*a*b*

Function: L*C*h->L*a*b* L*C*h

Function: CIEXYZ->sRGB xyz

Function: sRGB->CIEXYZ srgb

Function: CIEXYZ->xRGB xyz

Function: xRGB->CIEXYZ srgb

Function: sRGB->xRGB xyz

Function: xRGB->sRGB srgb

Function: CIEXYZ->e-sRGB n xyz

Function: e-sRGB->CIEXYZ n srgb

Function: sRGB->e-sRGB n srgb

Function: e-sRGB->sRGB n srgb

The integer n must be 10, 12, or 16. Because sRGB and e-sRGB use the same RGB709 chromaticities, conversion between them is simpler than conversion through CIEXYZ.

Do not convert e-sRGB precision through e-sRGB->sRGB then sRGB->e-sRGB -- values would be truncated to 8-bits!

Function: e-sRGB->e-sRGB n1 srgb n2

The integers n1 and n2 must be 10, 12, or 16. e-sRGB->e-sRGB converts srgb to e-sRGB of precision n2.

5.9.6 Color Names

<A NAME="Color_Names"></A> (require 'color-names)

Rather than ballast the color dictionaries with numbered grays, file->color-dictionary discards them. They are provided through the grey procedure:

Function: grey k

Returns (inexact->exact (round (* k 2.55))), the X11 color grey<k>.

The column names in a color dictionary are unimportant; the first field is the key, and the second is the color-string.

Function: color-name:canonicalize name

Returns a downcased copy of the string or symbol name with `_', `-', and whitespace removed.

Function: color-name->color name table1 table2 ...

table1, table2, ... must be color-dictionary tables. color-name->color searches for the canonical form of name in table1, table2, ... in order; returning the color-string of the first matching record; #f otherwise.

Function: color-dictionaries->lookup table1 table2 ...

table1, table2, ... must be color-dictionary tables. color-dictionaries->lookup returns a procedure which searches for the canonical form of its string argument in table1, table2, ...; returning the color-string of the first matching record; and #f otherwise.

Function: color-dictionary name rdb base-table-type

rdb must be a string naming a relational database file; and the symbol name a table therein. The database will be opened as base-table-type. color-dictionary returns the read-only table name in database name if it exists; #f otherwise.

Function: color-dictionary name rdb

rdb must be an open relational database or a string naming a relational database file; and the symbol name a table therein. color-dictionary returns the read-only table name in database name if it exists; #f otherwise.

Function: load-color-dictionary name rdb base-table-type

Function: load-color-dictionary name rdb

rdb must be a string naming a relational database file; and the symbol name a table therein. If the symbol base-table-type is provided, the database will be opened as base-table-type. load-color-dictionary creates a top-level definition of the symbol name to a lookup procedure for the color dictionary name in rdb.

The value returned by load-color-dictionary is unspecified.

Dictionary Creation

(require 'color-database)

Function: file->color-dictionary file table-name rdb base-table-type

Function: file->color-dictionary file table-name rdb

rdb must be an open relational database or a string naming a relational database file, table-name a symbol, and the string file must name an existing file with colornames and their corresponding xRGB (6-digit hex) values. file->color-dictionary creates a table table-name in rdb and enters the associations found in file into it.

Function: url->color-dictionary url table-name rdb base-table-type

Function: url->color-dictionary url table-name rdb

rdb must be an open relational database or a string naming a relational database file and table-name a symbol. url->color-dictionary retrieves the resource named by the string url using the wget program; then calls file->color-dictionary to enter its associations in table-name in url.

http://swissnet.ai.mit.edu/~jaffer/Color/Dictionaries.html

Describes and evaluates several color-name dictionaries on the web. The following procedure creates a database containing two of these dictionaries.

Function: make-slib-color-name-db

Creates an alist-table relational database in library-vicinity containing the Resene and saturate color-name dictionaries.

If the files `resenecolours.txt' and `saturate.txt' exist in the library-vicinity, then they used as the source of color-name data. Otherwise, make-slib-color-name-db calls url->color-dictionary with the URLs of appropriate source files.

The Short List

(require 'saturate)

Function: saturate name

Looks for name among the 19 saturated colors from Approximate Colors on CIE Chromaticity Diagram:

reddish orange	orange	yellowish orange	yellow
greenish yellow	yellow green	yellowish green	green
bluish green	blue green	greenish blue	blue
purplish blue	bluish purple	purple	reddish purple
red purple	purplish red	red

(http://swissnet.ai.mit.edu/~jaffer/Color/saturate.pdf). If name is found, the corresponding color is returned. Otherwise #f is returned. Use saturate only for light source colors.

Resene Paints Limited, New Zealand's largest privately-owned and operated paint manufacturing company, has generously made their Resene RGB Values List available.

(require 'resene)

Function: resene name

Looks for name among the 1300 entries in the Resene color-name dictionary (http://swissnet.ai.mit.edu/~jaffer/Color/resene.pdf). If name is found, the corresponding color is returned. Otherwise #f is returned. The Resene RGB Values List is an excellent source for surface colors.

If you include the Resene RGB Values List in binary form in a program, then you must include its license with your program:

Resene RGB Values List
For further information refer to http://www.resene.co.nz
Copyright Resene Paints Ltd 2001

Permission to copy this dictionary, to modify it, to redistribute it, to distribute modified versions, and to use it for any purpose is granted, subject to the following restrictions and understandings.

1. Any text copy made of this dictionary must include this copyright notice in full.

2. Any redistribution in binary form must reproduce this copyright notice in the documentation or other materials provided with the distribution.

3. Resene Paints Ltd makes no warranty or representation that this dictionary is error-free, and is under no obligation to provide any services, by way of maintenance, update, or otherwise.

4. There shall be no use of the name of Resene or Resene Paints Ltd in any advertising, promotional, or sales literature without prior written consent in each case.

5. These RGB colour formulations may not be used to the detriment of Resene Paints Ltd.

5.9.7 Daylight

<A NAME="Daylight"></A> (require 'daylight)

This package calculates the colors of sky as detailed in:
http://www.cs.utah.edu/vissim/papers/sunsky/sunsky.pdf
A Practical Analytic Model for Daylight
A. J. Preetham, Peter Shirley, Brian Smits

Function: solar-hour julian-day hour

Returns the solar-time in hours given the integer julian-day in the range 1 to 366, and the local time in hours.

To be meticulous, subtract 4 minutes for each degree of longitude west of the standard meridian of your time zone.

Function: solar-declination julian-day

Function: solar-polar declination latitude solar-hour

Returns a list of theta_s, the solar angle from the zenith, and phi_s, the solar azimuth. 0 <= theta_s measured in degrees. phi_s is measured in degrees from due south; west of south being positive.

Turbidity is a measure of the fraction of scattering due to haze as opposed to molecules. This is a convenient quantity because it can be estimated based on visibility of distant objects. This model fails for turbidity values less than 1.3.

    _______________________________________________________________
512|-:                                                             |
   | * pure-air                                                    |
256|-:**                                                           |
   | : ** exceptionally-clear                                      |
128|-:   *                                                         |
   | :    **                                                       |
 64|-:      *                                                      |
   | :       ** very-clear                                         |
 32|-:         **                                                  |
   | :           **                                                |
 16|-:             *** clear                                       |
   | :               ****                                          |
  8|-:                  ****                                       |
   | :                     **** light-haze                         |
  4|-:                         ****                                |
   | :                             ******                          |
  2|-:                                  ******** haze         thin-|
   | :                                          ***********    fog |
  1|-:----------------------------------------------------*******--|
   |_:____.____:____.____:____.____:____.____:____.____:____.____:_|
     1         2         4         8        16        32        64
              Meterorological range (km) versus Turbidity

Function: sunlight-spectrum turbidity theta_s

Returns a vector of 41 values, the spectrum of sunlight from 380.nm to 790.nm for a given turbidity and theta_s.

Function: sunlight-xyz turbidity theta_s

Returns (unnormalized) XYZ values for color of sunlight for a given turbidity and theta_s.

Function: sunlight-ciexyz turbidity theta_s

Given turbidity and theta_s, sunlight-ciexyz returns the CIEXYZ triple for color of sunlight scaled to be just inside the RGB709 gamut.

Function: zenith-xyy turbidity theta_s

Returns the xyY (chromaticity and luminance) at the zenith. The Luminance has units kcd/m^2.

Function: overcast-sky-color-xyy turbidity theta_s

turbidity is a positive real number expressing the amount of light scattering. The real number theta_s is the solar angle from the zenith in degrees.

overcast-sky-color-xyy returns a function of one angle theta, the angle from the zenith of the viewing direction (in degrees); and returning the xyY value for light coming from that elevation of the sky.

Function: clear-sky-color-xyy turbidity theta_s phi_s

Function: sky-color-xyy turbidity theta_s phi_s

turbidity is a positive real number expressing the amount of light scattering. The real number theta_s is the solar angle from the zenith in degrees. The real number phi_s is the solar angle from south.

clear-sky-color-xyy returns a function of two angles, theta and phi which specify the angles from the zenith and south meridian of the viewing direction (in degrees); returning the xyY value for light coming from that direction of the sky.

sky-color-xyY calls overcast-sky-color-xyY for turbidity <= 20; otherwise the clear-sky-color-xyy function.

5.10 Root Finding

(require 'root)

Function: newton: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.

Function: secant:find-root f x0 x1 prec

Function: secant:find-bracketed-root f x0 x1 prec

Given a real valued procedure f and two real valued starting points x0 and x1, 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 x0 and x1 are chosen such that they bracket a root, that is

(or (< (f x0) 0 (f x1)) (< (f x1) 0 (f x0)))

then the root returned will be between x0 and x1, and f will not be passed an argument outside of that interval.

secant:find-bracketed-root will return #f unless x0 and x1 bracket a root.

The secant method is used until a bracketing interval is found, at which point a modified regula falsi method is used.

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

If prec is a procedure it should accept 5 arguments: x0 f0 x1 f1 and count, where f0 will be (f x0), f1 (f x1), and count the number of iterations performed so far. prec should return non-false if the iteration should be stopped.

5.11 Minimizing

(require 'minimize)

The Golden Section Search (6) algorithm finds minima of functions which are expensive to compute or for which derivatives are not available. Although optimum for the general case, convergence is slow, requiring nearly 100 iterations for the example (x^3-2x-5).

If the derivative is available, Newton-Raphson is probably a better choice. If the function is inexpensive to compute, consider approximating the derivative.

Function: golden-section-search f x0 x1 prec

x_0 are x_1 real numbers. The (single argument) procedure f is unimodal over the open interval (x_0, x_1). That is, there is exactly one point in the interval for which the derivative of f is zero.

golden-section-search returns a pair (x . f(x)) where f(x) is the minimum. The prec parameter is the stop criterion. If prec is a positive number, then the iteration continues until x is within prec from the true value. If prec is a negative integer, then the procedure will iterate -prec times or until convergence. If prec is a procedure of seven arguments, x0, x1, a, b, fa, fb, and count, then the iterations will stop when the procedure returns #t.

Analytically, the minimum of x^3-2x-5 is 0.816497.

(define func (lambda (x) (+ (* x (+ (* x x) -2)) -5))) (golden-section-search func 0 1 (/ 10000)) ==> (816.4883855245578e-3 . -6.0886621077391165) (golden-section-search func 0 1 -5) ==> (819.6601125010515e-3 . -6.088637561916407) (golden-section-search func 0 1 (lambda (a b c d e f g ) (= g 500))) ==> (816.4965933140557e-3 . -6.088662107903635)

5.12 Commutative Rings

Scheme provides a consistent and capable set of numeric functions. Inexacts implement a field; integers a commutative ring (and Euclidean domain). This package allows one to use basic Scheme numeric functions with symbols and non-numeric elements of commutative rings.

(require 'commutative-ring)

The commutative-ring package makes the procedures +, -, *, /, and ^ careful in the sense that any non-numeric arguments they do not reduce appear in the expression output. In order to see what working with this package is like, self-set all the single letter identifiers (to their corresponding symbols).

(define a 'a)
...
(define z 'z)

Or just (require 'self-set). Now try some sample expressions:

(+ (+ a b) (- a b)) => (* a 2)
(* (+ a b) (+ a b)) => (^ (+ a b) 2)
(* (+ a b) (- a b)) => (* (+ a b) (- a b))
(* (- a b) (- a b)) => (^ (- a b) 2)
(* (- a b) (+ a b)) => (* (+ a b) (- a b))
(/ (+ a b) (+ c d)) => (/ (+ a b) (+ c d))
(^ (+ a b) 3) => (^ (+ a b) 3)
(^ (+ a 2) 3) => (^ (+ 2 a) 3)

Associative rules have been applied and repeated addition and multiplication converted to multiplication and exponentiation.

We can enable distributive rules, thus expanding to sum of products form:

(set! *ruleset* (combined-rulesets distribute* distribute/))

(* (+ a b) (+ a b)) => (+ (* 2 a b) (^ a 2) (^ b 2))
(* (+ a b) (- a b)) => (- (^ a 2) (^ b 2))
(* (- a b) (- a b)) => (- (+ (^ a 2) (^ b 2)) (* 2 a b))
(* (- a b) (+ a b)) => (- (^ a 2) (^ b 2))
(/ (+ a b) (+ c d)) => (+ (/ a (+ c d)) (/ b (+ c d)))
(/ (+ a b) (- c d)) => (+ (/ a (- c d)) (/ b (- c d)))
(/ (- a b) (- c d)) => (- (/ a (- c d)) (/ b (- c d)))
(/ (- a b) (+ c d)) => (- (/ a (+ c d)) (/ b (+ c d)))
(^ (+ a b) 3) => (+ (* 3 a (^ b 2)) (* 3 b (^ a 2)) (^ a 3) (^ b 3))
(^ (+ a 2) 3) => (+ 8 (* a 12) (* (^ a 2) 6) (^ a 3))

Use of this package is not restricted to simple arithmetic expressions:

(require 'determinant)

(determinant '((a b c) (d e f) (g h i))) =>
(- (+ (* a e i) (* b f g) (* c d h)) (* a f h) (* b d i) (* c e g))

Currently, only +, -, *, /, and ^ support non-numeric elements. Expressions with - are converted to equivalent expressions without -, so behavior for - is not defined separately. / expressions are handled similarly.

This list might be extended to include quotient, modulo, remainder, lcm, and gcd; but these work only for the more restrictive Euclidean (Unique Factorization) Domain.

5.13 Rules and Rulesets

The commutative-ring package allows control of ring properties through the use of rulesets.

Variable: *ruleset*

Contains the set of rules currently in effect. Rules defined by cring:define-rule are stored within the value of *ruleset* at the time cring:define-rule is called. If *ruleset* is #f, then no rules apply.

Function: make-ruleset rule1 ...

Function: make-ruleset name rule1 ...

Returns a new ruleset containing the rules formed by applying cring:define-rule to each 4-element list argument rule. If the first argument to make-ruleset is a symbol, then the database table created for the new ruleset will be named name. Calling make-ruleset with no rule arguments creates an empty ruleset.

Function: combined-rulesets ruleset1 ...

Function: combined-rulesets name ruleset1 ...

Returns a new ruleset containing the rules contained in each ruleset argument ruleset. If the first argument to combined-ruleset is a symbol, then the database table created for the new ruleset will be named name. Calling combined-ruleset with no ruleset arguments creates an empty ruleset.

Two rulesets are defined by this package.

Constant: distribute*

Contains the ruleset to distribute multiplication over addition and subtraction.

Constant: distribute/

Contains the ruleset to distribute division over addition and subtraction.

Take care when using both distribute* and distribute/ simultaneously. It is possible to put / into an infinite loop.

You can specify how sum and product expressions containing non-numeric elements simplify by specifying the rules for + or * for cases where expressions involving objects reduce to numbers or to expressions involving different non-numeric elements.

Function: cring:define-rule op sub-op1 sub-op2 reduction

Defines a rule for the case when the operation represented by symbol op is applied to lists whose cars are sub-op1 and sub-op2, respectively. The argument reduction is a procedure accepting 2 arguments which will be lists whose cars are sub-op1 and sub-op2.

Function: cring:define-rule op sub-op1 'identity reduction

Defines a rule for the case when the operation represented by symbol op is applied to a list whose car is sub-op1, and some other argument. Reduction will be called with the list whose car is sub-op1 and some other argument.

If reduction returns #f, the reduction has failed and other reductions will be tried. If reduction returns a non-false value, that value will replace the two arguments in arithmetic (+, -, and *) calculations involving non-numeric elements.

The operations + and * are assumed commutative; hence both orders of arguments to reduction will be tried if necessary.

The following rule is the definition for distributing * over +.

(cring:define-rule '* '+ 'identity (lambda (exp1 exp2) (apply + (map (lambda (trm) (* trm exp2)) (cdr exp1))))))

5.14 How to Create a Commutative Ring

The first step in creating your commutative ring is to write procedures to create elements of the ring. A non-numeric element of the ring must be represented as a list whose first element is a symbol or string. This first element identifies the type of the object. A convenient and clear convention is to make the type-identifying element be the same symbol whose top-level value is the procedure to create it.

(define (n . list1)
  (cond ((and (= 2 (length list1))
              (eq? (car list1) (cadr list1)))
         0)
        ((not (term< (first list1) (last1 list1)))
         (apply n (reverse list1)))
        (else (cons 'n list1))))

(define (s x y) (n x y))

(define (m . list1)
  (cond ((neq? (first list1) (term_min list1))
         (apply m (cyclicrotate list1)))
        ((term< (last1 list1) (cadr list1))
         (apply m (reverse (cyclicrotate list1))))
        (else (cons 'm list1))))

Define a procedure to multiply 2 non-numeric elements of the ring. Other multiplicatons are handled automatically. Objects for which rules have not been defined are not changed.

(define (n*n ni nj)
  (let ((list1 (cdr ni)) (list2 (cdr nj)))
    (cond ((null? (intersection list1 list2)) #f)
          ((and (eq? (last1 list1) (first list2))
                (neq? (first list1) (last1 list2)))
           (apply n (splice list1 list2)))
          ((and (eq? (first list1) (first list2))
                (neq? (last1 list1) (last1 list2)))
           (apply n (splice (reverse list1) list2)))
          ((and (eq? (last1 list1) (last1 list2))
                (neq? (first list1) (first list2)))
           (apply n (splice list1 (reverse list2))))
          ((and (eq? (last1 list1) (first list2))
                (eq? (first list1) (last1 list2)))
           (apply m (cyclicsplice list1 list2)))
          ((and (eq? (first list1) (first list2))
                (eq? (last1 list1) (last1 list2)))
           (apply m (cyclicsplice (reverse list1) list2)))
          (else #f))))

Test the procedures to see if they work.

;;; where cyclicrotate(list) is cyclic rotation of the list one step
;;; by putting the first element at the end
(define (cyclicrotate list1)
  (append (rest list1) (list (first list1))))
;;; and where term_min(list) is the element of the list which is
;;; first in the term ordering.
(define (term_min list1)
  (car (sort list1 term<)))
(define (term< sym1 sym2)
  (string<? (symbol->string sym1) (symbol->string sym2)))
(define first car)
(define rest cdr)
(define (last1 list1) (car (last-pair list1)))
(define (neq? obj1 obj2) (not (eq? obj1 obj2)))
;;; where splice is the concatenation of list1 and list2 except that their
;;; common element is not repeated.
(define (splice list1 list2)
  (cond ((eq? (last1 list1) (first list2))
         (append list1 (cdr list2)))
        (else (slib:error 'splice list1 list2))))
;;; where cyclicsplice is the result of leaving off the last element of
;;; splice(list1,list2).
(define (cyclicsplice list1 list2)
  (cond ((and (eq? (last1 list1) (first list2))
              (eq? (first list1) (last1 list2)))
         (butlast (splice list1 list2) 1))
        (else (slib:error 'cyclicsplice list1 list2))))

(N*N (S a b) (S a b)) => (m a b)

Then register the rule for multiplying type N objects by type N objects.

(cring:define-rule '* 'N 'N N*N))

Now we are ready to compute!

(define (t)
  (define detM
    (+ (* (S g b)
          (+ (* (S f d)
                (- (* (S a f) (S d g)) (* (S a g) (S d f))))
             (* (S f f)
                (- (* (S a g) (S d d)) (* (S a d) (S d g))))
             (* (S f g)
                (- (* (S a d) (S d f)) (* (S a f) (S d d))))))
       (* (S g d)
          (+ (* (S f b)
                (- (* (S a g) (S d f)) (* (S a f) (S d g))))
             (* (S f f)
                (- (* (S a b) (S d g)) (* (S a g) (S d b))))
             (* (S f g)
                (- (* (S a f) (S d b)) (* (S a b) (S d f))))))
       (* (S g f)
          (+ (* (S f b)
                (- (* (S a d) (S d g)) (* (S a g) (S d d))))
             (* (S f d)
                (- (* (S a g) (S d b)) (* (S a b) (S d g))))
             (* (S f g)
                (- (* (S a b) (S d d)) (* (S a d) (S d b))))))
       (* (S g g)
          (+ (* (S f b)
                (- (* (S a f) (S d d)) (* (S a d) (S d f))))
             (* (S f d)
                (- (* (S a b) (S d f)) (* (S a f) (S d b))))
             (* (S f f)
                (- (* (S a d) (S d b)) (* (S a b) (S d d))))))))
  (* (S b e) (S c a) (S e c)
     detM
     ))
(pretty-print (t))
-|
(- (+ (m a c e b d f g)
      (m a c e b d g f)
      (m a c e b f d g)
      (m a c e b f g d)
      (m a c e b g d f)
      (m a c e b g f d))
   (* 2 (m a b e c) (m d f g))
   (* (m a c e b d) (m f g))
   (* (m a c e b f) (m d g))
   (* (m a c e b g) (m d f)))

5.15 Matrix Algebra

(require 'determinant)

A Matrix can be either a list of lists (rows) or an array. As with linear-algebra texts, this package uses 1-based coordinates.

Function: matrix->lists matrix

Returns the list-of-lists form of matrix.

Function: matrix->array matrix

Returns the (ones-based) array form of matrix.

Function: determinant matrix

matrix must be a square matrix. determinant returns the determinant of matrix.

(require 'determinant) (determinant '((1 2) (3 4))) => -2 (determinant '((1 2 3) (4 5 6) (7 8 9))) => 0

Function: transpose matrix

Returns a copy of matrix flipped over the diagonal containing the 1,1 element.

Function: matrix:product m1 m2

Returns the product of matrices m1 and m2.

Function: matrix:inverse matrix

matrix must be a square matrix. If matrix is singlar, then matrix:inverse returns #f; otherwise matrix:inverse returns the matrix:product inverse of matrix.