package ICC::Support::spline;

use strict;
use Carp;

our $VERSION = 0.02;

# revised 2017-03-29
#
# Copyright © 2004-2018 by William B. Birkett

# add development directory
use lib 'lib';

# inherit from Shared
use parent qw(ICC::Shared);

# add complex math functions
use Math::Complex qw(root sqrt cbrt Im Re);

# enable static variables
use feature 'state';

# create new spline object
# hash keys are: 'input_range', 'output_values'
# parameters: ([ref_to_attribute_hash])
# returns: (object_reference)
sub new {

	# get object class
	my $class = shift();

	# create empty spline object
	my $self = [
		{}, # object header
		[], # input range
		[], # output values
		[], # derivative values
		[], # min/max output values
		[], # parametric derivative matrix
	];

	# if one parameter, a hash reference
	if (@_ == 1 && ref($_[0]) eq 'HASH') {
		
		# set object contents from hash
		_new_from_hash($self, $_[0]);
		
	} elsif (@_) {
		
		# error
		croak('invalid parameter(s)');
		
	}

	# return blessed object
	return(bless($self, $class));

}

# write spline object to ICC profile
# note: writes an equivalent curv object
# parameters: (ref_to_parent_object, file_handle, ref_to_tag_table_entry)
sub write_fh {

	# get tag reference
	my $self = shift();

	# verify 3 parameters
	(@_ == 3) || croak('wrong number of parameters');

	# write spline data to profile
	_writeICCspline($self, @_);

}

# get tag size (for writing to profile)
# note: writes an equivalent curv object
# returns: (tag_size)
sub size {

	# get parameters
	my ($self) = @_;

	# get count value (defaults to 4096, the maximum allowed in an 'mft2' tag)
	my $n = $self->[0]{'curv_points'} // 4096;

	# return size
	return(12 + $n * 2);

}

