1# math.tcl --
2#
3#	Collection of math functions.
4#
5# Copyright (c) 1998-2000 by Ajuba Solutions.
6#
7# See the file "license.terms" for information on usage and redistribution
8# of this file, and for a DISCLAIMER OF ALL WARRANTIES.
9#
10# RCS: @(#) $Id: misc.tcl,v 1.6 2005/10/10 14:02:47 arjenmarkus Exp $
11
12package require Tcl 8.2		;# uses [lindex $l end-$integer]
13namespace eval ::math {
14}
15
16# ::math::cov --
17#
18#	Return the coefficient of variation of three or more values
19#
20# Arguments:
21#	val1	first value
22#	val2	second value
23#	args	other values
24#
25# Results:
26#	cov	coefficient of variation expressed as percent value
27
28proc ::math::cov {val1 val2 args} {
29     set sum [ expr { $val1+$val2 } ]
30     set N [ expr { [ llength $args ] + 2 } ]
31     foreach val $args {
32        set sum [ expr { $sum+$val } ]
33     }
34     set mean [ expr { $sum/$N } ]
35     set sigma_sq 0
36     foreach val [ concat $val1 $val2 $args ] {
37        set sigma_sq [ expr { $sigma_sq+pow(($val-$mean),2) } ]
38     }
39     set sigma_sq [ expr { $sigma_sq/($N-1) } ]
40     set sigma [ expr { sqrt($sigma_sq) } ]
41     if { $mean != 0.0 } {
42        set cov [ expr { ($sigma/$mean)*100 } ]
43     } else {
44        return -code error -errorinfo "Cov undefined for data with zero mean" -errorcode {ARITH DOMAIN}
45     }
46     set cov
47}
48
49# ::math::fibonacci --
50#
51#	Return the n'th fibonacci number.
52#
53# Arguments:
54#	n	The index in the sequence to compute.
55#
56# Results:
57#	fib	The n'th fibonacci number.
58
59proc ::math::fibonacci {n} {
60    if { $n == 0 } {
61	return 0
62    } else {
63	set prev0 0
64	set prev1 1
65	for {set i 1} {$i < $n} {incr i} {
66	    set tmp $prev1
67	    incr prev1 $prev0
68	    set prev0 $tmp
69	}
70	return $prev1
71    }
72}
73
74# ::math::integrate --
75#
76#	calculate the area under a curve defined by a set of (x,y) data pairs.
77#	the x data must increase monotonically throughout the data set for the
78#	calculation to be meaningful, therefore the monotonic condition is
79#	tested, and an error is thrown if the x value is found to be
80#	decreasing.
81#
82# Arguments:
83#	xy_pairs	list of x y pairs (eg, 0 0 10 10 20 20 ...); at least 5
84#			data pairs are required, and if the number of data
85#			pairs is even, a padding value of (x0, 0) will be
86#			added.
87#
88# Results:
89#	result		A two-element list consisting of the area and error
90#			bound (calculation is "Simpson's rule")
91
92proc ::math::integrate { xy_pairs } {
93
94     set length [ llength $xy_pairs ]
95
96     if { $length < 10 } {
97        return -code error "at least 5 x,y pairs must be given"
98     }
99
100     ;## are we dealing with x,y pairs?
101     if { [ expr {$length % 2} ] } {
102        return -code error "unmatched xy pair in input"
103     }
104
105     ;## are there an even number of pairs?  Augment.
106     if { ! [ expr {$length % 4} ] } {
107        set xy_pairs [ concat [ lindex $xy_pairs 0 ] 0 $xy_pairs ]
108     }
109     set x0   [ lindex $xy_pairs 0     ]
110     set x1   [ lindex $xy_pairs 2     ]
111     set xn   [ lindex $xy_pairs end-1 ]
112     set xnminus1 [ lindex $xy_pairs end-3 ]
113
114     if { $x1 < $x0 } {
115        return -code error "monotonicity broken by x1"
116     }
117
118     if { $xn < $xnminus1 } {
119        return -code error "monotonicity broken by xn"
120     }
121
122     ;## handle the assymetrical elements 0, n, and n-1.
123     set sum [ expr {[ lindex $xy_pairs 1 ] + [ lindex $xy_pairs end ]} ]
124     set sum [ expr {$sum + (4*[ lindex $xy_pairs end-2 ])} ]
125
126     set data [ lrange $xy_pairs 2 end-4 ]
127
128     set xmax $x1
129     set i 1
130     foreach {x1 y1 x2 y2} $data {
131        incr i
132        if { $x1 < $xmax } {
133           return -code error "monotonicity broken by x$i"
134        }
135        set xmax $x1
136        incr i
137        if { $x2 < $xmax } {
138           return -code error "monotonicity broken by x$i"
139        }
140        set xmax $x2
141        set sum [ expr {$sum + (4*$y1) + (2*$y2)} ]
142     }
143
144     if { $xmax > $xnminus1 } {
145        return -code error "monotonicity broken by xn-1"
146     }
147
148     set h [ expr { ( $xn - $x0 ) / $i } ]
149     set area [ expr { ( $h / 3.0 ) * $sum } ]
150     set err_bound  [ expr { ( ( $xn - $x0 ) / 180.0 ) * pow($h,4) * $xn } ]
151     return [ list $area $err_bound ]
152}
153
154# ::math::max --
155#
156#	Return the maximum of two or more values
157#
158# Arguments:
159#	val	first value
160#	args	other values
161#
162# Results:
163#	max	maximum value
164
165proc ::math::max {val args} {
166    set max $val
167    foreach val $args {
168	if { $val > $max } {
169	    set max $val
170	}
171    }
172    set max
173}
174
175# ::math::mean --
176#
177#	Return the mean of two or more values
178#
179# Arguments:
180#	val	first value
181#	args	other values
182#
183# Results:
184#	mean	arithmetic mean value
185
186proc ::math::mean {val args} {
187    set sum $val
188    set N [ expr { [ llength $args ] + 1 } ]
189    foreach val $args {
190        set sum [ expr { $sum + $val } ]
191    }
192    set mean [expr { double($sum) / $N }]
193}
194
195# ::math::min --
196#
197#	Return the minimum of two or more values
198#
199# Arguments:
200#	val	first value
201#	args	other values
202#
203# Results:
204#	min	minimum value
205
206proc ::math::min {val args} {
207    set min $val
208    foreach val $args {
209	if { $val < $min } {
210	    set min $val
211	}
212    }
213    set min
214}
215
216# ::math::product --
217#
218#	Return the product of one or more values
219#
220# Arguments:
221#	val	first value
222#	args	other values
223#
224# Results:
225#	prod	 product of multiplying all values in the list
226
227proc ::math::product {val args} {
228    set prod $val
229    foreach val $args {
230        set prod [ expr { $prod*$val } ]
231    }
232    set prod
233}
234
235# ::math::random --
236#
237#	Return a random number in a given range.
238#
239# Arguments:
240#	args	optional arguments that specify the range within which to
241#		choose a number:
242#			(null)		choose a number between 0 and 1
243#			val		choose a number between 0 and val
244#			val1 val2	choose a number between val1 and val2
245#
246# Results:
247#	num	a random number in the range.
248
249proc ::math::random {args} {
250    set num [expr {rand()}]
251    if { [llength $args] == 0 } {
252	return $num
253    } elseif { [llength $args] == 1 } {
254	return [expr {int($num * [lindex $args 0])}]
255    } elseif { [llength $args] == 2 } {
256	foreach {lower upper} $args break
257	set range [expr {$upper - $lower}]
258	return [expr {int($num * $range) + $lower}]
259    } else {
260	set fn [lindex [info level 0] 0]
261	error "wrong # args: should be \"$fn ?value1? ?value2?\""
262    }
263}
264
265# ::math::sigma --
266#
267#	Return the standard deviation of three or more values
268#
269# Arguments:
270#	val1	first value
271#	val2	second value
272#	args	other values
273#
274# Results:
275#	sigma	population standard deviation value
276
277proc ::math::sigma {val1 val2 args} {
278     set sum [ expr { $val1+$val2 } ]
279     set N [ expr { [ llength $args ] + 2 } ]
280     foreach val $args {
281        set sum [ expr { $sum+$val } ]
282     }
283     set mean [ expr { $sum/$N } ]
284     set sigma_sq 0
285     foreach val [ concat $val1 $val2 $args ] {
286        set sigma_sq [ expr { $sigma_sq+pow(($val-$mean),2) } ]
287     }
288     set sigma_sq [ expr { $sigma_sq/($N-1) } ]
289     set sigma [ expr { sqrt($sigma_sq) } ]
290     set sigma
291}
292
293# ::math::stats --
294#
295#	Return the mean, standard deviation, and coefficient of variation as
296#	percent, as a list.
297#
298# Arguments:
299#	val1	first value
300#	val2	first value
301#	args	all other values
302#
303# Results:
304#	{mean stddev coefvar}
305
306proc ::math::stats {val1 val2 args} {
307     set sum [ expr { $val1+$val2 } ]
308     set N [ expr { [ llength $args ] + 2 } ]
309     foreach val $args {
310        set sum [ expr { $sum+$val } ]
311     }
312     set mean [ expr { $sum/double($N) } ]
313     set sigma_sq 0
314     foreach val [ concat $val1 $val2 $args ] {
315        set sigma_sq [ expr { $sigma_sq+pow(($val-$mean),2) } ]
316     }
317     set sigma_sq [ expr { $sigma_sq/double($N-1) } ]
318     set sigma [ expr { sqrt($sigma_sq) } ]
319     set cov [ expr { ($sigma/$mean)*100 } ]
320     return [ list $mean $sigma $cov ]
321}
322
323# ::math::sum --
324#
325#	Return the sum of one or more values
326#
327# Arguments:
328#	val	first value
329#	args	all other values
330#
331# Results:
332#	sum	arithmetic sum of all values in args
333
334proc ::math::sum {val args} {
335    set sum $val
336    foreach val $args {
337        set sum [ expr { $sum+$val } ]
338    }
339    set sum
340}
341
342#----------------------------------------------------------------------
343#
344# ::math::expectDouble --
345#
346#	Format an error message that an argument was expected to be
347#	double and wasn't
348#
349# Parameters:
350#	arg -- Misformatted argument
351#
352# Results:
353#	Returns an appropriate error message
354#
355# Side effects:
356#	None.
357#
358#----------------------------------------------------------------------
359
360proc ::math::expectDouble { arg } {
361    return [format "expected a floating-point number but found \"%.50s\"" $arg]
362}
363
364#----------------------------------------------------------------------
365#
366# ::math::expectInteger --
367#
368#	Format an error message that an argument was expected to be
369#	integer and wasn't
370#
371# Parameters:
372#	arg -- Misformatted argument
373#
374# Results:
375#	Returns an appropriate error message
376#
377# Side effects:
378#	None.
379#
380#----------------------------------------------------------------------
381
382proc ::math::expectInteger { arg } {
383    return [format "expected an integer but found \"%.50s\"" $arg]
384}
385
386