Tcl Library Source Code

[manpage_begin math::exact n 1.0]
[copyright "2015 Kevin B. Kenny <[email protected]>
Redistribution permitted under the terms of the Open\
Publication License <http://www.opencontent.org/openpub/>"]
[moddesc {Tcl Math Library}]
[titledesc {Exact Real Arithmetic}]
[category Mathematics]
[require Tcl 8.6]
[require grammar::aycock 1.0]
[require math::exact 1.0]
[description]
[para]
The [cmd exactexpr] command in the [cmd math::exact] package
allows for exact computations over the computable real numbers.
These are not arbitrary-precision calculations; rather they are
exact, with numbers represented by algorithms that produce successive
approximations. At the end of a calculation, the caller can
request a given precision for the end result, and intermediate results are
computed to whatever precision is necessary to satisfy the request.
[section "Procedures"]
The following procedure is the primary entry into the [cmd math::exact]
package.
[list_begin definitions]
[call [cmd ::math::exact::exactexpr] [arg expr]]

Accepts a mathematical expression in Tcl syntax, and returns an object
that represents the program to calculate successive approximations to
the expression's value. The result will be referred to as an
exact real number.

[call [arg number] [cmd ref]]

Increases the reference count of a given exact real number.

[call [arg number] [cmd unref]]

Decreases the reference count of a given exact real number, and destroys
the number if the reference count is zero.

[call [arg number] [cmd asPrint] [arg precision]]

Formats the given [arg number] for printing, with the specified [arg precision].
(See below for how [arg precision] is interpreted). Numbers that are known to
be rational are formatted as fractions.

[call [arg number] [cmd asFloat] [arg precision]]

Formats the given [arg number] for printing, with the specified [arg precision].
(See below for how [arg precision] is interpreted). All numbers are formatted
in floating-point E format.

[list_end]

[section Parameters]

[list_begin definitions]

[def [arg expr]]

Expression to evaluate. The syntax for expressions is the same as it is in Tcl,
but the set of operations is smaller. See [sectref Expressions] below
for details.

[def [arg number]]

The object returned by an earlier invocation of [cmd math::exact::exactexpr]

[def [arg precision]]

The requested 'precision' of the result. The precision is (approximately)
the absolute value of the binary exponent plus the number of bits of the
binary significand. For instance, to return results to IEEE-754 double
precision, 56 bits plus the exponent are required. Numbers between 1/2 and 2
will require a precision of 57; numbers between 1/4 and 1/2 or between 2 and 4
will require 58; numbers between 1/8 and 1/4 or between 4 and 8 will require
59; and so on.

[list_end]

[section Expressions]

The [cmd math::exact::exactexpr] command accepts expressions in a subset
of Tcl's syntax. The following components may be used in an expression.

[list_begin itemized]

[item]Decimal integers.
[item]Variable references with the dollar sign ([const \$]).
The value of the variable must be the result of another call to
[cmd math::exact::exactexpr]. The reference count of the value
will be increased by one for each position at which it appears
in the expression.
[item]The exponentiation operator ([const **]).
[item]Unary plus ([const +]) and minus ([const -]) operators.
[item]Multiplication ([const *]) and division ([const /]) operators.
[item]Parentheses used for grouping.
[item]Functions. See [sectref Functions] below for the functions that are
available.

[list_end]

[section Functions]

The following functions are available for use within exact real expressions.

[list_begin definitions]


[def [const acos(][arg x][const )]]
The inverse cosine of [arg x]. The result is expressed in radians. 
The absolute value of [arg x] must be less than 1.

[def [const acosh(][arg x][const )]]
The inverse hyperbolic cosine of [arg x]. 
[arg x] must be greater than 1.

[def [const asin(][arg x][const )]]
The inverse sine of [arg x]. The result is expressed in radians. 
The absolute value of [arg x] must be less than 1.

[def [const asinh(][arg x][const )]]
The inverse hyperbolic sine of [arg x].

[def [const atan(][arg x][const )]]
The inverse tangent of [arg x]. The result is expressed in radians.

[def [const atanh(][arg x][const )]]
The inverse hyperbolic tangent of [arg x].
The absolute value of [arg x] must be less than 1.

[def [const cos(][arg x][const )]]
The cosine of [arg x]. [arg x] is expressed in radians.

[def [const cosh(][arg x][const )]]
The hyperbolic cosine of [arg x].

[def [const e()]]
The base of the natural logarithms = [const 2.71828...]

[def [const exp(][arg x][const )]]
The exponential function of [arg x].

[def [const log(][arg x][const )]]
The natural logarithm of [arg x]. [arg x] must be positive.

[def [const pi()]]
The value of pi = [const 3.15159...]

[def [const sin(][arg x][const )]]
The sine of [arg x]. [arg x] is expressed in radians.

[def [const sinh(][arg x][const )]]
The hyperbolic sine of [arg x].

[def [const sqrt(][arg x][const )]]
The square root of [arg x]. [arg x] must be positive.

[def [const tan(][arg x][const )]]
The tangent of [arg x]. [arg x] is expressed in radians.

[def [const tanh(][arg x][const )]]
The hyperbolic tangent of [arg x].

[list_end]

[section Summary]

The [cmd math::exact::exactexpr] command provides a system that
performs exact arithmetic over computable real numbers, representing
the numbers as algorithms for successive approximation.

An example, which implements the high-school quadratic formula,
is shown below.

[example {
namespace import math::exact::exactexpr
proc exactquad {a b c} {
    set d [[exactexpr {sqrt($b*$b - 4*$a*$c)}] ref]
    set r0 [[exactexpr {(-$b - $d) / (2 * $a)}] ref]
    set r1 [[exactexpr {(-$b + $d) / (2 * $a)}] ref]
    $d unref
    return [list $r0 $r1]
}

set a [[exactexpr 1] ref]
set b [[exactexpr 200] ref]
set c [[exactexpr {(-3/2) * 10**-12}] ref]
lassign [exactquad $a $b $c] r0 r1
$a unref; $b unref; $c unref
puts [list [$r0 asFloat 70] [$r1 asFloat 110]]
$r0 unref; $r1 unref
}]

The program prints the result:
[example {
-2.000000000000000075e2 7.499999999999999719e-15
}]

Note that if IEEE-754 floating point had been used, a catastrophic
roundoff error would yield a smaller root that is a factor of two
too high:

[example {
-200.0 1.4210854715202004e-14
}]

The invocations of [cmd exactexpr] should be fairly self-explanatory.
The other commands of note are [cmd ref] and [cmd unref]. It is necessary
for the caller to keep track of references to exact expressions - to call
[cmd ref] every time an exact expression is stored in a variable and
[cmd unref] every time the variable goes out of scope or is overwritten.

The [cmd asFloat] method emits decimal digits as long as the requested
precision supports them. It terminates when the requested precision
yields an uncertainty of more than one unit in the least significant digit.

[vset CATEGORY mathematics]
[manpage_end]

# exact.tcl --
#
#	Tcl package for exact real arithmetic.
#
# Copyright (c) 2015 by Kevin B. Kenny
#
# See the file "license.terms" for information on usage and redistribution of
# this file, and for a DISCLAIMER OF ALL WARRANTIES.
#
# This package provides a library for performing exact
# computations over the computable real numbers. The algorithms
# are largely based on the ones described in:
#
# Potts, Peter John. _Exact Real Arithmetic using Möbius Transformations._
# PhD thesis, University of London, July 1998.
# http://www.doc.ic.ac.uk/~ae/papers/potts-phd.pdf
#
# Some of the algorithms for the elementary functions are found instead
# in:
#
# Menissier-Morain, Valérie. _Arbitrary Precision Real Arithmetic:
# Design and Algorithms. J. Symbolic Computation 11 (1996)
# http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.31.8983
# 
#-----------------------------------------------------------------------------

package require Tcl 8.6
package require grammar::aycock 1.0

namespace eval math::exact {
    
    namespace eval function {
	namespace path ::math::exact
    }
    namespace path ::tcl::mathop

    # math::exact::parser --
    #
    #	Grammar for parsing expressions in the exact real calculator
    #
    # The expression syntax is almost exactly that of Tcl expressions,
    # minus Tcl arrays, square-bracket substitution, and noncomputable
    # operations such as equality, comparisons, bit and Boolean operations,
    # and ?:.

    variable parser [grammar::aycock::parser {

	target ::= expression {
	    lindex $_ 0
	}

	expression ::= expression addop term {
	    {*}$_
	}
	expression ::= term {
	    lindex $_ 0
	}
	addop ::= + {
	    lindex $_ 0
	}
	addop ::= - {
	    lindex $_ 0
	}

	term ::= term mulop factor {
	    {*}$_
	}
	term ::= factor {
	    lindex $_ 0
	}
	mulop ::= * {
	    lindex $_ 0
	}
	mulop ::= / {
	    lindex $_ 0
	}

	factor ::= addop factor {
	    switch -exact -- [lindex $_ 0] {
		+ {
		    set result [lindex $_ 1]
		}
		- {
		    set result [[lindex $_ 1] U-]
		}
	    }
	    set result
	}
	factor ::= primary ** factor {
	    {*}$_
	}
	factor ::= primary {
	    lindex $_ 0
	}

	primary ::= {$} bareword {
	    uplevel [dict get $clientData caller] set [lindex $_ 1]
	}
	primary ::= number {
	    [dict get $clientData namespace]::V new [list [lindex $_ 0] 1]
	}
	primary ::= bareword ( ) {
	    [dict get $clientData namespace]::function::[lindex $_ 0]
	}	    
	primary ::= bareword ( arglist ) {
	    [dict get $clientData namespace]::function::[lindex $_ 0] \
		{*}[lindex $_ 2]
	}
	primary ::= ( expression ) {
	    lindex $_ 1
	}
	arglist ::= expression {
	    set _
	}
	arglist ::= arglist , expression {
	    linsert [lindex $_ 0] end [lindex $_ 2]
	}

    }]
}

# math::exact::Lexer --
#
#	Lexer for the arithmetic expressions that the 'math::exact' package
#	can evaluate.
#
# Results:
#	Returns a two element list. The first element is a list of the
#	lexical values of the tokens that were found in the expression;
#	the second is a list of the semantic values of the tokens. The
#	two sublists are the same length.

proc math::exact::Lexer {expression} {
    set start 0
    set tokens {}
    set values {}
    while {$expression ne {}} {
	if {[regexp {^\*\*(.*)} $expression -> rest]} {

	    # Exponentiation

	    lappend tokens **
	    lappend values **
	} elseif {[regexp {^([-+/*$(),])(.*)} $expression -> token rest]} {

	    # Single-character operators

	    lappend tokens $token
	    lappend values $token
	} elseif {[regexp {^([[:alpha:]][[:alnum:]_]*)(.*)} \
		       $expression -> token rest]} {

	    # Variable and function names

	    lappend tokens bareword
	    lappend values $token
	} elseif {[regexp -nocase {^([[:digit:]]+)(.*)} $expression -> \
		       token rest] } {

	    # Numbers

	    lappend tokens number
	    lappend values $token

	} elseif {[regexp {^[[:space:]]+(.*)} $expression -> rest]} {

	    # Whitespace

	} else {

	    # Anything else is an error

	    return -code error \
		-errorcode [list MATH EXACT EXPR INVCHAR \
				[string index $expression 0]] \
		[list invalid character [string index $expression 0]] \
	}
	set expression $rest
    }
    return [list $tokens $values]
}

# math::exact::K --
#
#	K combinator. Returns its first argumetn
#
# Parameters:
#	a - Return value
#	b - Value to discard
#
# Results:
#	Returns the first argument

proc math::exact::K {a b} {return $a}

# math::exact::exactexpr --
#
#	Evaluates an exact real expression.
#
# Parameters:
#	expr - Expression to evaluate. Variables in the expression are
#	       assumed to be reals, which are represented as Tcl objects.
#
# Results:
#	Returns a Tcl object representing the expression's value.
#
# The returned object must have its refcount incremented with [ref] if
# the caller retains a reference, and in general it is expected that a
# user of a real will [ref] the object when storing it in a variable and
# [unref] it again when the variable goes out of scope or is overwritten.

proc math::exact::exactexpr {expr} {
    variable parser
    set result [$parser parse {*}[Lexer $expr] \
		    [dict create \
			 caller "#[expr {[info level] - 1}]" \
			 namespace [namespace current]]]
}

# Basic data types

# A vector is a list {a b}. It can represent the rational number {a/b}

# A matrix is a list of its columns {{a b} {c d}}. In addition to
# the ordinary rules of linear algebra, it represents the linear
# transform (ax+b)/(cx+d). 

# If x is presumed to lie in the interval [0, Inf) then this transform 
# applied to x will lie in the interval [b/d, a/c), so the matrix 
# {{a b} {c d}} can represent that interval. The interval [0,Inf)
# can be represented by the identity matrix {{1 0} {0 1}}

# Moreover, if x  = {p/q} is a rational number, then 
#    (ax+b)/(cx+d) = (a(p/q)+b)/(c(p/q)+d)
#                  = ((ap+bq)/q)/(cp+dq)/q)
#                  = (ap+bq)/(cp+dq)
# which is the rational number represented by {{a c} {b d}} {p q}
# using the conventional rule of vector-matrix multiplication.

# Note that matrices used for this purpose are unique only up to scaling.
# If (ax+b)/(cx+d) is a rational number, then (eax+eb)/(ecx+ed) represents
# the same rational number. This means that matrix inversion may be replaced
# by matrix reversion: for {{a b} {c d}}, simply form the list of cofactors
# {{d -b} {-c a}}, without dividing by the determinant. The reverse of a matrix
# is well defined even if the matrix is singular.

# A tensor of the third degree is a list of its levels:
#  {{{a b} {c d}} {{e f} {g h}}}

# math::exact::gcd --
#
#	Greatest common divisor of a set of integers
#
# Parameters:
#	The integers whose gcd is to be found
#
# Results:
#	Returns the gcd

proc math::exact::gcd {a args} {
    foreach b $args {
	if {$a > $b} {
	    set t $b; set b $a; set a $t
	}
	while {$b > 0} {
	    set t $b
	    set b [expr {$a % $b}]
	    set a $t
	}
    }
    return $a
}

# math::exact::trans --
#
#	Transposes a 2x2 matrix or a 2x2x2 tensor
#
# Parameters:
#	x - Object to transpose
#
# Results:
#	Returns the transpose

proc math::exact::trans {x} {
    lassign $x ab cd
    lassign $ab a b
    lassign $cd c d
    tailcall list [list $a $c] [list $b $d]
}

# math::exact::determinant --
#
#	Calculates the determinant of a 2x2 matrix
#
# Parameters:
#	x - Matrix
#
# Results:
#	Returns the determinant.

proc math::exact::determinant {x} {
    lassign $x ab cd
    lassign $ab a b
    lassign $cd c d
    return [expr {$a*$d - $b*$c}]
}

# math::exact::reverse --
#
#	Calculates the reverse of a 2x2 matrix, which is its inverse times
#	its determinant.
#
# Parameters:
#	x - Matrix
#
# Results:
#	Returns reverse[x].
#
# Notes:
#	The reverse is well defined even for singular matrices.

proc math::exact::reverse {x} {
    lassign $x ab cd
    lassign $ab a b
    lassign $cd c d
    tailcall list [list $d [expr {-$b}]] [list [expr {-$c}] $a]
}

# math::exact::veven --
#
#	Tests if both components of a 2-vector are even.
#
# Parameters:
#	x - Vector to test
#
# Results:
#	Returns 1 if both components are even, 0 otherwise.

proc math::exact::veven {x} {
    lassign $x a b
    return [expr {($a % 2 == 0) && ($b % 2 == 0)}]
}

# math::exact::meven --
#
#	Tests if all components of a 2x2 matrix are even.
#
# Parameters:
#	x - Matrix to test
#
# Results:
#	Returns 1 if all components are even, 0 otherwise.

proc math::exact::meven {x} {
    lassign $x a b
    return [expr {[veven $a] && [veven $b]}]
}

# math::exact::teven --
#
#	Tests if all components of a 2x2x2 tensor are even
#
# Parameters:
#	x - Tensor to test
#
# Results:
#	Returns 1 if all components are even, 0 otherwise

proc math::exact::teven {x} {
    lassign $x a b
    return [expr {[meven $a] && [meven $b]}]
}

# math::exact::vhalf --
#
#	Divides both components of a 2-vector by 2
#
# Parameters:
#	x - Vector to scale
#
# Results:
#	Returns the scaled vector

proc math::exact::vhalf {x} {
    lassign $x a b
    tailcall list [expr {$a / 2}] [expr {$b / 2}]
}

# math::exact::mhalf --
#
#	Divides all components of a 2x2 matrix by 2
#
# Parameters:
#	x - Matrix to scale
#
# Results:
#	Returns the scaled matrix

proc math::exact::mhalf {x} {
    lassign $x a b
    tailcall list [vhalf $a] [vhalf $b]
}

# math::exact::thalf --
#
#	Divides all components of a 2x2x2 tensor by 2
#
# Parameters:
#	x - Tensor to scale
#
# Results:
#	Returns the scaled tensor

proc math::exact::thalf {x} {
    lassign $x a b
    tailcall list [mhalf $a] [mhalf $b]
}

# math::exact::vscale --
#
#	Removes all common factors of 2 from the two components of a 2-vector
#
# Paramters:
#	x - Vector to scale
#
# Results:
#	Returns the scaled vector

proc math::exact::vscale {x} {
    while {[veven $x]} {
	set x [vhalf $x]
    }
    return $x
}

# math::exact::mscale --
#
#	Removes all common factors of 2 from the two components of a
#	2x2 matrix
#
# Paramters:
#	x - Matrix to scale
#
# Results:
#	Returns the scaled matrix

proc math::exact::mscale {x} {
    while {[meven $x]} {
	set x [mhalf $x]
    }
    return $x
}

# math::exact::tscale --
#
#	Removes all common factors of 2 from the two components of a
#	2x2x2 tensor
#
# Paramters:
#	x - Tensor to scale
#
# Results:
#	Returns the scaled tensor

proc math::exact::tscale {x} {
    while {[teven $x]} {
	set x [thalf $x]
    }
    return $x
}

# math::exact::vreduce --
#
#	Reduces a vector (i.e., a rational number) to lowest terms
#
# Parameters:
#	x - Vector to scale
#
# Results:
#	Returns the scaled vector

proc math::exact::vreduce {x} {
    lassign $x a b
    set g [gcd $a $b]
    tailcall list [expr {$a / $g}] [expr {$b / $g}]
}

# math::exact::mreduce --
#
#	Removes all common factors from the two components of a
#	2x2 matrix
#
# Paramters:
#	x - Matrix to scale
#
# Results:
#	Returns the scaled matrix
#
# This procedure suffices to reduce the matrix to lowest terms if the matrix
# was constructed by pre- or post-multiplying a series of sign and digit
# matrices.

proc math::exact::mreduce {x} {
    lassign $x ab cd
    lassign $ab a b
    lassign $cd c d
    set g [gcd $a $b $c $d]
    tailcall list \
	[list [expr {$a / $g}] [expr {$b / $g}]] \
	[list [expr {$c / $g}] [expr {$d / $g}]]
}

# math::exact::treduce --
#
#	Removes all common factors from the components of a
#	2x2x2 tensor
#
# Paramters:
#	x - Tensor to scale
#
# Results:
#	Returns the scaled tensor
#
# This procedure suffices to reduce a tensor to lowest terms if it was 
# constructed by absorbing a digit matrix into a tensor that was already
# in lowest terms.

proc math::exact::treduce {x} {
    lassign $x abcd efgh
    lassign $abcd ab cd
    lassign $ab a b
    lassign $cd c d
    lassign $efgh ef gh
    lassign $ef e f
    lassign $gh g h
    set G [gcd $a $b $c $d $e $f $g $h]
    tailcall list \
	[list \
	     [list [expr {$a / $G}] [expr {$b / $G}]] \
	     [list [expr {$c / $G}] [expr {$d / $G}]]] \
	[list \
	     [list [expr {$e / $G}] [expr {$f / $G}]] \
	     [list [expr {$g / $G}] [expr {$h / $G}]]]
}

# math::exact::vadd --
#
#	Adds two 2-vectors
#
# Parameters:
#	x - First vector
#	y - Second vector
#
# Results:
#	Returns the vector sum

proc math::exact::vadd {x y} {
    lmap p $x q $y {expr {$p + $q}}
}

# math::exact::madd --
#
#	Adds two 2x2 matrices
#
# Parameters:
#	A - First matrix
#	B - Second matrix
#
# Results:
#	Returns the matrix sum

proc math::exact::madd {A B} {
    lmap x $A y $B {
	lmap p $x q $y {expr {$p + $q}}
    }
}

# math::exact::tadd --
#
#	Adds two 2x2x2 tensors
#
# Parameters:
#	U - First tensor
#	V - Second tensor
#
# Results:
#	Returns the tensor sum

proc math::exact::tadd {U V} {
    lmap A $U B $V {
	lmap x $A y $B {
	    lmap p $x q $y {expr {$p + $q}}
	}
    }
}

# math::exact::mdotv --
#
#	2x2 matrix times 2-vector
#
# Parameters;
#	A - Matrix
#	x - Vector
# 
# Results:
#	Returns the product vector

proc math::exact::mdotv {A x} {
    lassign $A ab cd
    lassign $ab a b
    lassign $cd c d
    lassign $x e f
    tailcall list [expr {$a*$e + $c*$f}] [expr {$b*$e + $d*$f}]
}

# math::exact::mdotm --
#
#	Product of two matrices
#
# Parameters:
#	A - Left matrix
#	B - Right matrix
#
# Results:
#	Returns the matrix product

proc math::exact::mdotm {A B} {
    lassign $B x y
    tailcall list [mdotv $A $x] [mdotv $A $y]
}

# math::exact::mdott --
#
#	Product of a matrix and a tensor
#
# Parameters:
#	A - Matrix
#	T - Tensor
#
# Results:
#	Returns the product tensor

proc math::exact::mdott {A T} {
    lassign $T B C
    tailcall list [mdotm $A $B] [mdotm $A $C]
}

# math::exact::trightv --
#
#	Right product of a tensor and a vector
#
# Parameters:
#	T - Tensor
#	v - Right-hand vector
#
# Results:
#	Returns the product matrix

proc math::exact::trightv {T v} {
    lassign $T m n
    tailcall list [mdotv $m $v] [mdotv $n $v]
}

# math::exact::trightm --
#
#	Right product of a tensor and a matrix
#
# Parameters:
#	T - Tensor
#	A - Right-hand matrix
#
# Results:
#	Returns the product tensor

proc math::exact::trightm {T A} {
    lassign $T m n
    tailcall list [mdotm $m $A] [mdotm $n $A]
}

# math::exact::tleftv --
#
#	Left product of a tensor and a vector
#
# Parameters:
#	T - Tensor
#	v - Left-hand vector
#
# Results:
#	Returns the product matrix

proc math::exact::tleftv {T v} {
    tailcall trightv [trans $T] $v
}

# math::exact::tleftm --
#
#	Left product of a tensor and a matrix
#
# Parameters:
#	T - Tensor
#	A - Left-hand matrix
#
# Results:
#	Returns the product tensor

proc math::exact::tleftm {T A} {
    tailcall trans [trightm [trans $T] $A]
}

# math::exact::vsign --
#
#	Computes the 'sign function' of a vector.
#
# Parameters:
#	v - Vector whose sign function is needed
#
# Results:
#	Returns the result of the sign function.
#
# a	b	sign
# -	-	 -1
# -	0	 -1
# -	+	  0
# 0	-	 -1
# 0	0	  0
# 0	+	  1
# +	-	  0
# +	0	  1
# +	+	  1
#
# If the quotient a/b is negative or indeterminate, the result is zero.
# If the quotient a/b is zero, the result is the sign of b.
# If the quotient a/b is positive, the result is the common sign of the
# operands, which are known to be of like sign
# If the quotient a/b is infinite, the result is the sign of a.

proc math::exact::sign {v} {
    lassign $v a b
    if {$a < 0} {
	if {$b <= 0} {
	    return -1
	} else {
	    return 0
	}
    } elseif {$a == 0} {
	if {$b < 0} {
	    return -1
	} elseif {$b == 0} {
	    return 0
	} else {
	    return 1
	}
    } else {
	if {$b < 0} {
	    return 0
	} else {
	    return 1
	}
    }
}

# math::exact::vrefines --
#
#	Test if a vector refines.
#
# Parameters:
#	v - Vector to test
#
# Results:
#	1 if the vector refines, 0 otherwise.

proc math::exact::vrefines {v} {
    return [expr {[sign $v] != 0}]
}

# math::exact::mrefines --
#
#	Test whether a matrix refines
#
# Parameters:
#	A - Matrix to test
#
# Results:
#	1 if the matrix refines, 0 otherwise.

proc math::exact::mrefines {A} {
    lassign $A v w
    set a [sign $v]
    set b [sign $w]
    return [expr {$a == $b && $b != 0}]
}

# math::exact::trefines --
#
#	Tests whether a tensor refines
#
# Parameters:
#	T - Tensor to test.
#
# Results:
#	1 if the tensor refines, 0 otherwise.

proc math::exact::trefines {T} {
    lassign $T vw xy
    lassign $vw v w
    lassign $xy x y
    set a [sign $v]
    set b [sign $w]
    set c [sign $x]
    set d [sign $y]
    return [expr {$a == $b && $b == $c && $c == $d && $d != 0}]
}

# math::exact::vlessv -
#
#	Test whether one rational is less than another
#
# Parameters:
#	v, w - Two rational numbers
#
# Returns:
#	The result of the comparison.

proc math::exact::vlessv {v w} {
    expr {[determinant [list $v $w]] < 0}
}

# math::exact::mlessv -
#
#	Tests whether a rational interval is less than a vector
#
# Parameters:
#	m - Matrix representing the interval
#	x - Rational to compare against
#
# Results:
#	Returns 1 if m < x, 0 otherwise

proc math::exact::mlessv {m x} {
    lassign $m v w
    expr {[vlessv $v $x] && [vlessv $w $x]}
}

# math::exact::mlessm -
#
#	Tests whether one rational interval is strictly less than another
#
# Parameters:
#	m - First interval
#	n - Second interval
#
# Results:
#	Returns 1 if m < n, 0 otherwise

proc math::exact::mlessm {m n} {
    lassign $n v w
    expr {[mlessv $m $v] && [mlessv $m $w]}
}

# math::exact::mdisjointm -
#
#	Tests whether two rational intervals are disjoint
#
# Parameters:
#	m - First interval
#	n - Second interval
#
# Results:
#	Returns 1 if the intervals are disjoint, 0 otherwise

proc math::exact::mdisjointm {m n} {
    expr {[mlessm $m $n] || [mlessm $n $m]}
}

# math::exact::mAsFloat
#
#	Formats a matrix that represents a rational interval as a floating 
#	point number, stopping as soon as a digit is not determined.
#
# Parameters:
#	m - Matrix to format
#
# Results:
#	Returns the floating point number in scientific notation, with no
#	digits to the left of the decimal point.

proc math::exact::mAsFloat {m} {

    # Special case: If a number is exact, the determinant is zero.

    set d [determinant $m]
    lassign [lindex $m 0] p q
    if {$d == 0} {
	if {$q < 0} {
	    set p [expr {-$p}]
	    set q [expr {-$q}]
	}
	if {$p == 0} {
	    if {$q == 0} {
		return NaN
	    } else {
		return 0
	    }
	} elseif {$q == 0} {
	    return Inf
	} elseif {$q == 1} {
	    return $p
	} else {
	    set G [gcd $p $q]
	    return [expr {$p/$G}]/[expr {$q/$G}]
	}
    } else {
	tailcall eFormat [scientificNotation $m]
    }
}

# math::exact::scientificNotation --
#
#	Takes a matrix representing a rational interval, and extracts as
#	many decimal digits as can be determined unambiguously
#
# Parameters:
#	m - Matrix to format
#
# Results:
#	Returns a list comprising the decimal exponent, followed by a series of
#	digits of the significand. The decimal point is to the left of the
#	leftmost digit of the significand.
#
#	Returns the empty string if a number is entirely undetermined.

proc math::exact::scientificNotation {m} {
    set n 0
    while {1} {
	if {[vrefines [mdotv [reverse $m] {1 0}]]} {
	    return {}
	} elseif {[mrefines [mdotm $math::exact::iszer $m]]} {
	    return [linsert [mantissa $m] 0 $n]
	} else {
	    set m [mdotm {{1 0} {0 10}} $m]
	    incr n
	}
    }
}

# math::exact::mantissa --
#
#	Given a matrix m that represents a rational interval whose
#	endpoints are in [0,1), formats as many digits of the represented
#	number as possible.
#
# Parameters:
#	m - Matrix to format
#
# Results:
#	Returns a list of digits

proc math::exact::mantissa {m} {
    set retval {}
    set done 0
    while {!$done} {
	set done 1

	# Brute force: try each digit in turn. This could no doubt be
	# improved on.

	for {set j -9} {$j <= 9} {incr j} {
	    set digitMatrix \
		[list [list [expr {$j+1}] 10] [list [expr {$j-1}] 10]]
	    if {[mrefines [mdotm [reverse $digitMatrix] $m]]} {
		lappend retval $j
		set nextdigit [list {10 0} [list [expr {-$j}] 1]]
		set m [mdotm $nextdigit $m]
		set done 0
		break
	    }
	}
    }
    return $retval
}

# math::exact::eFormat --
#
#	Formats a decimal exponent and significand in E format
#
# Parameters:
#	expAndDigits - List whose first element is the exponent and
#		       whose remaining elements are the digits of the
#		       significand.

proc math::exact::eFormat {expAndDigits} {

    # An empty sequence of digits is an indeterminate number

    if {[llength $expAndDigits] < 2} {
	return Undetermined
    }
    set significand [lassign $expAndDigits exponent]

    # Accumulate the digits
    set v 0
    foreach digit $significand {
	set v [expr {10 * $v + $digit}]
    }

    # Adjust the exponent if the significand has too few digits.

    set l [llength $significand]
    while {$l > 0 && abs($v) < 10**($l-1)} {
	incr l -1
	incr exponent -1
    }
    incr exponent -1

    # Put in the sign

    if {$v < 0} {
	set result -
	set v [expr {-$v}]
    } else {
	set result {}
    }

    # Put in the significand with the decimal point after the leading digit.

    if {$v == 0} {
	append result 0
    } else {
	append result [string index $v 0] . [string range $v 1 end]
    }

    # Put in the exponent

    append result e $exponent

    return $result
}

# math::exact::showRat --
#
#	Formats an exact rational for printing in E format.
#
# Parameters:
#	v - Two-element list of numerator and denominator.
#
# Results:
#	Returns the quotient in E format.  Nonzero/zero == Infinity,
#	0/0 == NaN.

proc math::exact::showRat {v} {
    lassign $v p q
    if {$p != 0 || $q != 0} {
	return [format %e [expr {double($p)/double($q)}]]
    } else {
	return NaN
    }
}

# math::exact::showInterval --
#
#	Formats a rational interval for printing
#
# Parameters:
#	m - Matrix representing the interval
#
# Results:
#	Returns a string representing the interval in E format.

proc math::exact::showInterval {m} {
    lassign $m v w
    return "\[[showRat $w] .. [showRat $v]\]"
}

# math::exact::showTensor --
#
#	Formats a tensor for printing
#
# Parameters:
#	t - Tensor to print
#
# Results:
#	Returns a string containing the left and right matrices of the
#	tensor, each represented as an interval.

proc math::exact::showTensor {t} {
    lassign $t m n
    return [list [showInterval $m] [showInterval $n]]
}

# math::exact::counted --
#
#	Reference counted object

oo::class create math::exact::counted {
    variable refcount_

    # Constructor builds an object with a zero refcount.
    constructor {} {
	if 0 {
	    puts {}
	    puts "construct: [self object] refcount now 0"
	    for {set i [info frame]} {$i > 0} {incr i -1} {
		set frame [info frame $i]
		if {[dict get $frame type] eq {source}} {
		    set line [dict get $frame line]
		    puts "\t[file tail [dict get $frame file]]:$line"
		    if {$line < 0} {
			if {[dict exists $frame proc]} {
			    puts "\t\t[dict get $frame proc]"
			}
			puts "\t\t\[[dict get $frame cmd]\]"
		    }
		} else {
		    puts $frame
		}
	    }
	}
	incr refcount_
	set refcount_ 0
    }

    # The 'ref' method adds a reference to this object, and returns this object
    method ref {} {
	if 0 {
	    puts {}
	    puts "ref: [self object] refcount now [expr {$refcount_ + 1}]"
	    if {$refcount_ == 0} {
		puts "\t[my dump]"
	    }
	    for {set i [info frame]} {$i > 0} {incr i -1} {
		set frame [info frame $i]
		if {[dict get $frame type] eq {source}} {
		    set line [dict get $frame line]
		    puts "\t[file tail [dict get $frame file]]:$line"
		    if {$line < 0} {
			if {[dict exists $frame proc]} {
			    puts "\t\t[dict get $frame proc]"
			}
			puts "\t\t\[[dict get $frame cmd]\]"
		    }
		} else {
		    puts $frame
		}
	    }
	}
	incr refcount_
	return [self]
    }

    # The 'unref' method removes a reference from this object, and destroys
    # this object if the refcount becomes nonpositive.
    method unref {} {
	if 0 {
	    puts {}
	    puts "unref: [self object] refcount now [expr {$refcount_ - 1}]"
	    for {set i [info frame]} {$i > 0} {incr i -1} {
		set frame [info frame $i]
		if {[dict get $frame type] eq {source}} {
		    set line [dict get $frame line]
		    puts "\t[file tail [dict get $frame file]]:$line"
		    if {$line < 0} {
			if {[dict exists $frame proc]} {
			    puts "\t\t[dict get $frame proc]"
			}
			puts "\t\t\[[dict get $frame cmd]\]"
		    }
		}
	    }
	}

	# Destroying this object can result in a long chain of object
	# destruction and eventually a stack overflow. Instead of destroying
	# immediately, list the objects to be destroyed in
	# math::exact::deleteStack, and destroy them only from the outermost
	# stack level that's running 'unref'.

	if {[incr refcount_ -1] <= 0} {
	    variable ::math::exact::deleteStack

	    # Is this the outermost level?
	    set queueActive [expr {[info exists deleteStack]}]

	    # Schedule this object's destruction
	    lappend deleteStack [self object]
	    if {!$queueActive} {

		# At outermost level, destroy all scheduled objects.
		# Destroying one may schedule another.
		while {[llength $deleteStack] != 0} {
		    set obj [lindex $deleteStack end]
		    set deleteStack \
			[lreplace $deleteStack[set deleteStack {}] end end]
		    $obj destroy
		}

		# Once everything quiesces, delete the list.
		unset deleteStack
	    }
	}
    }

    # The 'refcount' method returns the reference count of this object for
    # debugging.
    method refcount {} {
	return $refcount_
    }

    destructor {
    }
}

# An expression is a vector, a matrix applied to an expression,
# or a tensor applied to two expressions. The inner expressions
# may be constructed lazily.

oo::class create math::exact::Expression {
    superclass math::exact::counted

    # absorbed_, signAndMagnitude_, and leadingDigitAndRest_
    # memoize the return values of the 'absorb', 'getSignAndMagnitude',
    # and 'getLeadingDigitAndRest' methods.

    variable absorbed_ signAndMagnitude_ leadingDigitAndRest_

    # Constructor initializes refcount
    constructor {} {
	next
    }

    # Destructor releases memoized objects
    destructor {
	if {[info exists signAndMagnitude_]} {
	    [lindex $signAndMagnitude_ 1] unref
	}
	if {[info exists absorbed_]} {
	    $absorbed_ unref
	}
	if {[info exists leadingDigitAndRest_]} {
	    [lindex $leadingDigitAndRest_ 1] unref
	}
	next
    }

    # getSignAndMagnitude returns a two-element list:
    # the sign matrix, which is one of ispos, isneg, isinf, and iszer,
    # the magnitude, which is another exact real.
    method getSignAndMagnitude {} {
	if {![info exists signAndMagnitude_]} {
	    if {[my refinesM $::math::exact::ispos]} {
		set signAndMagnitude_ \
		    [list $::math::exact::spos \
			 [[my applyM $::math::exact::ispos] ref]]
	    } elseif {[my refinesM $::math::exact::isneg]} {
		set signAndMagnitude_ \
		    [list $::math::exact::sneg \
			 [[my applyM $::math::exact::isneg] ref]]
	    } elseif {[my refinesM $::math::exact::isinf]} {
		set signAndMagnitude_ \
		    [list $::math::exact::sinf \
			 [[my applyM $::math::exact::isinf] ref]]
	    } elseif {[my refinesM $::math::exact::iszer]} {
		set signAndMagnitude_ \
		    [list $::math::exact::szer \
			 [[my applyM $::math::exact::iszer] ref]]
	    } else {
		set absorbed_ [my absorb]
		set signAndMagnitude_ [$absorbed_ getSignAndMagnitude]
		[lindex $signAndMagnitude_ 1] ref
	    }
	}
	return $signAndMagnitude_
    }

    # The 'getLeadingDigitAndRest' method accepts a flag for whether
    # a digit must be extracted (1) or a rational number may be returned
    # directly (0). It returns a two-element list: a digit matrix, which
    # is one of $dpos, $dneg or $dzer, and an exact real representing
    # the number by which the given digit matrix must be postmultiplied.
    method getLeadingDigitAndRest {needDigit} {
	if {![info exists leadingDigitAndRest_]} {
	    if {[my refinesM $::math::exact::idpos]} {
		set leadingDigitAndRest_ \
		    [list $::math::exact::dpos \
			 [[my applyM $::math::exact::idpos] ref]]
	    } elseif {[my refinesM $::math::exact::idneg]} {
		set leadingDigitAndRest_ \
		    [list $::math::exact::dneg \
			 [[my applyM $::math::exact::idneg] ref]]
	    } elseif {[my refinesM $::math::exact::idzer]} {
		set leadingDigitAndRest_ \
		    [list $::math::exact::dzer \
			 [[my applyM $::math::exact::idzer] ref]]
	    } else {
		set absorbed_ [my absorb]
		set newval $absorbed_
		$newval ref
		set leadingDigitAndRest_ \
		    [$newval getLeadingDigitAndRest $needDigit]
		if {[llength $leadingDigitAndRest_] >= 2} {
		    [lindex $leadingDigitAndRest_ 1] ref
		}
		$newval unref
	    }
	}
	return $leadingDigitAndRest_
    }

    # getInterval --
    #	Accumulates 'nDigits' digit matrices, and returns their product,
    #	which is a matrix representing the interval that the digits represent.
    method getInterval {nDigits} {
	lassign [my getSignAndMagnitude] interval e
	$e ref
	lassign [$e getLeadingDigitAndRest 1] ld f
	set interval [math::exact::mdotm $interval $ld]
	$f ref; $e unref; set e $f
	set d $ld
	while {[incr nDigits -1] > 0} {
	    lassign [$e getLeadingDigitAndRest 1] d f
	    set interval [math::exact::mdotm $interval $d]
	    $f ref; $e unref; set e $f
	}
	$e unref
	return $interval
    }

    # asReal --
    #	Coerces an object from rational to real
    #
    # Parameters:
    #	None.
    #
    # Results:
    #	Returns this object
    method asReal {} {self object}

    # asFloat --
    #	Represents this number in E format, after accumulating 'relDigits'
    #	digit matrices.
    method asFloat {relDigits} {
	set v [[my asReal] ref]
	set result [math::exact::mAsFloat [$v getInterval $relDigits]]
	$v unref
	return $result
    }

    # asPrint --
    #	Represents this number for printing. Represents rationals exactly,
    #   others after accumulating 'relDigits' digit matrices.
    method asPrint {relDigits} {
	tailcall math::exact::mAsFloat [my getInterval $relDigits]
    }
    
    # Derived classes are expected to implement the following methods:
    # method dump {} {  
    #	# Formats the object for debugging 
    #	# Returns the formatted string
    # }
    method dump {} {
	error "[info object class [self object]] does not implement the 'dump' method."
    }

    # method refinesM {m} { 
    #	# Returns 1 if premultiplying by the matrix m refines this object
    #   # Returns 0 otherwise
    # }
    method refinesM {m} {
	error "[info object class [self object]] does not implement the 'refinesM' method."
    }
    
    # method applyM {m} { 
    #	# Premultiplies this object by the matrix m 
    # }
    method applyM {m} {
	error "[info object class [self object]] does not implement the 'applyM' method."
    }

    # method applyTLeft {t r} {
    # 	# Computes the left product of the tensor t with this object, and
    #	# applies the result to the right operand r.
    #	# Returns a new exact real representing the product
    # }
    method applyTLeft {t r} {
	error "[info object class [self object]] does not implement the 'applyTLeft' method."
    }
    
    # method applyTRight {t l} {
    # 	# Computes the right product of the tensor t with this object, and
    #	# applies the result to the left operand l.
    #	# Returns a new exact real representing the product
    # }
    method applyTRight {t l} {
	error "[info object class [self object]] does not implement the 'applyTRight' method."
    }
    
    # method absorb {} {
    #	# Absorbs the next subexpression or digit into this expression
    #	# Returns the result of absorption, which always represents a
    #	# smaller interval than this expression
    # }
    method absorb {} {
	error "[info object class [self object]] does not implement the 'absorb' method."
    }

    # U- --
    #
    #   Unary - applied to this object
    #
    # Results:
    #	Returns the negation.

    method U- {} {
	my ref
	lassign [my getSignAndMagnitude] sA mA
	set m [math::exact::mdotm {{-1 0} {0 1}} $sA]
	set result [math::exact::Mstrict new $m $mA]
	my unref
	return $result
    }; export U-

    # + --
    #	Adds this object to another.
    #
    # Parameters:
    #	r - Right addend
    #
    # Results:
    #	Returns the sum
    #
    # Either object may be rational (an instance of V) or real (any other
    # Expression).
    #
    # This method is a Consumer with respect to the current object and to r.
    # It is a Constructor with respect to its result, returning a zero-ref
    # object.
    
    method + {r} {
	return [$r E+ [self object]]
    }; export +

    # E+ --
    #	Adds two exact reals.
    #
    # Parameters:
    #	l - Left addend
    #
    # Results:
    #	Returns the sum.
    #
    # Neither object is an instance of V (that is, neither is a rational).
    #
    # This method is a Consumer with respect to the current object and to l.
    # It is a Constructor with respect to its result, returning a zero-ref
    # object.
    
    method E+ {l} {
	tailcall math::exact::+real $l [self object]
    }; export E+

    # V+ --
    #	Adds a rational and an exact real
    #
    # Parameters:
    #	l - Left addend
    #
    # Results:
    #	Returns the sum.
    #
    # This method is a Consumer with respect to the current object and to l.
    # It is a Constructor with respect to its result, returning a zero-ref
    # object.
    
    method V+ {l} {
	tailcall math::exact::+real $l [self object]
    }; export V+

    # - --
    #	Subtracts another object from this object
    #
    # Parameters:
    #	r - Subtrahend
    #
    # Results:
    #	Returns the difference
    #
    # Either object may be rational (an instance of V) or real (any other
    # Expression).
    #
    # This method is a Consumer with respect to the current object and to r.
    # It is a Constructor with respect to its result, returning a zero-ref
    # object.
    
    method - {r} {
	return [$r E- [self object]]
    }; export -

    # E- --
    #	Subtracts this exact real from another
    #
    # Parameters:
    #	l - Minuend
    #
    # Results:
    #	Returns the difference
    #
    # Neither object is an instance of V (that is, neither is a rational).
    #
    # This method is a Consumer with respect to the current object and to l.
    # It is a Constructor with respect to its result, returning a zero-ref
    # object.
    
    method E- {l} {
	tailcall math::exact::-real $l [self object]
    }; export E-

    # V- --
    #	Subtracts this exact real from a rational
    #
    # Parameters:
    #	l - Minuend
    #
    # Results:
    #	Returns the difference
    #
    # This method is a Consumer with respect to the current object and to l.
    # It is a Constructor with respect to its result, returning a zero-ref
    # object.
    
    method V- {l} {
	tailcall math::exact::-real $l [self object]
    }; export V-

    # * --
    #	Multiplies this object by another.
    #
    # Parameters:
    #	r - Right argument to the multiplication
    #
    # Results:
    #	Returns the product
    #
    # Either object may be rational (an instance of V) or real (any other
    # Expression).
    #
    # This method is a Consumer with respect to the current object and to r.
    # It is a Constructor with respect to its result, returning a zero-ref
    # object.
    
    method * {r} {
	return [$r E* [self object]]
    }; export *

    # E* --
    #	Multiplies two exact reals.
    #
    # Parameters:
    #	l - Left argument to the multiplication
    #
    # Results:
    #	Returns the product.
    #
    # Neither object is an instance of V (that is, neither is a rational).
    #
    # This method is a Consumer with respect to the current object and to l.
    # It is a Constructor with respect to its result, returning a zero-ref
    # object.
    
    method E* {l} {
	tailcall math::exact::*real $l [self object]
    }; export E*

    # V* --
    #	Multiplies a rational and an exact real
    #
    # Parameters:
    #	l - Left argument to the multiplication
    #
    # Results:
    #	Returns the product.
    #
    # This method is a Consumer with respect to the current object and to l.
    # It is a Constructor with respect to its result, returning a zero-ref
    # object.
    
    method V* {l} {
	tailcall math::exact::*real $l [self object]
    }; export V*

    # / --
    #	Divides this object by another.
    #
    # Parameters:
    #	r - Divisor
    #
    # Results:
    #	Returns the quotient
    #
    # Either object may be rational (an instance of V) or real (any other
    # Expression).
    #
    # This method is a Consumer with respect to the current object and to r.
    # It is a Constructor with respect to its result, returning a zero-ref
    # object.
    
    method / {r} {
	return [$r E/ [self object]]
    }; export /

    # E/ --
    #	Divides two exact reals.
    #
    # Parameters:
    #	l - Dividend
    #
    # Results:
    #	Returns the quotient.
    #
    # Neither object is an instance of V (that is, neither is a rational).
    #
    # This method is a Consumer with respect to the current object and to l.
    # It is a Constructor with respect to its result, returning a zero-ref
    # object.
    
    method E/ {l} {
	tailcall math::exact::/real $l [self object]
    }; export E/

    # V/ --
    #	Divides a rational by an exact real
    #
    # Parameters:
    #	l - Dividend
    #
    # Results:
    #	Returns the product.
    #
    # This method is a Consumer with respect to the current object and to l.
    # It is a Constructor with respect to its result, returning a zero-ref
    # object.
    
    method V/ {l} {
	tailcall math::exact::/real $l [self object]
    }; export V/

    # ** -
    #	Raises an exact real to a power
    #
    # Parameters:
    #	r - Exponent
    #
    # Results:
    #	Returns the power.
    #
    # This method is a Consumer with respect to the current object and to l.
    # It is a Constructor with respect to its result, returning a zero-ref
    # object.

    method ** {r} {
	tailcall $r E** [self object]
    }; export **

    # E** -
    #	Raises an exact real to the power of an exact real
    #
    # Parameters:
    #	l - Base to exponentiate
    #
    # Results:
    #	Returns the power
    #
    # This method is a Consumer with respect to the current object and to l.
    # It is a Constructor with respect to its result, returning a zero-ref
    # object.

    method E** {l} {
	# This doesn't work as a tailcall, because this object could have
	# been destroyed by the time we're trying to invoke the tailcall,
	# and that will keep command names from resolving because the
	# tailcall mechanism will try to find them in the destroyed namespace.
	return [math::exact::function::exp \
		    [my * [math::exact::function::log $l]]]
    }; export E**

    # V** -
    #	Raises a rational to the power of an exact real
    #
    # Parameters:
    #	l - Base to exponentiate
    #
    # Results:
    #	Returns the power
    #
    # This method is a Consumer with respect to the current object and to l.
    # It is a Constructor with respect to its result, returning a zero-ref
    # object.

    method V** {l} {
	# This doesn't work as a tailcall, because this object could have
	# been destroyed by the time we're trying to invoke the tailcall,
	# and that will keep command names from resolving because the
	# tailcall mechanism will try to find them in the destroyed namespace.
	return [math::exact::function::exp \
		    [my * [math::exact::function::log $l]]]
    }; export V**
    
    # sqrt --
    #
    #	Create an expression representing the square root of an exact
    #	real argument.
    #
    # Results:
    #	Returns the square root.
    #
    # This procedure is a Consumer with respect the the argument and a
    # Constructor with respect to the result, returning a zero-reference
    # result.

    method sqrt {} {
	variable ::math::exact::isneg
	variable ::math::exact::idzer
	variable ::math::exact::idneg
	variable ::math::exact::idpos
	
	# The algorithm is a modified Newton-Raphson from the Potts and
	# Menissier-Morain papers. The algorithm for sqrt(x) converges
	# rapidly only if x is close to 1, so we rescale to make sure that
	# x is between 1/3 and 3. Specifically:
	# - if x is known to be negative (that is, if $idneg refines it)
	#   then error.
	# - if x is close to 1, $idzer refines it, and we can calculate the
	#   square root directly.
	# - if x is less than 1, $idneg refines it, and we calculate sqrt(4*x)
	#   and multiply by 1/2.
	# - if x is greater than 1, $idpos refines it, and we calculate
	#   sqrt(x/4) and multiply by 2.
	# - if none of the above hold, we have insufficient information about
	#   the magnitude of x and perform a digit exchange.
	
	my ref
	if {[my refinesM $isneg]} {
	    # Negative argument is an error
	    return -code error -errorcode {MATH EXACT SQRTNEGATIVE} \
		"sqrt of negative argument"
	} elseif {[my refinesM $idzer]} {
	    # Argument close to 1
	    set res [::math::exact::SqrtWorker new [self object]]
	} elseif {[my refinesM $idneg]} {
	    # Small argument - multiply by 4 and halve the square root
	    set y [[my applyM {{4 0} {0 1}}] ref]
	    set z [[$y sqrt] ref]
	    set res [$z applyM {{1 0} {0 2}}]
	    $z unref
	    $y unref
	} elseif {[my refinesM $idpos]} {
	    # Large argument - divide by 4 and double the square root
	    set y [[my applyM {{1 0} {0 4}}] ref]
	    set z [[$y sqrt] ref]
	    set res [$z applyM {{2 0} {0 1}}]
	    $z unref
	    $y unref
	} else {
	    # Unclassified argyment. Perform a digit exchange and try again.
	    set y [[my absorb] ref]
	    set res [$y sqrt]
	    $y unref
	}
	my unref
	return $res
    }
}

# math::exact::V --
#	Vector object
#
# A vector object represents a rational number. It is always strict; no
# methods need to perform lazy evaluation.

oo::class create math::exact::V {
    superclass math::exact::Expression

    # v_ is the vector represented.
    variable v_

    # Constructor accepts the vector as a two-element list {n d}
    # where n is by convention the numerator and d the denominator.
    # It is expected that either n or d is nonzero, and that gcd(n,d) == 0.
    # It is also expected that the fraction will be in lowest terms.
    constructor {v} {
	next
	set v_ $v
    }

    # Destructor need only update reference counts
    destructor {
	next
    }
    
    # If a rational is acceptable, getLeadingDigitAndRest may simply return 
    # this object.
    method getLeadingDigitAndRest {needDigit} {
	if {$needDigit} {
	    return [next $needDigit]
	} else {
	    # Note that the result MUST NOT be memoized, as that would lead
	    # to a circular reference, breaking the refcount system.
	    return [self object]
	}
    }

    # Print this object
    method dump {} {
	return "V($v_)"
    }

    # Test if the vector refines when premultiplied by a matrix
    method refinesM {m} {
	return [math::exact::vrefines [math::exact::mdotv $m $v_]]
    }

    # Apply a matrix to the vector. 
    # Precondition: v is in lowest terms

    method applyM {m} {
	set d [math::exact::determinant $m]
	if {$d < 0} { set d [expr {-$d}] }
	if {($d & ($d-1)) != 0} {
	    return [math::exact::V new \
			[math::exact::vreduce [math::exact::mdotv $m $v_]]]
	} else {
	    return [math::exact::V new \
			[math::exact::vscale [math::exact::mdotv $m $v_]]]
	}
    }

    # Left-multiply a tensor t by the vector, and apply the result to
    # an expression 'r'
    method applyTLeft {t r} {
	set u [math::exact::mscale [math::exact::tleftv $t $v_]]
	set det [math::exact::determinant $u]
	if {$det < 0} { set det [expr {-$det}] }
	if {($det & ($det-1)) == 0} {
	    # determinant is a power of 2
	    set res [$r applyM $u]
	    return $res
	} else {
	    return [math::exact::Mstrict new $u $r]
	}
    }

    # Right-multiply a tensor t by the vector, and apply the result
    # to an expression 'l'
    method applyTRight {t l} {
	set u [math::exact::mscale [math::exact::trightv $t $v_]]
	set det [math::exact::determinant $u]
	if {$det < 0} { set det [expr {-$det}] }
	if {($det & ($det-1)) == 0} {
	    # determinant is a power of 2
	    set res [$l applyM $u]
	    return $res
	} else {
	    return [math::exact::Mstrict new $u $l]
	}
    }

    # Get the vector components
    method getV {} {
	return $v_
    }

    # Get the (zero-width) interval that the vector represents.
    method getInterval {nDigits} {
	return [list $v_ $v_]
    }

    # Absorb more information
    method absorb {} {
	# Nothing should ever call this, because a vector's information is
	# already complete.
	error "cannot absorb anything into a vector"
    }

    # asReal --
    #	Coerces an object from rational to real
    #
    # Parameters:
    #	None.
    #
    # Results:
    #	Returns this object converted to an exact real.
    method asReal {} {
	my ref
	lassign [my getSignAndMagnitude] s m
	set result [math::exact::Mstrict new $s $m]
	my unref
	return $result
    }

    # U- --
    #
    #   Unary - applied to this object
    #
    # Results:
    #	Returns the negation.

    method U- {} {
	my ref
	lassign $v_ p q
	set result [math::exact::V new [list [expr {-$p}] $q]]
	my unref
	return $result
    }; export U-

    # + --
    #	Adds a vector to another object
    #
    # Parameters:
    #	r - Right addend
    #
    # Results:
    #	Returns the sum
    #
    # The right-hand addend may be rational (an instance of V) or real
    # (any other Expression).
    #
    # This method is a Consumer with respect to the current object and to r.
    # It is a Constructor with respect to its result, returning a zero-ref
    # object.
    
    method + {r} {
	return [$r V+ [self object]]
    }; export +

    # E+ --
    #	Adds an exact real and a vector
    #
    # Parameters:
    #	l - Left addend
    #
    # Results:
    #	Returns the sim.
    #
    # This method is a Consumer with respect to the current object and to l.
    # It is a Constructor with respect to its result, returning a zero-ref
    # object.
    
    method E+ {l} {
	tailcall math::exact::+real $l [self object]
    }; export E+

    # V+ --
    #	Adds two rationals
    #
    # Parameters:
    #	l - Rational multiplicand
    #
    # Results:
    #	Returns the product.
    #
    # This method is a Consumer with respect to the current object and to l.
    # It is a Constructor with respect to its result, returning a zero-ref
    # object.
    method V+ {l} {
	my ref
	$l ref
	lassign [$l getV] a b
	lassign $v_ c d
	$l unref
	my unref
	return [math::exact::V new \
		    [math::exact::vreduce \
			 [list [expr {$a*$d+$b*$c}] [expr {$b*$d}]]]]
    }; export V+

    # - --
    #	Subtracts another object from a vector
    #
    # Parameters:
    #	r - Subtrahend
    #
    # Results:
    #	Returns the difference
    #
    # The right-hand operand may be rational (an instance of V) or real
    # (any other Expression).
    #
    # This method is a Consumer with respect to the current object and to r.
    # It is a Constructor with respect to its result, returning a zero-ref
    # object.
    
    method - {r} {
	return [$r V- [self object]]
    }; export -

    # E- --
    #	Subtracts this exact real from a rational
    #
    # Parameters:
    #	l - Left addend
    #
    # Results:
    #	Returns the difference.
    #
    # This method is a Consumer with respect to the current object and to l.
    # It is a Constructor with respect to its result, returning a zero-ref
    # object.
    
    method E- {l} {
	tailcall math::exact::-real $l [self object]
    }; export E-

    # V- --
    #	Subtracts this rational from another
    #
    # Parameters:
    #	l - Rational minuend
    #
    # Results:
    #	Returns the difference.
    #
    # This method is a Consumer with respect to the current object and to l.
    # It is a Constructor with respect to its result, returning a zero-ref
    # object.
    method V- {l} {
	my ref
	$l ref
	lassign [$l getV] a b
	lassign $v_ c d
	$l unref
	my unref
	return [math::exact::V new \
		    [math::exact::vreduce \
			 [list [expr {$a*$d-$b*$c}] [expr {$b*$d}]]]]
    }; export V-

    # * --
    #	Multiplies a rational by another object
    #
    # Parameters:
    #	r - Right-hand factor
    #
    # Results:
    #	Returns the difference
    #
    # The right-hand operand may be rational (an instance of V) or real
    # (any other Expression).
    #
    # This method is a Consumer with respect to the current object and to r.
    # It is a Constructor with respect to its result, returning a zero-ref
    # object.
    
    method * {r} {
	return [$r V* [self object]]
    }; export *

    # E* --
    #	Multiplies an exact real and a rational
    #
    # Parameters:
    #	l - Multiplicand
    #
    # Results:
    #	Returns the product.
    #
    # This method is a Consumer with respect to the current object and to l.
    # It is a Constructor with respect to its result, returning a zero-ref
    # object.
    
    method E* {l} {
	tailcall math::exact::*real $l [self object]
    }; export E*

    # V* --
    #	Multiplies two rationals
    #
    # Parameters:
    #	l - Rational multiplicand
    #
    # Results:
    #	Returns the product.
    #
    # This method is a Consumer with respect to the current object and to l.
    # It is a Constructor with respect to its result, returning a zero-ref
    # object.
    method V* {l} {
	my ref
	$l ref
	lassign [$l getV] a b
	lassign $v_ c d
	$l unref
	my unref
	return [math::exact::V new \
		    [math::exact::vreduce \
			 [list [expr {$a*$c}] [expr {$b*$d}]]]]
    }; export V*

    # / --
    #	Divides this object by another.
    #
    # Parameters:
    #	r - Divisor
    #
    # Results:
    #	Returns the quotient
    #
    # Either object may be rational (an instance of V) or real (any other
    # Expression).
    #
    # This method is a Consumer with respect to the current object and to r.
    # It is a Constructor with respect to its result, returning a zero-ref
    # object.
    
    method / {r} {
	return [$r V/ [self object]]
    }; export /
 
   # E/ --
    #	Divides an exact real and a rational
    #
    # Parameters:
    #	l - Dividend
    #
    # Results:
    #	Returns the quotient.
    #
    # The divisor is not a rationa.
    #
    # This method is a Consumer with respect to the current object and to l.
    # It is a Constructor with respect to its result, returning a zero-ref
    # object.
    
    method E/ {l} {
	tailcall math::exact::/real $l [self object]
    }; export E/

    # V/ --
    #	Divides two rationals
    #
    # Parameters:
    #	l - Dividend
    #
    # Results:
    #	Returns the quotient.
    #
    # This method is a Consumer with respect to the current object and to l.
    # It is a Constructor with respect to its result, returning a zero-ref
    # object.
    method V/ {l} {
	my ref
	$l ref
	lassign [$l getV] a b
	lassign $v_ c d
	set result [math::exact::V new \
			[math::exact::vreduce \
			     [list [expr {$a*$d}] [expr {$b*$c}]]]]
	$l unref
	my unref
	return $result
    }; export V/

    # ** -
    #	Raises a rational to a power
    #
    # Parameters:
    #	r - Exponent
    #
    # Results:
    #	Returns the power.
    #
    # This method is a Consumer with respect to the current object and to l.
    # It is a Constructor with respect to its result, returning a zero-ref
    # object.

    method ** {r} {
	tailcall $r V** [self object]
    }; export **

    # E** -
    #	Raises an exact real to a rational power
    #
    # Parameters:
    #	l - Base to exponentiate
    #
    # Results:
    #	Returns the power
    #
    # This method is a Consumer with respect to the current object and to l.
    # It is a Constructor with respect to its result, returning a zero-ref
    # object.

    method E** {l} {

	# Extract numerator and demominator of the exponent, and consume the
	# exponent.
	my ref
	lassign $v_ c d
	my unref
 
	# Normalize the sign of the exponent
	if {$d < 0} {
	    set c [expr {-$c}]
	    set d [expr {-$d}]
	}

	# Don't choke if somehow a 0/0 gets here.
	if {$c == 0 && $d == 0} {
	    $l unref
	    return -code error -errorcode "MATH EXACT ZERODIVZERO" \
		"zero divided by zero"
	}

	# Handle integer powers
	if {$d == 1} {
	    return [math::exact::real**int $l $c]
	}

	# Other rational powers come here.
	# We know that $d > 0, and we're not just doing
	# exponentiation by an integer
	    
	return [math::exact::real**rat $l $c $d]
    }; export E**

    # V** -
    #	Raises a rational base to a rational power
    #
    # Parameters:
    #	l - Base to exponentiate
    #
    # Results:
    #	Returns the power
    #
    # This method is a Consumer with respect to the current object and to l.
    # It is a Constructor with respect to its result, returning a zero-ref
    # object.

    method V** {l} {

	# Extract the numerator and denominator of the base and consume
	# the base.
	$l ref
	lassign [$l getV] a b
	$l unref

	# Extract numerator and demominator of the exponent, and consume the
	# exponent.
	my ref
	lassign $v_ c d
	my unref

	# Normalize the signs of the arguments
	if {$b < 0} {
	    set a [expr {-$a}]
	    set b [expr {-$b}]
	}
	if {$d < 0} {
	    set c [expr {-$c}]
	    set d [expr {-$d}]
	}

	# Don't choke if somehow a 0/0 gets here.
	if {$a == 0 && $b == 0 || $c == 0 && $d == 0} {
	    return -code error -errorcode "MATH EXACT ZERODIVZERO" \
		"zero divided by zero"
	}

	# b >= 0 and d >= 0

	if {$a == 0} {
	    if {$c == 0} {
		return -code error -errorcode "MATH EXACT ZEROPOWZERO" \
		    "zero to zero power"
	    } elseif {$d == 0} {
		return -code error -errorcode "MATH EXACT ZEROPOWINF" \
		    "zero to infinite power"
	    } else {
		return [math::exact::V new {0 1}]
	    }
	}

	# a != 0, b >= 0, d >= 0

	if {$b == 0} {
	    if {$c == 0} {
		return -code error -errorcode "MATH EXACT INFPOWZERO" \
		    "infinity to zero power"
	    } elseif {$c < 0} {
		return [math::exact::V new {0 1}]
	    } else {
		return [math::exact::V new {1 0}]
	    }
	}

	# a != 0, b > 0, d >= 0

	if {$c == 0} {
	    return [math::exact::V new {1 1}]
	}

	# handle integer exponents

	if {$d == 1} {
	    return [math::exact::rat**int $a $b $c]
	}

	# a != 0, b > 0, c != 0, d >= 0

	return [math::exact::rat**rat $a $b $c $d]
    }; export V**
    
    # sqrt --
    #
    #	Calculates the square root of this object
    #
    # Results:
    #	Returns the square root as an exact real.
    #
    # This method is a Consumer with respect to this object and a Constructor
    # with respect to the result, returning a zero-ref object.
    method sqrt {} {
	my ref
	if {([lindex $v_ 0] < 0) ^ ([lindex $v_ 1] < 0)} {
	    return -code error -errorCode "MATH EXACT SQRTNEGATIVE" \
		{square root of negative argument}
	}
	set result [::math::exact::Sqrtrat new {*}$v_]
	my unref
	return $result
    }
    
}

# math::exact::M --
#	Expression consisting of a matrix times another expression
#
# The matrix {a c} {b d} represents the homography (a*x + b) / (c*x + d).
#
# The inner expression may need to be evaluated lazily. Whether evaluation
# is strict or lazy, the 'e' method will return the expression.

oo::class create math::exact::M {
    superclass math::exact::Expression

    # m_ is the matrix; e_ the inner expression; absorbed_ a cache of the
    # result of the 'absorb' method.
    variable m_ e_ absorbed_

    # constructor accepts the matrix only. The expression is managed in
    # derived classes.
    constructor {m} {
	next
	set m_ $m
    }

    # destructor deletes the memoized expression if one has been stored.
    # The base class destructor handles cleaning up the result of 'absorb'
    destructor {
	if {[info exists e_]} {
	    $e_ unref
	}
	next
    }

    # Test if the matrix refines when premultiplied by another matrix n
    method refinesM {n} {
	return [math::exact::mrefines [math::exact::mdotm $n $m_]]
    }

    # Premultiply the matrix by another matrix n
    method applyM {n} {
	set d [math::exact::determinant $n]
	if {$d < 0} {set d [expr {-$d}]}
	if {($d & ($d-1)) != 0} {
	    return [math::exact::Mstrict new \
			[math::exact::mreduce [math::exact::mdotm $n $m_]] \
			[my e]]
	} else {
	    return [math::exact::Mstrict new \
			[math::exact::mscale [math::exact::mdotm $n $m_]] \
			[my e]]
	}
    }

    # Compute the left product of a tensor t with this matrix, and
    # apply the resulting tensor to the expression 'r'.
    method applyTLeft {t r} {
	return [math::exact::Tstrict new \
		    [math::exact::tscale [math::exact::tleftm $t $m_]] \
		    1 [my e] $r]
    }

    # Compute the right product of a tensor t with this matrix, and
    # apply the resulting tensor to the expression 'l'.
    method applyTRight {t l} {
	return [math::exact::Tstrict new \
		    [math::exact::tscale [math::exact::trightm $t $m_]] \
		    0 $l [my e]]
    }

    # Absorb a digit into this matrix.
    method absorb {} {
	if {![info exists absorbed_]} {
	    set absorbed_ [[[my e] applyM $m_] ref]
	}
	return $absorbed_
    }

    # Derived classes are expected to implement:
    # method e {} {
    #	# Returns the expression to which this matrix is applied.
    #	# Optionally memoizes the result in $e_.
    # }
    method e {} {
	error "[info object class [self object]] does not implement the 'e' method."
    }
}

# math::exact::Mstrict --
#
#	Expression representing the product of a matrix and another
#	expression.
#
# In this version of the class, the expression is known in advance - 
# evaluated strictly.

oo::class create math::exact::Mstrict {
    superclass math::exact::M

    # m_ is the matrix.
    # e_ is the expression
    # absorbed_ caches the result of the 'absorb' method.
    variable m_ e_ absorbed_

    # Constructor accepts the matrix and the expression to which
    # it applies.
    constructor {m e} {
	next $m 
	set e_ [$e ref]
    }

    # All the heavy lifting of destruction is performed in the base class.
    destructor {
	next
    }

    # The 'e' method returns the expression.
    method e {} {
	return $e_
    }

    # The 'dump' method formats this object for debugging.
    method dump {} {
	return "M($m_,[$e_ dump])"
    }
}

# math::exact::T --
#	Expression representing a 2x2x2 tensor of the third order,
#	applied to two subexpressions.

oo::class create math::exact::T {
    superclass math::exact::Expression

    # t_ - the tensor
    # i_ A flag indicating whether the next 'absorb' should come from the
    #    left (0) or the right (1).
    # l_ - the left subexpression
    # r_ - the right subexpression
    # absorbed_ - the result of an 'absorb' operation

    variable t_ i_ l_ r_ absorbed_

    # constructor accepts the tensor and the initial state for absorption
    constructor {t i} {
	next
	set t_ $t
	set i_ $i
    }

    # destructor removes cached items.
    destructor {
	if {[info exists l_]} {
	    $l_ unref
	}
	if {[info exists r_]} {
	    $r_ unref
	}
	next;			# The base class will clean up absorbed_
    }

    # refinesM --
    #
    #	Tests if this tensor refines when premultiplied by a matrix
    #
    # Parameters:
    #	m - matrix to test
    #
    # Results:
    #	Returns a Boolean indicator that is true if the product refines.

    method refinesM {m} {
	return [math::exact::trefines [math::exact::mdott $m $t_]]
    }

    # applyM --
    #
    #	Left multiplies this tensor by a matrix
    #
    # Parameters:
    #	m - Matrix to multiply
    #
    # Results:
    #	Returns the product
    #
    # This operation has the side effect of making the product strict at
    # the uppermost level, by calling [my l] [my r] to instantiate the
    # subexpressions.

    method applyM {m} {
	set d [math::exact::determinant $m]
	if {$d < 0} {set d [expr {-$d}]}
	if {($d & ($d-1)) != 0} {
	    return [math::exact::Tstrict new \
			[math::exact::treduce [math::exact::mdott $m $t_]] \
			0 [my l] [my r]]
	} else {
	    return [math::exact::Tstrict new \
			[math::exact::tscale [math::exact::mdott $m $t_]] \
			0 [my l] [my r]]
	}
    }

    # absorb --
    #
    #	Absorbs information from the subexpressions.
    #
    # Results:
    #	Returns a copy of the current object, with information from 
    #   at least one subexpression absorbed so that more information is
    #	immediately available.

    method absorb {} {
	if {![info exists absorbed_]} {
	    if {[math::exact::trefines $t_]} {
		lassign [math::exact::trans $t_] m n
		set side [math::exact::mdisjointm $m $n]
	    } else {
		set side $i_
	    }
	    if {$side} {
		set absorbed_ [[[my r] applyTRight $t_ [my l]] ref]
	    } else {
		set absorbed_ [[[my l] applyTLeft $t_ [my r]] ref]
	    }
	}
	return $absorbed_
    }

    # applyTRight --
    #
    #	Right-multiplies a tensor by this expression
    #
    # Results:
    #	Returns 't' left-product l right-product $r_.

    method applyTRight {t l} {
	# This is the hard case of digit exchange. We have to
	# get the leading digit from this tensor, absorbing as
	# necessary, right-multiply it into the tensor $t, and
	# compose the new object.
	#
	# Note that unless 'rest' is empty, 'ld' is a digit matrix,
	# so we need to check only for powers of 2 when reducing to
	# lowest terms
	lassign [my getLeadingDigitAndRest 0] ld rest
	if {$rest eq {}} {
	    set u [math::exact::mreduce [math::exact::trightv $t $ld]]
	    return [math::exact::Mstrict new $u $l]
	} else {
	    set u [math::exact::tscale [math::exact::trightm $t $ld]]
	    return [math::exact::Tstrict new $u 0 $l $rest]
	}
    }

    # applyTLeft --
    #
    #	Left-multiplies a tensor by this expression
    #
    # Results:
    #	Returns 't' left-product $l_ right-product 'r'
    method applyTLeft {t r} {
	# This is the hard case of digit exchange. We have to
	# get the leading digit from this tensor, absorbing as
	# necessary, left-multiply it into the tensor $t, and
	# compose the new object
	#
	# Note that unless 'rest' is empty, 'ld' is a digit matrix,
	# so we need to check only for powers of 2 when reducing to
	# lowest terms
	lassign [my getLeadingDigitAndRest 0] ld rest
	if {$rest eq {}} {
	    set u [math::exact::mreduce [math::exact::tleftv $t $ld]]
	    return [math::exact::Mstrict $u $r]
	} else {
	    set u [math::exact::tscale [math::exact::tleftm $t $ld]]
	    return [math::exact::Tstrict new $u 1 $rest $r]
	} 
    }

    # Derived classes are expected to implement the following:
    # l --
    #
    #	Returns the left operand
    method l {} {
	error "[info object class [self object]] does not implement the 'l' method"
    }

    # r --
    #
    #	Returns the right operand
    method r {} {
	error "[info object class [self object]] does not implement the 'r' method"
    }
    
}

# math::exact::Tstrict --
#
#	A strict tensor - one where the subexpressions are both known in
#	advance.

oo::class create math::exact::Tstrict {
    superclass math::exact::T

    # t_ - the tensor
    # i_ A flag indicating whether the next 'absorb' should come from the
    #    left (0) or the right (1).
    # l_ - the left subexpression
    # r_ - the right subexpression
    # absorbed_ - the result of an 'absorb' operation

    variable t_ i_ l_ r_ absorbed_

    # constructor accepts the tensor, the absorption state, and the
    # left and right operands.
    constructor {t i l r} {
	next $t $i
	set l_ [$l ref]
	set r_ [$r ref]
    }

    # base class handles all cleanup
    destructor {
	next
    }

    # l --
    #
    #	Returns the left operand
    method l {} {
	return $l_
    }

    # r --
    #
    #	Returns the right operand
    method r {} {
	return $r_
    }

    # dump --
    #
    #	Formats this object for debugging
    method dump {} {
	return T($t_,$i_\;[$l_ dump],[$r_ dump])
    }
}

# math::exact::opreal --
#
#	Applies a bihomography (bilinear fractional transformation)
#	to two expressions.
#
# Parameters:
#	op - Tensor {{{a b} {c d}} {{e f} {g h}}} representing the operation
#	x - left operand
#	y - right operand
#
# Results:
#	Returns an expression that represents the form:
#	(axy + cx + ey + g) / (bxy + dx + fy + h)
#
# Notes:
#	Note that the four basic arithmetic operations are included here.
#	In addition, this procedure may be used to craft other useful
#	transformations. For example, (1 - u**2) / (1 + u**2)
#	could be constructed as [opreal {{{-1 1} {0 0}} {{0 0} {1 1}}} $u $u]

proc math::exact::opreal {op x y {kludge {}}} {
    # split x and y into sign and magnitude
    $x ref; $y ref
    lassign [$x getSignAndMagnitude] sx mx
    lassign [$y getSignAndMagnitude] sy my
    $mx ref; $my ref
    $x unref; $y unref
    set t [tleftm [trightm $op $sy] $sx]
    set r [math::exact::Tstrict new $t 0 $mx $my]
    $mx unref; $my unref
    return $r
}

# math::exact::+real --
# math::exact::-real --
# math::exact::*real --
# math::exact::/real --
#
#	Sum, difference, product and quotient of exact reals
#
# Parameters:
#	x - First operand
#	y - Second operand
#
# Results:
#	Returns x+y, x-y, x*y or x/y as requested.

proc math::exact::+real {a b} { variable tadd; return [opreal $tadd $a $b] }
proc math::exact::-real {a b} { variable tsub; return [opreal $tsub $a $b] }
proc math::exact::*real {a b} { variable tmul; return [opreal $tmul $a $b] }
proc math::exact::/real {a b} { variable tdiv; return [opreal $tdiv $a $b] }

# real --
#
#	Coerce an argument to exact-real (possibly from rational)
#
# Parameters:
#	x - Argument to coerce.
#
# Results:
#	Returns the argument coerced to a real.
#
# This operation either does nothing and returns its argument, or is a
# Consumer with respect to its argument and a Constructor with respect to
# its result.

proc math::exact::function::real {x} {
    tailcall $x asReal
}

# SqrtWorker --
#
#	Class to calculate the square root of a real.


oo::class create math::exact::SqrtWorker {
    superclass math::exact::T
    variable l_ r_

    # e - The expression whose square root should be calculated.
    #     e should be between close to 1 for good performance. The
    #     'sqrtreal' procedure below handles the scaling.
    constructor {e} {
	next {{{1 0} {2 1}} {{1 2} {0 1}}} 0
	set l_ [$e ref]
    }
    method l {} {
	return $l_
    }
    method r {} {
	if {![info exists r_]} {
	    set r_ [[math::exact::SqrtWorker new $l_] ref]
	}
	return $r_
    }
    method dump {} {
	return "sqrt([$l_ dump])"
    }
}

# sqrt --
#
#	Returns the square root of a number
#
# Parameters:
#	x - Exact real number whose square root is needed.
#
# Results:
#	Returns the square root as an exact real.
#
# The number may be rational or real. There is a special optimization used
# if the number is rational

proc math::exact::function::sqrt {x} {
    tailcall $x sqrt
}

# ExpWorker --
#
#	Class that evaluates the exponential function for small exact reals

oo::class create math::exact::ExpWorker {
    superclass math::exact::T
    variable t_ l_ r_ n_

    # Constructor --
    #
    # Parameters:
    #	e - Argument whose exponential is to be computed. (What is
    #	    actually passed in is S0'(x) = (1+x)/(1-x))
    #	n - Number of the convergent of the continued fraction
    #
    # This class is implemented by expanding the continued fraction
    # as needed for precision. Each successive step becomes a new right
    # subexpression of the tensor product.

    constructor {e {n 0}} {
	next [list \
		  [list \
		       [list [expr {2*$n + 2}] [expr {2*$n + 1}]] \
		       [list [expr {2*$n + 1}] [expr {2*$n}]]] \
		  [list \
		       [list [expr {2*$n}] [expr {2*$n + 1}]] \
		       [list [expr {2*$n + 1}] [expr {2*$n + 2}]]]] 0
	set l_ [$e ref]
	set n_ [expr {$n + 1}]
    }

    # l --
    #
    #	Returns the left subexpression; that is, the argument to the
    #	exponential
    method l {} {
	return $l_
    }

    # r --
    #	Returns the right subexpresison - the next convergent, creating it
    #	if necessary
    method r {} {
	if {![info exists r_]} {
	    set r_ [[math::exact::ExpWorker new $l_ $n_] ref]
	}
	return $r_
    }

    # dump --
    #
    #	Displays this object for debugging
    method dump {} {
	return ExpWorker([$l_ dump],[expr {$n_-1}])
    }
}

# exp --
#
#	Evaluates the exponential function of an exact real
#
# Parameters:
#	x - Quantity to be exponentiated
#
# Results:
#	Returns the exact real function value.
#
# This procedure is a Consumer with respect to its argument and a
# Constructor with respect to its result, returning a zero-ref object.

proc math::exact::function::exp {x} {
    variable ::math::exact::iszer
    variable ::math::exact::tmul

    # The continued fraction converges only for arguments between -1 and 1.
    # If $iszer refines the argument, then it is in the correct range and
    # we launch ExpWorker to evaluate the continued fraction. If the argument
    # is outside the range [-1/2..1/2], then we evaluate exp(x/2) and square
    # the result. If neither of the above is true, then we perform a digit
    # exchange to get more information about the magnitude of the argument.

    $x ref
    if {[$x refinesM $iszer]} {
	# Argument's absolute value is small - evaluate the exponential
	set y [$x applyM $iszer]
	set result [ExpWorker new $y]
    } elseif {[$x refinesM {{2 2} {-1 1}}]} {
	# Argument's absolute value is large - evaluate exp(x/2)**2
	set xover2 [$x applyM {{1 0} {0 2}}]
	set expxover2 [exp $xover2]
	set result [*real $expxover2 $expxover2]
    } else {
	# Argument's absolute value is uncharacterized - perform a digit
	# exchange to get more information.
	set result [exp [$x absorb]]
    }
    $x unref
    return $result
}

# LogWorker --
#
#	Helper class for evaluating logarithm of an exact real argument.
#
# The algorithm used is a continued fraction representation from Peter Potts's
# paper. This worker evaluates the second and subsequent convergents. The
# first convergent is in the 'log' procedure below, and follows a different
# pattern from the rest of them.

oo::class create math::exact::LogWorker {
    superclass math::exact::T
    variable t_ l_ r_ n_

    # Constructor -
    #
    # Parameters:
    #	e - Argument whose log is to be extracted
    #   n - Number of the convergent.
    constructor {e {n 1}} {
	next [list \
		  [list \
		       [list $n 0] \
		       [list [expr {2*$n + 1}] [expr {$n+1}]]] \
		  [list \
		       [list [expr {$n + 1}] [expr {2*$n + 1}]] \
		       [list 0 $n]]] 0
	set l_ [$e ref]
	set n_ [expr {$n + 1}]
    }

    # l -
    #	Returns the argument whose log is to be extracted
    method l {} {
	return $l_
    }

    # r -
    #	Returns the next convergent, constructing it if necessary.
    method r {} {
	if {![info exists r_]} {
	    set r_ [[math::exact::LogWorker new $l_ $n_] ref]
	}
	return $r_
    }

    # dump -
    #	Dumps this object for debugging
    method dump {} {
	return LogWorker([$l_ dump],[expr {$n_-1}])
    }
}

# log -
#
#	Calculates the natural logarithm of an exact real argument.
#
# Parameters:
#	x - Quantity whose log is to be extracted.
#
# Results:
#	Returns the logarithm
#
# This procedure is a Consumer with respect to its argument and a Constructor
# with respect to its result, returning a zero-ref object.

proc math::exact::function::log {x} {
    variable ::math::exact::ispos
    variable ::math::exact::isneg
    variable ::math::exact::idpos
    variable ::math::exact::idneg
    variable ::math::exact::log2

    # If x is between 1/2 and 2, the continued fraction will converge. If
    # y = LogWorker(x), then log(x) = (xy + x - y - 1)/(x + y), and the
    # latter function is a bihomography that can be evaluated by 'opreal'
    # directly.
    #
    # If x is negative, that's an error.
    # If x > 1, idpos will refine it, and we compute log(x/2) + log(2)
    # If x < 1, idneg will refine it, and we compute log(2x) - log(2)
    # If none of the above can be proven, perform a digit exchange and
    # try again.

    $x ref
    if {[$x refinesM {{2 -1} {-1 2}}]} {
	# argument in bounds
	set result [math::exact::opreal {{{1 0} {1 1}} {{-1 1} {-1 0}}} \
			$x \
			[LogWorker new $x]]
    } elseif {[$x refinesM $isneg]} {
	# domain error
        return -code error -errorcode {MATH EXACT LOGNEGATIVE} \
	    "log of negative argument"
    } elseif {[$x refinesM $idpos]} {
	# large argument, reduce it and try again
	set result [+real [function::log [$x applyM {{1 0} {0 2}}]] $log2]
    } elseif {[$x refinesM $idneg]} {
	# small argument, increase it and try again
	set result [-real [function::log [$x applyM {{2 0} {0 1}}]] $log2]
    } else {
	# too little information, perform digit exchange.
	set result [function::log [$x absorb]]
    }
    $x unref
    return $result
}

# TanWorker --
#
#	Auxiliary function for tangent of an exact real argument
#
# This class develops the second and subsequent convergents of the continued
# fraction expansion in Potts's paper
oo::class create math::exact::TanWorker {
    superclass math::exact::T
    variable t_ l_ r_ n_

    # Constructor -
    #
    # Parameters:
    #	e - S0'(x) = (1+x)/(1-x), where we wish to evaluate tan(x).
    #   n - Ordinal position of the convergent
    constructor {e {n 1}} {
	next [list \
		  [list \
		       [list [expr {2*$n + 1}] [expr {2*$n + 3}]] \
		       [list [expr {2*$n - 1}] [expr {2*$n + 1}]]] \
		  [list \
		       [list [expr {2*$n + 1}] [expr {2*$n - 1}]] \
		       [list [expr {2*$n + 3}] [expr {2*$n + 1}]]]] 0
	set l_ [$e ref]
	set n_ [expr {$n + 1}]
    }

    # l -
    #  	Returns the argument S0'(x)
    method l {} {
	return $l_
    }

    # r -
    #	Returns the next convergent, constructing it if necessary
    method r {} {
	if {![info exists r_]} {
	    set r_ [[math::exact::TanWorker new $l_ $n_] ref]
	}
	return $r_
    }

    # dump -
    #	Displays this object for debugging
    method dump {} {
	return TanWorker([$l_ dump],[expr {$n_-1}])
    }
}

# tan --
#	Tangent of an exact real argument
#
# Parameters:
#	x - Quantity whose tangent is to be computed.
#
# Results:
#	Returns the tangent
#
# This procedure is a Consumer with respect to its argument and a Constructor
# with respect to its result, returning a zero-ref object.

proc math::exact::function::tan {x} {
    variable ::math::exact::iszer

    # If |x| < 1, then we use Potts's formula for the tangent.
    # If |x| > 1/2, then we compute y = tan(x/2) and then use the
    # trig identity tan(x) = 2*y/(1-y**2), recognizing that the latter
    # expression can be expressed as a bihomography applied to y and itself,
    # allowing opreal to do the job.
    # If neither can be proven, we perform a digit exchange to get more
    # information.
    # tan((2*n+1)*pi/2), for n an integer, is a well-behaved pole.
    # In particular, 1/tan(pi/2) will correctly return zero.

    $x ref
    if {[$x refinesM $iszer]} {
	set xx [$x applyM $iszer]
	set result [math::exact::Tstrict new {{{1 2} {1 0}} {{-1 0} {-1 2}}} 0 \
			$xx [TanWorker new $xx]]
    } elseif {[$x refinesM {{2 2} {-1 1}}]} {
	set xover2 [$x applyM {{1 0} {0 2}}]
	set tanxover2 [function::tan $xover2]
	set result [opreal {{{0 -1} {1 0}} {{1 0} {0 1}}} $tanxover2 $tanxover2]
    } else {
	set result [function::tan [$x absorb]]
    }
    $x unref
    return $result
}

# sin --
#	Sine of an exact real argument
#
# Parameters:
#	x - Quantity whose sine is to be computed.
#
# Results:
#	Returns the sine
#
# This procedure is a Consumer with respect to its argument and a Constructor
# with respect to its result, returning a zero-ref object.

proc math::exact::function::sin {x} {
    $x ref
    set tanxover2 [tan [$x applyM {{1 0} {0 2}}]]
    $x unref
    return [opreal {{{0 1} {1 0}} {{1 0} {0 1}}} $tanxover2 $tanxover2]
}

# cos --
#	Cosine of an exact real argument
#
# Parameters:
#	x - Quantity whose cosine is to be computed.
#
# Results:
#	Returns the cosine
#
# This procedure is a Consumer with respect to its argument and a Constructor
# with respect to its result, returning a zero-ref object.

proc math::exact::function::cos {x} {
    $x ref
    set tanxover2 [tan [$x applyM {{1 0} {0 2}}]]
    $x unref
    return [opreal {{{-1 1} {0 0}} {{0 0} {1 1}}} $tanxover2 $tanxover2]
}

# AtanWorker --
#
#	Auxiliary function for arctangent of an exact real argument
#
# This class develops the second and subsequent convergents of the continued
# fraction expansion in Potts's paper. The argument lies in [-1,1].

oo::class create math::exact::AtanWorker {
    superclass math::exact::T
    variable t_ l_ r_ n_
    # Constructor -
    #
    # Parameters:
    #	e - S0(x) = (x-1)/(x+1), where we wish to evaluate atan(x).
    #   n - Ordinal position of the convergent
    constructor {e {n 1}} {
	next [list \
		  [list \
		       [list [expr {2*$n + 1}] [expr {$n + 1}]] \
		       [list $n 0]] \
		  [list \
		       [list 0 $n] \
		       [list [expr {$n + 1}] [expr {2*$n + 1}]]]] 0
	set l_ [$e ref]
	set n_ [expr {$n + 1}]
    }

    # l -
    #  	Returns the argument S0(x)
    method l {} {
	return $l_
    }

    # r -
    #	Returns the next convergent, constructing it if necessary
    method r {} {
	if {![info exists r_]} {
	    set r_ [[math::exact::AtanWorker new $l_ $n_] ref]
	}
	return $r_
    }

    # dump -
    #	Displays this object for debugging
    method dump {} {
	return AtanWorker([$l_ dump],[expr {$n_-1}])
    }
}

# atanS0 -
#
#	Evaluates the arctangent of S0(x) = (x-1)/(x+1)
#
# Parameters:
#	x - Exact real argumetn
#
# Results:
#	Returns atan((x-1)/(x+1))
#
# This function is a Consumer with respect to its argument and a Constructor
# with respect to its result, returning a 0-reference object.

proc math::exact::atanS0 {x} {
    return [opreal {{{1 2} {1 0}} {{-1 0} {-1 2}}} $x [AtanWorker new $x]]
}

# atan -
#
#	Arctangent of an exact real
#
# Parameters:
#	x - Exact real argument
#
# Results:
#	Returns atan(x)
#
# This function is a Consumer with respect to its argument and a Constructor
# with respect to its result, returning a 0-reference object.
#
# atan(1/0) is undefined and may cause an infinite loop.

proc math::exact::function::atan {x} {

    # TODO - find p/q close to the real number x - can be done by
    #        getting a few digits - and do
    # arctan(p/q + eps) = arctan(p/q) + arctan(q**2*eps/(p*q*eps+p**q+q**2))
    # using [$eps applyM] to compute the argument of the second arctan

    variable ::math::exact::szer 
    variable ::math::exact::spos 
    variable ::math::exact::sinf 
    variable ::math::exact::sneg
    variable ::math::exact::pi

    # Four cases, depending on which octant the arctangent lies in.
    
    $x ref
    lassign [$x getSignAndMagnitude] signum mag
    $mag ref
    $x unref
    set aS0x [atanS0 $mag]
    $mag unref
    if {$signum eq $szer} {
	# -1 < x < 1
	return $aS0x
    } elseif {$signum eq $spos} {
	# x > 0
	return [opreal {{{0 0} {4 0}} {{1 0} {0 4}}} $aS0x $pi]
    } elseif {$signum eq $sinf} {
	# x < -1 or x > 1
	return [opreal {{{0 0} {2 0}} {{1 0} {0 2}}} $aS0x $pi]
    } elseif {$signum eq $sneg} {
	# x < 0
	return [opreal {{{0 0} {4 0}} {{-1 0} {0 4}}} $aS0x $pi]
    } else {
	# can't happen
	error "wrong sign: $signum"
    }
}

# asinreal -
#
#	Computes the arcsine of an exact real argument.
#
# The arcsine is computed from the arctangent by trigonometric identities
#
# This function is a Consumer with respect to its argument and a Constructor
# with respect to its result, returning a 0-reference object.
#
# The function is defined only over the open interval (-1,1). Outside
# that range INCLUDING AT THE ENDPOINTS, it may fail and give an infinite
# loop or stack overflow.

proc math::exact::asinreal {x} {
    variable iszer
    variable pi

    # Potts's formula doesn't work here - it's singular at zero,
    # and undefined over negative numbers. But some messing with the
    # algebra gives us:
    #     asin(S0*x) = 2*atan(sqrt(x)) - pi/2
    #                = (4*atan(sqrt(x)) - pi) / 2
    # which is continuous and computable over (-1..1)
    $x ref
    set y [$x applyM $iszer]
    $x unref
    return [opreal {{{0 0} {-1 0}} {{4 0} {0 2}}} \
		$pi \
		[function::atan [function::sqrt $y]]]
}

interp alias {} math::exact::function::asin {} math::exact::asinreal

# acosreal -
#
#	Computes the arccosine of an exact real argument.
#
# The arccosine is computed from the arctangent by trigonometric identities
#
# This function is a Consumer with respect to its argument and a Constructor
# with respect to its result, returning a 0-reference object.
#
# The function is defined only over the open interval (-1,1). Outside
# that range INCLUDING AT THE ENDPOINTS, it may fail and give an infinite
# loop or stack overflow.

proc math::exact::acosreal {x} {
    variable iszer
    variable pi
    # Potts's formula doesn't work here - it's singular at zero,
    # and undefined over negative numbers. But some messing with the
    # algebra gives us:
    # acos(S0*x) = pi - 2*atan(sqrt(x))
    $x ref
    set y [$x applyM $iszer]
    $x unref
    return [opreal {{{0 0} {1 0}} {{-2 0} {0 1}}} \
		$pi \
		[function::atan [function::sqrt $y]]]
}

interp alias {} math::exact::function::acos {} math::exact::acosreal

# sinhreal, coshreal, tanhreal --
#
#	Hyperbolic functions of exact real arguments
#
# Parameter:
#	x - Argument at which to evaluate the function
#
# Results:
#	Return sinh(x), cosh(x), tanh(x), respectively.
#
# These functions are all Consumers with respect to their arguments and
# Constructors with respect to their results, returning zero-ref objects.
#
# The three functions are well defined over all the finite reals, but
# are ill-behaved at infinity.

proc math::exact::sinhreal {x} {
    set expx [function::exp $x]
    return [opreal {{{1 0} {0 1}} {{0 1} {-1 0}}} $expx $expx]
}

interp alias {} math::exact::function::sinh {} math::exact::sinhreal

proc math::exact::coshreal {x} {
    set expx [function::exp $x]
    return [opreal {{{1 0} {0 1}} {{0 1} {1 0}}} $expx $expx]
}

interp alias {} math::exact::function::cosh {} math::exact::coshreal

proc math::exact::tanhreal {x} {
    set expx [function::exp $x]
    return [opreal {{{1 1} {0 0}} {{0 0} {-1 1}}} $expx $expx]
}

interp alias {} math::exact::function::tanh {} math::exact::tanhreal

# asinhreal, acoshreal, atanhreal --
#
#	Inverse hyperbolic functions of exact real arguments
#
# Parameter:
#	x - Argument at which to evaluate the function
#
# Results:
#	Return asinh(x), acosh(x), atanh(x), respectively.
#
# These functions are all Consumers with respect to their arguments and
# Constructors with respect to their results, returning zero-ref objects.
#
# asinh is defined over the entire real number line, with the exception
# of the point at infinity.  acosh is defined over x > 1 (NOT x=1, which
# is singular). atanh is defined over (-1..1) (NOT the endpoints of the
# interval.)

proc math::exact::asinhreal {x} {
    # domain (-Inf .. Inf)
    # asinh(x) = log(x + sqrt(x**2 + 1))
    $x ref
    set retval [function::log \
		    [+real $x \
			 [function::sqrt \
			      [opreal {{{1 0} {0 0}} {{0 0} {1 1}}} $x $x]]]]
    $x unref
    return $retval
}

interp alias {} math::exact::function::asinh {} math::exact::asinhreal

proc math::exact::acoshreal {x} {
    # domain (1 .. Inf)
    # asinh(x) = log(x + sqrt(x**2 - 1))
    $x ref
    set retval [function::log \
		    [+real $x \
			 [function::sqrt \
			      [opreal {{{1 0} {0 0}} {{0 0} {-1 1}}} $x $x]]]]
    $x unref
    return $retval
}

interp alias {} math::exact::function::acosh {} math::exact::acoshreal

proc math::exact::atanhreal {x} {
    # domain (-1 .. 1)
    variable sinf
    #atanh(x) = log(Sinf[x])/2

    $x ref
    set y [$x applyM $sinf]
    $y ref
    $x unref
    set z [function::log $y]
    $z ref
    $y unref
    set retval [$z applyM {{1 0} {0 2}}]
    $z unref
    return $retval
}

interp alias {} math::exact::function::atanh {} math::exact::atanhreal

# EWorker --
#
#	Evaluates the constant 'e' (the base of the natural logarithms
#
# This class is intended to be singleton. It returns 2.71828.... (the
# base of the natural logarithms) as an exact real.

oo::class create math::exact::EWorker {
    superclass math::exact::M
    variable m_ e_ n_

    # Constructor accepts the number of the continuant.

    constructor {{n 0}} {
	set n_ [expr {$n + 1}]
	next [list [list [expr {2*$n + 2}] [expr {2*$n + 1}]] \
		  [list [expr {2*$n + 1}] [expr {2*$n}]]]
    }
    destructor {
	next
    }

    # e -- Returns the next continuant after this one.

    method e {} {
	if {![info exists e_]} {
	    set e_ [[math::exact::EWorker new $n_] ref]
	}
	return $e_
    }

    # Formats this object for debugging
    
    method dump {} {
	return M($m_,EWorker($n_))
    }
}

# PiWorker --
#
#	Auxiliary object used in evaluating pi.
#
# This class evaluates the second and subsequent continuants in
# Ramanaujan's formula for sqrt(10005)/pi. The Potts paper presents
# the algorithm, almost without commentary.

oo::class create math::exact::PiWorker {
    superclass math::exact::M
    variable m_ e_ n_

    # Constructor accepts the number of the continuant

    constructor {{n 1}} {
	set n_ [expr {$n + 1}]
	set nsq [expr {$n * $n}]
	set n4 [expr {$nsq * $nsq}]
	set b [expr {(2*$n - 1) * (6*$n - 5) * (6*$n - 1)}]
	set c [expr {$b * (545140134 * $n + 13591409)}]
	set d [expr {$b * ($n + 1)}]
	set e [expr {10939058860032000 * $n4}]
	set p [list [expr {$e - $d - $c}] [expr {$e + $d + $c}]]
	set q [list [expr {$e + $d - $c}] [expr {$e - $d + $c}]]
	next [list $p $q]
    }
    destructor {
	next
    }

    # e --
    #
    #	Returns the next continuant after this one

    method e {} {
	if {![info exists e_]} {
	    set e_ [[math::exact::PiWorker new $n_] ref]
	}
	return $e_
    }

    # dump --
    #
    #	Formats this object for debugging
    method dump {} {
	return M($m_,PiWorker($n_))
    }
}

# Log2Worker --
#
#	Auxiliary class for evaluating log(2).
#
# This object represents the constant (1-2*log(2))/(log(2)-1), the
# product of the second, third, ... nth LFT's of the representation of log(2).

oo::class create math::exact::Log2Worker {
    superclass math::exact::M
    variable m_ e_ n_

    # Constructor accepts the number of the continuant
    constructor {{n 1}} {
	set n_ [expr {$n + 1}]
	set a [expr {3*$n + 1}]
	set b [expr {2*$n + 1}]
	set c [expr {4*$n + 2}]
	set d [expr {3*$n + 2}]
	next [list [list $a $b] [list $c $d]]
    }
    destructor {
	next
    }

    # e --
    #
    #	Returns the next continuant after this one.
    method e {} {
	if {![info exists e_]} {
	    set e_ [[math::exact::Log2Worker new $n_] ref]
	}
	return $e_
    }

    # dump --
    #
    #	Displays this object for debugging
    method dump {} {
	return M($m_,Log2Worker($n_))
    }
}

# Sqrtrat --
#
#	Class that evaluates the square root of a rational

oo::class create math::exact::Sqrtrat {
    superclass math::exact::M
    variable m_ e_ a_ b_ c_

    # Constructor accepts the numerator and denominator. The third argument
    # is an intermediate result for the second and later continuants.
    constructor {a b {c {}}} {
	if {$c eq {}} {
	    set c [expr {$a - $b}]
	}
	set d [expr {2*($b-$a) + $c}]
	if {$d >= 0} {
	    next $math::exact::dneg
	    set a_ [expr {4 * $a}]
	    set b_ $d
	    set c_ $c
	} else {
	    next $math::exact::dpos
	    set a_ [expr {-$d}]
	    set b_ [expr {4 * $b}]
	    set c_ $c
	}
    }
    destructor {
	next
    }

    # e --
    #
    #	Returns the next continuant after this one.
    method e {} {
	if {![info exists e_]} {
	    set e_ [[math::exact::Sqrtrat new $a_ $b_ $c_] ref]
	}
	return $e_
    }

    # dump --
    #	Formats this object for debugging.
    
    method dump {} {
	return "M($m_,Sqrtrat($a_,$b_,$c_))"
    }
}

# math::exact::rat**int --
#
#	Service procedure to raise a rational number to an integer power
#
# Parameters:
#	a - Numerator of the rational
#	b - Denominator of the rational
#	n - Power
#
# Preconditions:
#	n is not zero, a is not zero, b is positive.
#
# Results:
#	Returns the power
#
# This procedure is a Consumer with respect to its arguments and a
# Constructor with respect to its result, returning a zero-ref object.

proc math::exact::rat**int {a b n} {
    if {$n < 0} {
	return [V new [list [expr {$b**(-$n)}] [expr {$a**(-$n)}]]]
    } elseif {$n > 0} {
	return [V new [list [expr {$a**($n)}] [expr {$b**($n)}]]]
    } else { ;# zero power shouldn't get here
	return [V new {1 1}]
    }
}

# math::exact::rat**rat --
#
#	Service procedure to raise a rational number to a rational power
#
# Parameters:
#	a - Numerator of the base
#	b - Denominator of the base
#	m - Numerator of the exponent
#	n - Denominator of the exponent
#
# Results:
#	Returns the power as an exact real
#
# Preconditions:
#	a != 0, b > 0, m != 0, n > 0
#
# This procedure is a Constructor with respect to its result

proc math::exact::rat**rat {a b m n} {

    # It would be attractive to special case this, but the real mechanism
    # works as well for the moment.

    tailcall real**rat [V new [list $a $b]] $m $n
}

# PowWorker --
#
#	Auxiliary class to compute
#		((p/q)**n + b)**(m/n),
#	where 0<m<n are integers, p, q are integers, b is an exact real

oo::class create math::exact::PowWorker {
    superclass math::exact::T

    variable t_ l_ r_ delta_

    # Self-method: start
    #
    #	Sets up to find z**(m/n) (1 <= m < n), with
    #   z = (p/q)**n + y for integers p and q.
    #
    # Parameters:
    #	p - numerator of the estimated nth root
    #	q - denominator of the estimated nth root
    #	y - residual of the quantity whose root is being extracted
    #	m - numerator of the exponent
    #	n - denominator of the exponent (1 <= m < n)
    #
    # Results:
    #	Returns the power, as an exact real.
    
    self method start {p q y m n} {
	set pm [expr {$p ** $m}]
	set pnmm [expr {$p ** ($n-$m)}]
	set pn [expr {$pm * $pnmm}]
	set qm [expr {$q ** $m}]
	set qnmm [expr {$q ** ($n-$m)}]
	set qn [expr {$qm * $qnmm}]

	set t0 \
	    [list \
		 [list \
		      [list [expr {$m * $qn}] [expr {$n*$pnmm*$qm}]] \
		      [list 0 [expr {($n-$m) * $qn}]]] \
		 [list \
		      [list [expr {2 * $n * $pn}] 0] \
		      [list [expr {2 * ($n-$m) * $pm * $qnmm}] 0]]]
	set t1 \
	    [list \
		 [list \
		      [list [expr {$n * $qn}] [expr {2*$n * $pnmm*$qm}]] \
		      [list 0 [expr {$n * $qn}]]] \
		 [list \
		      [list [expr {4 * $n * $pn}] 0] \
		      [list [expr {2 * $n * $pm * $qnmm}] 0]]]
    
	set tinit \
	    [list \
		 [list \
		      [list [expr {$m * $qn}] 0] \
		      [list 0 0]] \
		 [list \
		      [list [expr {$n * $pn}] [expr {$n * $pnmm * $qm}]] \
		      [list \
			   [expr {($n-$m) * $pm * $qnmm}] \
			   [expr {($n-$m) * $qn}]]]]
	$y ref
	set result [$y applyTLeft $tinit [my new $t0 $t1 $y]]
	$y unref
	return $result
    }

    # Constructor --
    #
    # Parameters:
    #	t0 - Tensor from the previous iteration
    #	delta - Increment to use
    #	y - Residual
    #
    # The constructor should not be called directly. Instead, the 'start'
    # method should be called to initialize the iteration

    constructor {t0 delta y} {
	set t [math::exact::tadd $t0 $delta]
	next $t 0
	set l_ [$y ref]
	set delta_ $delta
    }

    # l --
    #
    #	Returns the left subexpression: that is, the 'y' parameter
    method l {} {
	return $l_
    }

    # r --
    #
    #	Returns the right subexpression: that is, the next continuant,
    #	creating it if necessary
    method r {} {
	if {![info exists r_]} {
	    set r_ [[math::exact::PowWorker new $t_ $delta_ $l_] ref]
	}
	return $r_
    }

    method dump {} {
	set res "PowWorker($t_,$delta_,[$l_ dump],"
	if {[info exists r_]} {
	    append res [$r_ dump]
	} else {
	    append res ...
	}
	append res ")"
	return $res
    }
    
}

# math::exact::real**int --
#
#	Service procedure to raise a real number to an integer power.
#
# Parameters:
#	b - Number to exponentiate
#	e - Power to raise b to.
#
# Results:
#	Returns the power.
#
# This procedure is a Consumer with respect to its arguments and a
# Constructor with respect to its result, returning a zero-ref object.

proc math::exact::real**int {b e} {

    # Handle a negative power by raising the reciprocal of the base to
    # a positive power
    if {$e < 0} {
	set e [expr {-$e}]
	set b [K [[$b ref] applyM {{0 1} {1 0}}] [$b unref]]
    }
    
    # Reduce using square-and-add
    $b ref
    set result [V new {1 1}]
    while {$e != 0} {
	if {$e & 1} {
	    set result [$b * $result]
	    set e [expr {$e & ~1}]
	}
	if {$e == 0} break
	set b [K [[$b * $b] ref] [$b unref]]
	set e [expr {$e>>1}]
    }
    $b unref
    return $result
}

# math::exact::real**rat --
#
#	Service procedure to raise a real number to a rational power.
#
# Parameters -
#
#	b - The base to be exponentiated
#	m - The numerator of the power
#	n - The denominator of the power
#
# Preconditions:
#	n > 0
#
# Results:
#	Returns the power.
#
# This procedure is a Consumer with respect to its arguments and a
# Constructor with respect to its result, returning a zero-ref object.

proc math::exact::real**rat {b m n} {

    variable isneg
    variable ispos

    # At this point we need to know the sign of b. Try to determine it.
    # (This can be an infinite loop if b is zero or infinite)
    while {1} {
	if {[$b refinesM $ispos]} {
	    break
	} elseif {[$b refinesM $isneg]} {
	    # negative number to rational power. The denominator must be
	    # odd.
	    if {$n % 2 == 0} {
		return -code error -errorCode {MATH EXACT NEGATIVEPOWREAL} \
		    "negative number to real power"
	    } else {
		set b [K [[$b ref] U-] [$b unref]]
		tailcall [math::exact::real**rat $b $m $n] U-
	    }
	} else {
	    # can't determine positive or negative yet
	    $b ref
	    set nextb [$b absorb]
	    set result [math::exact::real**rat $nextb $m $n]
	    $b unref
	    return $result
	}
    }
	    
    # Handle b(-m/n) by taking (1/b)(m/n)
    if {$m < 0} {
	set m [expr {-$m}]
	set b [K [[$b ref] applyM {{0 1} {1 0}}] [$b unref]]
    }

    # Break m/n apart into integer and fractional parts
    set i [expr {$m / $n}]
    set m [expr {$m % $n}]

    # Do the integer part
    $b ref
    set result [real**int $b $i]
    if {$m == 0} {
	# We really shouldn't get here if m/n is an integer, but don't choke
	$b unref
	return $result
    }

    # Come up with a rational approximation for b**(1/n)
    # real: exp(log(b)/n)
    set approx [[math::exact::function::exp \
		     [[math::exact::function::log $b] \
			  * [math::exact::V new [list 1 $n]]]] ref]
    lassign [$approx getSignAndMagnitude] partial rest
    $rest ref
    $approx unref
    while {1} {
	lassign [$rest getLeadingDigitAndRest 0] digit y
	$y ref
	$rest unref
	set partial [math::exact::mscale [math::exact::mdotm $partial $digit]]
	set rest $y
	lassign $partial pq rs
	lassign $pq p q
	lassign $rs r s
	set qrn [expr {($q*$r)**$n}]
	set t1 [expr {$qrn}]
	set t2 [expr {2 * ($p*$s)**$n}]
	set t3 [expr {4 * $qrn}]
	if {$t1 < $t2 && $t2 < $t3} break
    }
    $y unref
    
    # Get the residual

    lassign [math::exact::vscale [list $r $s]] p q
    set xn [math::exact::V new [list [expr {$p**$n}] [expr {$q**$n}]]]
    set y [$b - $xn]; $b unref

    # Launch a worker process to perform quasi-Newton iteration to refine
    # the result
    
    set retval [$result * [math::exact::PowWorker start $p $q $y $m $n]]
    return $retval
}

# pi --
#
#	Returns pi as an exact real

proc math::exact::function::pi {} {
    variable ::math::exact::pi
    return $pi
}

# e --
#
#	Returns e as an exact real

proc math::exact::function::e {} {
    variable ::math::exact::e
    return $e
}

# math::exact::signum1 --
#
#	Tests an argument's sign.
#
# Parameters:
#	x - Exact real number to test.
#
# Results:
#	Returns -1 if x < -1. Returns 1 if x > 1. May return -1, 0 or 1 if
#	-1 <= x <= 1.
#
# Equality of exact reals is not decidable, so a weaker version of comparison
# testing is needed. This function provides the guts of such a thing. It
# returns an approximation to the signum function that is exact for
# |x| > 1, and arbitrary for |x| < 1.
#
# A typical use would be to replace a test p < q with a test that
# looks like signum1((p-q) / epsilon) == -1. This test is decidable,
# and becomes a test that is true if p < q - epsilon, false if p > q+epsilon,
# and indeterminate if p lies within epsilon of q.  This test is enough for
# most checks for convergence or for selecting a branch of a function.
#
# This function is not decidable if it is not decidable whether x is finite.

proc math::exact::signum1 {x} {
    variable ispos
    variable isneg
    variable iszer
    while {1} {
	if {[$x refinesM $ispos]} {
	    return 1
	} elseif {[$x refinesM $isneg]} {
	    return -1
	} elseif {[$x refinesM $iszer]} {
	    return 0
	} else {
	    set x [$x absorb]
	}
    }
}

# math::exact::abs1 -
#
#	Test whether an exact real is 'small' in absolute value.
#
# Parameters:
#	x - Exact real number to test
#
# Results:
#	Returns 0 if |x| is 'close to zero', 1 if |x| is 'far from zero'
#	and either 0, or 1 if |x| is close to 1.
#
# This function is another useful comparator for convergence testing.
# It returns a three-way indication:
#	|x| < 1/2 : 0
#	|x| > 1 : 1
#	1/2 <= |x| <= 2 : May return -1, 0, 1
#
# This function is useful for convergence testing, where it is desired
# to know whether a given value has an absolute value less than a given
# tolerance.

proc math::exact::abs1 {x} {
    variable iszer
    while 1 {
	if {[$x refinesM $iszer]} {
	    return 0
	} elseif {[$x refinesM {{2 1} {-2 1}}]} {
	    return 1
	} else {
	    set x [$x absorb]
	}
    }
}

namespace eval math::exact {

    # Constant vectors, matrices and tensors

    ;				# the identity matrix
    variable identity		{{ 1  0} { 0  1}}
    ;				# sign matrices for exact floating point
    variable spos		$identity
    variable sinf		{{ 1 -1} { 1  1}}
    variable sneg		{{ 0  1} {-1  0}}
    variable szer		{{ 1  1} {-1  1}}

    ;				# inverses of the sign matrices
    variable ispos		[reverse $spos]
    variable isinf		[reverse $sinf]
    variable isneg		[reverse $sneg]
    variable iszer		[reverse $szer]

    ;				# digit matrices for exact floating point
    variable dneg		{{ 1  1} { 0  2}}
    variable dzer		{{ 3  1} { 1  3}}
    variable dpos		{{ 2  0} { 1  1}}

    ;				# inverses of the digit matrices
    variable idneg 		[reverse $dneg]
    variable idzer		[reverse $dzer]
    variable idpos		[reverse $dpos]

    ;				# aritmetic operators as tensors
    variable tadd		{{{ 0  0} { 1  0}} {{ 1  0} { 0  1}}}
    variable tsub		{{{ 0  0} { 1  0}} {{-1  0} { 0  1}}}
    variable tmul		{{{ 1  0} { 0  0}} {{ 0  0} { 0  1}}}
    variable tdiv		{{{ 0  0} { 1  0}} {{ 0  1} { 0  0}}}

    proc init {} {

	# Variables for fundamental constants e, pi, log2

	variable e [[EWorker new] ref]

	set worker \
	    [[math::exact::Mstrict new {{6795705 213440} {6795704 213440}} \
		  [math::exact::PiWorker new]] ref]
	variable pi [[/real [function::sqrt [V new {10005 1}]] $worker] ref]
	$worker unref

	set worker [[Log2Worker new] ref]
	variable log2 [[$worker applyM {{1 1} {1 2}}] ref]
	$worker unref

    }
    init
    rename init {}

    namespace export exactexpr abs1 signum1
}

package provide math::exact 1.0

#-----------------------------------------------------------------------

# exact.test --
#
#    Test cases for the math::exact package
#
# Copyright (c) 2015 by Kevin B. Kenny
#
# See the file "license.terms" for information on usage and redistribution of
# this file, and for a DISCLAIMER OF ALL WARRANTIES.
#
#-----------------------------------------------------------------------------

source [file join \
	    [file dirname [file dirname [file join [pwd] [info script]]]] \
	    devtools testutilities.tcl]

testsNeedTcl     8.6
testsNeedTcltest 2.3

support {
    use     grammar_aycock/aycock-runtime.tcl grammar::aycock::runtime grammar::aycock
    useKeep grammar_aycock/aycock-debug.tcl   grammar::aycock::debug   grammar::aycock
    useKeep grammar_aycock/aycock-build.tcl   grammar::aycock          grammar::aycock
}
testing {
    useLocal exact.tcl             math::exact
}

package require Tcl             8.6
package require grammar::aycock 1.0
package require math::exact     1.0

#-----------------------------------------------------------------------------

namespace eval math::exact::test {

    namespace import ::math::exact::exactexpr

    proc signum {x} {expr {($x > 0) - ($x < 0)}}

    proc leakBaseline {} {
	variable leakBaseline
	foreach o [info commands ::oo::Obj*] {
	    dict set leakBaseline $o {}
	}
	return
    }

    proc leakCheck {} {
	variable leakBaseline
	set trouble {}
	set sep {}
	foreach o [lsort -dictionary [info commands ::oo::Obj*]] {
	    if {![dict exists $leakBaseline $o]} {
		if {[info object isa typeof $o math::exact::counted]} {
		    append trouble $sep "Leaked counted object " \
			$o ": " [$o dump] \n
		} else {
		    append trouble $sep "Leaked object " $o \n
		}
	    }
	}
	if {$trouble ne {}} {
	    return -code error -errorcode {LEAKCHECK} $trouble
	}
	return
    }
			  
    namespace import ::tcltest::test

    test math::exact-1.0 {unit test gcd} {
	math::exact::gcd 2
    } 2
    test math::exact-1.1 {unit test gcd} {
	math::exact::gcd 2 0
    } 2
    test math::exact-1.2 {unit test gcd} {
	math::exact::gcd 0 2
    } 2
    test math::exact-1.3 {unit test gcd} {
	math::exact::gcd 2 3
    } 1
    test math::exact-1.4 {unit test gcd} {
	math::exact::gcd 3 2
    } 1
    test math::exact-1.5 {unit test gcd} {
	math::exact::gcd 21 12
     } 3
    test math::exact-1.6 {unit test gcd} {
	math::exact::gcd 12 21
    } 3
    test math::exact-1.5 {unit test gcd} {
	math::exact::gcd 21 12
     } 3
    test math::exact-1.6 {unit test gcd} {
	math::exact::gcd 12 21
    } 3
    test math::exact-1.7 {unit test gcd} {
	math::exact::gcd 108 66
     } 6
    test math::exact-1.8 {unit test gcd} {
	math::exact::gcd 66 108
    } 6
    test math::exact-1.9 {unit test gcd} {
	math::exact::gcd 66 108 88
    } 2

    test math::exact-2.0 {unit test transpose matrix} {
	math::exact::trans {{0 1} {2 3}}
    } {{0 2} {1 3}}
    test math::exact-2.1 {unit test transpose 2x2x2} {
	math::exact::trans {{{0 1} {2 3}} {{4 5} {6 7}}}
    } {{{0 1} {4 5}} {{2 3} {6 7}}}

    test math::exact-3.1 {unit test determinant} {
	math::exact::determinant {{2 3} {5 7}}
    } -1

    test math::exact-4.1 {unit test reverse} {
	math::exact::reverse {{2 3} {5 7}}
    } {{7 -3} {-5 2}}

    test math::exact-5.1 {unit test veven} {
	math::exact::veven {2 4}
    } 1
    test math::exact-5.2 {unit test veven} {
	math::exact::veven {2 3}
    } 0

    test math::exact-6.1 {unit test meven} {
	math::exact::meven {{2 4} {6 8}}
    } 1
    test math::exact-6.2 {unit test meven} {
	math::exact::meven {{2 3} {6 8}}
    } 0

    test math::exact-7.1 {unit test teven} {
	math::exact::teven {{{2 4} {6 8}} {{10 12} {14 16}}}
    } 1
    test math::exact-7.2 {unit test teven} {
	math::exact::teven {{{2 4} {6 8}} {{10 13} {14 16}}}
    } 0

    test math::exact-8.1 {unit test vhalf} {
	math::exact::vhalf {6 8}
    } {3 4}

    test math::exact-9.1 {unit test mhalf} {
	math::exact::mhalf {{6 8} {10 12}}
    } {{3 4} {5 6}}

    test math::exact-10.1 {unit test thalf} {
	math::exact::thalf {{{6 8} {10 12}} {{14 16} {18 20}}}
    } {{{3 4} {5 6}} {{7 8} {9 10}}}

    test math::exact-11.1 {unit test sign} {
	set trouble {}
	set sep \n
	for {set a -1} {$a <= 1} {incr a} {
	    for {set b -1} {$b <= 1} {incr b} {
		if {$a ==0 && $b == 0} {
		    set sb 0
		} elseif {$a == 0} {
		    set sb [signum $b]
		} elseif {$b == 0} {
		    set sb [signum $a]
		} elseif {$a/$b < 0} {
		    set sb 0
		} else {
		    set sb [signum $a]
		}
		set is [math::exact::sign [list $a $b]]
		if {$is != $sb} {
		    append trouble "sign(" $a "," $b ") is " $is \
			", should be " $sb "\n"
		}
	    }
	}
	set trouble
    } {}

    test math::exact-12.1 {unit test vrefines} {
	set trouble {}
	set sep {}
	for {set a -1} {$a <= 1} {incr a} {
	    for {set b -1} {$b <= 1} {incr b} {
		if {$a ==0 && $b == 0} {
		    set sb 0
		} elseif {$a == 0} {
		    set sb 1
		} elseif {$b == 0} {
		    set sb 1
		} elseif {$a/$b < 0} {
		    set sb 0
		} else {
		    set sb 1
		}
		set is [math::exact::vrefines [list $a $b]]
		if {$is != $sb} {
		    append trouble $sep "vrefines(" $a "," $b ") is " $is \
			", should be " $sb
		    set sep \n
		}
	    }
	}
	set trouble
    } {}

    test math::exact-13.1 {unit test mrefines} {
	math::exact::mrefines {{1 2} {3 4}}
    } 1
    test math::exact-13.2 {unit test mrefines} {
	math::exact::mrefines {{1 2} {-3 -4}}
    } 0
    test math::exact-13.3 {unit test mrefines} {
	math::exact::mrefines {{-1 -2} {-3 -4}}
    } 1
    test math::exact-13.4 {unit test mrefines} {
	math::exact::mrefines {{-1 2} {-3 4}}
    } 0

    test math::exact-14.1 {unit test trefines} {
	math::exact::trefines {{{1 2} {3 4}} {{5 6} {7 8}}}
    } 1
    test math::exact-14.2 {unit test trefines} {
	math::exact::trefines {{{-1 -2} {-3 -4}} {{-5 -6} {-7 -8}}}
    } 1
    test math::exact-14.3 {unit test trefines} {
	math::exact::trefines {{{-1 2} {-3 4}} {{-5 6} {-7 8}}}
    } 0
    test math::exact-14.4 {unit test trefines} {
	math::exact::trefines {{{1 2} {3 4}} {{5 6}} {{-7 -8}}}
    } 0
    test math::exact-14.5 {unit test trefines} {
	math::exact::trefines {{{1 2} {3 4}} {{-5 -6}} {{-7 -8}}}
    } 0
    test math::exact-14.6 {unit test trefines} {
	math::exact::trefines {{{1 2} {-3 -4}} {{-5 -6}} {{-7 -8}}}
    } 0

    test math::exact-15.1 {unit test vlessv} {
	set intervals {
	    {-1 0} {-2 1} {-1 1} {-1 2} {0 1} {2 4} {3 3} {14 7} {1 0}
	}
	set trouble {}
	set sep {}
	set i 0
	foreach a $intervals {
	    set j 0
	    foreach b $intervals {
		set is [math::exact::vlessv $a $b]
		if {[lindex $a 1] == 0 && [lindex $b 1] == 0} {
		    set sb 0
		} else {
		    set sb [expr {$i < $j}]
		}
		if {$is != $sb} {
		    append trouble $sep "vlessv(" $a ";" $b ") is " $is \
			" should be " $sb
		    set sep \n
		}
		incr j
	    }
	    incr i
	}
	set trouble
    } {}
    
    test math::exact-16.1 {unit test mlessm - also tests mlessv} {
	set intervals {
	    {-2 1} {-1 1} {-1 2} {0 1} {2 4} {3 3} {14 7} {1 0}
	}
	set trouble {}
	set sep {}
	set i 0
	foreach a $intervals {
	    set j $i
	    foreach b [lrange $intervals $i end] {
		set k 0
		foreach c $intervals {
		    set l $k
		    foreach d [lrange $intervals $k end] {
			if {[lindex $b 1] == 0 && [lindex $c 1] == 0} {
			    set sb 0
			} else {
			    set sb [expr {$j < $k}]
			}
			set is [math::exact::mlessm [list $b $a] [list $d $c]]
			if {$is != $sb} {
			    append trouble $sep "mlessm(" $a "," $b ";" \
				$c "," $d ") is " $is \
				" should be " $sb " -- " \
				[list i $i j $j k $k l $k]
			    set sep \n
			}
			incr l
		    }
		    incr k
		}
		incr j
	    }
	    incr i
	}
	set trouble
    } {}

    test math::exact-17.1 {unit test vscale} {
	math::exact::vscale {2 3}
    } {2 3}
    test math::exact-17.2 {unit test vscale} {
	math::exact::vscale {4 6}
    } {2 3}
    test math::exact-17.1 {unit test vscale} {
	math::exact::vscale {8 12}
    } {2 3}

    test math::exact-18.1 {unit test mscale} {
	math::exact::mscale {{2 3} {4 5}}
    } {{2 3} {4 5}}
    test math::exact-18.2 {unit test mscale} {
	math::exact::mscale {{4 6} {8 10}}
    } {{2 3} {4 5}}
    test math::exact-18.3 {unit test mscale} {
	math::exact::mscale {{8 12} {16 20}}
    } {{2 3} {4 5}}
    
    test math::exact-19.1 {unit test tscale} {
	math::exact::tscale {{{2 3} {4 5}} {{6 7} {8 9}}}
    } {{{2 3} {4 5}} {{6 7} {8 9}}}
    test math::exact-19.2 {unit test tscale} {
	math::exact::tscale {{{4 6} {8 10}} {{12 14} {16 18}}}
    } {{{2 3} {4 5}} {{6 7} {8 9}}}
    test math::exact-10.3 {unit test tscale} {
	math::exact::tscale {{{8 12} {16 20}} {{24 28} {32 36}}}
    } {{{2 3} {4 5}} {{6 7} {8 9}}}
    
    test math::exact-20.1 {unit test vreduce} {
	math::exact::vreduce {2 3}
    } {2 3}
    test math::exact-20.2 {unit test vreduce} {
	math::exact::vreduce {4 6}
    } {2 3}
    test math::exact-20.1 {unit test vreduce} {
	math::exact::vreduce {8 12}
    } {2 3}

    test math::exact-21.1 {unit test mreduce} {
	math::exact::mreduce {{2 3} {4 5}}
    } {{2 3} {4 5}}
    test math::exact-21.2 {unit test mreduce} {
	math::exact::mreduce {{4 6} {8 10}}
    } {{2 3} {4 5}}
    test math::exact-21.3 {unit test mreduce} {
	math::exact::mreduce {{8 12} {16 20}}
    } {{2 3} {4 5}}
    
    test math::exact-22.1 {unit test treduce} {
	math::exact::treduce {{{2 3} {4 5}} {{6 7} {8 9}}}
    } {{{2 3} {4 5}} {{6 7} {8 9}}}
    test math::exact-22.2 {unit test treduce} {
	math::exact::treduce {{{4 6} {8 10}} {{12 14} {16 18}}}
    } {{{2 3} {4 5}} {{6 7} {8 9}}}
    test math::exact-22.3 {unit test treduce} {
	math::exact::treduce {{{8 12} {16 20}} {{24 28} {32 36}}}
    } {{{2 3} {4 5}} {{6 7} {8 9}}}

    test math::exact-23.1 {unit test mdotv} {
	math::exact::mdotv {{2 3} {4 5}} {10 1}
    } {24 35}

    test math::exact-24.1 {unit test mdotm} {
	math::exact::mdotm {{2 3} {4 5}} {{1000 10} {100 1}}
    } {{2040 3050} {204 305}}

    test math::exact-25.1 {unit test mdott} {
	math::exact::mdott {{1000 10} {100 1}} {{{2 3} {4 5}} {{6 7} {8 9}}}
    } {{{2300 23} {4500 45}} {{6700 67} {8900 89}}}

    test math::exact-26.1 {unit test tleftv} {
	math::exact::tleftv {{{2 3} {4 5}} {{6 7} {8 9}}} {10 1}
    } {{26 37} {48 59}}

    test math::exact-27.1 {unit test trightv} {
	math::exact::trightv {{{2 3} {4 5}} {{6 7} {8 9}}} {10 1}
    } {{24 35} {68 79}}

    test math::exact-28.1 {unit test tleftm} {
	math::exact::tleftm {{{2 3} {4 5}} {{6 7} {8 9}}} {{1000 10} {100 1}}
    } {{{2060 3070} {4080 5090}} {{206 307} {408 509}}}

    test math::exact-29.1 {unit test trightm} {
	math::exact::trightm {{{2 3} {4 5}} {{6 7} {8 9}}} {{1000 10} {100 1}}
    } {{{2040 3050} {204 305}} {{6080 7090} {608 709}}}

    test math::exact-30.1 {unit test mdisjointm} {
	set intervals {
	    {-2 1} {-1 1} {-1 2} {0 1} {2 4} {3 3} {14 7} {1 0}
	}
	set trouble {}
	set sep {}
	set i 0
	foreach a $intervals {
	    set j $i
	    foreach b [lrange $intervals $i end] {
		set k 0
		foreach c $intervals {
		    set l $k
		    foreach d [lrange $intervals $k end] {
			set sb [expr {$j < $k || $l < $i}]
			set is [math::exact::mdisjointm \
				    [list $b $a] [list $d $c]]
			if {$is != $sb} {
			    append trouble $sep "mdisjointm(" $a "," $b ";" \
				$c "," $d ") is " $is \
				" should be " $sb " -- " \
				[list i $i j $j k $k l $k]
			    set sep \n
			}
			incr l
		    }
		    incr k
		}
		incr j
	    }
	    incr i
	}
	set trouble
    } {}

    test math::exact-31.0 {mAsFloat, rational} {
	math::exact::mAsFloat {{1 3} {1 3}}
    } 1/3

    test math::exact-31.1 {mAsFloat, scientificNotation, mantissa, eFormat} {
	set p 0
	set q 1
	set res {}
	for {set i 0} {$i < 16} {incr i} {
	    set r [expr {$p + $q}]
	    if {$q * $q > $p * $r} {
		set m [list [list $q $p] \
			   [list $r $q]]
	    } else {
		set m [list [list $r $q] \
			   [list $q $p]]
	    }
	    lappend res [math::exact::mAsFloat $m]
	    set p $q
	    set q $r
	}
	set res
    } [list \
	   Undetermined 1.e0 1.e0 1.6e0 1.6e0 1.6e0 1.62e0 1.61e0 1.61e0 \
	   1.618e0 1.618e0 1.6180e0 1.6180e0 1.61803e0 1.61803e0 1.61803e0]
    
    test math::exact-31.2 {mAsFloat, scientificNotation, mantissa, eFormat} {
	set p 0
	set q 1
	set res {}
	for {set i 0} {$i < 16} {incr i} {
	    set r [expr {$p + $q}]
	    if {$q * $q > $p * $r} {
		set m [list [list [expr {1000*$q}] $p] \
			   [list [expr {1000*$r}] $q]]
	    } else {
		set m [list [list [expr {1000*$r}] $q] \
			   [list [expr {1000*$q}] $p]]
	    }
	    lappend res [math::exact::mAsFloat $m]
	    set p $q
	    set q $r
	}
	set res
    } [list \
	   Undetermined 1.e3 1.e3 1.6e3 1.6e3 1.6e3 1.62e3 1.61e3 1.61e3 \
	   1.618e3 1.618e3 1.6180e3 1.6180e3 1.61803e3 1.61803e3 1.61803e3]

    test math::exact-31.3 {mAsFloat, scientificNotation, mantissa, eFormat} {
	set p 0
	set q 1
	set res {}
	for {set i 0} {$i < 16} {incr i} {
	    set r [expr {$p + $q}]
	    if {$q * $q > $p * $r} {
		set m [list [list $q [expr {1000*$p}]] \
			   [list $r [expr {1000*$q}]]]
	    } else {
		set m [list [list $r [expr {1000*$q}]] \
			   [list $q [expr {1000*$p}]]]
	    }
	    lappend res [math::exact::mAsFloat $m]
	    set p $q
	    set q $r
	}
	set res
    } [list \
	   Undetermined 1.e-3 1.e-3 1.6e-3 1.6e-3 \
	   1.6e-3 1.62e-3 1.61e-3 1.61e-3 \
	   1.618e-3 1.618e-3 1.6180e-3 1.6180e-3 \
	   1.61803e-3 1.61803e-3 1.61803e-3]
    
    test math::exact-31.4 {mAsFloat, scientificNotation, mantissa, eFormat} {
	set p 0
	set q 1
	set res {}
	for {set i 0} {$i < 16} {incr i} {
	    set r [expr {$p + $q}]
	    set mq [expr {-$q}]
	    set mr [expr {-$r}]
	    if {$q * $q > $p * $r} {
		set m [list [list $mq $p] \
			   [list $mr $q]]
	    } else {
		set m [list [list $mr $q] \
			   [list $mq $p]]
	    }
	    lappend res [math::exact::mAsFloat $m]
	    set p $q
	    set q $r
	}
	set res
    } [list \
	   Undetermined -2.e0 -2.e0 -1.6e0 \
	   -1.7e0 -1.7e0 -1.62e0 -1.62e0 \
	   -1.62e0 -1.618e0 -1.618e0 -1.6180e0 \
	   -1.6181e0 -1.61803e0 -1.61804e0 -1.61804e0]

    test math::exact-31.5 {mAsFloat, 0/0} {
	math::exact::mAsFloat {{0 0} {0 0}}
    } NaN

    test math::exact-31.6 {mAsFloat, infinity} {
	math::exact::mAsFloat {{1 0} {1 0}}
    } Inf

    test math::exact-31.7 {mAsFloat, zero} {
	math::exact::mAsFloat {{0 1} {0 1}}
    } 0

    test math::exact-31.8 {mAsFloat, integer} {
	math::exact::mAsFloat {{2 1} {2 1}}
    } 2

    test math::exact-31.9 {mAsFloat, reverse signs} {
	list [math::exact::mAsFloat {{2 -1} {2 -1}}] \
	    [math::exact::mAsFloat {{-2 -1} {-2 -1}}]
    } {-2 2}

    test math::exact-40.1 {simple expr} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {1}] ref]
	    set result [list [$v asPrint 57] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 1.0000000000000000e0}
    }
    
    test math::exact-40.2 {unary plus} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {+ 1}] ref]
	    set result [list [$v asPrint 57] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 1.0000000000000000e0}
    }
    
    test math::exact-40.3 {unary minus} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {- 1}] ref]
	    set result [list [$v asPrint 57] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {-1 -9.999999999999999e-1}
    }

    test math::exact-40.4 {product} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {2 * 3}] ref]
	    set result [list [$v asPrint 57] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {6 6.000000000000000e0}
    }

    test math::exact-40.5 {quotient} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {2 / 3}] ref]
	    set result [list [$v asPrint 57] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {2/3 6.666666666666666e-1}
    }

    test math::exact-40.6 {associativity of /} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {2 / 3 / 4}] ref]
	    set result [list [$v asPrint 57] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1/6 1.6666666666666667e-1}
    }

    test math::exact-40.7 {associativity of */} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {2 / 3 * 4}] ref]
	    set result [list [$v asPrint 57] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {8/3 2.6666666666666667e0}
    }

    test math::exact-40.8 {associativity of */} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {2 * 3 / 4}] ref]
	    set result [list [$v asPrint 57] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {3/2 1.5000000000000000e0}
    }

    test math::exact-40.9 {sum} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {2 + 3}] ref]
	    set result [list [$v asPrint 57] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {5 5.000000000000000e0}
    }

    test math::exact-40.10 {difference} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {2 - 3}] ref]
	    set result [list [$v asPrint 57] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {-1 -9.999999999999999e-1}
    }

    test math::exact-40.11 {associativity of -} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {2 - 3 - 4}] ref]
	    set result [list [$v asPrint 57] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {-5 -5.000000000000000e0}
    }

    test math::exact-40.12 {associativity of +-} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {2 - 3 + 4}] ref]
	    set result [list [$v asPrint 57] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {3 3.0000000000000001e0}
    }

    test math::exact-40.13 {associativity of +-} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {2 + 3 - 4}] ref]
	    set result [list [$v asPrint 57] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 1.0000000000000000e0}
    }

    test math::exact-40.14 {precedence of +*} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {3 + 5 * 7}] ref]
	    set result [list [$v asPrint 57] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {38 3.800000000000000e1}
    }

    test math::exact-40.15 {precedence of +*} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {3 * 5 + 7}] ref]
	    set result [list [$v asPrint 57] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {22 2.200000000000000e1}
    }

    test math::exact-40.16 {parentheses} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {(2 + 3) * 5}] ref]
	    set result [list [$v asPrint 57] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {25 2.500000000000000e1}
    }

    test math::exact-40.17 {V + E} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {2 + real(-3/5)}] ref]
	    set result [list [$v asPrint 57] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1.4000000000000000e0 1.4000000000000000e0}
    }

    test math::exact-40.18 {V - E} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {2 - real(3/5)}] ref]
	    set result [list [$v asPrint 57] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1.4000000000000000e0 1.4000000000000000e0}
    }

    test math::exact-40.19 {E / E} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {real(3/5)/real(2/5)}] ref]
	    set result [list [$v asPrint 57] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1.5000000000000000e0 1.5000000000000000e0}
    }

    test math::exact-40.20 {E + V} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {real(3/2) + (2/5)}] ref]
	    set result [list [$v asPrint 57] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1.9000000000000000e0 1.9000000000000000e0}
    }
	
    test math::exact-40.21 {E - V} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {real(3/2) - (2/5)}] ref]
	    set result [list [$v asPrint 57] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1.1000000000000000e0 1.1000000000000000e0}
    }

    test math::exact-40.22 {E * V} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {real(3/2) * (2/5)}] ref]
	    set result [list [$v asPrint 57] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {6.0000000000000001e-1 6.0000000000000001e-1}
    }
	
    test math::exact-40.23 {E / V} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {real(3/2) / (5/2)}] ref]
	    set result [list [$v asPrint 57] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {6.0000000000000001e-1 6.0000000000000001e-1}
    }

    test math::exact-40.24 {lexical error} {
	-setup leakBaseline
	-body {
	    set result [list [catch {exactexpr {2 ! 1}} m] $m]
	    leakCheck
	    set result
	}
	-match glob
	-result {1 {invalid character*}}
    }
	
    test math::exact-40.25 {syntax error} {
	-setup leakBaseline
	-body {
	    set result [list [catch {exactexpr {2 $ 1}} m] $m]
	    leakCheck
	    set result
	}
	-match glob
	-result {1 {syntax error*}}
    }
	
    test math::exact-41.1 {square root} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {sqrt(25/16)}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    #leakCheck
	    set result
	}
	-result {1 1.2500000000000000e0}
    }

    test math::exact-41.2 {square root} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {sqrt(2)}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 1.4142135623730950e0}
    }

    test math::exact-41.3 {square root} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {sqrt(2000000)}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 1.4142135623731e3}
    }

    test math::exact-41.4 {square root} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {sqrt(2 / 1000000)}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 1.41421356237309e-3}
    }

    test math::exact-41.5 {square root} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {sqrt(sqrt(1/81))}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 3.3333333333333333e-1}
    }

    test math::exact-41.6 {square root of negative rational} {
	-setup {
	    leakBaseline
	    catch {unset v}
	}
	-body {
	    set v [[exactexpr {sqrt(-1)}] ref]
	    set result [$v asPrint 57]
	    $v unref
	    unset v
	    set result
	}
	-cleanup {
	    if {[info exists v]} {$v unref}
	}
	-match glob
	-returnCodes error
	-result {*negative argument*}
    }

    test math::exact-41.7 {square root of negative real} {
	-setup {
	    leakBaseline
	    catch {unset v}
	}
	-body {
	    set v [[exactexpr {sqrt(-sqrt(81))}] ref]
	    set result [$v asPrint 57]
	    $v unref
	    unset v
	    set result
	}
	-cleanup {
	    if {[info exists v]} {$v unref}
	}
	-match glob
	-returnCodes error
	-result {*negative argument*}
    }

    test math::exact-41.8 {square root, cached result} {
	-setup leakBaseline
	-body {
	    set x [[exactexpr {sqrt(2)}] ref]
	    set y [[exactexpr {$x * $x}] ref]
	    $x unref
	    set result [list [$y asFloat 57] [$y asFloat 57]]
	    $y unref
	    leakCheck
	    set result
	}
	-result {2.0000000000000000e0 2.0000000000000000e0}
    }
    
    test math::exact-42.1 {exponential} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {exp(1)}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 2.7182818284590452e0}
    }

    test math::exact-42.2 {exponential} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {exp(4)}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 5.45981500331442e1}
    }

    test math::exact-42.3 {exponential} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {exp(0)}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 1.0000000000000000e0}
    }

    test math::exact-43.1 {logarithm} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {log(1)}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 0e-18}
    }

    test math::exact-43.2 {logarithm} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {log(2)}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 6.931471805599453e-1}
    }

    test math::exact-43.3 {logarithm} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {log(1/2)}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 -6.931471805599454e-1}
    }

    test math::exact-43.4 {logarithm} {
	-setup {
	    # Consume digits from math::exact::log2 to avoid appearance of
	    # a leak in its cache
	    $math::exact::log2 asFloat 100
	    leakBaseline
	}
	-body {
	    set v [[exactexpr {log(4)}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 1.3862943611198906e0}
    }

    test math::exact-43.5 {logarithm} {
	-setup {
	    # Consume digits from math::exact::log2 to avoid appearance of
	    # a leak in its cache
	    $math::exact::log2 asFloat 100
	    leakBaseline
	}
	-body {
	    set v [[exactexpr {log(1/4)}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 -1.3862943611198907e0}
    }

    test math::exact-43.6 {logarithm} {
	-setup {
	    # Consume digits from math::exact::log2 to avoid appearance of
	    # a leak in its cache
	    $math::exact::log2 asFloat 100
	    leakBaseline
	}
	-body {
	    set v [[exactexpr {log(exp(10))}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 1.0000000000000000e1}
    }

    test math::exact-43.7 {logarithm} {
	-setup {
	    # Consume digits from math::exact::log2 to avoid appearance of
	    # a leak in its cache
	    $math::exact::log2 asFloat 100
	    leakBaseline
	}
	-body {
	    set v [[exactexpr {log(exp(1/10))}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 1.0000000000000000e-1}
    }

    test math::exact-43.8 {logarithm of negative argument} {
	-setup {
	    leakBaseline
	    catch {unset v}
	}
	-body {
	    set v [[exactexpr {log(-sqrt(81))}] ref]
	    set result [$v asPrint 57]
	    $v unref
	    unset v
	    set result
	}
	-cleanup {
	    if {[info exists v]} {$v unref}
	}
	-match glob
	-returnCodes error
	-result {*negative argument*}
    }

    test math::exact-44.1 {pi} {
	-setup {
	    # Consume digits from math::exact::pi to avoid appearance of
	    # a leak in its cache
	    $math::exact::pi asFloat 3000
	    leakBaseline
	}
	-body {
	    set v [[exactexpr {pi()}] ref]
	    set result [$v asFloat 3000]
	    $v unref
	    leakCheck
	    list [string range $result 0 4] \
		[string first 999999 $result] \
		[string range $result end-1 end]
	}
	-result {3.141 763 e0}
    }
    
    test math::exact-44.2 {Ramanujan constant} {
	-setup {
	    # Consume digits from math::exact::pi to avoid appearance of
	    # a leak in its cache
	    $math::exact::pi asFloat 100
	    leakBaseline
	}
	-body {
	    set v [[exactexpr {exp(pi()*sqrt(163))}] ref]
	    set result [$v asFloat 160]
	    $v unref
	    leakCheck
	    set result
	}
	-result 2.625374126407687439999999999992e17
    }
	

    test math::exact-45.1 {tangent} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {tan(0)}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 0e-18}
    }

    test math::exact-45.2 {tangent} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {tan(pi()/4)}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 1.0000000000000000e0}
    }

    test math::exact-45.3 {tangent} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {tan(pi()/-4)}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 -1.0000000000000000e0}
    }

    test math::exact-45.4 {tangent} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {tan(pi()/3)-sqrt(3)}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 0e-18}
    }

    test math::exact-45.5 {tangent} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {tan(-pi()/3)+sqrt(3)}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 0e-18}
    }

    test math::exact-45.6 {tangent} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {1/tan(pi()/2)}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 0e-18}
    }

    test math::exact-45.7 {tangent} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {tan(pi()/2)}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 Undetermined}
    }

    test math::exact-46.1 {sine} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {sin(0)}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 0e-18}
    }

    test math::exact-46.2 {sine} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {sin(pi()/6)}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 5.000000000000000e-1}
    }

    test math::exact-46.3 {sine} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {sin(-pi()/6)}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 -5.000000000000000e-1}
    }

    test math::exact-46.4 {sine} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {sin(pi()/2)}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 1.0000000000000000e0}
    }

    test math::exact-46.5 {sine} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {sin(-pi()/2)}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 -1.0000000000000000e0}
    }

    test math::exact-46.6 {sine} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {sin(pi())}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 0e-18}
    }

    test math::exact-46.7 {sine} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {sin(-pi())}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 0e-18}
    }

    test math::exact-46.8 {sine} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {sin(13*pi()/6)}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 5.000000000000000e-1}
    }

    test math::exact-46.9 {sine} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {sin(-13*pi()/6)}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 -5.000000000000000e-1}
    }
    
    test math::exact-47.1 {cosine} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {cos(0)}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 9.999999999999999e-1}
    }

    test math::exact-47.2 {cosine} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {cos(pi()/3)}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 5.000000000000000e-1}
    }

    test math::exact-47.3 {cosine} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {cos(-pi()/3)}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 5.000000000000000e-1}
    }

    test math::exact-47.4 {cosine} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {cos(pi()/2)}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 0e-18}
    }

    test math::exact-47.5 {cosine} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {cos(-pi()/2)}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 0e-18}
    }

    test math::exact-47.6 {cosine} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {cos(pi())}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 -1.0000000000000000e0}
    }

    test math::exact-47.7 {cosine} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {cos(-pi())}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 -1.0000000000000000e0}
    }

    test math::exact-47.8 {cosine} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {cos(7*pi()/3)}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 5.000000000000000e-1}
    }

    test math::exact-47.9 {cosine} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {cos(-7*pi()/3)}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 5.000000000000000e-1}
    }
    
    test math::exact-45.1 {arctangent} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {atan(0)}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 0e-18}
    }

    test math::exact-45.2 {arctangent} {
	-setup leakBaseline
	-body {
	    # Hack to get $szer as a sign matrix
	    set v [[exactexpr {atan(pi()-pi())}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 0e-18}
    }

    test math::exact-45.3 {arctangent} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {4*atan(1)-pi()}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 0e-18}
    }
	
    test math::exact-45.4 {arctangent} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {4*atan(-1)+pi()}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 0e-18}
    }

    test math::exact-45.5 {arctangent} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {atan(tan(157/100))}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 1.5700000000000000e0}
    }
	    
    test math::exact-45.6 {arctangent, cached} {
	-setup leakBaseline
	-body {
	    set u [[exactexpr {atan(1)}] ref]
	    set v [[exactexpr {$u + $u + $u + $u}] ref]
	    $u unref
	    set result [$v asFloat 57]
	    $v unref
	    leakCheck
	    set result
	}
	-result 3.1415926535897933e0
    }

    test math::exact-46.1 {arcsine} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {asin(0)}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 0e-18}
    }

    test math::exact-46.2 {arcsine} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {asin(1/2)-pi()/6}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 0e-18}
    }

    test math::exact-46.3 {arcsine} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {asin(-1/2)+pi()/6}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 0e-18}
    }

    test math::exact-47.1 {arccosine} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {acos(0)-pi()/2}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 0e-18}
    }

    test math::exact-47.2 {arccosine} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {acos(1/2)-pi()/3}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 0e-18}
    }

    test math::exact-47.3 {arccosine} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {acos(-1/2)-2*pi()/3}] ref]
	    set result [list [$v refcount] [$v asFloat 57]]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1 0e-18}
    }
    
    test math::exact-48.1 {hyperbolic functions} {
	-setup leakBaseline
	-body {
	    set result {}
	    set v [[exactexpr {sinh(0)}] ref]
	    lappend result [$v asFloat 57]
	    $v unref
	    set v [[exactexpr {cosh(0)}] ref]
	    lappend result [$v asFloat 57]
	    $v unref
	    set v [[exactexpr {tanh(0)}] ref]
	    lappend result [$v asFloat 57]
	    $v unref
	    leakCheck
	    set result
	}
	-result {0e-18 1.0000000000000000e0 0e-18}
    }

    test math::exact-48.2 {hyperbolic functions} {
	-setup leakBaseline
	-body {
	    set result {}
	    set v [[exactexpr {sinh(1)}] ref]
	    lappend result [$v asFloat 57]
	    $v unref
	    set v [[exactexpr {cosh(1)}] ref]
	    lappend result [$v asFloat 57]
	    $v unref
	    set v [[exactexpr {tanh(1)}] ref]
	    lappend result [$v asFloat 57]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1.1752011936438014e0 1.5430806348152437e0 7.615941559557649e-1}
    }

    test math::exact-48.3 {hyperbolic functions} {
	-setup leakBaseline
	-body {
	    set result {}
	    set v [[exactexpr {sinh(-1)}] ref]
	    lappend result [$v asFloat 57]
	    $v unref
	    set v [[exactexpr {cosh(-1)}] ref]
	    lappend result [$v asFloat 57]
	    $v unref
	    set v [[exactexpr {tanh(-1)}] ref]
	    lappend result [$v asFloat 57]
	    $v unref
	    leakCheck
	    set result
	}
	-result {-1.1752011936438015e0 1.5430806348152437e0 -7.615941559557649e-1}
    }

    test math::exact-49.1 {asinh} {
	-setup leakBaseline
	-body {
	    set result {}
	    set v [[exactexpr {asinh(-1)}] ref]
	    lappend result [$v asFloat 57]
	    $v unref
	    set v [[exactexpr {asinh(0)}] ref]
	    lappend result [$v asFloat 57]
	    $v unref
	    set v [[exactexpr {asinh(1)}] ref]
	    lappend result [$v asFloat 57]
	    $v unref
	    leakCheck
	    set result
	}
	-result {-8.813735870195431e-1 0e-18 8.813735870195430e-1}
    }
	
    test math::exact-50.1 {acosh} {
	-setup leakBaseline
	-body {
	    set result {}
	    set v [[exactexpr {acosh(3/2)}] ref]
	    lappend result [$v asFloat 57]
	    $v unref
	    set v [[exactexpr {acosh(2)}] ref]
	    lappend result [$v asFloat 57]
	    $v unref
	    set v [[exactexpr {acosh(3)}] ref]
	    lappend result [$v asFloat 57]
	    $v unref
	    leakCheck
	    set result
	}
	-result {9.624236501192069e-1 1.3169578969248167e0 1.7627471740390860e0}
    }

    test math::exact-51.1 {atanh} {
	-setup leakBaseline
	-body {
	    set result {}
	    set v [[exactexpr {atanh(-1/2)}] ref]
	    lappend result [$v asFloat 57]
	    $v unref
	    set v [[exactexpr {atanh(0)}] ref]
	    lappend result [$v asFloat 57]
	    $v unref
	    set v [[exactexpr {atanh(1/2)}] ref]
	    lappend result [$v asFloat 57]
	    $v unref
	    leakCheck
	    set result
	}
	-result {-5.4930614433405485e-1 0e-18 5.4930614433405485e-1}
    }
    
    test math::exact-52.1 {e} {
	-setup {
	    # don't report cached digits of e as a leak
	    $math::exact::e asPrint 100;
	    leakBaseline
	}
	-body {
	    set v [[exactexpr {e()}] ref]
	    set result [$v asPrint 57]
	    $v unref
	    leakCheck
	    set result
	}
	-result 2.7182818284590452e0
    }
    
    test math::exact-52.2 {e} {
	-setup {
	    # don't report cached digits of e as a leak
	    $math::exact::e asPrint 100;
	    leakBaseline
	}
	-body {
	    set v [[exactexpr {log(e())}] ref]
	    set result [$v asPrint 57]
	    $v unref
	    leakCheck
	    set result
	}
	-result 1.0000000000000000e0
    }

    test math::exact-52.2 {e} {
	-setup {
	    # don't report cached digits of e as a leak
	    $math::exact::e asPrint 100;
	    leakBaseline
	}
	-body {
	    set v [[exactexpr {asinh((e() - 1/e()) / 2)}] ref]
	    set result [$v asPrint 57]
	    $v unref
	    leakCheck
	    set result
	}
	-result 1.0000000000000000e0
    }

    test math::exact-53.1 {real**real} {
	-setup {
	    # Consume digits from math::exact::e and math::exact::log2
	    # to avoid appearance of a leak in the cache
	    $math::exact::e asFloat 100
	    $math::exact::log2 asFloat 100
	    leakBaseline
	}
	-body {
	    set v [[exactexpr {e() ** log(2)}] ref]
	    set result [$v asPrint 57]
	    $v unref
	    leakCheck
	    set result
	}
	-result 2.0000000000000000e0
    }

    test math::exact-53.2 {rational**real} {
	-setup {
	    # Consume digits from math::exact::e
	    # to avoid appearance of a leak in the cache
	    $math::exact::e asFloat 100
	    leakBaseline
	}
	-body {
	    set v [[exactexpr {2 ** e()}] ref]
	    set result [$v asPrint 57]
	    $v unref
	    leakCheck
	    set result
	}
	-result 6.580885991017921e0
    }

    test math::exact-53.3 {real**1} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {sqrt(4) ** 1}] ref]
	    set result [$v asPrint 57]
	    $v unref
	    leakCheck
	    set result
	}
	-result 2.0000000000000000e0
    }

    test math::exact-53.4 {real**-1} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {sqrt(4) ** (-1)}] ref]
	    set result [$v asPrint 57]
	    $v unref
	    leakCheck
	    set result
	}
	-result 5.000000000000000e-1
    }

    test math::exact-53.5 {real**0} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {sqrt(4) ** 0}] ref]
	    set result [$v asPrint 57]
	    $v unref
	    leakCheck
	    set result
	}
	-result 1
    }

    test math::exact-53.6 {real**+int} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {sqrt(4)**2}] ref]
	    set result [$v asPrint 57]
	    $v unref
	    leakCheck
	    set result
	}
	-result 4.000000000000000e0
    }

    test math::exact-53.7 {real**+int} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {sqrt(4)**5}] ref]
	    set result [$v asPrint 57]
	    $v unref
	    leakCheck
	    set result
	}
	-result 3.200000000000000e1
    }

    test math::exact-53.6 {real**-int} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {sqrt(4)**-2}] ref]
	    set result [$v asPrint 57]
	    $v unref
	    leakCheck
	    set result
	}
	-result 2.5000000000000000e-1
    }

    test math::exact-53.7 {real**+int} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {sqrt(4)**-5}] ref]
	    set result [$v asPrint 57]
	    $v unref
	    leakCheck
	    set result
	}
	-result 3.125000000000000e-2
    }

    test math::exact-53.8 {real**rational} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {sqrt(64)**(10/3)}] ref]
	    set result [$v asPrint 57]
	    $v unref
	    leakCheck
	    set result
	}
	-result 1.02400000000000e3
    }

    test math::exact-53.9 {real**rational} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {sqrt(64)**(1/-3)}] ref]
	    set result [$v asPrint 57]
	    $v unref
	    leakCheck
	    set result
	}
	-result 5.000000000000000e-1
    }
	
    test math::exact-53.10 {real**integer, accidental} {
	-setup leakBaseline
	-body {
	    set v [[math::exact::real**rat [exactexpr {3}] 2 1] ref]
	    set result [$v asPrint 57]
	    $v unref
	    leakCheck
	    set result
	}
	-result 9
    }
	
    test math::exact-53.11 {zero to zero power} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {0**0}] ref]
	    set result [$v asPrint 57]
	    $v unref
	    leakCheck
	    set result
	}
	-returnCodes error
	-result "zero to zero power"
    }

    test math::exact-53.12 {zero to infinite power} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {0**(1/0)}] ref]
	    set result [$v asPrint 57]
	    $v unref
	    leakCheck
	    set result
	}
	-returnCodes error
	-result "zero to infinite power"
    }
	
    test math::exact-53.13 {zero to rational power} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {0**(1/2)}] ref]
	    set result [$v asPrint 57]
	    $v unref
	    leakCheck
	    set result
	}
	-result 0
    }
	
    test math::exact-53.14 {infinity to zero power} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {(1/0)**0}] ref]
	    set result [$v asPrint 57]
	    $v unref
	    leakCheck
	    set result
	}
	-returnCodes error
	-result "infinity to zero power"
    }
	
    test math::exact-53.15 {infinity to negative power} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {(1/0)**-1}] ref]
	    set result [$v asPrint 57]
	    $v unref
	    leakCheck
	    set result
	}
	-result 0
    }
	
    test math::exact-53.15 {infinity to positive power} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {(1/0)**1}] ref]
	    set result [$v asPrint 57]
	    $v unref
	    leakCheck
	    set result
	}
	-result Inf
    }
	
    test math::exact-53.16 {rational to zero power} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {(2/3)**0}] ref]
	    set result [$v asPrint 57]
	    $v unref
	    leakCheck
	    set result
	}
	-result 1
    }

    test math::exact-53.17 {rational power of negative real argument} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {(-sqrt(64))**(1/3)}] ref]
	    set result [$v asPrint 57]
	    $v unref
	    leakCheck
	    set result
	}
	-result -2.0000000000000000e0
    }

    test math::exact-53.18 {rational power of argument near zero} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {log(exp(1/8))**(1/3)}] ref]
	    set result [$v asPrint 57]
	    $v unref
	    leakCheck
	    set result
	}
	-result 5.000000000000000e-1
    }

    test math::exact-53.19 {negative real to real power} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {(-sqrt(4))**(1/2)}] ref]
	    set result [$v asPrint 57]
	    $v unref
	    leakCheck
	    set result
	}
	-returnCodes error
	-result "negative number to real power"
    }
	
    test math::exact-53.20 {rational to zero power} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {(2/3)**0}] ref]
	    set result [$v asPrint 57]
	    $v unref
	    leakCheck
	    set result
	}
	-result 1
    }

    test math::exact-53.21 {rational to positive integer power} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {(2/3)**2}] ref]
	    set result [$v asPrint 57]
	    $v unref
	    leakCheck
	    set result
	}
	-result 4/9
    }
	
    test math::exact-53.22 {rational to negative integer power} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {(2/3)**-2}] ref]
	    set result [$v asPrint 57]
	    $v unref
	    leakCheck
	    set result
	}
	-result 9/4
    }

    test math::exact-53.23 {rational to rational} {
	-setup leakBaseline
	-body {
	    set v [[exactexpr {(-8)**(1/3)}] ref]
	    set result [$v asPrint 57]
	    $v unref
	    leakCheck
	    set result
	}
	-result -2.0000000000000000e0
    }

    test math::exact-53.24 {real to 0/0} {
	-setup leakBaseline
	-body {
	    set bad [[math::exact::V new {0 0}] ref]
	    set v [[exactexpr {sqrt(2)**$bad}] ref]
	    $bad unref
	    set result [$v asPrint 57]
	    $v unref
	    leakCheck
	    set result
	}
	-returnCodes error
	-result {zero divided by zero}
    }
	
    test math::exact-53.24 {rational to 0/0} {
	-setup leakBaseline
	-body {
	    set bad [[math::exact::V new {0 0}] ref]
	    set v [[exactexpr {(1/2)**$bad}] ref]
	    $bad unref
	    set result [$v asPrint 57]
	    $v unref
	    leakCheck
	    set result
	}
	-returnCodes error
	-result {zero divided by zero}
    }
	
    test math::exact-53.26 {unit test - rat**int (2/3)**0} {
	-setup leakBaseline
	-body {
	    set v [[math::exact::rat**int 2 3 0] ref]
	    set result [$v asPrint 57]
	    $v unref
	    leakCheck
	    set result
	}
	-result {1}
    }

    test math::exact-53.27 {rational powers - normalize base and exponent} {
	-setup leakBaseline
	-body {
	    set p [[math::exact::V new {-2 -1}] ref]
	    set q [[math::exact::V new {-3 -1}] ref]
	    set r [[exactexpr {$p ** $q}] ref]
	    $p unref
	    $q unref
	    set result [$r asPrint 57]
	    $r unref
	    leakCheck
	    set result
	}
	-result 8
    }

    test math::exact-54.1 {abs1, signum1} {
	-setup leakBaseline
	-body {
	    set p [[exactexpr {0}] ref]
	    set q [[exactexpr {2}] ref]
	    while 1 {
		set t [[exactexpr {($q-$p) * 10**36}] ref]
		set f [math::exact::abs1 $t]; $t unref
		if {!$f} break
		set x [[exactexpr {($p+$q)/2}] ref]
		set resid [[exactexpr {$x*$x-2}] ref]
		set t [[exactexpr {$resid * 10**36}] ref]
		if {[math::exact::signum1 $t] > 0} {
		    $q unref; set q $x
		} else {
		    $p unref; set p $x
		}
		$t unref; $resid unref
	    }
	    set result [$p asFloat 100]
	    $p unref
	    $q unref
	    leakCheck
	    set result
	}
	-result 1.41421356237309504880168872421e0
    }

    # following are demos that I don't know where to put, yet

    if 0 {
    set p 1
	for {set i 0} {$i < 20} {incr i} {
	    set f [expr {sin(0.01 * $p* acos(-1))}]
	    set v [[exactexpr "sin($p * pi() / 100)"] ref]
	    set a [$v asPrint 57]
	    set r [expr {$f - $a}]
	    puts "i: $i p: $p float: $f exact: $a difference: $r"
	    $v unref
	    set p [expr {11 * $p}]
    }
    }

    if 0 {
	for {set x 100} {$x <= 12200} {incr x 100} {
	    set ex [[exactexpr $x] ref]
	    puts "x $x ex [$ex asPrint 57]"
	    set fa [expr {-(double($x)**-4)}]
	    set ea [[exactexpr {-($ex**-4)}] ref]
	    puts "fa $fa ea [$ea asPrint 57]"
	    set fb [expr {exp($fa)}]
	    set eb [[exactexpr {exp($ea)}] ref]
	    puts "fb $fb eb [$eb asPrint 120]"
	    set fc [expr {log($fb)}]
	    set ec [[exactexpr {log($eb)}] ref]
	    puts "fc $fc ec [$ec asPrint 120]"
	    catch {expr {(-$fc) ** -0.25}} ff
	    set ef [[exactexpr {(-$ec)**(-1/4)}] ref]
	    puts [format "kahan's function: %s %g" $ff [$ef asFloat 28]]
	    $ef unref
	    $ec unref
	    $eb unref
	    $ea unref
	    $ex unref
	}
    }

    if 0 {
	set x0 4.0
	set x1 4.25
	set ex0 [[exactexpr 4] ref]
	set ex1 [[exactexpr 4+25/100] ref]
	for {set i 1} {$i < 100} {incr i} {
	    set x2 [expr {108. - (815. - 1500. / $x0) / $x1}]
	    set x0 $x1
	    set x1 $x2
	    set ex2 [[exactexpr {108 - (815 - 1500 / $ex0) / $ex1}] ref]
	    $ex0 unref
	    set ex0 $ex1
	    set ex1 $ex2
	    puts "$i $x2 [$ex2 asFloat 57]"
	}
	$ex0 unref
	$ex1 unref
    }

    testsuiteCleanup

}

#-----------------------------------------------------------------------------

# End of test cases

testsuiteCleanup

# Exit if running this test standalone, to allow for Nagelfar coverage
if {$::argv0 eq [info script]} {
    exit
}

# Local Variables:
# mode: tcl
# End:

if {![package vsatisfies [package provide Tcl] 8.5]} {return}
package ifneeded math::calculus::symdiff 1.0   [list source [file join $dir symdiff.tcl]]
package ifneeded math::bigfloat          2.0.2 [list source [file join $dir bigfloat2.tcl]]
package ifneeded math::numtheory         1.0   [list source [file join $dir numtheory.tcl]]
package ifneeded math::decimal           1.0.3 [list source [file join $dir decimal.tcl]]

if {![package vsatisfies [package require Tcl] 8.6]} {return}
package ifneeded math::exact             1.0   [list source [file join $dir exact.tcl]]