# compute curve function
# parameters: (input_value)
# returns: (output_value)
sub transform {

	# get parameters
	my ($self, $in) = @_;

	# if an extrapolated solution (s < 0)
	if (($self->[1][1] > $self->[1][0] && $in < $self->[1][0]) || ($self->[1][1] < $self->[1][0] && $in > $self->[1][0])) {
		
		# return extrapolated value
		return($self->[2][0] + $self->[3][0] * $#{$self->[2]} * ($in - $self->[1][0])/($self->[1][1] - $self->[1][0]));
		
	# if an extrapolated solution (s >= 1)
	} elsif (($self->[1][1] > $self->[1][0] && $in >= $self->[1][1]) || ($self->[1][1] < $self->[1][0] && $in <= $self->[1][1])) {
		
		# return extrapolated value
		return($self->[2][-1] + $self->[3][-1] * $#{$self->[2]} * ($self->[1][1] - $in)/($self->[1][0] - $self->[1][1]));
		
	} else {
		
		# return interpolated value
		return(_fwd($self, POSIX::modf($#{$self->[2]} * ($in - $self->[1][0])/($self->[1][1] - $self->[1][0]))));
		
	}
	
}

# compute inverse curve function
# note: there may be multiple solutions
# parameters: (input_value)
# returns: (output_value -or- array_of_output_values)
sub inverse {

	# get parameters
	my ($self, $in) = @_;

	# local variables
	my ($n, $sf, @sol, @t);

	# get number of segments
	$n = $#{$self->[2]};

	# compute scale factor
	$sf = ($self->[1][1] - $self->[1][0])/$n;

	# if an extrapolated solution (s < 0)
	if (($in < $self->[2][0] && $self->[3][0] > 0) || ($in > $self->[2][0] && $self->[3][0] < 0)) {
		
		# add solution
		push(@sol, $self->[1][0] + $sf * ($in - $self->[2][0])/$self->[3][0]);
		
	}

	# for each segment (0 <= s < 1)
	for my $i (0 .. $n - 1) {
		
		# if input is lower knot
		if ($in == $self->[2][$i]) {
			
			# add solution
			push(@sol, $self->[1][0] + $i * $sf);
			
		}
		
		# if input lies within the y-value range
		if ($in >= $self->[4][$i][0] && $in <= $self->[4][$i][1]) {
			
			# compute inverse values (t-values)
			@t = _rev($self, $in, $i);
			
			# add solution(s)
			push(@sol, map {$self->[1][0] + ($i + $_) * $sf} @t);
			
		}
		
	}

	# if input is last knot (s == 1)
	if ($in == $self->[2][-1]) {
		
		# add solution
		push(@sol, $self->[1][1]);
		
	}

	# if an extrapolated solution (s > 1)
	if (($in > $self->[2][-1] && $self->[3][-1] > 0) || ($in < $self->[2][-1] && $self->[3][-1] < 0)) {
		
		# add solution
		push(@sol, $self->[1][1] + $sf * ($in - $self->[2][-1])/$self->[3][-1]);
		
	}

	# return result (array or first solution)
	return(wantarray ? @sol : $sol[0]);

}

# compute curve derivative
# parameters: (input_value)
# returns: (derivative_value)
sub derivative {

	# get parameters
	my ($self, $in) = @_;

	# if an extrapolated solution (s < 0)
	if (($self->[1][1] > $self->[1][0] && $in < $self->[1][0]) || ($self->[1][1] < $self->[1][0] && $in > $self->[1][0])) {
		
		# return endpoint derivative
		return($self->[3][0] * $#{$self->[2]}/($self->[1][1] - $self->[1][0]));
		
	# if an extrapolated solution (s >= 1)
	} elsif (($self->[1][1] > $self->[1][0] && $in >= $self->[1][1]) || ($self->[1][1] < $self->[1][0] && $in <= $self->[1][1])) {
		
		# return endpoint derivative
		return($self->[3][-1] * $#{$self->[2]}/($self->[1][1] - $self->[1][0]));
		
	} else {
		
		# return interpolated derivative value
		return(_derv($self, POSIX::modf($#{$self->[2]} * ($in - $self->[1][0])/($self->[1][1] - $self->[1][0]))) * $#{$self->[2]}/($self->[1][1] - $self->[1][0]));
		
	}
	
}

# compute parametric partial derivatives
# parameters: (input_value)
# returns: (partial_derivative_array)
sub parametric {

	# get parameters
	my ($self, $in) = @_;

	# local variables
	my ($low, $t, $tc, $h00, $h01, $h10, $h11, @pj);

	# update parametric derivative matrix, if necessary
	_objpara($self) if (! defined($self->[5]) || $#{$self->[5]} != $#{$self->[2]});

	# if an extrapolated solution (s < 0)
	if (($self->[1][1] > $self->[1][0] && $in < $self->[1][0]) || ($self->[1][1] < $self->[1][0] && $in > $self->[1][0])) {
		
		# for each output value (parameter)
		for my $i (0 .. $#{$self->[2]}) {
			
			# compute partial derivative
			$pj[$i] = ($i == 0) + $self->[5][0][$i] * $#{$self->[2]} * ($in - $self->[1][0])/($self->[1][1] - $self->[1][0]);
			
		}
		
	# if an extrapolated solution (s >= 1)
	} elsif (($self->[1][1] > $self->[1][0] && $in >= $self->[1][1]) || ($self->[1][1] < $self->[1][0] && $in <= $self->[1][1])) {
		
		# for each output value (parameter)
		for my $i (0 .. $#{$self->[2]}) {
			
			# compute partial derivative
			$pj[$i] = ($i == $#{$self->[2]}) + $self->[5][-1][$i] * $#{$self->[2]} * ($in - $self->[1][1])/($self->[1][1] - $self->[1][0]);
			
		}
		
	} else {
		
		# compute position (0 - 1) and lower knot index
		($t, $low) = POSIX::modf($#{$self->[2]} * ($in - $self->[1][0])/($self->[1][1] - $self->[1][0]));
		
		# compute Hermite coefficients
		$tc = 1 - $t;
		$h00 = (1 + 2 * $t) * $tc * $tc;
		$h01 = 1 - $h00;
		$h10 = $t * $tc * $tc;
		$h11 = -$t * $t * $tc;
		
		# for each output value (parameter)
		for my $i (0 .. $#{$self->[2]}) {
			
			# compute partial derivative
			$pj[$i] = $h00 * ($low == $i) + $h01 * ($low + 1 == $i) + $h10 * $self->[5][$low][$i] + $h11 * $self->[5][$low + 1][$i];
			
		}
		
	}

	# return
	return(@pj);

}

# get/set reference to header hash
# parameters: ([ref_to_new_hash])
# returns: (ref_to_hash)
sub header {

	# get object reference
	my $self = shift();

	# if there are parameters
	if (@_) {
	
		# if one parameter, a hash reference
		if (@_ == 1 && ref($_[0]) eq 'HASH') {
		
			# set header to copy of hash
			$self->[0] = Storable::dclone(shift());
		
		} else {
		
			# error
			croak('parameter must be a hash reference');
		
		}
	
	}

	# return reference
	return($self->[0]);

}

# get/set input range array reference
# parameters: ([ref_to_array])
# returns: (ref_to_array)
sub range {

	# get object reference
	my $self = shift();

	# if one parameter supplied
	if (@_ == 1) {
		
		# verify parameter is an array reference
		(ref($_[0]) eq 'ARRAY') || croak('parameter not an array reference');
		
		# verify array contents
		(@{$_[0]} == 2 && (@{$_[0]} == grep {Scalar::Util::looks_like_number($_)} @{$_[0]}) && $_[0]->[0] != $_[0]->[1]) || croak('array must contain 2 unequal numeric values');
		
		# store input range
		$self->[1] = [@{$_[0]}];
		
	} elsif (@_) {
		
		# error
		croak('too many parameters');
		
	}
	
	# return array reference
	return($self->[1]);

}

# get/set output value array reference
# parameters: ([ref_to_array])
# returns: (ref_to_array)
sub array {

	# get object reference
	my $self = shift();

	# if one parameter supplied
	if (@_ == 1) {
		
		# verify parameter is an array reference
		(ref($_[0]) eq 'ARRAY') || croak('parameter not an array reference');
		
		# verify array contents
		(@{$_[0]} >= 2 && (@{$_[0]} == grep {Scalar::Util::looks_like_number($_)} @{$_[0]})) || croak('array must contain 2 or more numeric values');
		
		# store output values
		$self->[2] = [@{$_[0]}];
		
		# recompute knot derivatives
		_objderv($self);
		
		# recompute segment min/max y-values
		_minmax($self);
		
	} elsif (@_) {
		
		# error
		croak('too many parameters');
		
	}
	
	# return array reference
	return($self->[2]);

}

# check if monotonic
# returns min/max values
# returns s-values by default
# returns input-values if format is 'x'
# returns output-values if format is 'y'
# curve is monotonic if no min/max values
# parameters: ([format])
# returns: (array_of_values)
sub monotonic {

	# get object reference
	my ($self, $fmt) = @_;

	# if format undefined or 's'
	if (! defined($fmt) || $fmt eq 's') {
		
		# return s-values
		return(_minmax($self));
		
	# if format 'x'
	} elsif ($fmt eq 'x') {
		
		# return x-values
		return(map {(1 - $_/$#{$self->[2]}) * $self->[1][0] + $_/$#{$self->[2]} * $self->[1][1]} _minmax($self));
		
	# if format 'y'
	} elsif ($fmt eq 'y') {
		
		# return y-values
		return(map {_fwd($self, POSIX::modf($_))} _minmax($self));
		
	} else {
		
		# error
		croak("unsupported format for min/max values");
		
	}
	
}

# make table for 'curv' objects
# assumes curve domain/range is (0 - 1)
# parameters: (number_of_table_entries, [direction])
# returns: (ref_to_table_array)
sub table {

	# get parameters
	my ($self, $n, $dir) = @_;

	# local variables
	my ($up, $table);

	# validate number of table entries
	($n == int($n) && $n >= 2) || carp('invalid number of table entries');

	# purify direction flag
	$dir = ($dir) ? 1 : 0;

	# array upper index
	$up = $n - 1;

	# for each table entry
	for my $i (0 .. $up) {
		
		# compute table value
		$table->[$i] = _transform($self, $dir, $i/$up);
		
	}

	# return table reference
	return($table);

}

# make 'curv' object
# assumes curve domain/range is (0 - 1)
# parameters: (number_of_table_entries, [direction])
# returns: (ref_to_curv_object)
sub curv {

	# return 'curv' object reference
	return(ICC::Profile::curv->new(table(@_)));

}

# normalize transform
# adjusts object values linearly
# so that endpoint values are 0 or 1
# adjusts input-range and output-values by default
# adjusts input-values if format is 'x'
# adjusts output-values if format is 'y'
# parameters: ([format])
sub normalize {

	# get parameters
	my ($self, $fmt) = @_;

	# local variables
	my ($off, $range);

	# verify flag parameter
	(! defined($fmt) || $fmt eq 'x' || $fmt eq 'y') || croak('invalid normalize format parameter');

	# if range selected
	if (! defined($fmt) || $fmt eq 'x') {
		
		# if range is increasing
		if ($self->[1][1] > $self->[1][0]) {
			
			# set range (0 - 1)
			$self->[1][0] = 0;
			$self->[1][1] = 1;
			
		} else {
			
			# set range (1 - 0)
			$self->[1][0] = 1;
			$self->[1][1] = 0;
			
		}
		
	}

	# if output values selected
	if (! defined($fmt) || $fmt eq 'y') {
		
		# if output values are increasing
		if ($self->[2][-1] > $self->[2][0]) {
			
			# set offset
			$off = $self->[2][0];
			
			# compute range
			$range = ($self->[2][-1] - $self->[2][0]);
			
		} else {
			
			# set offset
			$off = $self->[2][-1];
			
			# compute range
			$range = ($self->[2][0] - $self->[2][-1]);
			
		}
		
		# verify range
		($range != 0) || croak('cannot normalize output values');
		
		# for each output value
		for my $i (0 .. $#{$self->[2]}) {
			
			# adjust using offset and scale factor
			$self->[2][$i] = ($self->[2][$i] - $off)/$range;
			
		}
		
		# recompute knot derivatives
		_objderv($self);
		
		# recompute segment min/max y-values
		_minmax($self);
		
		# clear parametric derivative matrix
		$self->[5] = [];
		
	}
	
}

# print object contents to string
# format is an array structure
# parameter: ([format])
# returns: (string)
sub sdump {

	# get parameters
	my ($self, $p) = @_;

	# local variables
	my ($s, $fmt);

	# resolve parameter to an array reference
	$p = defined($p) ? ref($p) eq 'ARRAY' ? $p : [$p] : [];

	# get format string
	$fmt = defined($p->[0]) && ! ref($p->[0]) ? $p->[0] : 'undef';

	# set string to object ID
	$s = sprintf("'%s' object, (0x%x)\n", ref($self), $self);

	# return
	return($s);

}

# compute parametric partial derivatives
# direction: 0 - normal, 1 - inverse
# parameters: (object_reference, direction, input_value)
# returns: (partial_derivative_array)
sub _parametric {

	# get parameters
	my ($self, $dir, $in) = @_;

	# if inverse direction
	if ($dir) {
		
		# unimplemented function error
		croak('unimplemented function');
		
	} else {
		
		# return array of partial derivatives
		return(parametric($self, $in));
		
	}
	
}

# combined derivative
# direction: 0 - normal, 1 - inverse
# parameters: (object_reference, direction, input_value)
# returns: (derivative_value)
sub _derivative {

	# get parameters
	my ($self, $dir, $in) = @_;

	# if inverse direction
	if ($dir) {
		
		# unimplemented function error
		croak('unimplemented function');
		
	} else {
		
		# return derivative
		return(derivative($self, $in));
		
	}
	
}

# combined transform
# direction: 0 - normal, 1 - inverse
# parameters: (object_reference, direction, input_value)
# returns: (output_value)
sub _transform {

	# get parameters
	my ($self, $dir, $in) = @_;

	# if inverse direction
	if ($dir) {
		
		# return inverse
		return(scalar(inverse($self, $in)));
		
	} else {
		
		# return transform
		return(transform($self, $in));
		
	}
	
}

# compute knot derivatives for natural spline
# parameters: (ref_to_object)
sub _objderv {

	# get parameter
	my $self = shift();

	# local variables
	my ($ix, $rhs, $info, $derv);

	# verify object has two or more knots
	(($ix = $#{$self->[2]}) > 0) || croak('object must have two or more knots');

	# check if ICC::Support::Lapack module is loaded
	state $lapack = defined($INC{'ICC/Support/Lapack.pm'});

	# make empty Math::Matrix object
	$rhs = bless([], 'Math::Matrix');

	# for each input element
	for my $i (1 .. $ix - 1) {
		
		# compute rhs (3 * (y[i + 1] - y[i - 1]))
		$rhs->[$i][0] = 3 * ($self->[2][$i + 1] - $self->[2][$i - 1]);
		
	}

	# set rhs endpoint values
	$rhs->[0][0] = 3 * ($self->[2][1] - $self->[2][0]);
	$rhs->[$ix][0] = 3 * ($self->[2][$ix] - $self->[2][$ix - 1]);

	# if ICC::Support::Lapack module is loaded
	if ($lapack) {
		
		# solve for derivative matrix
		($info, $derv) = ICC::Support::Lapack::trisolve([(1) x $ix], [2, (4) x ($ix - 1), 2], [(1) x $ix], $rhs);
		
	# otherwise, use Math::Matrix module
	} else {
		
		# solve for derivative matrix
		$derv = Math::Matrix->tridiagonal([2, (4) x ($ix - 1), 2])->concat($rhs)->solve();
		
	}

	# for each knot
	for my $i (0 .. $ix) {
		
		# set derivative value
		$self->[3][$i] = $derv->[$i][0];
		
	}

}

# compute partial derivative matrix for natural spline
# parameters: (ref_to_object)
sub _objpara {

	# get parameter
	my $self = shift();

	# local variables
	my ($ix, $rhs, $info, $derv);

	# verify object has two or more knots
	(($ix = $#{$self->[2]}) > 0) || croak('object must have two or more knots');

	# check if ICC::Support::Lapack module is loaded
	state $lapack = defined($INC{'ICC/Support/Lapack.pm'});

	# make rhs matrix (filled with zeros)
	$rhs = $lapack ? ICC::Support::Lapack::zeros($ix + 1) : bless([map {[(0) x ($ix + 1)]} (0 .. $ix)], 'Math::Matrix');

	# for each row
	for my $i (1 .. $ix - 1) {
		
		# set rhs diagonal values
		$rhs->[$i - 1][$i] = 3;
		$rhs->[$i + 1][$i] = -3;
		
	}

	# set rhs endpoint values
	$rhs->[0][0] = -3;
	$rhs->[1][0] = -3;
	$rhs->[$ix][$ix] = 3;
	$rhs->[$ix - 1][$ix] = 3;

	# if ICC::Support::Lapack module is loaded
	if ($lapack) {
		
		# solve for derivative matrix
		($info, $derv) = ICC::Support::Lapack::trisolve([(1) x $ix], [2, (4) x ($ix - 1), 2], [(1) x $ix], $rhs);
		
		# set object
		$self->[5] = bless($derv, 'Math::Matrix');
		
	# otherwise, use Math::Matrix module
	} else {
		
		# solve for derivative matrix
		$self->[5] = Math::Matrix->tridiagonal([2, (4) x ($ix - 1), 2])->concat($rhs)->solve();
		
	}
	
}

# add local min/max values
# parameters: (ref_to_object)
# returns: (s-value_array)
sub _minmax {

	# get object reference
	my $self = shift();

	# local variables
	my (@t, @y, @s);

	# for each interval
	for my $i (0 .. $#{$self->[2]} - 1) {
		
		# get min/max location(s)
		@t = _local($self, $i);
		
		# compute local min/max y-values
		@y = map {_fwd($self, $_, $i)} @t;
		
		# add the knot values
		push(@y, ($self->[2][$i], $self->[2][$i + 1]));
		
		# sort the values
		@y = sort {$a <=> $b} @y;
		
		# save y-value min/max span
		$self->[4][$i] = [$y[0], $y[-1]];
		
		# add min/max s-values
		push(@s, map {$_ + $i} @t);
		
	}

	# return min/max s-values
	return(sort {$a <=> $b} @s);

}

# compute local min/max value(s)
# values are within the segment (t)
# parameters: (ref_to_object, segment)
# returns: (t-value_array)
sub _local {

	# get parameters
	my ($self, $low) = @_;

	# local variables
	my ($y1, $y2, $m1, $m2);
	my ($a, $b, $c, $dscr, @t);

	# get endpoint values
	$y1 = $self->[2][$low];
	$y2 = $self->[2][$low + 1];
	$m1 = $self->[3][$low];
	$m2 = $self->[3][$low + 1];

	# compute coefficients of quadratic equation  (at^2 + bt + c = 0)
	$a = 6 * ($y1 - $y2) + 3 * ($m1 + $m2);
	$b = -6 * ($y1 - $y2) - 2 * $m2 - 4 * $m1;
	$c = $m1;

	# return if constant
	return() if (abs($a) < 1E-15 && abs($b) < 1E-15);

	# if linear equation (a is zero)
	if (abs($a) < 1E-15) {
		
		# add solution
		push(@t, -$c/$b);
		
	# if quadratic equation
	} else {
		
		# compute discriminant
		$dscr = $b**2 - 4 * $a * $c;
		
		# if discriminant > zero
		if ($dscr > 0) {
			
			# add solutions (critical points)
			push(@t, -($b + sqrt($dscr))/(2 * $a));
			push(@t, -($b - sqrt($dscr))/(2 * $a));
			
		}
		
	}

	# return solution(s) within interval (0 < t < 0)
	return (grep {$_ > 0 && $_ < 1} @t);

}

# compute second derivative for unit spline segment
# parameters: (ref_to_object, t-value, segment)
# returns: (second_derivative)
sub _derv2 {

	# get parameters
	my ($self, $t, $low) = @_;

	# local variables
	my ($h00, $h01, $h10, $h11);

	# compute Hermite derivative coefficients
	$h00 = 12 * $t - 6;
	$h01 = - $h00;
	$h10 = 6 * $t - 4;
	$h11 = 6 * $t - 2;

	# interpolate value
	return($h00 * $self->[2][$low] + $h01 * $self->[2][$low + 1] + $h10 * $self->[3][$low] + $h11 * $self->[3][$low + 1]);

}

# compute derivative for unit spline segment
# parameters: (ref_to_object, t-value, segment)
# returns: (derivative)
sub _derv {

	# get parameters
	my ($self, $t, $low) = @_;

	# local variables
	my ($tc, $ttc, $h00, $h01, $h10, $h11);

	# if position is non-zero
	if ($t) {
		
		# compute Hermite derivative coefficients
		$tc = (1 - $t);
		$ttc = -3 * $t * $tc;
		$h00 = 2 * $ttc;
		$h01 = -2 * $ttc;
		$h10 = $ttc + $tc;
		$h11 = $ttc + $t;
		
		# interpolate value
		return($h00 * $self->[2][$low] + $h01 * $self->[2][$low + 1] + $h10 * $self->[3][$low] + $h11 * $self->[3][$low + 1]);
		
	} else {
		
		# use lower knot derivative value
		return($self->[3][$low]);
		
	}
	
}

# compute transform value
# parameters: (ref_to_object, t-value, segment)
# returns: (y-value)
sub _fwd {

	# get parameters
	my ($self, $t, $low) = @_;

	# local variables
	my ($tc, $h00, $h01, $h10, $h11);

	# check if ICC::Support::Lapack module is loaded
	state $lapack = defined($INC{'ICC/Support/Lapack.pm'});

	# if position is non-zero
	if ($t) {
		
		# if ICC::Support::Lapack module is loaded
		if ($lapack) {
			
			# interpolate value (Lapack XS)
			return(ICC::Support::Lapack::hermite($t, $self->[2][$low], $self->[2][$low + 1], $self->[3][$low], $self->[3][$low + 1]));
			
		} else {
			
			# compute Hermite coefficients
			$tc = 1 - $t;
			$h00 = (1 + 2 * $t) * $tc * $tc;
			$h01 = 1 - $h00;
			$h10 = $t * $tc * $tc;
			$h11 = -$t * $t * $tc;
			
			# interpolate value (no Lapack XS)
			return($h00 * $self->[2][$low] + $h01 * $self->[2][$low + 1] + $h10 * $self->[3][$low] + $h11 * $self->[3][$low + 1]);
			
		}
		
	} else {
		
		# use lower knot output value
		return($self->[2][$low]);
		
	}
	
}

# compute inverse value
# parameters: (ref_to_object, y-value, segment)
# there are 0 to 3 solutions in the range 0 < t < 1
# returns: (t-value_array)
sub _rev {

	# get parameters
	my ($self, $in, $low) = @_;

	# local variables
	my ($y1, $y2, $m1, $m2);
	my ($a, $b, $c, $d, $dscr, @t);
	my ($d0, $d1, $cs, $cc, @r, $ccr, $sol, $lim0, $lim1);

	# get endpoint values
	$y1 = $self->[2][$low];
	$y2 = $self->[2][$low + 1];
	$m1 = $self->[3][$low];
	$m2 = $self->[3][$low + 1];

	# compute coefficients of cubic equation (at^3 + bt^2 + ct + d = 0)
	$a = 2 * ($y1 - $y2) + $m1 + $m2;
	$b = -3 * ($y1 - $y2) -2 * $m1 - $m2; 
	$c = $m1;
	$d = $y1 - $in;

	# constant equation (a, b, c are zero)
	if (abs($a) < 5E-15 && abs($b) < 5E-15 && abs($c) < 5E-15) {
		
		# add solution
		push(@t, 0.5) if ($d == 0);
		
	# linear equation (a, b are zero)
	} elsif (abs($a) < 5E-15 && abs($b) < 5E-15) {
		
		# add solution
		push(@t, -$d/$c);
		
	# quadratic equation (a is zero)
	} elsif (abs($a) < 5E-15) {
		
		# compute discriminant: > 0 two real, == 0 one real, < 0 two complex
		$dscr = $c**2 - 4 * $b * $d;
		
		# if discriminant is zero
		if ($dscr == 0) {
			
			# add solution (double root)
			push(@t, -$c/(2 * $b));
			
		# if discriminant > zero
		} elsif ($dscr > 0) {
			
			# add solutions
			push(@t, -($c + sqrt($dscr))/(2 * $b));
			push(@t, -($c - sqrt($dscr))/(2 * $b));
			
		}
		
	# cubic equation
	} else {
		
		# compute discriminant: > 0 three real, == 0 one or two real, < 0 one real and two complex
		$dscr = 18 * $a * $b * $c * $d - 4 * $b**3 * $d + $b**2 * $c**2 - 4 * $a * $c**3 - 27 * $a**2 * $d**2;
		
		# compute ∆0
		$d0 = $b**2 - 3 * $a * $c;
		
		# if discriminant is zero
		if ($dscr == 0) {
			
			# if ∆0 is zero
			if ($d0 == 0) {
				
				# add solution (triple root)
				push(@t, -$b/(3 * $a));
				
			} else {
				
				# add solutions (double root and single root)
				push(@t, (9 * $a * $d - $b * $c)/(2 * $d0));
				push(@t, (4 * $a * $b * $c - 9 * $a**2 * $d - $b**3)/($a * $d0));
				
			}
			
		} else {
			
			# compute ∆1
			$d1 = 2 * $b**3 - 9 * $a * $b * $c + 27 * $a**2 * $d;
			
			# compute C (a complex number)
			$cs = sqrt($d1**2 - 4 * $d0**3);
			$cc = cbrt($d1 == $cs ? ($d1 + $cs)/2 : ($d1 - $cs)/2);
			
			# compute cube roots of 1 (three complex numbers)
			@r = root(1, 3);
			
			# for each root
			for my $i (0 .. 2) {
				
				# multiply by cube root of 1
				$ccr = $r[$i] * $cc;
				
				# compute solution (a complex number)
				$sol = ($b + $ccr + $d0/$ccr)/(-3 * $a);
				
				# add solution, if real
				push(@t, Re($sol)) if ($dscr > 0 || abs(Im($sol)) < 1E-15);
				
			}
			
		}
		
	}

	# set test limits to avoid duplicates at knots
	$lim0 = ($in == $y1) ? 1E-14 : 0;
	$lim1 = ($in == $y2) ? (1 - 1E-14) : 1;

	# return valid solutions
	return(sort {$ICC::Support::spline::a <=> $ICC::Support::spline::b} grep {$_ > $lim0 && $_ < $lim1} @t);

}

# set object contents from parameter hash
# parameters: (object_reference, ref_to_parameter_hash)
sub _new_from_hash {

	# get parameters
	my ($self, $hash) = @_;

	# local variables
	my ($output, $range);

	# if 'output_values' is defined
	if (defined($output = $hash->{'output_values'})) {
		
		# verify 'output_values' (two or more numeric values)
		(ref($output) eq 'ARRAY' && 2 <= @{$output} && (@{$output} == grep {Scalar::Util::looks_like_number($_)} @{$output})) || croak('invalid output values');
		
		# set 'output_values'
		$self->[2] = [@{$output}];
		
	} else {
		
		# use default
		$self->[2] = [0, 1];
		
	}

	# if 'input_range' is defined
	if (defined($range = $hash->{'input_range'})) {
		
		# verify 'input_range' (two unequal numeric values)
		(ref($range) eq 'ARRAY' && 2 == @{$range} && (2 == grep {Scalar::Util::looks_like_number($_)} @{$range}) && $range->[0] != $range->[1]) || croak('invalid input range');
		
		# set 'input_range'
		$self->[1] = [@{$range}];
		
	} else {
		
		# use default
		$self->[1] = [0, 1];
		
	}

	# compute knot derivatives
	_objderv($self);

	# add segment min/max y-values
	_minmax($self);

}

# write spline object to ICC profile
# note: writes an equivalent curv object
# parameters: (ref_to_object, ref_to_parent_object, file_handle, ref_to_tag_table_entry)
sub _writeICCspline {

	# get parameters
	my ($self, $parent, $fh, $tag) = @_;

	# local variables
	my ($n, $up, @table);

	# get count value (defaults to 4096, the maximum allowed in an 'mft2' tag)
	$n = $self->[0]{'curv_points'} // 4096;

	# validate number of table entries
	($n == int($n) && $n >= 2) || carp('invalid number of table entries');

	# array upper index
	$up = $n - 1;

	# for each table entry
	for my $i (0 .. $up) {
		
		# transform value
		$table[$i] = transform($self, $i/$up);
		
	}

	# seek start of tag
	seek($fh, $tag->[1], 0);

	# write tag type signature and count
	print $fh pack('a4 x4 N', 'curv', $n);

	# write array
	print $fh pack('n*', map {$_ < 0 ? 0 : ($_ > 1 ? 65535 : $_ * 65535 + 0.5)} @table);

}

1;